workshop: add the FIR case to the work-in-progress gaussian-blur

Things missing are:

- OpenCL handles only 4-components babl formats,
- The IIR case doesn't make the extent grow

4 files changed, 676 insertions, 61 deletions

diff --git a/opencl/gblur-1d.cl b/opencl/gblur-1d.cl
new file mode 100644
index 00000000..834cc858
--- /dev/null
+++ b/opencl/gblur-1d.cl
@@ -0,0 +1,73 @@
+/* This file is an image processing operation for GEGL
+ *
+ * GEGL is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 3 of the License, or (at your option) any later version.
+ *
+ * GEGL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GEGL; if not, see <http://www.gnu.org/licenses/>.
+ *
+ * Copyright 2013 Victor Oliveira <victormatheus@gmail.com>
+ * Copyright 2013 Téo Mazars      <teomazars@gmail.com>
+ */
+
+__kernel void fir_ver_blur(const global float4 *src_buf,
+                                 global float4 *dst_buf,
+                           const global float  *cmatrix,
+                           const        int     clen)
+{
+    const int gidx          = get_global_id (0);
+    const int gidy          = get_global_id (1);
+    const int src_rowstride = get_global_size (0);
+    const int dst_rowstride = get_global_size (0);
+
+    const int half_clen     = clen / 2;
+
+    const int src_offset    = gidx + (gidy + half_clen) * src_rowstride;
+    const int dst_offset    = gidx +  gidy              * dst_rowstride;
+
+    const int src_start_ind = src_offset - half_clen * src_rowstride;
+
+    float4 v = 0.0f;
+
+    for (int i = 0; i < clen; i++)
+      {
+        v += src_buf[src_start_ind + i * src_rowstride] * cmatrix[i];
+      }
+
+    dst_buf[dst_offset] = v;
+}
+
+
+__kernel void fir_hor_blur(const global float4 *src_buf,
+                                 global float4 *dst_buf,
+                           const global float  *cmatrix,
+                           const        int     clen)
+{
+    const int gidx          = get_global_id (0);
+    const int gidy          = get_global_id (1);
+    const int src_rowstride = get_global_size (0) + clen - 1; /*== 2*(clen/2) */
+    const int dst_rowstride = get_global_size (0);
+
+    const int half_clen     = clen / 2;
+
+    const int src_offset    = gidx + gidy * src_rowstride + half_clen;
+    const int dst_offset    = gidx + gidy * dst_rowstride;
+
+    const int src_start_ind = src_offset - half_clen;
+
+    float4 v = 0.0f;
+
+    for (int i = 0; i < clen; i++)
+      {
+        v += src_buf[src_start_ind + i] * cmatrix[i];
+      }
+
+    dst_buf[dst_offset] = v;
+}
diff --git a/opencl/gblur-1d.cl.h b/opencl/gblur-1d.cl.h
new file mode 100644
index 00000000..a3391f9d
--- /dev/null
+++ b/opencl/gblur-1d.cl.h
@@ -0,0 +1,75 @@
+static const char* gblur_1d_cl_source =
+"/* This file is an image processing operation for GEGL                        \n"
+" *                                                                            \n"
+" * GEGL is free software; you can redistribute it and/or                      \n"
+" * modify it under the terms of the GNU Lesser General Public                 \n"
+" * License as published by the Free Software Foundation; either               \n"
+" * version 3 of the License, or (at your option) any later version.           \n"
+" *                                                                            \n"
+" * GEGL is distributed in the hope that it will be useful,                    \n"
+" * but WITHOUT ANY WARRANTY; without even the implied warranty of             \n"
+" * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU          \n"
+" * Lesser General Public License for more details.                            \n"
+" *                                                                            \n"
+" * You should have received a copy of the GNU Lesser General Public           \n"
+" * License along with GEGL; if not, see <http://www.gnu.org/licenses/>.       \n"
+" *                                                                            \n"
+" * Copyright 2013 Victor Oliveira <victormatheus@gmail.com>                   \n"
+" * Copyright 2013 T\303\251o Mazars      <teomazars@gmail.com>                \n"
+" */                                                                           \n"
+"                                                                              \n"
+"__kernel void fir_ver_blur(const global float4 *src_buf,                      \n"
+"                                 global float4 *dst_buf,                      \n"
+"                           const global float  *cmatrix,                      \n"
+"                           const        int     clen)                         \n"
+"{                                                                             \n"
+"    const int gidx          = get_global_id (0);                              \n"
+"    const int gidy          = get_global_id (1);                              \n"
+"    const int src_rowstride = get_global_size (0);                            \n"
+"    const int dst_rowstride = get_global_size (0);                            \n"
+"                                                                              \n"
+"    const int half_clen     = clen / 2;                                       \n"
+"                                                                              \n"
+"    const int src_offset    = gidx + (gidy + half_clen) * src_rowstride;      \n"
+"    const int dst_offset    = gidx +  gidy              * dst_rowstride;      \n"
+"                                                                              \n"
+"    const int src_start_ind = src_offset - half_clen * src_rowstride;         \n"
+"                                                                              \n"
+"    float4 v = 0.0f;                                                          \n"
+"                                                                              \n"
+"    for (int i = 0; i < clen; i++)                                            \n"
+"      {                                                                       \n"
+"        v += src_buf[src_start_ind + i * src_rowstride] * cmatrix[i];         \n"
+"      }                                                                       \n"
+"                                                                              \n"
+"    dst_buf[dst_offset] = v;                                                  \n"
+"}                                                                             \n"
+"                                                                              \n"
+"                                                                              \n"
+"__kernel void fir_hor_blur(const global float4 *src_buf,                      \n"
+"                                 global float4 *dst_buf,                      \n"
+"                           const global float  *cmatrix,                      \n"
+"                           const        int     clen)                         \n"
+"{                                                                             \n"
+"    const int gidx          = get_global_id (0);                              \n"
+"    const int gidy          = get_global_id (1);                              \n"
+"    const int src_rowstride = get_global_size (0) + clen - 1; /*== 2*(clen/2) */\n"
+"    const int dst_rowstride = get_global_size (0);                            \n"
+"                                                                              \n"
+"    const int half_clen     = clen / 2;                                       \n"
+"                                                                              \n"
+"    const int src_offset    = gidx + gidy * src_rowstride + half_clen;        \n"
+"    const int dst_offset    = gidx + gidy * dst_rowstride;                    \n"
+"                                                                              \n"
+"    const int src_start_ind = src_offset - half_clen;                         \n"
+"                                                                              \n"
+"    float4 v = 0.0f;                                                          \n"
+"                                                                              \n"
+"    for (int i = 0; i < clen; i++)                                            \n"
+"      {                                                                       \n"
+"        v += src_buf[src_start_ind + i] * cmatrix[i];                         \n"
+"      }                                                                       \n"
+"                                                                              \n"
+"    dst_buf[dst_offset] = v;                                                  \n"
+"}                                                                             \n"
+;
diff --git a/operations/workshop/gaussian-blur-iir.c b/operations/workshop/gaussian-blur-iir.c
index 52cf36e4..016335a6 100644
--- a/operations/workshop/gaussian-blur-iir.c
+++ b/operations/workshop/gaussian-blur-iir.c
@@ -21,14 +21,33 @@
 
 #ifdef GEGL_CHANT_PROPERTIES
 
