[gegl] operations: attenpt to fix lens-distortion.c



commit 1962d2daa6bea331a474e69ec2f768551432aea4
Author: Téo Mazars <teo mazars ensimag fr>
Date:   Thu Oct 10 16:05:22 2013 +0200

    operations: attenpt to fix lens-distortion.c
    
    - Never allocate a too large buffer
    - Give it GeglOperationAreaFilter's flavors
    - Improve the use of the allocated buffer
    - Format fixes and some renaming for clarity
    
    Still not perfect, get_required_for_output may not
    be 100% bullet-proof with corner values

 operations/common/lens-distortion.c |  328 +++++++++++++++++++++++++----------
 1 files changed, 234 insertions(+), 94 deletions(-)
---
diff --git a/operations/common/lens-distortion.c b/operations/common/lens-distortion.c
index ee2ef47..db614e8 100644
--- a/operations/common/lens-distortion.c
+++ b/operations/common/lens-distortion.c
@@ -13,15 +13,15 @@
  * You should have received a copy of the GNU Lesser General Public
  * License along with GEGL; if not, see <http://www.gnu.org/licenses/>.
  *
+ * Lens plug-in - adjust for lens distortion
+ *
  * Copyright (C) 1995 Spencer Kimball and Peter Mattis
+ * Copyright (C) 2001-2005 David Hodson <hodsond acm org>
+ * Copyright (C) 2011 Robert Sasu <sasu robert gmail com>
+ * Copyright (C) 2013 Téo Mazars <teo mazars ensimag fr>
  *
- * Lens plug-in - adjust for lens distortion
- * Copyright (C) 2001-2005 David Hodson hodsond acm org
  * Many thanks for Lars Clausen for the original inspiration,
  *   useful discussion, optimisation and improvements.
- * Framework borrowed from many similar plugins...
- *
- * Copyright (C) 2011 Robert Sasu <sasu robert gmail com>
  */
 
 #include "config.h"
