Blame operations/common/gaussian-blur.c

Packit Service 2781ba
/* This file is an image processing operation for GEGL
Packit Service 2781ba
 * GEGL is free software; you can redistribute it and/or
Packit Service 2781ba
 * modify it under the terms of the GNU Lesser General Public
Packit Service 2781ba
 * License as published by the Free Software Foundation; either
Packit Service 2781ba
 * version 3 of the License, or (at your option) any later version.
Packit Service 2781ba
 *
Packit Service 2781ba
 * GEGL is distributed in the hope that it will be useful,
Packit Service 2781ba
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Packit Service 2781ba
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Packit Service 2781ba
 * Lesser General Public License for more details.
Packit Service 2781ba
 *
Packit Service 2781ba
 * You should have received a copy of the GNU Lesser General Public
Packit Service 2781ba
 * License along with GEGL; if not, see <http://www.gnu.org/licenses/>.
Packit Service 2781ba
 *
Packit Service 2781ba
 * Copyright 2006 Dominik Ernst <dernst@gmx.de>
Packit Service 2781ba
 *
Packit Service 2781ba
 * Recursive Gauss IIR Filter as described by Young / van Vliet
Packit Service 2781ba
 * in "Signal Processing 44 (1995) 139 - 151"
Packit Service 2781ba
 *
Packit Service 2781ba
 * NOTE: The IIR filter should not be used for radius < 0.5, since it
Packit Service 2781ba
 *       becomes very inaccurate.
Packit Service 2781ba
 */
Packit Service 2781ba
Packit Service 2781ba
#include "config.h"
Packit Service 2781ba
#include <glib/gi18n-lib.h>
Packit Service 2781ba
Packit Service 2781ba
Packit Service 2781ba
#ifdef GEGL_CHANT_PROPERTIES
Packit Service 2781ba
Packit Service 2781ba
gegl_chant_double_ui (std_dev_x, _("Size X"), 0.0, 10000.0, 4.0, 0.0, 1000.0, 1.5,
Packit Service 2781ba
   _("Standard deviation for the horizontal axis. (multiply by ~2 to get radius)"))
Packit Service 2781ba
gegl_chant_double_ui (std_dev_y, _("Size Y"), 0.0, 10000.0, 4.0, 0.0, 1000.0, 1.5,
Packit Service 2781ba
   _("Standard deviation for the vertical axis. (multiply by ~2 to get radius.)"))
Packit Service 2781ba
gegl_chant_string (filter, _("Filter"), "auto",
Packit Service 2781ba
   _("Optional parameter to override the automatic selection of blur filter. "
Packit Service 2781ba
     "Choices are fir, iir, auto"))
Packit Service 2781ba
Packit Service 2781ba
#else
Packit Service 2781ba
Packit Service 2781ba
#define GEGL_CHANT_TYPE_AREA_FILTER
Packit Service 2781ba
#define GEGL_CHANT_C_FILE       "gaussian-blur.c"
Packit Service 2781ba
Packit Service 2781ba
#include "gegl-chant.h"
Packit Service 2781ba
#include <math.h>
Packit Service 2781ba
#include <stdio.h>
Packit Service 2781ba
Packit Service 2781ba
#define RADIUS_SCALE   4
Packit Service 2781ba
Packit Service 2781ba
Packit Service 2781ba
static void
Packit Service 2781ba
iir_young_find_constants (gfloat   radius,
Packit Service 2781ba
                          gdouble *B,
Packit Service 2781ba
                          gdouble *b);
Packit Service 2781ba
Packit Service 2781ba
static gint
Packit Service 2781ba
fir_gen_convolve_matrix (gdouble   sigma,
Packit Service 2781ba
                         gdouble **cmatrix_p);
Packit Service 2781ba
Packit Service 2781ba
Packit Service 2781ba
static void
Packit Service 2781ba
iir_young_find_constants (gfloat   sigma,
Packit Service 2781ba
                          gdouble *B,
Packit Service 2781ba
                          gdouble *b)
