[gegl] buffer: add box filtering downsampling paths to linear and cubic resamplers

[Øyvind Kolås <ok src gnome org>]

commit ce7383e705c9d8106400671d5c92faba134a2a6d
Author: Øyvind Kolås <pippin gimp org>
Date:   Sat May 7 13:33:22 2016 +0200

    buffer: add box filtering downsampling paths to linear and cubic resamplers
    
    Having resamplers that do interpolation for scaling down is just wrong,
    use some simple 2x2 samples and 4x4 samples box filters for linear and cubic samplers to improve default 
results.

 gegl/buffer/gegl-sampler-cubic.c  |  151 +++++++++++++------
 gegl/buffer/gegl-sampler-linear.c |  292 ++++++++++++++++++++++---------------
 2 files changed, 276 insertions(+), 167 deletions(-)
---
diff --git a/gegl/buffer/gegl-sampler-cubic.c b/gegl/buffer/gegl-sampler-cubic.c
index ad4f4a7..401e80e 100644
--- a/gegl/buffer/gegl-sampler-cubic.c
+++ b/gegl/buffer/gegl-sampler-cubic.c
@@ -152,6 +152,51 @@ gegl_sampler_cubic_init (GeglSamplerCubic *self)
     }
 }
 
+static void inline
+gegl_sampler_box_get (GeglSampler*    restrict  self,
+                      const gdouble             absolute_x,
+                      const gdouble             absolute_y,
+                      GeglMatrix2              *scale,
+                      void*           restrict  output,
+                      GeglAbyssPolicy           repeat_mode)
+{
+  gfloat result[4] = {0,0,0,0};
+  const float iabsolute_x = (float) absolute_x - 0.5;
+  const float iabsolute_y = (float) absolute_y - 0.5;
+  const gint ix = floorf (iabsolute_x - scale->coeff[0][0]/2);
+  const gint iy = floorf (iabsolute_y - scale->coeff[1][1]/2);
+  const gint xx = ceilf (iabsolute_x + scale->coeff[0][0]/2);
+  const gint yy = ceilf (iabsolute_y + scale->coeff[1][1]/2);
+  int u, v;
+  int count = 0;
+  int hskip = scale->coeff[0][0] / 5;
+  int vskip = scale->coeff[1][1] / 5;
+
+  if (hskip <= 0)
+    hskip = 1;
+  if (vskip <= 0)
+    vskip = 1;
+
+  for (v = iy; v < yy; v += vskip)
+  {
+    for (u = ix; u < xx; u += hskip)
+    {
+      int c;
+      gfloat input[4];
+      GeglRectangle rect = {u, v, 1, 1};
+      gegl_buffer_get (self->buffer, &rect, 1.0, self->interpolate_format, input, GEGL_AUTO_ROWSTRIDE, 
repeat_mode);
+      for (c = 0; c < 4; c++)
+        result[c] += input[c];
+      count ++;
+    }
+  }
+  result[0] /= count;
+  result[1] /= count;
+  result[2] /= count;
+  result[3] /= count;
+  babl_process (self->fish, result, output, 1);
+}
+
 void
 gegl_sampler_cubic_get (      GeglSampler     *self,
                         const gdouble          absolute_x,
@@ -160,62 +205,72 @@ gegl_sampler_cubic_get (      GeglSampler     *self,
                               void            *output,
                               GeglAbyssPolicy  repeat_mode)
 {
-  GeglSamplerCubic *cubic       = (GeglSamplerCubic*)(self);
-  const gint        offsets[16] = {
-                                    -4-GEGL_SAMPLER_MAXIMUM_WIDTH   *4, 4, 4, 4,
-                                      (GEGL_SAMPLER_MAXIMUM_WIDTH-3)*4, 4, 4, 4,
-                                      (GEGL_SAMPLER_MAXIMUM_WIDTH-3)*4, 4, 4, 4,
-                                      (GEGL_SAMPLER_MAXIMUM_WIDTH-3)*4, 4, 4, 4
-                                  };
-  gfloat           *sampler_bptr;
-  gfloat            factor;
-  gfloat            newval[4]   = {0, 0, 0, 0};
-  gint              i,
-                    j,
-                    k           = 0;
-
-  /*
-   * The "-1/2"s are there because we want the index of the pixel
-   * center to the left and top of the location, and with GIMP's
-   * convention the top left of the top left pixel is located at
-   * (0,0), and its center is at (1/2,1/2), so that anything less than
-   * 1/2 needs to go negative. Another way to look at this is that we
-   * are converting from a coordinate system in which the origin is at
-   * the top left corner of the pixel with index (0,0), to a
-   * coordinate system in which the origin is at the center of the
-   * same pixel.
-   */
-  const double iabsolute_x = (double) absolute_x - 0.5;
-  const double iabsolute_y = (double) absolute_y - 0.5;
-
-  const gint ix = floorf (iabsolute_x);
-  const gint iy = floorf (iabsolute_y);
-
-  /*
-   * x is the x-coordinate of the sampling point relative to the
-   * position of the center of the top left pixel. Similarly for
-   * y. Range of values: [0,1].
-   */
-  const gfloat x = iabsolute_x - ix;
-  const gfloat y = iabsolute_y - iy;
-
-  sampler_bptr = gegl_sampler_get_ptr (self, ix, iy, repeat_mode);
-
-  for (j=-1; j<3; j++)
-    for (i=-1; i<3; i++)
-      {
-        sampler_bptr += offsets[k++];
-
-        factor = cubicKernel (y - j, cubic->b, cubic->c) *
-                 cubicKernel (x - i, cubic->b, cubic->c);
-
-        newval[0] += factor * sampler_bptr[0];
-        newval[1] += factor * sampler_bptr[1];
-        newval[2] += factor * sampler_bptr[2];
-        newval[3] += factor * sampler_bptr[3];
-      }
-
-  babl_process (self->fish, newval, output, 1);
+  if (scale && 
+      (scale->coeff[0][0] * scale->coeff[0][0] +
+      scale->coeff[1][1] * scale->coeff[1][1])
+    > 2.5)
+  {
+    gegl_sampler_box_get (self, absolute_x, absolute_y, scale, output, repeat_mode);
+  }
+  else
+  {
+    GeglSamplerCubic *cubic       = (GeglSamplerCubic*)(self);
+    const gint        offsets[16] = {
+                                      -4-GEGL_SAMPLER_MAXIMUM_WIDTH   *4, 4, 4, 4,
+                                        (GEGL_SAMPLER_MAXIMUM_WIDTH-3)*4, 4, 4, 4,
+                                        (GEGL_SAMPLER_MAXIMUM_WIDTH-3)*4, 4, 4, 4,
+                                        (GEGL_SAMPLER_MAXIMUM_WIDTH-3)*4, 4, 4, 4
+                                    };
+    gfloat           *sampler_bptr;
+    gfloat            factor;
+    gfloat            newval[4]   = {0, 0, 0, 0};
+    gint              i,
+                      j,
+                      k           = 0;
+
+    /*
+     * The "-1/2"s are there because we want the index of the pixel
+     * center to the left and top of the location, and with GIMP's
+     * convention the top left of the top left pixel is located at
+     * (0,0), and its center is at (1/2,1/2), so that anything less than
+     * 1/2 needs to go negative. Another way to look at this is that we
+     * are converting from a coordinate system in which the origin is at
+     * the top left corner of the pixel with index (0,0), to a
+     * coordinate system in which the origin is at the center of the
+     * same pixel.
+     */
+    const double iabsolute_x = (double) absolute_x - 0.5;
+    const double iabsolute_y = (double) absolute_y - 0.5;
+
+    const gint ix = floorf (iabsolute_x);
+    const gint iy = floorf (iabsolute_y);
+
+    /*
+     * x is the x-coordinate of the sampling point relative to the
+     * position of the center of the top left pixel. Similarly for
+     * y. Range of values: [0,1].
+     */
+    const gfloat x = iabsolute_x - ix;
+    const gfloat y = iabsolute_y - iy;
+
+    sampler_bptr = gegl_sampler_get_ptr (self, ix, iy, repeat_mode);
+
+    for (j=-1; j<3; j++)
+      for (i=-1; i<3; i++)
+        {
+          sampler_bptr += offsets[k++];
+
+          factor = cubicKernel (y - j, cubic->b, cubic->c) *
+                   cubicKernel (x - i, cubic->b, cubic->c);
+
+          newval[0] += factor * sampler_bptr[0];
+          newval[1] += factor * sampler_bptr[1];
+          newval[2] += factor * sampler_bptr[2];
+          newval[3] += factor * sampler_bptr[3];
+        }
+
+    babl_process (self->fish, newval, output, 1);
+  }
 }
 
 static void
diff --git a/gegl/buffer/gegl-sampler-linear.c b/gegl/buffer/gegl-sampler-linear.c
index add9f04..47a54d1 100644
--- a/gegl/buffer/gegl-sampler-linear.c
+++ b/gegl/buffer/gegl-sampler-linear.c
@@ -73,130 +73,184 @@ gegl_sampler_linear_init (GeglSamplerLinear *self)
   GEGL_SAMPLER (self)->interpolate_format = babl_format ("RaGaBaA float");
 }
 