@@ -36,6 +36,7 @@ gegl_chant_double (main, _("Main"),
 gegl_chant_double (zoom, _("Zoom"),
                    -100.0, 100.0, 0.0,
                    _("Main value of distortion"))
+
 gegl_chant_double (edge, _("Edge"),
                    -100.0, 100.0, 0.0,
                    _("Edge value of distortion"))
@@ -59,6 +60,17 @@ gegl_chant_double (y_shift, _("Y shift"),
 
 #include "gegl-chant.h"
 #include <math.h>
+#include <stdio.h>
+
+#define SQR(x) ((x)*(x))
+
+#define MIN3(x,y,z) (MIN (MIN ((x),(y)),(z)))
+
+#define MAX3(x,y,z) (MAX (MAX ((x),(y)),(z)))
+
+#define MAX_WH     1024
+#define CHUNK_SIZE 512
+
 
 typedef struct
 {
@@ -69,76 +81,186 @@ typedef struct
   gdouble rescale;
   gdouble brighten;
   gdouble norm;
-} LensDistortion;
+} LensValues;
+
 
 static void
-lens_setup_calc (GeglChantO     *o,
-                 GeglRectangle   boundary,
-                 LensDistortion *old)
+reorder (gdouble *low,
+         gdouble *high)
 {
-  old->norm = 4.0 / (boundary.width * boundary.width +
-                     boundary.height * boundary.height);
-
-  old->centre_x = boundary.width * (100.0 + o->x_shift) / 200.0;
-  old->centre_y = boundary.height * (100.0 + o->y_shift) / 200.0;
-  old->mult_sq = o->main / 200.0;
-  old->mult_qd = o->edge / 200.0;
-  old->rescale = pow (2.0, -o->zoom / 100.0);
-  old->brighten = -o->brighten / 10.0;
+  gdouble temp;
+
+  if (*low < *high) return;
+
+  temp = *low;
+  *low = *high;
+  *high = temp;
 }
 
-static GeglRectangle
-get_bounding_box (GeglOperation *operation)
+static LensValues
+lens_setup_calc (GeglChantO    *o,
+                 GeglRectangle  boundary)
 {
-  GeglRectangle  result = {0,0,0,0};
-  GeglRectangle *in_rect;
+  LensValues lens;
 
-  in_rect = gegl_operation_source_get_bounding_box (operation, "input");
-  if (!in_rect)
-    return result;
+  lens.norm     = 4.0 / (SQR (boundary.width) + SQR (boundary.height));
+  lens.centre_x = boundary.width  * (100.0 + o->x_shift) / 200.0;
+  lens.centre_y = boundary.height * (100.0 + o->y_shift) / 200.0;
+  lens.mult_sq  = o->main / 200.0;
+  lens.mult_qd  = o->edge / 200.0;
+  lens.rescale  = pow (2.0, - o->zoom / 100.0);
+  lens.brighten = - o->brighten / 10.0;
 
-  return *in_rect;
+  return lens;
 }
 
-static GeglRectangle
-get_required_for_output (GeglOperation       *operation,
-                         const gchar         *input_pad,
-                         const GeglRectangle *roi)
+static void
+lens_get_source_coord (gdouble     i,
+                       gdouble     j,
+                       gdouble    *x,
+                       gdouble    *y,
+                       gdouble    *mag,
+                       LensValues *lens)
 {
-  return get_bounding_box (operation);
+  gdouble radius_sq, off_x, off_y, radius_mult;
+
+  off_x = i - lens->centre_x;
+  off_y = j - lens->centre_y;
+
+  radius_sq = SQR (off_x) + SQR (off_y);
+
+  radius_sq *= lens->norm;
+
+  radius_mult = radius_sq * lens->mult_sq + SQR (radius_sq) * lens->mult_qd;
+
+  *mag = radius_mult;
+
+  radius_mult = lens->rescale * (1.0 + radius_mult);
+
+  *x = lens->centre_x + radius_mult * off_x;
+  *y = lens->centre_y + radius_mult * off_y;
 }
 
-static void
-prepare (GeglOperation *operation)
+/* FIXME: not 100% bullet proof */
+static GeglRectangle
+get_required (GeglRectangle       *boundary,
+              const GeglRectangle *roi,
+              GeglOperation       *operation)
 {
-  gegl_operation_set_format (operation, "input", babl_format ("RGBA float"));
-  gegl_operation_set_format (operation, "output", babl_format ("RGBA float"));
+
+  GeglChantO    *o;
+  GeglRectangle  area;
+  LensValues     lens;
+  gdouble        x1, y1, x2, y2, x3, y3, x4, y4, mag;
+  gint           x, y, width, height;
+
+  o = GEGL_CHANT_PROPERTIES (operation);
+
+  lens = lens_setup_calc (o, *boundary);
+
+  x = roi->x;
+  y = roi->y;
+  width = roi->width;
+  height = roi->height;
+
+  lens_get_source_coord (x,         y,          &x1, &y1, &mag, &lens);
+  lens_get_source_coord (x + width, y,          &x2, &y2, &mag, &lens);
+  lens_get_source_coord (x,         y + height, &x3, &y3, &mag, &lens);
+  lens_get_source_coord (x + width, y + height, &x4, &y4, &mag, &lens);
+
+  /* This is ugly, and happens
+   * with a crazy set of parameters */
+  reorder (&x1, &x2);
+  reorder (&x3, &x4);
+
+  reorder (&y1, &y3);
+  reorder (&y2, &y4);
+
+  if (lens.centre_y > y && lens.centre_y < y + height)
+    {
+      gdouble x5, y5, x6, y6;
+
+      lens_get_source_coord (x,         lens.centre_y, &x5, &y5, &mag, &lens);
+      lens_get_source_coord (x + width, lens.centre_y, &x6, &y6, &mag, &lens);
+
+      reorder (&x5, &x6);
+
+      area.x = floor (MIN3 (x1, x3, x5)) - 1;
+      area.width = ceil (MAX3 (x2, x4, x6)) + 3 - area.x;
+    }
+  else
+    {
+      area.x = floor (MIN (x1, x3)) - 1;
+      area.width = ceil (MAX (x2, x4)) + 3 - area.x;
+    }
+
+  if (lens.centre_x > x && lens.centre_x < x + width)
+    {
+      gdouble x5, y5, x6, y6;
+
+      lens_get_source_coord (lens.centre_x, y,          &x5, &y5, &mag, &lens);
+      lens_get_source_coord (lens.centre_x, y + height, &x6, &y6, &mag, &lens);
+
+      reorder (&y5, &y6);
+
+      area.y = floor (MIN3 (y1, y2, y5)) - 1;
+      area.height = ceil (MAX3 (y3, y4, y6)) + 3 - area.y;
+    }
+  else
+    {
+      area.y = floor (MIN (y1, y2)) - 1;
+      area.height = ceil (MAX (y3, y4)) + 3 - area.y;
+    }
+
+  return area;
 }
 
 static void
-lens_get_source_coord (gdouble         i,
-                       gdouble         j,
-                       gdouble        *x,
-                       gdouble        *y,
-                       gdouble        *mag,
-                       LensDistortion *o)
+clamp_area (GeglRectangle *area,
+            gdouble        center_x,
+            gdouble        center_y)
 {
-  gdouble radius_sq, off_x, off_y, radius_mult;
+  if (center_x <= area->x || area->width < 1)
+    {
+      area->width = CLAMP (area->width, 1, MAX_WH);
+    }
+  else
+    {
+      area->x += area->width - CLAMP (area->width, 1, MAX_WH);
+      area->width = CLAMP (area->width, 1, MAX_WH);
+    }
 
-  off_x = i - o->centre_x;
-  off_y = j - o->centre_y;
-  radius_sq = (off_x * off_x) + (off_y * off_y);
+  if (center_y <= area->y || area->height < 1)
+    {
+      area->height = CLAMP (area->height, 1, MAX_WH);
+    }
+  else
+    {
+      area->y += area->height - CLAMP (area->height, 1, MAX_WH);
+      area->height = CLAMP (area->height, 1, MAX_WH);
+    }
+}
 
-  radius_sq *= o->norm;
+static GeglRectangle
+get_required_for_output (GeglOperation       *operation,
+                         const gchar         *input_pad,
+                         const GeglRectangle *roi)
+{
+  GeglRectangle *boundary;
 
-  radius_mult = radius_sq * o->mult_sq + radius_sq * radius_sq * o->mult_qd;
+  boundary = gegl_operation_source_get_bounding_box (operation, "input");
 
-  *mag = radius_mult;
-  radius_mult = o->rescale * (1.0 + radius_mult);
+  if (strcmp (input_pad, "input") || !boundary)
+    return *GEGL_RECTANGLE (0, 0, 0, 0);
 
-  *x = o->centre_x + radius_mult * off_x;
-  *y = o->centre_y + radius_mult * off_y;
+  return get_required (boundary, roi, operation);
+}
 
+static void
+prepare (GeglOperation *operation)
+{
+  gegl_operation_set_format (operation, "input", babl_format ("RGBA float"));
+  gegl_operation_set_format (operation, "output", babl_format ("RGBA float"));
 }
 
 /*
@@ -176,26 +298,24 @@ lens_cubic_interpolate (gfloat  *src,
   vp1 = ((-1.5 * dy + 2.0) * dy + 0.5) * dy;
   vp2 = (0.5 * dy - 0.5) * dy * dy;
 
-  /* Note: if dst_depth < src_depth, we calculate unneeded pixels here */
-  /* later - select or create index array */
   for (c = 0; c < 4 * 4; ++c)
     {
-      verts[c] = vm1 * src[c] + v * src[c+row_stride] +
-        vp1 * src[c+row_stride*2] + vp2 * src[c+row_stride*3];
+      verts[c] = vm1 * src[c] + v * src[c + row_stride] +
+        vp1 * src[c + row_stride * 2] + vp2 * src[ c + row_stride * 3];
     }
 
   for (c = 0; c < 4; ++c)
     {
       gfloat result;
 
-      result = um1 * verts[c] + u * verts[c+4] +
-        up1 * verts[c+4*2] + up2 * verts[c+4*3];
+      result = um1 * verts[c] + u * verts[c + 4] +
+        up1 * verts[c + 4 * 2] + up2 * verts[c + 4 * 3];
 
-      result *= brighten;
+      if (c != 3)
+        result *= brighten;
 
       dst[c] = CLAMP (result, 0.0, 1.0);
     }
-
 }
 
 static void