-gegl_chant_enum      (abyss_policy, _("Abyss policy"), GeglAbyssPolicy,
-                      gegl_abyss_policy, GEGL_ABYSS_CLAMP, _(""))
+gegl_chant_register_enum (gegl_gaussian_blur_filter)
+  enum_value (GEGL_GAUSSIAN_BLUR_FILTER_AUTO, "Auto")
+  enum_value (GEGL_GAUSSIAN_BLUR_FILTER_FIR,  "FIR")
+  enum_value (GEGL_GAUSSIAN_BLUR_FILTER_IIR,  "IIR")
+gegl_chant_register_enum_end (GeglGaussianBlurFilter)
+
+gegl_chant_register_enum (gegl_gaussian_blur_policy)
+   enum_value (GEGL_GAUSSIAN_BLUR_ABYSS_NONE,  "None")
+   enum_value (GEGL_GAUSSIAN_BLUR_ABYSS_CLAMP, "Clamp")
+   enum_value (GEGL_GAUSSIAN_BLUR_ABYSS_BLACK, "Black")
+   enum_value (GEGL_GAUSSIAN_BLUR_ABYSS_WHITE, "White")
+gegl_chant_register_enum_end (GeglGaussianBlurPolicy)
+
+
 gegl_chant_double_ui (std_dev_x, _("Horizontal Std. Dev."),
-                      0.5, 1500.0, 1.5, 0.5, 100.0, 3.0,
+                      0.0, 1500.0, 1.5, 0.0, 100.0, 3.0,
                       _("Standard deviation (spatial scale factor)"))
 gegl_chant_double_ui (std_dev_y, _("Vertical Std. Dev."),
-                      0.5, 1500.0, 1.5, 0.5, 100.0, 3.0,
+                      0.0, 1500.0, 1.5, 0.0, 100.0, 3.0,
                       _("Standard deviation (spatial scale factor)"))