-static void
-gegl_sampler_linear_get (      GeglSampler*    restrict  self,
-                         const gdouble                   absolute_x,
-                         const gdouble                   absolute_y,
-                               GeglMatrix2              *scale,
-                               void*           restrict  output,
-                               GeglAbyssPolicy           repeat_mode)
+static inline void
+gegl_sampler_box_get (GeglSampler*    restrict  self,
+                      const gdouble             absolute_x,
+                      const gdouble             absolute_y,
+                      GeglMatrix2              *scale,
+                      void*           restrict  output,
+                      GeglAbyssPolicy           repeat_mode)
 {
-  const gint pixels_per_buffer_row = GEGL_SAMPLER_MAXIMUM_WIDTH;
-  const gint channels = 4;
-
-  /*
-   * The "-1/2"s are there because we want the index of the pixel to
-   * the left and top of the location, and with GIMP's convention the
-   * top left of the top left pixel is located at
-   * (1/2,1/2). Basically, we are converting from a coordinate system
-   * in which the origin is at the top left pixel of the pixel with
-   * index (0,0), to a coordinate system in which the origin is at the
-   * center of the same pixel.
-   */
+  gfloat result[4] = {0,0,0,0};
   const float iabsolute_x = (float) absolute_x - 0.5;
   const float iabsolute_y = (float) absolute_y - 0.5;
+  const gint ix = floorf (iabsolute_x - scale->coeff[0][0]/2);
+  const gint iy = floorf (iabsolute_y - scale->coeff[1][1]/2);
+  const gint xx = ceilf (iabsolute_x + scale->coeff[0][0]/2);
+  const gint yy = ceilf (iabsolute_y + scale->coeff[1][1]/2);
+  int u, v;
+  int count = 0;
+  int hskip = scale->coeff[0][0] / 3;
+  int vskip = scale->coeff[1][1] / 3;
+
+  if (hskip <= 0)
+    hskip = 1;
+  if (vskip <= 0)
+    vskip = 1;
 
-  const gint ix = floorf (iabsolute_x);
-  const gint iy = floorf (iabsolute_y);
-
-  /*
-   * Point the data tile pointer to the first channel of the top_left
-   * pixel value:
-   */
-  const gfloat* restrict in_bptr =
-    gegl_sampler_get_ptr (self, ix, iy, repeat_mode);
-
-  /*
-   * x is the x-coordinate of the sampling point relative to the
-   * position of the center of the top left pixel. Similarly for
-   * y. Range of values: [0,1].
-   */
-  const gfloat x = iabsolute_x - ix;
-  const gfloat y = iabsolute_y - iy;
-
-  /*
-   * First bilinear weight:
-   */
-  const gfloat x_times_y = x * y;
-
-  /*
-   * Load top row:
-   */
-  const gfloat top_left_0 = *in_bptr++;
-  const gfloat top_left_1 = *in_bptr++;
-  const gfloat top_left_2 = *in_bptr++;
-  const gfloat top_left_3 = *in_bptr++;
-  const gfloat top_rite_0 = *in_bptr++;
-  const gfloat top_rite_1 = *in_bptr++;
-  const gfloat top_rite_2 = *in_bptr++;
-  const gfloat top_rite_3 = *in_bptr;
-
-  in_bptr += 1 + ( pixels_per_buffer_row - 2 ) * channels;
+  for (v = iy; v < yy; v += vskip)
+  {
+    for (u = ix; u < xx; u += hskip)
+    {
+      int c;
+      gfloat input[4];
+      GeglRectangle rect = {u, v, 1, 1};
+      gegl_buffer_get (self->buffer, &rect, 1.0, self->interpolate_format, input, GEGL_AUTO_ROWSTRIDE, 
repeat_mode);
+      for (c = 0; c < 4; c++)
+        result[c] += input[c];
+      count ++;
+    }
+  }
+  result[0] /= count;
+  result[1] /= count;
+  result[2] /= count;
+  result[3] /= count;
+  babl_process (self->fish, result, output, 1);
+}
 