@@ -204,7 +324,7 @@ lens_distort_func (gfloat              *src_buf,
                    const GeglRectangle *extended,
                    const GeglRectangle *result,
                    const GeglRectangle *boundary,
-                   LensDistortion       old,
+                   LensValues          *lens,
                    gint                 xx,
                    gint                 yy,
                    GeglBuffer          *input)
@@ -217,9 +337,10 @@ lens_distort_func (gfloat              *src_buf,
 
   temp[0] = temp[1] = temp[2] = temp[3] = 0.0;
 
-  lens_get_source_coord ((gdouble)xx, (gdouble)yy, &sx, &sy, &mag, &old);
+  lens_get_source_coord ((gdouble) xx, (gdouble) yy, &sx, &sy, &mag, lens);
 
-  brighten = 1.0 + mag * old.brighten;
+  /* pseudo gamma transformation, since the input is scRGB */
+  brighten = pow (MAX (1.0 + mag * lens->brighten, 0.0), 2.4);
 
   x_int = floor (sx);
   dx = sx - x_int;
@@ -233,31 +354,31 @@ lens_distort_func (gfloat              *src_buf,
         {
           gint b;
 
-          if (x >= extended->x && x<(extended->x + extended->width)
+          if (x >= extended->x && x< (extended->x + extended->width)
               && y >= extended->y && y < (extended->y + extended->height))
             {
               gint src_off;
               src_off = (y - extended->y) * extended->width * 4 +
                 (x - extended->x) * 4;
-              for (b=0; b<4; b++)
-                temp[b] = src_buf[src_off++];
 
+              for (b = 0; b < 4; b++)
+                temp[b] = src_buf[src_off++];
             }
-          else if (x >= boundary->x && x < boundary->x + boundary->width &&
-                   y >= boundary->y && y < boundary->y + boundary->height)
+          else
             {
               gegl_buffer_sample (input, x, y, NULL, temp,
                                   babl_format ("RGBA float"),
-                                  GEGL_SAMPLER_CUBIC,
-                                  GEGL_ABYSS_NONE);
+                                  GEGL_SAMPLER_LINEAR,
+                                  GEGL_ABYSS_CLAMP);
             }
-          else
+
+          if (x < boundary->x || x >= (boundary->x + boundary->width) ||
+              y < boundary->y || y >= (boundary->y + boundary->height))
             {
-              for (b=0; b<4; b++)
-                temp[b] = 0.0;
+              temp[3] = 0.0;
             }
 
-          for (b=0; b<4; b++)
+          for (b = 0; b < 4; b++)
             pixel_buffer[offset++] = temp[b];
         }
     }
@@ -265,7 +386,8 @@ lens_distort_func (gfloat              *src_buf,
   lens_cubic_interpolate (pixel_buffer, temp, dx, dy, brighten);
 
   offset = (yy - result->y) * result->width * 4 + (xx - result->x) * 4;
-  for (x=0; x<4; x++)
+
+  for (x = 0; x < 4; x++)
     dst_buf[offset++] = temp[x];
 }
 
@@ -276,32 +398,51 @@ process (GeglOperation       *operation,
          const GeglRectangle *result,
          gint                 level)
 {
-  GeglChantO          *o = GEGL_CHANT_PROPERTIES (operation);
-  LensDistortion       old_lens;
-  GeglRectangle        boundary = *gegl_operation_source_get_bounding_box
-    (operation, "input");
+  GeglChantO    *o = GEGL_CHANT_PROPERTIES (operation);
+  LensValues     lens;
+  GeglRectangle  boundary;
+  gint           i, j;
+  gfloat        *src_buf, *dst_buf;
+
+  boundary = *gegl_operation_source_get_bounding_box (operation, "input");
+  lens     =  lens_setup_calc (o, boundary);
 
-  gint     x, y;
-  gfloat  *src_buf, *dst_buf;
+  src_buf = g_new0 (gfloat, SQR (MAX_WH) * 4);
+  dst_buf = g_new0 (gfloat, SQR (CHUNK_SIZE) * 4);
 
-  src_buf    = g_new0 (gfloat, result->width * result->height * 4);
-  dst_buf    = g_new0 (gfloat, result->width * result->height * 4);
+  for (j = 0; (j-1) * CHUNK_SIZE < result->height; j++)
+    for (i = 0; (i-1) * CHUNK_SIZE < result->width; i++)
+      {
+        GeglRectangle chunked_result;
+        GeglRectangle area;
+        gint          x, y;
 
-  lens_setup_calc (o, boundary, &old_lens);
+        chunked_result = *GEGL_RECTANGLE (result->x + i * CHUNK_SIZE,
+                                          result->y + j * CHUNK_SIZE,
+                                          CHUNK_SIZE, CHUNK_SIZE);
 
-  gegl_buffer_get (input, result, 1.0, babl_format ("RGBA float"), src_buf,
-                   GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
+        gegl_rectangle_intersect (&chunked_result, &chunked_result, result);
 
-  for (y = result->y; y < result->y + result->height; y++)
-    for (x = result->x; x < result->x + result->width; x++)
-      {
-        lens_distort_func (src_buf, dst_buf, result, result, &boundary,
-                           old_lens, x, y, input);
-      }
+        if (chunked_result.width < 1  || chunked_result.height < 1)
+          continue;
 
+        area = get_required (&boundary, &chunked_result, operation);
 
-  gegl_buffer_set (output, result, 0, babl_format ("RGBA float"),
-                   dst_buf, GEGL_AUTO_ROWSTRIDE);
+        clamp_area (&area, lens.centre_x, lens.centre_y);
+
+        gegl_buffer_get (input, &area, 1.0, babl_format ("RGBA float"), src_buf,
+                         GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_CLAMP);
+
+        for (y = chunked_result.y; y < chunked_result.y + chunked_result.height; y++)
+          for (x = chunked_result.x; x < chunked_result.x + chunked_result.width; x++)
+            {
+              lens_distort_func (src_buf, dst_buf, &area, &chunked_result, &boundary,
+                                 &lens, x, y, input);
+            }
+
+        gegl_buffer_set (output, &chunked_result, 0, babl_format ("RGBA float"),
+                         dst_buf, GEGL_AUTO_ROWSTRIDE);
+      }
 
   g_free (dst_buf);
   g_free (src_buf);
@@ -337,7 +478,6 @@ gegl_chant_class_init (GeglChantClass *klass)
   filter_class    = GEGL_OPERATION_FILTER_CLASS (klass);
 
   operation_class->prepare                 = prepare;
-  operation_class->get_bounding_box        = get_bounding_box;
   operation_class->get_required_for_output = get_required_for_output;
 
   filter_class->process                    = process;


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