[gegl] convolution-matrix: refactor conditionals and common expressions out of inner loops



commit 1b365c85b8737f960e6d42dce5de86abbffc9d13
Author: Øyvind Kolås <pippin gimp org>
Date:   Sun Feb 14 19:31:54 2016 +0100

    convolution-matrix: refactor conditionals and common expressions out of inner loops

 operations/common/convolution-matrix.c |  470 +++++++++++++++++++++++++-------
 1 files changed, 368 insertions(+), 102 deletions(-)
---
diff --git a/operations/common/convolution-matrix.c b/operations/common/convolution-matrix.c
index 5a66cbf..b8cd27e 100644
--- a/operations/common/convolution-matrix.c
+++ b/operations/common/convolution-matrix.c
@@ -75,16 +75,43 @@ property_enum (border, _("Border"),
 #include <math.h>
 #include <stdio.h>
 
-#define MATRIX_SIZE   (5)
-#define HALF_WINDOW   (MATRIX_SIZE/2)
-#define CHANNELS      (5)
+#define MAX_MATRIX_SIZE 5
+
+static gboolean
+enough_with_3x3 (GeglProperties *o)
+{
+  if (o->a1 == 0.0 &&
+      o->a2 == 0.0 &&
+      o->a3 == 0.0 &&
+      o->a4 == 0.0 &&
+      o->a5 == 0.0 &&
+      o->b1 == 0.0 &&
+      o->b5 == 0.0 &&
+      o->c1 == 0.0 &&
+      o->c5 == 0.0 &&
+      o->d1 == 0.0 &&
+      o->d5 == 0.0 &&
+      o->e1 == 0.0 &&
+      o->e2 == 0.0 &&
+      o->e3 == 0.0 &&
+      o->e4 == 0.0 &&
+      o->e5 == 0.0)
+  {
+    return TRUE;
+  }
+  return FALSE;
+}
 
 static void
 prepare (GeglOperation *operation)
 {
   GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
 
-  op_area->left = op_area->right = op_area->top = op_area->bottom = HALF_WINDOW;
+  GeglProperties          *o       = GEGL_PROPERTIES (operation);
+  if (enough_with_3x3 (o))
+    op_area->left = op_area->right = op_area->top = op_area->bottom = 1; /* 3 */
+  else
+    op_area->left = op_area->right = op_area->top = op_area->bottom = 2; /* 5 */
 
   gegl_operation_set_format (operation, "output",
                              babl_format ("RGBA float"));
@@ -92,49 +119,68 @@ prepare (GeglOperation *operation)
 
 static void
 make_matrix (GeglProperties  *o,
-             gdouble        **matrix)
+             gfloat           matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
+             gint             matrix_size)
 {
-  matrix[0][0] = o->a1;
-  matrix[0][1] = o->a2;
-  matrix[0][2] = o->a3;
-  matrix[0][3] = o->a4;
-  matrix[0][4] = o->a5;
-
-  matrix[1][0] = o->b1;
-  matrix[1][1] = o->b2;
-  matrix[1][2] = o->b3;
-  matrix[1][3] = o->b4;
-  matrix[1][4] = o->b5;
-
-  matrix[2][0] = o->c1;
-  matrix[2][1] = o->c2;
-  matrix[2][2] = o->c3;
-  matrix[2][3] = o->c4;
-  matrix[2][4] = o->c5;
-
-  matrix[3][0] = o->d1;
-  matrix[3][1] = o->d2;
-  matrix[3][2] = o->d3;
-  matrix[3][3] = o->d4;
-  matrix[3][4] = o->d5;
-
-  matrix[4][0] = o->e1;
-  matrix[4][1] = o->e2;
-  matrix[4][2] = o->e3;
-  matrix[4][3] = o->e4;
-  matrix[4][4] = o->e5;
+  if (matrix_size == 3)
+  {
+    matrix[0][0] = o->b2;
+    matrix[0][1] = o->b3;
+    matrix[0][2] = o->b4;
+
+    matrix[1][0] = o->c2;
+    matrix[1][1] = o->c3;
+    matrix[1][2] = o->c4;
+
+    matrix[2][0] = o->d2;
+    matrix[2][1] = o->d3;
+    matrix[2][2] = o->d4;
+  }
+  else
+  {
+    matrix[0][0] = o->a1;
+    matrix[0][1] = o->a2;
+    matrix[0][2] = o->a3;
+    matrix[0][3] = o->a4;
+    matrix[0][4] = o->a5;
+
+    matrix[1][0] = o->b1;
+    matrix[1][1] = o->b2;
+    matrix[1][2] = o->b3;
+    matrix[1][3] = o->b4;
+    matrix[1][4] = o->b5;
+
+    matrix[2][0] = o->c1;
+    matrix[2][1] = o->c2;
+    matrix[2][2] = o->c3;
+    matrix[2][3] = o->c4;
+    matrix[2][4] = o->c5;
+
+    matrix[3][0] = o->d1;
+    matrix[3][1] = o->d2;
+    matrix[3][2] = o->d3;
+    matrix[3][3] = o->d4;
+    matrix[3][4] = o->d5;
+
+    matrix[4][0] = o->e1;
+    matrix[4][1] = o->e2;
+    matrix[4][2] = o->e3;
+    matrix[4][3] = o->e4;
+    matrix[4][4] = o->e5;
+  }
 }
 
 static gboolean
 normalize_o (GeglProperties  *o,
-             gdouble        **matrix)
+             gfloat           matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
+             gint             matrix_size)
 {
   gint      x, y;
   gboolean  valid = FALSE;
   gfloat    sum   = 0.0;
 
-  for (y = 0; y < MATRIX_SIZE; y++)
-    for (x = 0; x < MATRIX_SIZE; x++)
+  for (y = 0; y < matrix_size; y++)
+    for (x = 0; x < matrix_size; x++)
       {
         sum += matrix[x][y];
         if (matrix[x][y] != 0.0)
@@ -160,29 +206,27 @@ normalize_o (GeglProperties  *o,
   return valid;
 }
 
-static void
-convolve_pixel (GeglProperties       *o,
+static void inline
+convolve_pixel_componentwise (GeglProperties       *o,
                 gfloat               *src_buf,
                 gfloat               *dst_buf,
                 const GeglRectangle  *result,
                 const GeglRectangle  *extended,
-                gdouble             **matrix,
+                gfloat                matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
+                gint                  matrix_size,
+                gint                  d_offset,
+                gint                  ss_offset,
                 gint                  xx,
                 gint                  yy,
-                gdouble               matrixsum)
+                gfloat                matrixsum,
+                gfloat                inv_divisor)
 {
-  gfloat  color[4];
-  gint    d_offset;
-  gint    s_offset;
-  gint    i;
-
-  d_offset = (yy - result->y) * result->width * 4 +
-             (xx - result->x) * 4;
-
+  gint   i;
+  gint   s_stride = (extended->width - matrix_size) * 4;
   for (i = 0; i < 4; i++)
     {
-      gdouble sum      = 0.0;
-      gdouble alphasum = 0.0;
+      gfloat sum = 0.0;
+      gint s_offset = ss_offset;
 
       if ((i == 0 && o->red)   ||
           (i == 1 && o->blue)  ||
@@ -190,50 +234,208 @@ convolve_pixel (GeglProperties       *o,
           (i == 3 && o->alpha))
         {
           gint x, y;
-
-          for (x = 0; x < MATRIX_SIZE; x++)
-            for (y = 0; y < MATRIX_SIZE; y++)
+          for (y = 0; y < matrix_size; y++)
+          {
+            for (x = 0; x < matrix_size; x++)
               {
-                gint s_x = x + xx - HALF_WINDOW;
-                gint s_y = y + yy - HALF_WINDOW;
+                sum += matrix[x][y] * src_buf[s_offset + i];
+                s_offset += 4;
+              }
+            s_offset += s_stride;
+          }
+          sum = sum * inv_divisor;
+          sum += o->offset;
+          dst_buf[d_offset + i] = sum;
+        }
+      else
+        {
+          dst_buf[d_offset + i] = src_buf[s_offset + i];
+        }
+    }
+}
 
-                s_offset = (s_y - extended->y) * extended->width * 4 +
-                           (s_x - extended->x) * 4;
+static void inline
+convolve_pixel (GeglProperties       *o,
+                gfloat               *src_buf,
+                gfloat               *dst_buf,
+                const GeglRectangle  *result,
+                const GeglRectangle  *extended,
+                gfloat                matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
+                gint                  matrix_size,
+                gint                  d_offset,
+                gint                  ss_offset,
+                gint                  xx,
+                gint                  yy,
+                gfloat                matrixsum,
+                gfloat                inv_divisor)
+{
+  gint   i;
+  gint   s_stride = (extended->width - matrix_size) * 4;
+  for (i = 0; i < 4; i++)
+    {
+      gfloat sum = 0.0;
+      gint s_offset = ss_offset;
+      gint x, y;
+      for (y = 0; y < matrix_size; y++)
+      {
+        for (x = 0; x < matrix_size; x++)
+          {
+            sum += matrix[x][y] * src_buf[s_offset + i];
+            s_offset += 4;
+          }
+        s_offset += s_stride;
+      }
+      sum = sum * inv_divisor;
+      sum += o->offset;
+      dst_buf[d_offset + i] = sum;
+    }
+}
 
-                if (i != 3 && o->alpha_weight)
-                  sum += matrix[x][y] * src_buf[s_offset + i] * src_buf[s_offset + 3];
-                else
-                  sum += matrix[x][y] * src_buf[s_offset + i];
+static void inline
+convolve_pixel_alpha_weight (GeglProperties       *o,
+                             gfloat               *src_buf,
+                             gfloat               *dst_buf,
+                             const GeglRectangle  *result,
+                             const GeglRectangle  *extended,
+                             gfloat                matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
+                             gint                  matrix_size,
+                             gint                  d_offset,
+                             gint                  ss_offset,
+                             gint                  xx,
+                             gint                  yy,
+                             gfloat                matrixsum,
+                             gfloat                inv_divisor)
+{
+  gint   i;
+  gint   s_stride = (extended->width - matrix_size) * 4;
 
-                if (i == 3)
-                  alphasum += fabs (matrix[x][y] * src_buf[s_offset + i]);
-              }
+  for (i = 0; i < 3; i++)
+    {
+      gfloat sum      = 0.0;
+      gint   s_offset = ss_offset;
+      gint x, y;
+      for (y = 0; y < matrix_size; y++)
+      {
+        for (x = 0; x < matrix_size; x++)
+          {
+            sum += matrix[x][y] * src_buf[s_offset + i] * src_buf[s_offset + 3];
+            s_offset += 4;
+          }
+        s_offset += s_stride;
+      }
+      sum = sum * inv_divisor + o->offset;
+      dst_buf[d_offset + i] = sum;
+    }
+    {
+      gfloat sum      = 0.0;
+      gfloat alphasum = 0.0;
+      gint   s_offset = ss_offset;
+      gint x, y;
 
-          sum = sum / o->divisor;
+      for (y = 0; y < matrix_size; y++)
+      {
+        for (x = 0; x < matrix_size; x++)
+          {
+            float val = matrix[x][y] * src_buf[s_offset + i];
+            sum += val;
+            alphasum += fabsf (val);
+            s_offset += 4;
+          }
+        s_offset += s_stride;
+      }
 
-          if (i == 3 && o->alpha_weight)
-            {
-              if (alphasum != 0)
-                sum = sum * matrixsum / alphasum;
-              else
-                sum = 0.0;
-            }
+      if (alphasum > 0.0)
+      {
+        sum = sum * inv_divisor;
+        sum = sum * matrixsum / alphasum + o->offset;
+      }
+      else
+        sum = o->offset;
+      dst_buf[d_offset + i] = sum;
+    }
+}
 
-          sum += o->offset;
+static void inline
+convolve_pixel_alpha_weight_componentwise (GeglProperties       *o,
+                             gfloat               *src_buf,
+                             gfloat               *dst_buf,
+                             const GeglRectangle  *result,
+                             const GeglRectangle  *extended,
+                             gfloat                matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
+                             gint                  matrix_size,
+                             gint                  d_offset,
+                             gint                  ss_offset,
+                             gint                  xx,
+                             gint                  yy,
+                             gfloat                matrixsum,
+                             gfloat                inv_divisor)
+{
+  gint   i;
+  gint   s_stride = (extended->width - matrix_size) * 4;
+
+  for (i = 0; i < 3; i++)
+    {
+      gfloat sum      = 0.0;
+      gint   s_offset = ss_offset;
 
-          color[i] = sum;
+      if ((i == 0 && o->red)   ||
+          (i == 1 && o->blue)  ||
+          (i == 2 && o->green))
+        {
+          gint x, y;
+          for (y = 0; y < matrix_size; y++)
+          {
+            for (x = 0; x < matrix_size; x++)
+              {
+                sum += matrix[x][y] * src_buf[s_offset + i] * src_buf[s_offset + 3];
+                s_offset += 4;
+              }
+            s_offset += s_stride;
+          }
+          sum = sum * inv_divisor + o->offset;
         }
       else
         {
-          s_offset = (yy - result->y + HALF_WINDOW) * extended->width * 4 +
-                     (xx - result->x + HALF_WINDOW) * 4;
-
-          color[i] = src_buf[s_offset + i];
+          sum = src_buf[s_offset + i];
         }
+      dst_buf[d_offset + i] = sum;
     }
+    {
+      gfloat sum      = 0.0;
+      gfloat alphasum = 0.0;
+      gint   s_offset = ss_offset;
 
-  for (i = 0; i < 4; i++)
-    dst_buf[d_offset + i] = color[i];
+      if (o->alpha)
+        {
+          gint x, y;
+
+          for (y = 0; y < matrix_size; y++)
+          {
+            for (x = 0; x < matrix_size; x++)
+              {
+                float val = matrix[x][y] * src_buf[s_offset + i];
+                sum += val;
+                alphasum += fabsf (val);
+                s_offset += 4;
+              }
+            s_offset += s_stride;
+          }
+
+
+          if (alphasum != 0)
+          {
+            sum = sum * inv_divisor;
+            sum = sum * matrixsum / alphasum + o->offset;
+          }
+          else
+            sum = o->offset;
+        }
+      else
+        {
+          sum = src_buf[s_offset + i];
+        }
+      dst_buf[d_offset + i] = sum;
+    }
 }
 
 static gboolean
@@ -249,44 +451,109 @@ process (GeglOperation       *operation,
   GeglRectangle            rect;
   gfloat                  *src_buf;
   gfloat                  *dst_buf;
-  gdouble                **matrix;
-  gdouble                  matrixsum = 0.0;
+  gfloat                   matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE]={{0,}};
+  gfloat                   matrixsum = 0.0;
+  gfloat                   inv_divisor;
   gint                     x, y;
+  gint                     matrix_size = MAX_MATRIX_SIZE;
 
-  matrix = g_new0 (gdouble *, MATRIX_SIZE);
+  if (enough_with_3x3 (o))
+    matrix_size = 3;
 
-  for (x = 0; x < MATRIX_SIZE ;x++)
-    matrix[x] = g_new0 (gdouble, MATRIX_SIZE);
-
-  make_matrix (o, matrix);
+  make_matrix (o, matrix, matrix_size);
 
   if (o->normalize)
-    normalize_o (o, matrix);
+    normalize_o (o, matrix, matrix_size);
+  inv_divisor = 1.0 / o->divisor;
 
-  for (x = 0; x < MATRIX_SIZE; x++)
-    for (y = 0; y < MATRIX_SIZE; y++)
-      matrixsum += fabs (matrix[x][y]);
+  for (x = 0; x < matrix_size; x++)
+    for (y = 0; y < matrix_size; y++)
+      matrixsum += fabsf (matrix[x][y]);
 
   rect.x      = result->x - op_area->left;
   rect.width  = result->width + op_area->left + op_area->right;
   rect.y      = result->y - op_area->top;
   rect.height = result->height + op_area->top + op_area->bottom;
 
-  src_buf = g_new0 (gfloat, rect.width * rect.height * 4);
-  dst_buf = g_new0 (gfloat, result->width * result->height * 4);
+  src_buf = g_new (gfloat, rect.width * rect.height * 4);
+  dst_buf = g_new (gfloat, result->width * result->height * 4);
 
   gegl_buffer_get (input, &rect, 1.0, format, src_buf,
                    GEGL_AUTO_ROWSTRIDE, o->border);
 
   if (o->divisor != 0)
     {
-      for (y = result->y; y < result->height + result->y; y++)
-        for (x = result->x; x < result->width + result->x; x++)
-          convolve_pixel (o, src_buf, dst_buf, result, &rect,
-                          matrix, x, y, matrixsum);
-
-      gegl_buffer_set (output, result, 0, format,
-                       dst_buf, GEGL_AUTO_ROWSTRIDE);
+      gint ss_offset = (result->y - matrix_size/2 - rect.y) * rect.width * 4 +
+                       (result->x - matrix_size/2 - rect.x) * 4;
+      int d_offset = 0;
+      if (o->alpha_weight)
+        {
+          if (o->red == FALSE || o->green == FALSE || o->blue == FALSE || o->alpha == FALSE)
+          {
+            gint x, y;
+            for (y = result->y; y < result->height + result->y; y++)
+            {
+              for (x = result->x; x < result->width + result->x; x++)
+              {
+                convolve_pixel_alpha_weight_componentwise (o, src_buf, dst_buf, result, &rect,
+                                           matrix, matrix_size, d_offset, ss_offset, x, y, matrixsum, 
inv_divisor);
+                d_offset += 4;
+                ss_offset += 4;
+              }
+              ss_offset += (rect.width - result->width) * 4;
+            }
+          }
+          else
+          {
+            gint x, y;
+            for (y = result->y; y < result->height + result->y; y++)
+            {
+              for (x = result->x; x < result->width + result->x; x++)
+              {
+                convolve_pixel_alpha_weight (o, src_buf, dst_buf, result, &rect,
+                                           matrix, matrix_size, d_offset, ss_offset, x, y, matrixsum, 
inv_divisor);
+                d_offset += 4;
+                ss_offset += 4;
+              }
+              ss_offset += (rect.width - result->width) * 4;
+            }
+          }
+        }
+      else
+        {
+          if (o->red == FALSE || o->green == FALSE || o->blue == FALSE || o->alpha == FALSE)
+          {
+            gint x, y;
+            for (y = result->y; y < result->height + result->y; y++)
+            {
+              for (x = result->x; x < result->width + result->x; x++)
+              {
+                convolve_pixel_componentwise (o, src_buf, dst_buf, result, &rect,
+                                matrix, matrix_size, d_offset, ss_offset, x, y, matrixsum, inv_divisor);
+                d_offset += 4;
+                ss_offset += 4;
+              }
+              ss_offset += (rect.width - result->width) * 4;
+            }
+          }
+          else
+          {
+            gint x, y;
+            for (y = result->y; y < result->height + result->y; y++)
+            {
+              for (x = result->x; x < result->width + result->x; x++)
+              {
+                convolve_pixel (o, src_buf, dst_buf, result, &rect,
+                                matrix, matrix_size, d_offset, ss_offset, x, y, matrixsum, inv_divisor);
+                d_offset += 4;
+                ss_offset += 4;
+              }
+              ss_offset += (rect.width - result->width) * 4;
+            }
+          }
+        }
+        gegl_buffer_set (output, result, 0, format,
+                         dst_buf, GEGL_AUTO_ROWSTRIDE);
     }
   else
     {
@@ -329,7 +596,6 @@ gegl_op_class_init (GeglOpClass *klass)
 
   operation_class = GEGL_OPERATION_CLASS (klass);
   filter_class    = GEGL_OPERATION_FILTER_CLASS (klass);
-
   filter_class->process                    = process;
   operation_class->prepare                 = prepare;
   operation_class->get_bounding_box        = get_bounding_box;


[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]