/* 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 2005 Øyvind Kolås <pippin@gimp.org>,

 *           2007 Øyvind Kolås <oeyvindk@hig.no>

 */

#include "config.h"

#include <glib/gi18n-lib.h>

#ifdef GEGL_CHANT_PROPERTIES

gegl_chant_double_ui (blur_radius, _("Blur radius"), 0.0, 1000.0, 4.0, 0.0, 100.0, 1.5,

  _("Radius of square pixel region, (width and height will be radius*2+1)."))

gegl_chant_double (edge_preservation, _("Edge preservation"), 0.0, 100.0, 8.0,

  _("Amount of edge preservation"))

#else

#define GEGL_CHANT_TYPE_AREA_FILTER

#define GEGL_CHANT_C_FILE       "bilateral-filter.c"

#include "gegl-chant.h"

#include <math.h>

static void

bilateral_filter (GeglBuffer          *src,

                  const GeglRectangle *src_rect,

                  GeglBuffer          *dst,

                  const GeglRectangle *dst_rect,

                  gdouble              radius,

                  gdouble              preserve);

#include <stdio.h>

static void prepare (GeglOperation *operation)

{

  GeglOperationAreaFilter *area = GEGL_OPERATION_AREA_FILTER (operation);

  GeglChantO              *o = GEGL_CHANT_PROPERTIES (operation);

  area->left = area->right = area->top = area->bottom = ceil (o->blur_radius);

  gegl_operation_set_format (operation, "input", babl_format ("RGBA float"));

  gegl_operation_set_format (operation, "output", babl_format ("RGBA float"));

}

#include "opencl/gegl-cl.h"

#include "buffer/gegl-buffer-cl-iterator.h"

static const char* kernel_source =

"#define POW2(a) ((a) * (a))                                           \n"

"kernel void bilateral_filter(global float4 *in,                       \n"

"                             global float4 *out,                      \n"

"                             const  float radius,                     \n"

"                             const  float preserve)                   \n"

"{                                                                     \n"

"    int gidx       = get_global_id(0);                                \n"

"    int gidy       = get_global_id(1);                                \n"

"    int n_radius   = ceil(radius);                                    \n"

"    int dst_width  = get_global_size(0);                              \n"

"    int src_width  = dst_width + n_radius * 2;                        \n"

"                                                                      \n"

"    int u, v, i, j;                                                   \n"

"    float4 center_pix =                                               \n"

"        in[(gidy + n_radius) * src_width + gidx + n_radius];          \n"

"    float4 accumulated = 0.0f;                                        \n"

"    float4 tempf       = 0.0f;                                        \n"

"    float  count       = 0.0f;                                        \n"

"    float  diff_map, gaussian_weight, weight;                         \n"

"                                                                      \n"

"    for (v = -n_radius;v <= n_radius; ++v)                            \n"

"    {                                                                 \n"

"        for (u = -n_radius;u <= n_radius; ++u)                        \n"

"        {                                                             \n"

"            i = gidx + n_radius + u;                                  \n"

"            j = gidy + n_radius + v;                                  \n"

"                                                                      \n"

"            int gid1d = i + j * src_width;                            \n"

"            tempf = in[gid1d];                                        \n"

"                                                                      \n"

"            diff_map = exp (                                          \n"

"                - (   POW2(center_pix.x - tempf.x)                    \n"

"                    + POW2(center_pix.y - tempf.y)                    \n"

"                    + POW2(center_pix.z - tempf.z))                   \n"

"                * preserve);                                          \n"

"                                                                      \n"

"            gaussian_weight =                                         \n"

"                exp( - 0.5f * (POW2(u) + POW2(v)) / radius);          \n"

"                                                                      \n"

"            weight = diff_map * gaussian_weight;                      \n"

"                                                                      \n"

"            accumulated += tempf * weight;                            \n"

"            count += weight;                                          \n"

"        }                                                             \n"

"    }                                                                 \n"

"    out[gidx + gidy * dst_width] = accumulated / count;               \n"

"}                                                                     \n";

static gegl_cl_run_data *cl_data = NULL;

static cl_int

cl_bilateral_filter (cl_mem                in_tex,

                     cl_mem                out_tex,

                     size_t                global_worksize,

                     const GeglRectangle  *roi,

                     gfloat                radius,

                     gfloat                preserve)

{

  cl_int cl_err = 0;

  size_t global_ws[2];

  if (!cl_data)

  {

    const char *kernel_name[] = {"bilateral_filter", NULL};

    cl_data = gegl_cl_compile_and_build (kernel_source, kernel_name);

  }

  if (!cl_data) return 1;

  global_ws[0] = roi->width;

  global_ws[1] = roi->height;

  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 0, sizeof(cl_mem),   (void*)&in_tex);

  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 1, sizeof(cl_mem),   (void*)&out_tex);

  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 2, sizeof(cl_float), (void*)&radius);

  cl_err |= gegl_clSetKernelArg(cl_data->kernel[0], 3, sizeof(cl_float), (void*)&preserve);

  if (cl_err != CL_SUCCESS) return cl_err;

  cl_err = gegl_clEnqueueNDRangeKernel(gegl_cl_get_command_queue (),

                                       cl_data->kernel[0], 2,

                                       NULL, global_ws, NULL,

                                       0, NULL, NULL);

  if (cl_err != CL_SUCCESS) return cl_err;

  return cl_err;

}

static gboolean

cl_process (GeglOperation       *operation,

            GeglBuffer          *input,

            GeglBuffer          *output,

            const GeglRectangle *result)