+gegl_chant_enum      (filter, _("Filter"),
+                      GeglGaussianBlurFilter, gegl_gaussian_blur_filter,
+                      GEGL_GAUSSIAN_BLUR_FILTER_AUTO,
+                      _("How the gaussian kernel is discretized"))
+gegl_chant_enum      (abyss_policy, _("Abyss policy"), GeglGaussianBlurPolicy,
+                      gegl_gaussian_blur_policy, GEGL_GAUSSIAN_BLUR_ABYSS_NONE,
+                      _("How image edges are handled"))
 
 #else
 
@@ -40,17 +59,30 @@ gegl_chant_double_ui (std_dev_y, _("Vertical Std. Dev."),
 static void
 attach (GeglOperation *operation)
 {
-  GeglNode *gegl = operation->node;
+  GeglNode *gegl   = operation->node;
   GeglNode *output = gegl_node_get_output_proxy (gegl, "output");
-  GeglNode *vblur = gegl_node_new_child (gegl, "operation", "gegl:gblur-1d", "orientation", 1, "abyss-policy", 0, NULL);
-  GeglNode *hblur = gegl_node_new_child (gegl, "operation", "gegl:gblur-1d", "abyss-policy", 0, NULL);
-  GeglNode *input = gegl_node_get_input_proxy (gegl, "input");
+
+  GeglNode *vblur  = gegl_node_new_child (gegl,
+                                          "operation", "gegl:gblur-1d",
+                                          "orientation", 1,
+                                          NULL);
+
+  GeglNode *hblur  = gegl_node_new_child (gegl,
+                                          "operation", "gegl:gblur-1d",
+                                          "orientation", 0,
+                                          NULL);
+
+  GeglNode *input  = gegl_node_get_input_proxy (gegl, "input");
 
   gegl_node_link_many (input, hblur, vblur, output, NULL);
-  gegl_operation_meta_redirect (operation, "std-dev-x" , hblur, "std-dev");
+
+  gegl_operation_meta_redirect (operation, "std-dev-x",    hblur, "std-dev");
   gegl_operation_meta_redirect (operation, "abyss-policy", hblur, "abyss-policy");
-  gegl_operation_meta_redirect (operation, "std-dev-y" , vblur, "std-dev");
+  gegl_operation_meta_redirect (operation, "filter",       hblur, "filter");
+
+  gegl_operation_meta_redirect (operation, "std-dev-y",    vblur, "std-dev");
   gegl_operation_meta_redirect (operation, "abyss-policy", vblur, "abyss-policy");
+  gegl_operation_meta_redirect (operation, "filter",       vblur, "filter");
 }
 
 static void
diff --git a/operations/workshop/gblur-1d.c b/operations/workshop/gblur-1d.c
index 92db86fb..77d0801a 100644
--- a/operations/workshop/gblur-1d.c
+++ b/operations/workshop/gblur-1d.c
@@ -25,32 +25,40 @@
 #include "config.h"
 #include <glib/gi18n-lib.h>
 
-
 #ifdef GEGL_CHANT_PROPERTIES
-#if 0
+
 gegl_chant_register_enum (gegl_gblur_1d_policy)
-   enum_value (GEGL_BLUR_1D_ABYSS_NONE, "None")
-   enum_value (GEGL_BLUR_1D_ABYSS_CLAMP, "Clamp")
-//   enum_value (GEGL_BLUR_1D_ABYSS_LOOP,
-   enum_value (GEGL_BLUR_1D_ABYSS_BLACK, "Black")
-   enum_value (GEGL_BLUR_1D_ABYSS_WHITE, "White")
+   enum_value (GEGL_GBLUR_1D_ABYSS_NONE,  "None")
+   enum_value (GEGL_GBLUR_1D_ABYSS_CLAMP, "Clamp")
+   enum_value (GEGL_GBLUR_1D_ABYSS_BLACK, "Black")
+   enum_value (GEGL_GBLUR_1D_ABYSS_WHITE, "White")
 gegl_chant_register_enum_end (GeglGblur1dPolicy)
-#endif
+
 gegl_chant_register_enum (gegl_gblur_1d_orientation)
-  enum_value (GEGL_BLUR_1D_HORIZONTAL, "Horizontal")
-  enum_value (GEGL_BLUR_1D_VERTICAL,   "Vertical")
-gegl_chant_register_enum_end (GeglGblur1dOrientation)
+  enum_value (GEGL_GBLUR_1D_HORIZONTAL, "Horizontal")
+  enum_value (GEGL_GBLUR_1D_VERTICAL,   "Vertical")
+  gegl_chant_register_enum_end (GeglGblur1dOrientation)
+
+gegl_chant_register_enum (gegl_gblur_1d_filter)
+  enum_value (GEGL_GBLUR_1D_AUTO, "Auto")
+  enum_value (GEGL_GBLUR_1D_FIR,  "FIR")
+  enum_value (GEGL_GBLUR_1D_IIR,  "IIR")
+gegl_chant_register_enum_end (GeglGblur1dFilter)
 
 gegl_chant_double_ui (std_dev, _("Size"),
-                      0.5, 1500.0, 1.5, 0.5, 100.0, 3.0,
-                      _("Standard deviation "
-                        "(multiply by ~2 to get radius)"))
+                      0.0, 1500.0, 1.5, 0.0, 100.0, 3.0,
+                      _("Standard deviation (spatial scale factor)"))
 gegl_chant_enum      (orientation, _("Orientation"),
                       GeglGblur1dOrientation, gegl_gblur_1d_orientation,
-                      GEGL_BLUR_1D_HORIZONTAL,
+                      GEGL_GBLUR_1D_HORIZONTAL,
                       _("The orientation of the blur - hor/ver"))