Packit Service 2781ba
{
Packit Service 2781ba
  gdouble q;
Packit Service 2781ba
Packit Service 2781ba
  if (sigma == 0.0)
Packit Service 2781ba
    {
Packit Service 2781ba
      /* to prevent unexpected ringing at tile boundaries,
Packit Service 2781ba
         we define an (expensive) copy operation here */
Packit Service 2781ba
      *B = 1.0;
Packit Service 2781ba
      b[0] = 1.0;
Packit Service 2781ba
      b[1] = b[2] = b[3] = 0.0;
Packit Service 2781ba
      return;
Packit Service 2781ba
    }
Packit Service 2781ba
Packit Service 2781ba
  if(sigma >= 2.5)
Packit Service 2781ba
    q = 0.98711*sigma - 0.96330;
Packit Service 2781ba
  else
Packit Service 2781ba
    q = 3.97156 - 4.14554*sqrt(1-0.26891*sigma);
Packit Service 2781ba
Packit Service 2781ba
  b[0] = 1.57825 + (2.44413*q) + (1.4281*q*q) + (0.422205*q*q*q);
Packit Service 2781ba
  b[1] = (2.44413*q) + (2.85619*q*q) + (1.26661*q*q*q);
Packit Service 2781ba
  b[2] = -((1.4281*q*q) + (1.26661*q*q*q));
Packit Service 2781ba
  b[3] = 0.422205*q*q*q;
Packit Service 2781ba
Packit Service 2781ba
  *B = 1 - ( (b[1]+b[2]+b[3])/b[0] );
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
static inline void
Packit Service 2781ba
iir_young_blur_1D (gfloat  * buf,
Packit Service 2781ba
                   gint      offset,
Packit Service 2781ba
                   gint      delta_offset,
Packit Service 2781ba
                   gdouble   B,
Packit Service 2781ba
                   gdouble * b,
Packit Service 2781ba
                   gfloat  * w,
Packit Service 2781ba
                   gint      w_len)
Packit Service 2781ba
{
Packit Service 2781ba
  gint wcount, i;
Packit Service 2781ba
  gdouble tmp;
Packit Service 2781ba
  gdouble recip = 1.0 / b[0];
Packit Service 2781ba
Packit Service 2781ba
  /* forward filter */
Packit Service 2781ba
  wcount = 0;
Packit Service 2781ba
Packit Service 2781ba
  while (wcount < w_len)
Packit Service 2781ba
    {
Packit Service 2781ba
      tmp = 0;
Packit Service 2781ba
Packit Service 2781ba
      for (i=1; i<4; i++)
Packit Service 2781ba
        {
Packit Service 2781ba
          if (wcount-i >= 0)
Packit Service 2781ba
            tmp += b[i]*w[wcount-i];
Packit Service 2781ba
        }
Packit Service 2781ba
Packit Service 2781ba
      tmp *= recip;
Packit Service 2781ba
      tmp += B*buf[offset];
Packit Service 2781ba
      w[wcount] = tmp;
Packit Service 2781ba
Packit Service 2781ba
      wcount++;
Packit Service 2781ba
      offset += delta_offset;
Packit Service 2781ba
    }
Packit Service 2781ba
Packit Service 2781ba
  /* backward filter */
Packit Service 2781ba
  wcount = w_len - 1;
Packit Service 2781ba
  offset -= delta_offset;
Packit Service 2781ba
Packit Service 2781ba
  while (wcount >= 0)
Packit Service 2781ba
    {
Packit Service 2781ba
      tmp = 0;
Packit Service 2781ba
Packit Service 2781ba
      for (i=1; i<4; i++)
Packit Service 2781ba
        {
Packit Service 2781ba
          if (wcount+i < w_len)
Packit Service 2781ba
            tmp += b[i]*buf[offset+delta_offset*i];
Packit Service 2781ba
        }
Packit Service 2781ba
Packit Service 2781ba
      tmp *= recip;
Packit Service 2781ba
      tmp += B*w[wcount];
Packit Service 2781ba
      buf[offset] = tmp;
Packit Service 2781ba
Packit Service 2781ba
      offset -= delta_offset;
Packit Service 2781ba
      wcount--;
Packit Service 2781ba
    }
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
/* expects src and dst buf to have the same height and no y-offset */
Packit Service 2781ba
static void
Packit Service 2781ba
iir_young_hor_blur (GeglBuffer          *src,
Packit Service 2781ba
                    const GeglRectangle *src_rect,
Packit Service 2781ba
                    GeglBuffer          *dst,
Packit Service 2781ba
                    const GeglRectangle *dst_rect,
Packit Service 2781ba
                    gdouble              B,
Packit Service 2781ba
                    gdouble             *b)
Packit Service 2781ba
{
Packit Service 2781ba
  gint v;
Packit Service 2781ba
  gint c;
Packit Service 2781ba
  gint w_len;
Packit Service 2781ba
  gfloat *buf;
Packit Service 2781ba
  gfloat *w;
Packit Service 2781ba
Packit Service 2781ba
  buf = g_new0 (gfloat, src_rect->height * src_rect->width * 4);
Packit Service 2781ba
  w   = g_new0 (gfloat, src_rect->width);
Packit Service 2781ba
Packit Service 2781ba
  gegl_buffer_get (src, src_rect, 1.0, babl_format ("RaGaBaA float"),
Packit Service 2781ba
                   buf, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
Packit Service 2781ba
Packit Service 2781ba
  w_len = src_rect->width;
Packit Service 2781ba
  for (v=0; v<src_rect->height; v++)
Packit Service 2781ba
    {
Packit Service 2781ba
      for (c = 0; c < 4; c++)
Packit Service 2781ba
        {
Packit Service 2781ba
          iir_young_blur_1D (buf,
Packit Service 2781ba
                             v*src_rect->width*4+c,
Packit Service 2781ba
                             4,
Packit Service 2781ba
                             B,
Packit Service 2781ba
                             b,
Packit Service 2781ba
                             w,
Packit Service 2781ba
                             w_len);
Packit Service 2781ba
        }
Packit Service 2781ba
    }
Packit Service 2781ba
Packit Service 2781ba
  gegl_buffer_set (dst, src_rect, 0.0, babl_format ("RaGaBaA float"),
Packit Service 2781ba
                   buf, GEGL_AUTO_ROWSTRIDE);
Packit Service 2781ba
  g_free (buf);
Packit Service 2781ba
  g_free (w);
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
/* expects src and dst buf to have the same width and no x-offset */
Packit Service 2781ba
static void
Packit Service 2781ba
iir_young_ver_blur (GeglBuffer          *src,
Packit Service 2781ba
                    const GeglRectangle *src_rect,
Packit Service 2781ba
                    GeglBuffer          *dst,
Packit Service 2781ba
                    const GeglRectangle *dst_rect,
Packit Service 2781ba
                    gdouble              B,
Packit Service 2781ba
                    gdouble             *b)
Packit Service 2781ba
{
Packit Service 2781ba
  gint u;
Packit Service 2781ba
  gint c;
Packit Service 2781ba
  gint w_len;
Packit Service 2781ba
  gfloat *buf;
Packit Service 2781ba
  gfloat *w;
Packit Service 2781ba
Packit Service 2781ba
  buf = g_new0 (gfloat, src_rect->height * src_rect->width * 4);
Packit Service 2781ba
  w   = g_new0 (gfloat, src_rect->height);
Packit Service 2781ba
Packit Service 2781ba
  gegl_buffer_get (src, src_rect, 1.0, babl_format ("RaGaBaA float"),
Packit Service 2781ba
                   buf, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
Packit Service 2781ba
Packit Service 2781ba
  w_len = src_rect->height;
Packit Service 2781ba
  for (u=0; u<dst_rect->width; u++)
Packit Service 2781ba
    {
Packit Service 2781ba
      for (c = 0; c < 4; c++)
Packit Service 2781ba
        {
Packit Service 2781ba
          iir_young_blur_1D (buf,
Packit Service 2781ba
                             u*4 + c,
Packit Service 2781ba
                             src_rect->width*4,
Packit Service 2781ba
                             B,
Packit Service 2781ba
                             b,
Packit Service 2781ba
                             w,
Packit Service 2781ba
                             w_len);
Packit Service 2781ba
        }
Packit Service 2781ba
    }
Packit Service 2781ba
Packit Service 2781ba
  gegl_buffer_set (dst, src_rect, 0,
Packit Service 2781ba
                   babl_format ("RaGaBaA float"), buf, GEGL_AUTO_ROWSTRIDE);
Packit Service 2781ba
  g_free (buf);
Packit Service 2781ba
  g_free (w);
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
Packit Service 2781ba
static gint
Packit Service 2781ba
fir_calc_convolve_matrix_length (gdouble sigma)
Packit Service 2781ba
{
Packit Service 2781ba
  return sigma ? ceil (sigma)*6.0+1.0 : 1;
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
static gint
Packit Service 2781ba
fir_gen_convolve_matrix (gdouble   sigma,
Packit Service 2781ba
                         gdouble **cmatrix_p)
Packit Service 2781ba
{
Packit Service 2781ba
  gint     matrix_length;
Packit Service 2781ba
  gdouble *cmatrix;
Packit Service 2781ba
Packit Service 2781ba
  matrix_length = fir_calc_convolve_matrix_length (sigma);
Packit Service 2781ba
  cmatrix = g_new (gdouble, matrix_length);
Packit Service 2781ba
  if (!cmatrix)
Packit Service 2781ba
    return 0;
Packit Service 2781ba
Packit Service 2781ba
  if (matrix_length == 1)
Packit Service 2781ba
    {
Packit Service 2781ba
      cmatrix[0] = 1;
Packit Service 2781ba
    }
Packit Service 2781ba
  else
Packit Service 2781ba
    {
Packit Service 2781ba
      gint i,x;
Packit Service 2781ba
      gdouble sum = 0;
Packit Service 2781ba
Packit Service 2781ba
      for (i=0; i
Packit Service 2781ba
        {
Packit Service 2781ba
          gdouble y;
Packit Service 2781ba
Packit Service 2781ba
          x = i - (matrix_length/2);
Packit Service 2781ba
          y = (1.0/(sigma*sqrt(2.0*G_PI))) * exp(-(x*x) / (2.0*sigma*sigma));
Packit Service 2781ba
Packit Service 2781ba
          cmatrix[i] = y;
Packit Service 2781ba
          sum += cmatrix[i];
Packit Service 2781ba
        }
Packit Service 2781ba
Packit Service 2781ba
      for (i=matrix_length/2 + 1; i
Packit Service 2781ba
        {
Packit Service 2781ba
          cmatrix[i] = cmatrix[matrix_length-i-1];
Packit Service 2781ba
          sum += cmatrix[i];
Packit Service 2781ba
        }
Packit Service 2781ba
Packit Service 2781ba
      for (i=0; i
Packit Service 2781ba
      {
Packit Service 2781ba
        cmatrix[i] /= sum;
Packit Service 2781ba
      }
Packit Service 2781ba
    }
Packit Service 2781ba
Packit Service 2781ba
  *cmatrix_p = cmatrix;
Packit Service 2781ba
  return matrix_length;
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
static inline float
Packit Service 2781ba
fir_get_mean_component_1D (gfloat  * buf,
Packit Service 2781ba
                           gint      offset,
Packit Service 2781ba
                           gint      delta_offset,
Packit Service 2781ba
                           gdouble * cmatrix,
Packit Service 2781ba
                           gint      matrix_length)
Packit Service 2781ba
{
Packit Service 2781ba
  gint    i;
Packit Service 2781ba
  gdouble acc=0;
Packit Service 2781ba
Packit Service 2781ba
  for (i=0; i < matrix_length; i++)
Packit Service 2781ba
    {
Packit Service 2781ba
      acc += buf[offset] * cmatrix[i];
Packit Service 2781ba
      offset += delta_offset;
Packit Service 2781ba
    }
Packit Service 2781ba
Packit Service 2781ba
  return acc;
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
/* expects src and dst buf to have the same height and no y-offset */
Packit Service 2781ba
static void
Packit Service 2781ba
fir_hor_blur (GeglBuffer          *src,
Packit Service 2781ba
              const GeglRectangle *src_rect,
Packit Service 2781ba
              GeglBuffer          *dst,
Packit Service 2781ba
              const GeglRectangle *dst_rect,
Packit Service 2781ba
              gdouble             *cmatrix,
Packit Service 2781ba
              gint                 matrix_length,
Packit Service 2781ba
              gint                 xoff) /* offset between src and dst */
Packit Service 2781ba
{
Packit Service 2781ba
  gint        u,v;
Packit Service 2781ba
  gint        offset;
Packit Service 2781ba
  gfloat     *src_buf;
Packit Service 2781ba
  gfloat     *dst_buf;
Packit Service 2781ba
  const gint  radius = matrix_length/2;
Packit Service 2781ba
  const gint  src_width = src_rect->width;
Packit Service 2781ba
Packit Service 2781ba
  g_assert (xoff >= radius);
Packit Service 2781ba
Packit Service 2781ba
  src_buf = g_new0 (gfloat, src_rect->height * src_rect->width * 4);
Packit Service 2781ba
  dst_buf = g_new0 (gfloat, dst_rect->height * dst_rect->width * 4);
Packit Service 2781ba
Packit Service 2781ba
  gegl_buffer_get (src, src_rect, 1.0, babl_format ("RaGaBaA float"),
Packit Service 2781ba
                   src_buf, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
Packit Service 2781ba
Packit Service 2781ba
  offset = 0;
Packit Service 2781ba
  for (v=0; v<dst_rect->height; v++)
Packit Service 2781ba
    for (u=0; u<dst_rect->width; u++)
Packit Service 2781ba
      {
Packit Service 2781ba
        gint src_offset = (u-radius+xoff + v*src_width) * 4;
Packit Service 2781ba
        gint c;
Packit Service 2781ba
        for (c=0; c<4; c++)
Packit Service 2781ba
          dst_buf [offset++] = fir_get_mean_component_1D (src_buf,
Packit Service 2781ba
                                                          src_offset + c,
Packit Service 2781ba
                                                          4,
Packit Service 2781ba
                                                          cmatrix,
Packit Service 2781ba
                                                          matrix_length);
Packit Service 2781ba
      }
Packit Service 2781ba
Packit Service 2781ba
  gegl_buffer_set (dst, dst_rect, 0.0, babl_format ("RaGaBaA float"),
Packit Service 2781ba
                   dst_buf, GEGL_AUTO_ROWSTRIDE);
Packit Service 2781ba
  g_free (src_buf);
Packit Service 2781ba
  g_free (dst_buf);
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
/* expects src and dst buf to have the same width and no x-offset */
Packit Service 2781ba
static void
Packit Service 2781ba
fir_ver_blur (GeglBuffer          *src,
Packit Service 2781ba
              const GeglRectangle *src_rect,
Packit Service 2781ba
              GeglBuffer          *dst,
Packit Service 2781ba
              const GeglRectangle *dst_rect,
Packit Service 2781ba
              gdouble             *cmatrix,
Packit Service 2781ba
              gint                 matrix_length,
Packit Service 2781ba
              gint                 yoff) /* offset between src and dst */
Packit Service 2781ba
{
Packit Service 2781ba
  gint        u,v;
Packit Service 2781ba
  gint        offset;
Packit Service 2781ba
  gfloat     *src_buf;
Packit Service 2781ba
  gfloat     *dst_buf;
Packit Service 2781ba
  const gint  radius = matrix_length/2;
Packit Service 2781ba
  const gint  src_width = src_rect->width;
Packit Service 2781ba
Packit Service 2781ba
  g_assert (yoff >= radius);
Packit Service 2781ba
Packit Service 2781ba
  src_buf = g_new0 (gfloat, src_rect->width * src_rect->height * 4);
Packit Service 2781ba
  dst_buf = g_new0 (gfloat, dst_rect->width * dst_rect->height * 4);
Packit Service 2781ba
Packit Service 2781ba
  gegl_buffer_get (src, src_rect, 1.0, babl_format ("RaGaBaA float"),
Packit Service 2781ba
                   src_buf, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
Packit Service 2781ba
Packit Service 2781ba
  offset=0;
Packit Service 2781ba
  for (v=0; v< dst_rect->height; v++)
Packit Service 2781ba
    for (u=0; u< dst_rect->width; u++)
Packit Service 2781ba
      {
Packit Service 2781ba
        gint src_offset = (u + (v-radius+yoff)*src_width) * 4;
Packit Service 2781ba
        gint c;
Packit Service 2781ba
        for (c=0; c<4; c++)
Packit Service 2781ba
          dst_buf [offset++] = fir_get_mean_component_1D (src_buf,
Packit Service 2781ba
                                                          src_offset + c,
Packit Service 2781ba
                                                          src_width * 4,
Packit Service 2781ba
                                                          cmatrix,
Packit Service 2781ba
                                                          matrix_length);
Packit Service 2781ba
      }
Packit Service 2781ba
Packit Service 2781ba
  gegl_buffer_set (dst, dst_rect, 0, babl_format ("RaGaBaA float"),
Packit Service 2781ba
                   dst_buf, GEGL_AUTO_ROWSTRIDE);
Packit Service 2781ba
  g_free (src_buf);
Packit Service 2781ba
  g_free (dst_buf);
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
static void prepare (GeglOperation *operation)
Packit Service 2781ba
{
Packit Service 2781ba
#define max(A,B) ((A) > (B) ? (A) : (B))
Packit Service 2781ba
  GeglOperationAreaFilter *area = GEGL_OPERATION_AREA_FILTER (operation);
Packit Service 2781ba
  GeglChantO              *o    = GEGL_CHANT_PROPERTIES (operation);
Packit Service 2781ba
Packit Service 2781ba
  gfloat fir_radius_x = fir_calc_convolve_matrix_length (o->std_dev_x) / 2;
Packit Service 2781ba
  gfloat fir_radius_y = fir_calc_convolve_matrix_length (o->std_dev_y) / 2;
Packit Service 2781ba
  gfloat iir_radius_x = o->std_dev_x * RADIUS_SCALE;
Packit Service 2781ba
  gfloat iir_radius_y = o->std_dev_y * RADIUS_SCALE;
Packit Service 2781ba
Packit Service 2781ba
  /* XXX: these should be calculated exactly considering o->filter, but we just
Packit Service 2781ba
   * make sure there is enough space */
Packit Service 2781ba
  area->left = area->right = ceil ( max (fir_radius_x, iir_radius_x));
Packit Service 2781ba
  area->top = area->bottom = ceil ( max (fir_radius_y, iir_radius_y));
Packit Service 2781ba
Packit Service 2781ba
  gegl_operation_set_format (operation, "input",
Packit Service 2781ba
                             babl_format ("RaGaBaA float"));
Packit Service 2781ba
  gegl_operation_set_format (operation, "output",
Packit Service 2781ba
                             babl_format ("RaGaBaA float"));
Packit Service 2781ba
#undef max
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
#include "opencl/gegl-cl.h"
Packit Service 2781ba
#include "buffer/gegl-buffer-cl-iterator.h"
Packit Service 2781ba
Packit Service 2781ba
static const char* kernel_source =
Packit Service 2781ba
"float4 fir_get_mean_component_1D_CL(const global float4 *buf,     \n"
Packit Service 2781ba
"                                    int offset,                   \n"
Packit Service 2781ba
"                                    const int delta_offset,       \n"
Packit Service 2781ba
"                                    constant float *cmatrix,      \n"
Packit Service 2781ba
"                                    const int matrix_length)      \n"
Packit Service 2781ba
"{                                                                 \n"
Packit Service 2781ba
"    float4 acc = 0.0f;                                            \n"
Packit Service 2781ba
"    int i;                                                        \n"
Packit Service 2781ba
"                                                                  \n"
Packit Service 2781ba
"    for(i=0; i
Packit Service 2781ba
"      {                                                           \n"
Packit Service 2781ba
"        acc    += buf[offset] * cmatrix[i];                       \n"
Packit Service 2781ba
"        offset += delta_offset;                                   \n"
Packit Service 2781ba
"      }                                                           \n"
Packit Service 2781ba
"    return acc;                                                   \n"
Packit Service 2781ba
"}                                                                 \n"
Packit Service 2781ba
"                                                                  \n"
Packit Service 2781ba
"__kernel void fir_ver_blur_CL(const global float4 *src_buf,       \n"
Packit Service 2781ba
"                              const int src_width,                \n"
Packit Service 2781ba
"                              global float4 *dst_buf,             \n"
Packit Service 2781ba
"                              constant float *cmatrix,            \n"
Packit Service 2781ba
"                              const int matrix_length,            \n"
Packit Service 2781ba
"                              const int yoff)                     \n"
Packit Service 2781ba
"{                                                                 \n"
Packit Service 2781ba
"    int gidx = get_global_id(0);                                  \n"
Packit Service 2781ba
"    int gidy = get_global_id(1);                                  \n"
Packit Service 2781ba
"    int gid  = gidx + gidy * get_global_size(0);                  \n"
Packit Service 2781ba
"                                                                  \n"
Packit Service 2781ba
"    int radius = matrix_length / 2;                               \n"
Packit Service 2781ba
"    int src_offset = gidx + (gidy - radius + yoff) * src_width;   \n"
Packit Service 2781ba
"                                                                  \n"
Packit Service 2781ba
"    dst_buf[gid] = fir_get_mean_component_1D_CL(                  \n"
Packit Service 2781ba
"        src_buf, src_offset, src_width, cmatrix, matrix_length);  \n"
Packit Service 2781ba
"}                                                                 \n"
Packit Service 2781ba
"                                                                  \n"
Packit Service 2781ba
"__kernel void fir_hor_blur_CL(const global float4 *src_buf,       \n"
Packit Service 2781ba
"                              const int src_width,                \n"
Packit Service 2781ba
"                              global float4 *dst_buf,             \n"
Packit Service 2781ba
"                              constant float *cmatrix,            \n"
Packit Service 2781ba
"                              const int matrix_length,            \n"
Packit Service 2781ba
"                              const int yoff)                     \n"
Packit Service 2781ba
"{                                                                 \n"
Packit Service 2781ba
"    int gidx = get_global_id(0);                                  \n"
Packit Service 2781ba
"    int gidy = get_global_id(1);                                  \n"
Packit Service 2781ba
"    int gid  = gidx + gidy * get_global_size(0);                  \n"
Packit Service 2781ba
"                                                                  \n"
Packit Service 2781ba
"    int radius = matrix_length / 2;                               \n"
Packit Service 2781ba
"    int src_offset = gidy * src_width + (gidx - radius + yoff);   \n"
Packit Service 2781ba
"                                                                  \n"
Packit Service 2781ba
"    dst_buf[gid] = fir_get_mean_component_1D_CL(                  \n"
Packit Service 2781ba
"        src_buf, src_offset, 1, cmatrix, matrix_length);          \n"
Packit Service 2781ba
"}                                                                 \n";
Packit Service 2781ba
Packit Service 2781ba
static gegl_cl_run_data *cl_data = NULL;
Packit Service 2781ba
Packit Service 2781ba
static cl_int
Packit Service 2781ba
cl_gaussian_blur (cl_mem                in_tex,
Packit Service 2781ba
                  cl_mem                out_tex,
Packit Service 2781ba
                  cl_mem                aux_tex,
Packit Service 2781ba
                  size_t                global_worksize,
Packit Service 2781ba
                  const GeglRectangle  *roi,
Packit Service 2781ba
                  const GeglRectangle  *src_rect,
Packit Service 2781ba
                  const GeglRectangle  *aux_rect,
Packit Service 2781ba
                  gfloat               *dmatrix_x,
Packit Service 2781ba
                  gint                  matrix_length_x,
Packit Service 2781ba
                  gint                  xoff,
Packit Service 2781ba
                  gfloat               *dmatrix_y,
Packit Service 2781ba
                  gint                  matrix_length_y,
Packit Service 2781ba
                  gint                  yoff)
Packit Service 2781ba
{
Packit Service 2781ba
  cl_int cl_err = 0;
Packit Service 2781ba
Packit Service 2781ba
  size_t global_ws[2];
Packit Service 2781ba
Packit Service 2781ba
  cl_mem cl_matrix_x, cl_matrix_y;
Packit Service 2781ba
Packit Service 2781ba
  if (!cl_data)
Packit Service 2781ba
    {
Packit Service 2781ba
      const char *kernel_name[] = {"fir_ver_blur_CL", "fir_hor_blur_CL", NULL};
Packit Service 2781ba
      cl_data = gegl_cl_compile_and_build (kernel_source, kernel_name);
Packit Service 2781ba
    }
Packit Service 2781ba
  if (!cl_data) return 1;
Packit Service 2781ba
Packit Service 2781ba
  cl_matrix_x = gegl_clCreateBuffer(gegl_cl_get_context(),
Packit Service 2781ba
                                    CL_MEM_ALLOC_HOST_PTR | CL_MEM_READ_ONLY,
Packit Service 2781ba
                                    matrix_length_x * sizeof(cl_float), NULL, &cl_err);
Packit Service 2781ba
  if (cl_err != CL_SUCCESS) return cl_err;
Packit Service 2781ba
Packit Service 2781ba
  cl_err = gegl_clEnqueueWriteBuffer(gegl_cl_get_command_queue(), cl_matrix_x,
Packit Service 2781ba
                                     CL_TRUE, 0, matrix_length_x * sizeof(cl_float), dmatrix_x,
Packit Service 2781ba
                                     0, NULL, NULL);
Packit Service 2781ba
  if (cl_err != CL_SUCCESS) return cl_err;
Packit Service 2781ba
Packit Service 2781ba
  cl_matrix_y = gegl_clCreateBuffer(gegl_cl_get_context(),
Packit Service 2781ba
                                    CL_MEM_ALLOC_HOST_PTR | CL_MEM_READ_ONLY,
Packit Service 2781ba
                                    matrix_length_y * sizeof(cl_float), NULL, &cl_err);
Packit Service 2781ba
  if (cl_err != CL_SUCCESS) return cl_err;
Packit Service 2781ba
Packit Service 2781ba
  cl_err = gegl_clEnqueueWriteBuffer(gegl_cl_get_command_queue(), cl_matrix_y,
Packit Service 2781ba
                                     CL_TRUE, 0, matrix_length_y * sizeof(cl_float), dmatrix_y,
Packit Service 2781ba
                                     0, NULL, NULL);
Packit Service 2781ba
  if (cl_err != CL_SUCCESS) return cl_err;
Packit Service 2781ba
Packit Service 2781ba
  {
Packit Service 2781ba
  global_ws[0] = aux_rect->width;
Packit Service 2781ba
  global_ws[1] = aux_rect->height;
Packit Service 2781ba
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[1], 0, sizeof(cl_mem), (void*)&in_tex);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[1], 1, sizeof(cl_int), (void*)&src_rect->width);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[1], 2, sizeof(cl_mem), (void*)&aux_tex);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[1], 3, sizeof(cl_mem), (void*)&cl_matrix_x);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[1], 4, sizeof(cl_int), (void*)&matrix_length_x);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[1], 5, sizeof(cl_int), (void*)&xoff);
Packit Service 2781ba
  if (cl_err != CL_SUCCESS) return cl_err;
Packit Service 2781ba
Packit Service 2781ba
  cl_err = gegl_clEnqueueNDRangeKernel(gegl_cl_get_command_queue (),
Packit Service 2781ba
                                       cl_data->kernel[1], 2,
Packit Service 2781ba
                                       NULL, global_ws, NULL,
Packit Service 2781ba
                                       0, NULL, NULL);
Packit Service 2781ba
  if (cl_err != CL_SUCCESS) return cl_err;
Packit Service 2781ba
  }
Packit Service 2781ba
Packit Service 2781ba
  {
Packit Service 2781ba
  global_ws[0] = roi->width;
Packit Service 2781ba
  global_ws[1] = roi->height;
Packit Service 2781ba
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 0, sizeof(cl_mem), (void*)&aux_tex);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 1, sizeof(cl_int), (void*)&aux_rect->width);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 2, sizeof(cl_mem), (void*)&out_tex);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 3, sizeof(cl_mem), (void*)&cl_matrix_y);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 4, sizeof(cl_int), (void*)&matrix_length_y);
Packit Service 2781ba
  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 5, sizeof(cl_int), (void*)&yoff);
Packit Service 2781ba
  if (cl_err != CL_SUCCESS) return cl_err;
Packit Service 2781ba
Packit Service 2781ba
  cl_err = gegl_clEnqueueNDRangeKernel(gegl_cl_get_command_queue (),
Packit Service 2781ba
                                       cl_data->kernel[0], 2,
Packit Service 2781ba
                                       NULL, global_ws, NULL,
Packit Service 2781ba
                                       0, NULL, NULL);
Packit Service 2781ba
  if (cl_err != CL_SUCCESS) return cl_err;
Packit Service 2781ba
  }
Packit Service 2781ba
Packit Service 2781ba
  gegl_clFinish(gegl_cl_get_command_queue ());
Packit Service 2781ba
Packit Service 2781ba
  gegl_clReleaseMemObject(cl_matrix_x);
Packit Service 2781ba
  gegl_clReleaseMemObject(cl_matrix_y);
Packit Service 2781ba
Packit Service 2781ba
  return CL_SUCCESS;
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
static gboolean
Packit Service 2781ba
cl_process (GeglOperation       *operation,
Packit Service 2781ba
            GeglBuffer          *input,
Packit Service 2781ba
            GeglBuffer          *output,
Packit Service 2781ba
            const GeglRectangle *result)
Packit Service 2781ba
{
Packit Service 2781ba
  const Babl *in_format  = gegl_operation_get_format (operation, "input");
Packit Service 2781ba
  const Babl *out_format = gegl_operation_get_format (operation, "output");
Packit Service 2781ba
  gint err;
Packit Service 2781ba
  gint j;
Packit Service 2781ba
  cl_int cl_err;
Packit Service 2781ba
Packit Service 2781ba
  GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
Packit Service 2781ba
  GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);
Packit Service 2781ba
Packit Service 2781ba
  gdouble *cmatrix_x, *cmatrix_y;
Packit Service 2781ba
  gint cmatrix_len_x, cmatrix_len_y;
Packit Service 2781ba
Packit Service 2781ba
  gfloat *fmatrix_x, *fmatrix_y;
Packit Service 2781ba
Packit Service 2781ba
  cmatrix_len_x = fir_gen_convolve_matrix (o->std_dev_x, &cmatrix_x);
Packit Service 2781ba
  cmatrix_len_y = fir_gen_convolve_matrix (o->std_dev_y, &cmatrix_y);
Packit Service 2781ba
Packit Service 2781ba
  fmatrix_x = g_new (gfloat, cmatrix_len_x);
Packit Service 2781ba
  fmatrix_y = g_new (gfloat, cmatrix_len_y);
Packit Service 2781ba
Packit Service 2781ba
  for(j=0; j
Packit Service 2781ba
    fmatrix_x[j] = (gfloat) cmatrix_x[j];
Packit Service 2781ba
Packit Service 2781ba
  for(j=0; j
Packit Service 2781ba
    fmatrix_y[j] = (gfloat) cmatrix_y[j];
Packit Service 2781ba
Packit Service 2781ba
  {
Packit Service 2781ba
  GeglBufferClIterator *i = gegl_buffer_cl_iterator_new (output, result, out_format, GEGL_CL_BUFFER_WRITE, GEGL_ABYSS_NONE);
Packit Service 2781ba
  gint read = gegl_buffer_cl_iterator_add_2 (i, input, result, in_format, GEGL_CL_BUFFER_READ,
Packit Service 2781ba
                                             op_area->left, op_area->right, op_area->top, op_area->bottom, GEGL_ABYSS_NONE);
Packit Service 2781ba
  gint aux  = gegl_buffer_cl_iterator_add_2 (i, NULL, result, in_format,  GEGL_CL_BUFFER_AUX,
Packit Service 2781ba
                                             0, 0, op_area->top, op_area->bottom, GEGL_ABYSS_NONE);
Packit Service 2781ba
Packit Service 2781ba
  while (gegl_buffer_cl_iterator_next (i, &err))
Packit Service 2781ba
    {
Packit Service 2781ba
      if (err) return FALSE;
Packit Service 2781ba
      for (j=0; j < i->n; j++)
Packit Service 2781ba
        {
Packit Service 2781ba
           cl_err = cl_gaussian_blur(i->tex[read][j],
Packit Service 2781ba
                                     i->tex[0][j],
Packit Service 2781ba
                                     i->tex[aux][j],
Packit Service 2781ba
                                     i->size[0][j],
Packit Service 2781ba
                                     &i->roi[0][j],
Packit Service 2781ba
                                     &i->roi[read][j],
Packit Service 2781ba
                                     &i->roi[aux][j],
Packit Service 2781ba
                                     fmatrix_x,
Packit Service 2781ba
                                     cmatrix_len_x,
Packit Service 2781ba
                                     op_area->left,
Packit Service 2781ba
                                     fmatrix_y,
Packit Service 2781ba
                                     cmatrix_len_y,
Packit Service 2781ba
                                     op_area->top);
Packit Service 2781ba
          if (cl_err != CL_SUCCESS)
Packit Service 2781ba
            {
Packit Service 2781ba
              g_warning("[OpenCL] Error in gegl:gaussian-blur: %s", gegl_cl_errstring(cl_err));
Packit Service 2781ba
              return FALSE;
Packit Service 2781ba
            }
Packit Service 2781ba
        }
Packit Service 2781ba
    }
Packit Service 2781ba
  }
Packit Service 2781ba
Packit Service 2781ba
  g_free (fmatrix_x);
Packit Service 2781ba
  g_free (fmatrix_y);
Packit Service 2781ba
Packit Service 2781ba
  g_free (cmatrix_x);
Packit Service 2781ba
  g_free (cmatrix_y);
Packit Service 2781ba
  return TRUE;
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
static gboolean
Packit Service 2781ba
process (GeglOperation       *operation,
Packit Service 2781ba
         GeglBuffer          *input,
Packit Service 2781ba
         GeglBuffer          *output,
Packit Service 2781ba
         const GeglRectangle *result,
Packit Service 2781ba
         gint                 level)
Packit Service 2781ba
{
Packit Service 2781ba
  GeglRectangle rect;
Packit Service 2781ba
  GeglBuffer *temp;
Packit Service 2781ba
  GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
Packit Service 2781ba
  GeglChantO              *o       = GEGL_CHANT_PROPERTIES (operation);
Packit Service 2781ba
Packit Service 2781ba
  GeglRectangle temp_extend;
Packit Service 2781ba
  gdouble       B, b[4];
Packit Service 2781ba
  gdouble      *cmatrix;
Packit Service 2781ba
  gint          cmatrix_len;
Packit Service 2781ba
  gboolean      force_iir;
Packit Service 2781ba
  gboolean      force_fir;
Packit Service 2781ba
Packit Service 2781ba
  rect.x      = result->x - op_area->left;
Packit Service 2781ba
  rect.width  = result->width + op_area->left + op_area->right;
Packit Service 2781ba
  rect.y      = result->y - op_area->top;
Packit Service 2781ba
  rect.height = result->height + op_area->top + op_area->bottom;
Packit Service 2781ba
Packit Service 2781ba
  force_iir = o->filter && !strcmp (o->filter, "iir");
Packit Service 2781ba
  force_fir = o->filter && !strcmp (o->filter, "fir");
Packit Service 2781ba
Packit Service 2781ba
  if (gegl_cl_is_accelerated () && !force_iir)
Packit Service 2781ba
    if (cl_process(operation, input, output, result))
Packit Service 2781ba
      return TRUE;
Packit Service 2781ba
Packit Service 2781ba
  temp_extend = rect;
Packit Service 2781ba
  temp_extend.x      = result->x;
Packit Service 2781ba
  temp_extend.width  = result->width;
Packit Service 2781ba
  temp = gegl_buffer_new (&temp_extend, babl_format ("RaGaBaA float"));
Packit Service 2781ba
Packit Service 2781ba
  if ((force_iir || o->std_dev_x > 1.0) && !force_fir)
Packit Service 2781ba
    {
Packit Service 2781ba
      iir_young_find_constants (o->std_dev_x, &B, b);
Packit Service 2781ba
      iir_young_hor_blur (input, &rect, temp, &temp_extend,  B, b);
Packit Service 2781ba
    }
Packit Service 2781ba
  else
Packit Service 2781ba
    {
Packit Service 2781ba
      cmatrix_len = fir_gen_convolve_matrix (o->std_dev_x, &cmatrix);
Packit Service 2781ba
      fir_hor_blur (input, &rect, temp, &temp_extend,
Packit Service 2781ba
                    cmatrix, cmatrix_len, op_area->left);
Packit Service 2781ba
      g_free (cmatrix);
Packit Service 2781ba
    }
Packit Service 2781ba
Packit Service 2781ba
Packit Service 2781ba
  if ((force_iir || o->std_dev_y > 1.0) && !force_fir)
Packit Service 2781ba
    {
Packit Service 2781ba
      iir_young_find_constants (o->std_dev_y, &B, b);
Packit Service 2781ba
      iir_young_ver_blur (temp, &temp_extend, output, result, B, b);
Packit Service 2781ba
    }
Packit Service 2781ba
  else
Packit Service 2781ba
    {
Packit Service 2781ba
      cmatrix_len = fir_gen_convolve_matrix (o->std_dev_y, &cmatrix);
Packit Service 2781ba
      fir_ver_blur (temp, &temp_extend, output, result, cmatrix, cmatrix_len,
Packit Service 2781ba
                    op_area->top);
Packit Service 2781ba
      g_free (cmatrix);
Packit Service 2781ba
    }
Packit Service 2781ba
Packit Service 2781ba
  g_object_unref (temp);
Packit Service 2781ba
  return  TRUE;
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
Packit Service 2781ba
static void
Packit Service 2781ba
gegl_chant_class_init (GeglChantClass *klass)
Packit Service 2781ba
{
Packit Service 2781ba
  GeglOperationClass       *operation_class;
Packit Service 2781ba
  GeglOperationFilterClass *filter_class;
Packit Service 2781ba
Packit Service 2781ba
  operation_class = GEGL_OPERATION_CLASS (klass);
Packit Service 2781ba
  filter_class    = GEGL_OPERATION_FILTER_CLASS (klass);
Packit Service 2781ba
Packit Service 2781ba
  filter_class->process    = process;
Packit Service 2781ba
  operation_class->prepare = prepare;
Packit Service 2781ba
Packit Service 2781ba
  operation_class->opencl_support = TRUE;
Packit Service 2781ba
Packit Service 2781ba
  gegl_operation_class_set_keys (operation_class,
Packit Service 2781ba
    "name",       "gegl:gaussian-blur",
Packit Service 2781ba
    "categories", "blur",
Packit Service 2781ba
    "description",
Packit Service 2781ba
        _("Performs an averaging of neighboring pixels with the "
Packit Service 2781ba
          "normal distribution as weighting"),
Packit Service 2781ba
        NULL);
Packit Service 2781ba
}
Packit Service 2781ba
Packit Service 2781ba
#endif