{

  const Babl *in_format  = gegl_operation_get_format (operation, "input");

  const Babl *out_format = gegl_operation_get_format (operation, "output");

  gint err;

  gint j;

  cl_int cl_err;

  GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);

  GeglChantO *o = GEGL_CHANT_PROPERTIES (operation);

  GeglBufferClIterator *i = gegl_buffer_cl_iterator_new (output,   result, out_format, GEGL_CL_BUFFER_WRITE, GEGL_ABYSS_NONE);

                gint read = gegl_buffer_cl_iterator_add_2 (i, input, result, in_format, GEGL_CL_BUFFER_READ, op_area->left, op_area->right, op_area->top, op_area->bottom, GEGL_ABYSS_NONE);

  while (gegl_buffer_cl_iterator_next (i, &err))

  {

    if (err) return FALSE;

    for (j=0; j < i->n; j++)

    {

      cl_err = cl_bilateral_filter(i->tex[read][j], i->tex[0][j], i->size[0][j], &i->roi[0][j], ceil(o->blur_radius), o->edge_preservation);

      if (cl_err != CL_SUCCESS)

      {

        g_warning("[OpenCL] Error in gegl:bilateral-filter: %s", gegl_cl_errstring(cl_err));

        return FALSE;

      }

    }

  }

  return TRUE;

}

static gboolean

process (GeglOperation       *operation,

         GeglBuffer          *input,

         GeglBuffer          *output,

         const GeglRectangle *result,

         gint                 level)

{

  GeglChantO   *o = GEGL_CHANT_PROPERTIES (operation);

  GeglRectangle compute;

  if (o->blur_radius >= 1.0 && gegl_cl_is_accelerated ())

    if (cl_process (operation, input, output, result))

      return TRUE;

  compute = gegl_operation_get_required_for_output (operation, "input",result);

  if (o->blur_radius < 1.0)

    {

      output = g_object_ref (input);

    }

  else

    {

      bilateral_filter (input, &compute, output, result, o->blur_radius, o->edge_preservation);

    }

  return  TRUE;

}

static void

bilateral_filter (GeglBuffer          *src,

                  const GeglRectangle *src_rect,

                  GeglBuffer          *dst,

                  const GeglRectangle *dst_rect,

                  gdouble              radius,

                  gdouble              preserve)

{

  gfloat *gauss;

  gint x,y;

  gint offset;

  gfloat *src_buf;

  gfloat *dst_buf;

  gint width = (gint) radius * 2 + 1;

  gint iradius = radius;

  gint src_width = src_rect->width;

  gint src_height = src_rect->height;

  gauss = g_newa (gfloat, width * width);

  src_buf = g_new0 (gfloat, src_rect->width * src_rect->height * 4);

  dst_buf = g_new0 (gfloat, dst_rect->width * dst_rect->height * 4);

  gegl_buffer_get (src, src_rect, 1.0, babl_format ("RGBA float"), src_buf, GEGL_AUTO_ROWSTRIDE,

                   GEGL_ABYSS_NONE);

  offset = 0;

#define POW2(a) ((a)*(a))

  for (y=-iradius;y<=iradius;y++)

    for (x=-iradius;x<=iradius;x++)

      {

        gauss[x+(int)radius + (y+(int)radius)*width] = exp(- 0.5*(POW2(x)+POW2(y))/radius   );

      }

  for (y=0; y<dst_rect->height; y++)

    for (x=0; x<dst_rect->width; x++)

      {

        gint u,v;

        gfloat *center_pix = src_buf + ((x+iradius)+((y+iradius) * src_width)) * 4;

        gfloat  accumulated[4]={0,0,0,0};

        gfloat  count=0.0;

        for (v=-iradius;v<=iradius;v++)

          for (u=-iradius;u<=iradius;u++)

            {

              gint i,j;

              i = x + radius + u;

              j = y + radius + v;

              if (i >= 0 && i < src_width &&

                  j >= 0 && j < src_height)

                {

                  gint c;

                  gfloat *src_pix = src_buf + (i + j * src_width) * 4;

                  gfloat diff_map   = exp (- (POW2(center_pix[0] - src_pix[0])+

                                              POW2(center_pix[1] - src_pix[1])+

                                              POW2(center_pix[2] - src_pix[2])) * preserve

                                          );

                  gfloat gaussian_weight;

                  gfloat weight;

                  gaussian_weight = gauss[u+(int)radius+(v+(int)radius)*width];

                  weight = diff_map * gaussian_weight;

                  for (c=0;c<4;c++)

                    {

                      accumulated[c] += src_pix[c] * weight;

                    }

                  count += weight;

                }

            }

        for (u=0; u<4;u++)

          dst_buf[offset*4+u] = accumulated[u]/count;

        offset++;

      }

  gegl_buffer_set (dst, dst_rect, 0, babl_format ("RGBA float"), dst_buf,

                   GEGL_AUTO_ROWSTRIDE);

  g_free (src_buf);

  g_free (dst_buf);

}

static void

gegl_chant_class_init (GeglChantClass *klass)

{

  GeglOperationClass       *operation_class;

  GeglOperationFilterClass *filter_class;

  operation_class  = GEGL_OPERATION_CLASS (klass);

  filter_class     = GEGL_OPERATION_FILTER_CLASS (klass);

  filter_class->process   = process;

  operation_class->prepare = prepare;

  operation_class->opencl_support = TRUE;

  gegl_operation_class_set_keys (operation_class,

           "name", "gegl:bilateral-filter",

           "categories", "misc",

           "description",

           _("An edge preserving blur filter that can be used for noise reduction. "

          "It is a gaussian blur where the contribution of neighbourhood pixels "

          "are weighted by the color difference from the center pixel."),

           NULL);

}

#endif