+void
+gegl_sampler_linear_get (     GeglSampler     *self,
+                        const gdouble          absolute_x,
+                        const gdouble          absolute_y,
+                              GeglMatrix2     *scale,
+                              void            *output,
+                              GeglAbyssPolicy  repeat_mode)
+{
+  if (scale && (scale->coeff[0][0] * scale->coeff[0][0] +
+      scale->coeff[1][1] * scale->coeff[1][1])
+    > 2.5)
   {
-  /*
-   * More bilinear weights:
-   *
-   * (Note: w = 1-x and z = 1-y.)
-   */
-  const gfloat w_times_y = y - x_times_y;
-  const gfloat x_times_z = x - x_times_y;
-
-  /*
-   * Load bottom row:
-   */
-  const gfloat bot_left_0 = *in_bptr++;
-  const gfloat bot_left_1 = *in_bptr++;
-  const gfloat bot_left_2 = *in_bptr++;
-  const gfloat bot_left_3 = *in_bptr++;
-  const gfloat bot_rite_0 = *in_bptr++;
-  const gfloat bot_rite_1 = *in_bptr++;
-  const gfloat bot_rite_2 = *in_bptr++;
-  const gfloat bot_rite_3 = *in_bptr;
-
-  /*
-   * Last bilinear weight:
-   */
-  const gfloat w_times_z = (gfloat) 1. - ( x + w_times_y );
-
-  gfloat newval[4];
-
-  newval[0] =
-    x_times_y * bot_rite_0
-    +
-    w_times_y * bot_left_0
-    +
-    x_times_z * top_rite_0
-    +
-    w_times_z * top_left_0;
-
-  newval[1] =
-    x_times_y * bot_rite_1
-    +
-    w_times_y * bot_left_1
-    +
-    x_times_z * top_rite_1
-    +
-    w_times_z * top_left_1;
-
-  newval[2] =
-    x_times_y * bot_rite_2
-    +
-    w_times_y * bot_left_2
-    +
-    x_times_z * top_rite_2
-    +
-    w_times_z * top_left_2;
-
-  newval[3] =
-    x_times_y * bot_rite_3
-    +
-    w_times_y * bot_left_3
-    +
-    x_times_z * top_rite_3
-    +
-    w_times_z * top_left_3;
-
-    babl_process (self->fish, newval, output, 1);
+    gegl_sampler_box_get (self, absolute_x, absolute_y, scale, output, repeat_mode);
+  }
+  else
+  {
+    const gint pixels_per_buffer_row = GEGL_SAMPLER_MAXIMUM_WIDTH;
+    const gint channels = 4;
+
+    /*
+     * The "-1/2"s are there because we want the index of the pixel to
+     * the left and top of the location, and with GIMP's convention the
+     * top left of the top left pixel is located at
+     * (1/2,1/2). Basically, we are converting from a coordinate system
+     * in which the origin is at the top left pixel of the pixel with
+     * index (0,0), to a coordinate system in which the origin is at the
+     * center of the same pixel.
+     */
+    const float iabsolute_x = (float) absolute_x - 0.5;
+    const float iabsolute_y = (float) absolute_y - 0.5;
+
+    const gint ix = floorf (iabsolute_x);
+    const gint iy = floorf (iabsolute_y);
+
+    /*
+     * Point the data tile pointer to the first channel of the top_left
+     * pixel value:
+     */
+    const gfloat* restrict in_bptr =
+      gegl_sampler_get_ptr (self, ix, iy, repeat_mode);
+
+    /*
+     * x is the x-coordinate of the sampling point relative to the
+     * position of the center of the top left pixel. Similarly for
+     * y. Range of values: [0,1].
+     */
+    const gfloat x = iabsolute_x - ix;
+    const gfloat y = iabsolute_y - iy;
+
+    /*
+     * First bilinear weight:
+     */
+    const gfloat x_times_y = x * y;
+
+    /*
+     * Load top row:
+     */
+    const gfloat top_left_0 = *in_bptr++;
+    const gfloat top_left_1 = *in_bptr++;
+    const gfloat top_left_2 = *in_bptr++;
+    const gfloat top_left_3 = *in_bptr++;
+    const gfloat top_rite_0 = *in_bptr++;
+    const gfloat top_rite_1 = *in_bptr++;
+    const gfloat top_rite_2 = *in_bptr++;
+    const gfloat top_rite_3 = *in_bptr;
+
+    in_bptr += 1 + ( pixels_per_buffer_row - 2 ) * channels;
+
+    {
+    /*
+     * More bilinear weights:
+     *
+     * (Note: w = 1-x and z = 1-y.)
+     */
+    const gfloat w_times_y = y - x_times_y;
+    const gfloat x_times_z = x - x_times_y;
+
+    /*
+     * Load bottom row:
+     */
+    const gfloat bot_left_0 = *in_bptr++;
+    const gfloat bot_left_1 = *in_bptr++;
+    const gfloat bot_left_2 = *in_bptr++;
+    const gfloat bot_left_3 = *in_bptr++;
+    const gfloat bot_rite_0 = *in_bptr++;
+    const gfloat bot_rite_1 = *in_bptr++;
+    const gfloat bot_rite_2 = *in_bptr++;
+    const gfloat bot_rite_3 = *in_bptr;
+
+    /*
+     * Last bilinear weight:
+     */
+    const gfloat w_times_z = (gfloat) 1. - ( x + w_times_y );
+
+    gfloat newval[4];
+
+    newval[0] =
+      x_times_y * bot_rite_0
+      +
+      w_times_y * bot_left_0
+      +
+      x_times_z * top_rite_0
+      +
+      w_times_z * top_left_0;
+
+    newval[1] =
+      x_times_y * bot_rite_1
+      +
+      w_times_y * bot_left_1
+      +
+      x_times_z * top_rite_1
+      +
+      w_times_z * top_left_1;
+
+    newval[2] =
+      x_times_y * bot_rite_2
+      +
+      w_times_y * bot_left_2
+      +
+      x_times_z * top_rite_2
+      +
+      w_times_z * top_left_2;
+
+    newval[3] =
+      x_times_y * bot_rite_3
+      +
+      w_times_y * bot_left_3
+      +
+      x_times_z * top_rite_3
+      +
+      w_times_z * top_left_3;
+
+      babl_process (self->fish, newval, output, 1);
+    }
   }
 }