-gegl_chant_enum      (abyss_policy, _("Abyss policy"), GeglAbyssPolicy,
-                      gegl_abyss_policy, GEGL_ABYSS_CLAMP, _(""))
+gegl_chant_enum      (filter, _("Filter"),
+                      GeglGblur1dFilter, gegl_gblur_1d_filter,
+                      GEGL_GBLUR_1D_AUTO,
+                      _("How the gaussian kernel is discretized"))
+gegl_chant_enum      (abyss_policy, _("Abyss policy"), GeglGblur1dPolicy,
+                      gegl_gblur_1d_policy, GEGL_GBLUR_1D_ABYSS_NONE,
+                      _("How image edges are handled"))
 #else
 
 #define GEGL_CHANT_TYPE_FILTER
@@ -60,6 +68,13 @@ gegl_chant_enum      (abyss_policy, _("Abyss policy"), GeglAbyssPolicy,
 #include <math.h>
 #include <stdio.h>
 
+
+/**********************************************
+ *
+ * Infinite Impulse Response (IIR)
+ *
+ **********************************************/
+
 static void
 iir_young_find_constants (gfloat   sigma,
                           gdouble *b,
@@ -74,10 +89,10 @@ iir_young_find_constants (gfloat   sigma,
   const gdouble b2 =         - q*           q*(1.4281  + q*1.26661);
   const gdouble b3 =           q*           q*           q*0.422205;
 
-  const double a1 = b1 / b0;
-  const double a2 = b2 / b0;
-  const double a3 = b3 / b0;
-  const gdouble c = 1. / ((1+a1-a2+a3) * (1+a2+(a1-a3)*a3));
+  const gdouble a1 = b1 / b0;
+  const gdouble a2 = b2 / b0;
+  const gdouble a3 = b3 / b0;
+  const gdouble c  = 1. / ((1+a1-a2+a3) * (1+a2+(a1-a3)*a3));
 
   m[0][0] = c * (-a3*(a1+a3)-a2 + 1);
   m[0][1] = c * (a3+a1)*(a2+a3*a1);
@@ -127,7 +142,7 @@ iir_young_blur_1D (gfloat        *buf,
 {
   gfloat  white[4] = { 1, 1, 1, 1 };
   gfloat  black[4] = { 0, 0, 0, 1 };
-  gfloat  none[4]  = { 0, };
+  gfloat  none[4]  = { 0, 0, 0, 0 };
   gfloat *uplus, *iminus;
   gint    i, j, c;
 
@@ -241,6 +256,340 @@ iir_young_ver_blur (GeglBuffer          *src,
   g_free (col);
 }
 
+
+/**********************************************
+ *
+ * Finite Impulse Response (FIR)
+ *
+ **********************************************/
+
+static inline void
+fir_blur_1D (const gfloat *input,
+                   gfloat *output,
+             const gfloat *cmatrix,
+             const gint    clen,
+             const gint    len,
+             const gint    nc)
+{
+  gint i;
+  for (i = 0; i < len; i++)
+    {
+      gint c;
+      for (c = 0; c < nc; c++)
+        {
+          gint   index = i * nc + c;
+          gfloat acc   = 0.0f;
+          gint   m;
+
+          for (m = 0; m < clen; m++)
+            acc += input [index + m * nc] * cmatrix [m];
+
+          output [index] = acc;
+        }
+    }
+}
+
+static void
+fir_hor_blur (GeglBuffer          *src,
+              const GeglRectangle *rect,
+              GeglBuffer          *dst,
+              gfloat              *cmatrix,
+              gint                 clen,
+              GeglAbyssPolicy      policy,
+              const Babl          *format)
+{
+  GeglRectangle  cur_row = *rect;
+  GeglRectangle  in_row;
+  const gint     nc = babl_format_get_n_components (format);
+  gfloat        *row;
+  gfloat        *out;
+  gint           v;
+
+  cur_row.height = 1;
+
+  in_row         = cur_row;
+  in_row.width  += clen - 1;
+  in_row.x      -= clen / 2;
+
+  row = gegl_malloc (sizeof (gfloat) * in_row.width  * nc);
+  out = gegl_malloc (sizeof (gfloat) * cur_row.width * nc);
+
+  for (v = 0; v < rect->height; v++)
+    {
+      cur_row.y = in_row.y = rect->y + v;
+
+      gegl_buffer_get (src, &in_row, 1.0, format, row, GEGL_AUTO_ROWSTRIDE, policy);
+
+      fir_blur_1D (row, out, cmatrix, clen, rect->width, nc);
+
+      gegl_buffer_set (dst, &cur_row, 0, format, out, GEGL_AUTO_ROWSTRIDE);
+    }
+
+  gegl_free (out);
+  gegl_free (row);
+}
+
+static void
+fir_ver_blur (GeglBuffer          *src,
+              const GeglRectangle *rect,
+              GeglBuffer          *dst,
+              gfloat              *cmatrix,
+              gint                 clen,
+              GeglAbyssPolicy      policy,
+              const Babl          *format)
+{
+  GeglRectangle  cur_col = *rect;
+  GeglRectangle  in_col;
+  const gint     nc = babl_format_get_n_components (format);
+  gfloat        *col;
+  gfloat        *out;
+  gint           v;
+
+  cur_col.width = 1;
+
+  in_col         = cur_col;
+  in_col.height += clen - 1;
+  in_col.y      -= clen / 2;
+
+  col = gegl_malloc (sizeof (gfloat) * in_col.height  * nc);
+  out = gegl_malloc (sizeof (gfloat) * cur_col.height * nc);
+
+  for (v = 0; v < rect->width; v++)
+    {
+      cur_col.x = in_col.x = rect->x + v;
+
+      gegl_buffer_get (src, &in_col, 1.0, format, col, GEGL_AUTO_ROWSTRIDE, policy);
+
+      fir_blur_1D (col, out, cmatrix, clen, rect->height, nc);
+
+      gegl_buffer_set (dst, &cur_col, 0, format, out, GEGL_AUTO_ROWSTRIDE);
+    }
+
+  gegl_free (out);
+  gegl_free (col);
+}
+
+
+#include "opencl/gegl-cl.h"
+#include "buffer/gegl-buffer-cl-iterator.h"
+
+#include "opencl/gblur-1d.cl.h"
+
+static GeglClRunData *cl_data = NULL;
+
+
+static gboolean
+cl_gaussian_blur (cl_mem                 in_tex,
+                  cl_mem                 out_tex,
+                  const GeglRectangle   *roi,
+                  cl_mem                 cl_cmatrix,
+                  gint                   clen,
+                  GeglGblur1dOrientation orientation)
+{
+  cl_int cl_err = 0;
+  size_t global_ws[2];
+  gint   kernel_num;
+
+  if (!cl_data)
+    {
+      const char *kernel_name[] = {"fir_ver_blur", "fir_hor_blur", NULL};
+      cl_data = gegl_cl_compile_and_build (gblur_1d_cl_source, kernel_name);
+    }
+
+  if (!cl_data)
+    return TRUE;
+
+  if (orientation == GEGL_GBLUR_1D_VERTICAL)
+    kernel_num = 0;
+  else
+    kernel_num = 1;
+
+  global_ws[0] = roi->width;
+  global_ws[1] = roi->height;
+
+  cl_err = gegl_cl_set_kernel_args (cl_data->kernel[kernel_num],
+                                    sizeof(cl_mem), (void*)&in_tex,
+                                    sizeof(cl_mem), (void*)&out_tex,
+                                    sizeof(cl_mem), (void*)&cl_cmatrix,
+                                    sizeof(cl_int), (void*)&clen,
+                                    NULL);
+  CL_CHECK;
+
+  cl_err = gegl_clEnqueueNDRangeKernel (gegl_cl_get_command_queue (),
+                                        cl_data->kernel[kernel_num], 2,
+                                        NULL, global_ws, NULL,
+                                        0, NULL, NULL);
+  CL_CHECK;
+
+  cl_err = gegl_clFinish (gegl_cl_get_command_queue ());
+  CL_CHECK;
+
+  return FALSE;
+
+error:
+  return TRUE;
+}
+
+static gboolean
+fir_cl_process (GeglBuffer            *input,
+                GeglBuffer            *output,
+                const GeglRectangle   *result,
+                const Babl            *format,
+                gfloat                *cmatrix,
+                gint                   clen,
+                GeglGblur1dOrientation orientation,
+                GeglAbyssPolicy        abyss)
+{
+  gboolean              err = FALSE;
+  cl_int                cl_err;
+  cl_mem                cl_cmatrix = NULL;
+  GeglBufferClIterator *i;
+  gint                  read;
+  gint                  left, right, top, bottom;
+
+  if (orientation == GEGL_GBLUR_1D_HORIZONTAL)
+    {
+      right = left = clen / 2;
+      top = bottom = 0;
+    }
+  else
+    {
+      right = left = 0;
+      top = bottom = clen / 2;
+    }
+
+  i = gegl_buffer_cl_iterator_new (output,
+                                   result,
+                                   format,
+                                   GEGL_CL_BUFFER_WRITE);
+
+  read = gegl_buffer_cl_iterator_add_2 (i,
+                                        input,
+                                        result,
+                                        format,
+                                        GEGL_CL_BUFFER_READ,
+                                        left, right,
+                                        top, bottom,
+                                        abyss);
+
+  cl_cmatrix = gegl_clCreateBuffer (gegl_cl_get_context(),
+                                    CL_MEM_COPY_HOST_PTR | CL_MEM_READ_ONLY,
+                                    clen * sizeof(cl_float), cmatrix, &cl_err);
+  CL_CHECK;
+
+  while (gegl_buffer_cl_iterator_next (i, &err) && !err)
+    {
+      err = cl_gaussian_blur(i->tex[read],
+                             i->tex[0],
+                             &i->roi[0],
+                             cl_cmatrix,
+                             clen,
+                             orientation);
+
+      if (err)
+        {
+          gegl_buffer_cl_iterator_stop (i);
+          break;
+        }
+    }
+
+  cl_err = gegl_clReleaseMemObject (cl_cmatrix);
+  CL_CHECK;
+
+  cl_cmatrix = NULL;
+
+  return !err;
+
+error:
+  if (cl_cmatrix)
+    gegl_clReleaseMemObject (cl_cmatrix);
+
+  return FALSE;
+}
+
+static gfloat
+gaussian_func_1d (gfloat x,
+                  gfloat sigma)
+{
+  return (1.0 / (sigma * sqrt (2.0 * G_PI))) *
+          exp (-(x * x) / (2.0 * sigma * sigma));
+}
+
+static gint
+fir_calc_convolve_matrix_length (gfloat sigma)
+{
+#if 1
+  return sigma > GEGL_FLOAT_EPSILON ? ceil (sigma) * 6 + 1 : 1;
+#else
+  if (sigma > GEGL_FLOAT_EPSILON)
+    {
+      gint x = 0;
+      while (gaussian_func_1d (x, sigma) > 1e-3)
+        x++;
+      return 2 * x + 1;
+    }
+  else return 1;
+#endif
+}
+
+static gint
+fir_gen_convolve_matrix (gfloat   sigma,
+                         gfloat **cmatrix)
+{
+  gint    clen;
+  gfloat *cmatrix_p;
+
+  clen = fir_calc_convolve_matrix_length (sigma);
+
+  *cmatrix  = gegl_malloc (sizeof (gfloat) * clen);
+  cmatrix_p = *cmatrix;
+
+  if (clen == 1)
+    {
+      cmatrix_p [0] = 1;
+    }
+  else
+    {
+      gint    i;
+      gdouble sum = 0;
+      gint    half_clen = clen / 2;
+
+      for (i = 0; i < clen; i++)
+        {
+          cmatrix_p [i] = gaussian_func_1d (i - half_clen, sigma);
+          sum += cmatrix_p [i];
+        }
+
+      for (i = 0; i < clen; i++)
+        {
+          cmatrix_p [i] /= sum;
+        }
+    }
+
+  return clen;
+}
+
+static GeglGblur1dFilter
+filter_disambiguation (GeglGblur1dFilter filter,
+                       gfloat            std_dev)
+{
+  if (filter == GEGL_GBLUR_1D_AUTO)
+    {
+      if (std_dev < 1.0)
+        filter = GEGL_GBLUR_1D_FIR;
+      else
+        filter = GEGL_GBLUR_1D_IIR;
+    }
+  return filter;
+}
+
+
+/**********************************************
+ *
+ * GEGL Operation API
+ *
+ **********************************************/
+
 static void
 gegl_gblur_1d_prepare (GeglOperation *operation)
 {
@@ -269,23 +618,45 @@ gegl_gblur_1d_get_required_for_output (GeglOperation       *operation,
                                        const GeglRectangle *output_roi)
 {
   GeglRectangle        required_for_output = { 0, };
-  const GeglRectangle *in_rect = gegl_operation_source_get_bounding_box (operation, input_pad);
+  GeglChantO          *o       = GEGL_CHANT_PROPERTIES (operation);
+  GeglGblur1dFilter    filter  = filter_disambiguation (o->filter, o->std_dev);
+
+  if (filter == GEGL_GBLUR_1D_IIR)
+    {
+      const GeglRectangle *in_rect =
+        gegl_operation_source_get_bounding_box (operation, input_pad);
 
-  if (in_rect && ! gegl_rectangle_is_infinite_plane (in_rect))
+      if (in_rect && ! gegl_rectangle_is_infinite_plane (in_rect))
+        {
+          required_for_output = *output_roi;
+
+          if (o->orientation == GEGL_GBLUR_1D_HORIZONTAL)
+            {
+              required_for_output.x     = in_rect->x;
+              required_for_output.width = in_rect->width;
+            }
+          else
+            {
+              required_for_output.y      = in_rect->y;
+              required_for_output.height = in_rect->height;
+            }
+        }
+    }
+  else
     {
-      GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);
+      gint clen = fir_calc_convolve_matrix_length (o->std_dev);
 
       required_for_output = *output_roi;
 
-      if (o->orientation == GEGL_BLUR_1D_HORIZONTAL)
+      if (o->orientation == GEGL_GBLUR_1D_HORIZONTAL)
         {
-          required_for_output.x     = in_rect->x;
-          required_for_output.width = in_rect->width;
+          required_for_output.x     -= clen / 2;
+          required_for_output.width += clen - 1;
         }
       else
         {
-          required_for_output.y      = in_rect->y;
-          required_for_output.height = in_rect->height;
+          required_for_output.y      -= clen / 2;
+          required_for_output.height += clen - 1;
         }
     }
 
@@ -296,31 +667,64 @@ static GeglRectangle
 gegl_gblur_1d_get_cached_region (GeglOperation       *operation,
                                  const GeglRectangle *output_roi)
 {
-  GeglRectangle        cached_region = { 0, };
-  const GeglRectangle *in_rect = gegl_operation_source_get_bounding_box (operation,
-                                                                         "input");
+  GeglRectangle      cached_region = { 0, };
+  GeglChantO        *o       = GEGL_CHANT_PROPERTIES (operation);
+  GeglGblur1dFilter  filter  = filter_disambiguation (o->filter, o->std_dev);
 
-  if (in_rect && ! gegl_rectangle_is_infinite_plane (in_rect))
+  if (filter == GEGL_GBLUR_1D_IIR)
     {
-      GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);
-
-      cached_region = *output_roi;
+      const GeglRectangle *in_rect =
+        gegl_operation_source_get_bounding_box (operation, "input");
 
-      if (o->orientation == GEGL_BLUR_1D_HORIZONTAL)
+      if (in_rect && ! gegl_rectangle_is_infinite_plane (in_rect))
         {
-          cached_region.x     = in_rect->x;
-          cached_region.width = in_rect->width;
-        }
-      else
-        {
-          cached_region.y      = in_rect->y;
-          cached_region.height = in_rect->height;
+          GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);
+
+          cached_region = *output_roi;
+
+          if (o->orientation == GEGL_GBLUR_1D_HORIZONTAL)
+            {
+              cached_region.x     = in_rect->x;
+              cached_region.width = in_rect->width;
+            }
+          else
+            {
+              cached_region.y      = in_rect->y;
+              cached_region.height = in_rect->height;
+            }
         }
     }
+  else
+    {
+      cached_region = *output_roi;
+    }
 
   return cached_region;
 }
 
+static GeglAbyssPolicy
+to_gegl_policy (GeglGblur1dPolicy policy)
+{
+  switch (policy)
+    {
+    case (GEGL_GBLUR_1D_ABYSS_NONE):
+      return GEGL_ABYSS_NONE;
+      break;
+    case (GEGL_GBLUR_1D_ABYSS_CLAMP):
+      return GEGL_ABYSS_CLAMP;
+      break;
+    case (GEGL_GBLUR_1D_ABYSS_WHITE):
+      return GEGL_ABYSS_WHITE;
+      break;
+    case (GEGL_GBLUR_1D_ABYSS_BLACK):
+      return GEGL_ABYSS_BLACK;
+      break;
+    default:
+      g_warning ("gblur-1d: unsupported abyss policy");
+      return GEGL_ABYSS_NONE;
+    }
+}
+
 static gboolean
 gegl_gblur_1d_process (GeglOperation       *operation,
                        GeglBuffer          *input,
@@ -328,16 +732,46 @@ gegl_gblur_1d_process (GeglOperation       *operation,
                        const GeglRectangle *result,
                        gint                 level)
 {
-  GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);
+  GeglChantO *o      = GEGL_CHANT_PROPERTIES (operation);
   const Babl *format = gegl_operation_get_format (operation, "output");
-  gdouble b[4], m[3][3];
 
-  iir_young_find_constants (o->std_dev, b, m);
+  GeglGblur1dFilter filter;
+  GeglAbyssPolicy   abyss_policy = to_gegl_policy (o->abyss_policy);
+
+  filter = filter_disambiguation (o->filter, o->std_dev);
+
+  if (filter == GEGL_GBLUR_1D_IIR)
+    {
+      gdouble b[4], m[3][3];
 
-  if (o->orientation == GEGL_BLUR_1D_HORIZONTAL)
-    iir_young_hor_blur (input, result, output, b, m, o->abyss_policy, format);
+      iir_young_find_constants (o->std_dev, b, m);
+
+      if (o->orientation == GEGL_GBLUR_1D_HORIZONTAL)
+        iir_young_hor_blur (input, result, output, b, m, abyss_policy, format);
+      else
+        iir_young_ver_blur (input, result, output, b, m, abyss_policy, format);
+    }
   else
-    iir_young_ver_blur (input, result, output, b, m, o->abyss_policy, format);
+    {
+      gfloat *cmatrix;
+      gint    clen;
+
+      clen = fir_gen_convolve_matrix (o->std_dev, &cmatrix);
+
+      /* FIXME: implement others format cases */
+      if (gegl_operation_use_opencl (operation) &&
+          format == babl_format ("RaGaBaA float"))
+        if (fir_cl_process(input, output, result, format,
+                           cmatrix, clen, o->orientation, abyss_policy))
+          return TRUE;
+
+      if (o->orientation == GEGL_GBLUR_1D_HORIZONTAL)
+        fir_hor_blur (input, result, output, cmatrix, clen, abyss_policy, format);
+      else
+        fir_ver_blur (input, result, output, cmatrix, clen, abyss_policy, format);
+
+      gegl_free (cmatrix);
+    }
 
   return  TRUE;
 }
@@ -355,6 +789,7 @@ gegl_chant_class_init (GeglChantClass *klass)
   operation_class->prepare                 = gegl_gblur_1d_prepare;
   operation_class->get_required_for_output = gegl_gblur_1d_get_required_for_output;
   operation_class->get_cached_region       = gegl_gblur_1d_get_cached_region;
+  operation_class->opencl_support          = TRUE;
   gegl_operation_class_set_keys (operation_class,
     "name",       "gegl:gblur-1d",
     "categories", "hidden:blur",