[gegl] warp: rework the operation, now matches more closely the iwarp behaviour

[Simon Budig <simon src gnome org>]

commit 67553d4f5dfdaf17124f9413b4d03e0dbafe2e89
Author: Simon Budig <simon budig de>
Date:   Sun Sep 25 01:50:11 2016 +0200

    warp: rework the operation, now matches more closely the iwarp behaviour

 operations/common/warp.c |  133 ++++++++++++++++++++++++++++++++-------------
 1 files changed, 94 insertions(+), 39 deletions(-)
---
diff --git a/operations/common/warp.c b/operations/common/warp.c
index bd34443..03417a7 100644
--- a/operations/common/warp.c
+++ b/operations/common/warp.c
@@ -185,7 +185,7 @@ get_stamp_force (GeglProperties *o,
       calc_lut (o);
     }
 
-  radius = sqrt(x*x+y*y);
+  radius = hypot (x, y);
 
   if (radius < 0.5 * o->size + 1)
     {
@@ -219,10 +219,17 @@ stamp (GeglProperties          *o,
   gdouble              x_mean = 0.0;
   gdouble              y_mean = 0.0;
   gint                 x_iter, y_iter;
-  GeglRectangle        area = {x - o->size / 2.0,
-                               y - o->size / 2.0,
-                               o->size,
-                               o->size};
+  GeglRectangle        area;
+  const GeglRectangle *src_extent;
+  gfloat              *srcbuf, *stampbuf;
+  gint                 buf_rowstride = 0;
+
+  area.x = floor (x - o->size / 2.0);
+  area.y = floor (y - o->size / 2.0);
+  area.width   = ceil (x + o->size / 2.0);
+  area.height  = ceil (y + o->size / 2.0);
+  area.width  -= area.x;
+  area.height -= area.y;
 
   /* first point of the stroke */
   if (!priv->last_point_set)
@@ -264,68 +271,116 @@ stamp (GeglProperties          *o,
       y_mean /= pixel_count;
     }
 
-  it = gegl_buffer_iterator_new (priv->buffer, &area, 0, format,
-                                 GEGL_ACCESS_READWRITE, GEGL_ABYSS_NONE);
+  srcbuf = gegl_buffer_linear_open (priv->buffer, NULL, &buf_rowstride, NULL);
+  buf_rowstride /= sizeof (gfloat);
+  src_extent = gegl_buffer_get_extent (priv->buffer);
 
-  while (gegl_buffer_iterator_next (it))
+  stampbuf = g_new0 (gfloat, 2 * area.height * area.width);
+
+  for (y_iter = 0; y_iter < area.height; y_iter++)
     {
-      /* iterate inside the stamp roi */
-      gint    n_pixels = it->length;
-      gfloat *coords   = it->data[0];
+      for (x_iter = 0; x_iter < area.width; x_iter++)
+        {
+          gfloat nvx, nvy;
+          gfloat xi, yi;
+          gfloat *vals;
+          gint dx, dy;
+          gfloat weight_x, weight_y;
+          gfloat *srcptr;
 
-      x_iter = it->roi->x; /* initial x         */
-      y_iter = it->roi->y; /* and y coordinates */
+          xi = area.x + x_iter;
+          xi += -x + 0.5;
+          yi = area.y + y_iter;
+          yi += -y + 0.5;
 
-      while (n_pixels--)
-        {
-          influence = 0.01 * o->strength * get_stamp_force (o,
-                                                            x_iter - x,
-                                                            y_iter - y);
+          influence = 0.01 * o->strength * get_stamp_force (o, xi, yi);
 
           switch (o->behavior)
             {
               case GEGL_WARP_BEHAVIOR_MOVE:
-                coords[0] += influence * (priv->last_x - x);
-                coords[1] += influence * (priv->last_y - y);
+                nvx =  influence * (priv->last_x - x);
+                nvy =  influence * (priv->last_y - y);
                 break;
               case GEGL_WARP_BEHAVIOR_GROW:
-                coords[0] -= influence * 0.1 * (x_iter - x);
-                coords[1] -= influence * 0.1 * (y_iter - y);
+                nvx = -influence * 0.1 * xi;
+                nvy = -influence * 0.1 * yi;
                 break;
               case GEGL_WARP_BEHAVIOR_SHRINK:
-                coords[0] += influence * 0.1 * (x_iter - x);
-                coords[1] += influence * 0.1 * (y_iter - y);
+                nvx =  influence * 0.1 * xi;
+                nvy =  influence * 0.1 * yi;
                 break;
               case GEGL_WARP_BEHAVIOR_SWIRL_CW:
-                coords[0] += 3.0 * influence * 0.1 * (y_iter - y);
-                coords[1] -= 5.0 * influence * 0.1 * (x_iter - x);
+                nvx =  influence * 0.1 * yi;
+                nvy = -influence * 0.1 * xi;
                 break;
               case GEGL_WARP_BEHAVIOR_SWIRL_CCW:
-                coords[0] -= 3.0 * influence * 0.1 * (y_iter - y);
-                coords[1] += 5.0 * influence * 0.1 * (x_iter - x);
+                nvx = -influence * 0.1 * yi;
+                nvy =  influence * 0.1 * xi;
                 break;
               case GEGL_WARP_BEHAVIOR_ERASE:
-                coords[0] *= 1.0 - MIN (influence, 1.0);
-                coords[1] *= 1.0 - MIN (influence, 1.0);
-                break;
               case GEGL_WARP_BEHAVIOR_SMOOTH:
-                coords[0] -= influence * (coords[0] - x_mean);
-                coords[1] -= influence * (coords[1] - y_mean);
+              default:
+                nvx = 0.0;
+                nvy = 0.0;
                 break;
             }
 
-          coords += 2;
+          vals = stampbuf + (y_iter * area.width + x_iter) * 2;
+
+          dx = floorf (nvx);
+          dy = floorf (nvy);
+
+          if (area.x + x_iter + dx     <  src_extent->x                     ||
+              area.x + x_iter + dx + 1 >= src_extent->x + src_extent->width ||
+              area.y + y_iter + dy     <  src_extent->y                     ||
+              area.y + y_iter + dy + 1 >= src_extent->y + src_extent->height)
+            {
+              continue;
+            }
+
+          srcptr = srcbuf + (area.y + y_iter + dy) * buf_rowstride +
+                            (area.x + x_iter + dx) * 2;
 
-          /* update x and y coordinates */
-          x_iter++;
-          if (x_iter >= (it->roi->x + it->roi->width))
+          if (o->behavior == GEGL_WARP_BEHAVIOR_ERASE)
+            {
+              vals[0] = srcptr[0] * (1.0 - MIN (influence, 1.0));
+              vals[1] = srcptr[1] * (1.0 - MIN (influence, 1.0));
+            }
+          else if (o->behavior == GEGL_WARP_BEHAVIOR_SMOOTH)
+            {
+              vals[0] = (1.0 - influence) * srcptr[0] + influence * x_mean;
+              vals[1] = (1.0 - influence) * srcptr[1] + influence * y_mean;
+            }
+          else
             {
-              x_iter = it->roi->x;
-              y_iter++;
+              weight_x = nvx - dx;
+              weight_y = nvy - dy;
+
+              /* bilinear interpolation of the vectors */
+
+              vals[0]  = srcptr[0] * (1.0 - weight_x) * (1.0 - weight_y);
+              vals[1]  = srcptr[1] * (1.0 - weight_x) * (1.0 - weight_y);
+
+              vals[0] += srcptr[2] * weight_x * (1.0 - weight_y);
+              vals[1] += srcptr[3] * weight_x * (1.0 - weight_y);
+
+              vals[0] += srcptr[buf_rowstride + 0] * (1.0 - weight_x) * weight_y;
+              vals[1] += srcptr[buf_rowstride + 1] * (1.0 - weight_x) * weight_y;
+
+              vals[0] += srcptr[buf_rowstride + 2] * weight_x * weight_y;
+              vals[1] += srcptr[buf_rowstride + 3] * weight_x * weight_y;
+
+              vals[0] += nvx;
+              vals[1] += nvy;
             }
         }
     }
 
+  gegl_buffer_linear_close (priv->buffer, srcbuf);
+  gegl_buffer_set (priv->buffer, &area, 0, format,
+                   stampbuf, GEGL_AUTO_ROWSTRIDE);
+  g_free (stampbuf);
+
   /* Memorize the stamp location for movement dependant behavior like move */
   priv->last_x = x;
   priv->last_y = y;