operations/common/bilateral-filter.c

   1 /* This file is an image processing operation for GEGL
   2  *
   3  * GEGL is free software; you can redistribute it and/or
   4  * modify it under the terms of the GNU Lesser General Public
   5  * License as published by the Free Software Foundation; either
   6  * version 3 of the License, or (at your option) any later version.
   7  *
   8  * GEGL is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11  * Lesser General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU Lesser General Public
  14  * License along with GEGL; if not, see <https://www.gnu.org/licenses/>.
  15  *
  16  * Copyright 2005 Øyvind Kolås <pippin@gimp.org>,
  17  *           2007 Øyvind Kolås <oeyvindk@hig.no>
  18  */
  19
  20 #include "config.h"
  21 #include <glib/gi18n-lib.h>
  22
  23
  24 #ifdef GEGL_PROPERTIES
  25
  26 property_double (blur_radius, _("Blur radius"), 4.0)
  27   description(_("Radius of square pixel region, (width and height will be radius*2+1)."))
  28   value_range   (0.0, 1000.0)
  29   ui_range      (0.0, 100.0)
  30   ui_gamma      (1.5)
  31
  32 property_double (edge_preservation, _("Edge preservation"), 8.0)
  33   description   (_("Amount of edge preservation"))
  34   value_range   (0.0, 100.0)
  35
  36 #else
  37
  38 #define GEGL_OP_AREA_FILTER
  39 #define GEGL_OP_NAME     bilateral_filter
  40 #define GEGL_OP_C_SOURCE bilateral-filter.c
  41
  42 #include "gegl-op.h"
  43
  44 static void
  45 bilateral_filter (GeglBuffer          *src,
  46                   const GeglRectangle *src_rect,
  47                   GeglBuffer          *dst,
  48                   const GeglRectangle *dst_rect,
  49                   gdouble              radius,
  50                   gdouble              preserve,
  51                   const Babl          *format);
  52
  53 #include <stdio.h>
  54
  55 static void prepare (GeglOperation *operation)
  56 {
  57   const Babl *space = gegl_operation_get_source_space (operation, "input");
  58   const Babl *format = babl_format_with_space ("RGBA float", space);
  59   GeglOperationAreaFilter *area = GEGL_OPERATION_AREA_FILTER (operation);
  60   GeglProperties              *o = GEGL_PROPERTIES (operation);
  61
  62   area->left = area->right = area->top = area->bottom = ceil (o->blur_radius);
  63   gegl_operation_set_format (operation, "input", format);
  64   gegl_operation_set_format (operation, "output", format);
  65 }
  66
  67 #include "opencl/gegl-cl.h"
  68 #include "gegl-buffer-cl-iterator.h"
  69
  70 #include "opencl/bilateral-filter.cl.h"
  71
  72 static GeglClRunData *cl_data = NULL;
  73
  74 static gboolean
  75 cl_bilateral_filter (cl_mem                in_tex,
  76                      cl_mem                out_tex,
  77                      size_t                global_worksize,
  78                      const GeglRectangle  *roi,
  79                      gfloat                radius,
  80                      gfloat                preserve)
  81 {
  82   cl_int cl_err = 0;
  83   size_t global_ws[2];
  84
  85   if (!cl_data)
  86     {
  87       const char *kernel_name[] = {"bilateral_filter", NULL};
  88       cl_data = gegl_cl_compile_and_build (bilateral_filter_cl_source, kernel_name);
  89     }
  90   if (!cl_data) return TRUE;
  91
  92   global_ws[0] = roi->width;
  93   global_ws[1] = roi->height;
  94
  95   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 0, sizeof(cl_mem),   (void*)&in_tex);
  96   CL_CHECK;
  97   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 1, sizeof(cl_mem),   (void*)&out_tex);
  98   CL_CHECK;
  99   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 2, sizeof(cl_float), (void*)&radius);
 100   CL_CHECK;
 101   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 3, sizeof(cl_float), (void*)&preserve);
 102   CL_CHECK;
 103
 104   cl_err = gegl_clEnqueueNDRangeKernel(gegl_cl_get_command_queue (),
 105                                        cl_data->kernel[0], 2,
 106                                        NULL, global_ws, NULL,
 107                                        0, NULL, NULL);
 108   CL_CHECK;
 109
 110   return FALSE;
 111
 112 error:
 113   return TRUE;
 114 }
 115
 116 static gboolean
 117 cl_process (GeglOperation       *operation,
 118             GeglBuffer          *input,
 119             GeglBuffer          *output,
 120             const GeglRectangle *result)
 121 {
 122   const Babl *in_format  = gegl_operation_get_format (operation, "input");
 123   const Babl *out_format = gegl_operation_get_format (operation, "output");
 124   gint err;
 125
 126   GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
 127   GeglProperties *o = GEGL_PROPERTIES (operation);
 128
 129   GeglBufferClIterator *i = gegl_buffer_cl_iterator_new (output,
 130                                                          result,
 131                                                          out_format,
 132                                                          GEGL_CL_BUFFER_WRITE);
 133
 134   gint read = gegl_buffer_cl_iterator_add_2 (i,
 135                                              input,
 136                                              result,
 137                                              in_format,
 138                                              GEGL_CL_BUFFER_READ,
 139                                              op_area->left,
 140                                              op_area->right,
 141                                              op_area->top,
 142                                              op_area->bottom,
 143                                              GEGL_ABYSS_NONE);
 144
 145   while (gegl_buffer_cl_iterator_next (i, &err))
 146     {
 147       if (err) return FALSE;
 148
 149       err = cl_bilateral_filter(i->tex[read],
 150                                 i->tex[0],
 151                                 i->size[0],
 152                                 &i->roi[0],
 153                                 ceil(o->blur_radius),
 154                                 o->edge_preservation);
 155
 156       if (err) return FALSE;
 157     }
 158
 159   return TRUE;
 160 }
 161
 162 static gboolean
 163 process (GeglOperation       *operation,
 164          GeglBuffer          *input,
 165          GeglBuffer          *output,
 166          const GeglRectangle *result,
 167          gint                 level)
 168 {
 169   GeglProperties *o = GEGL_PROPERTIES (operation);
 170   GeglRectangle compute;
 171   const Babl *format = gegl_operation_get_format (operation, "output");
 172
 173   if (o->blur_radius >= 1.0 && gegl_operation_use_opencl (operation))
 174     if (cl_process (operation, input, output, result))
 175       return TRUE;
 176
 177   compute = gegl_operation_get_required_for_output (operation, "input",result);
 178
 179   if (o->blur_radius < 1.0)
 180     {
 181       gegl_buffer_copy (input, result, GEGL_ABYSS_NONE,
 182                         output, result);
 183     }
 184   else
 185     {
 186       bilateral_filter (input, &compute, output, result, o->blur_radius, o->edge_preservation, format);
 187     }
 188
 189   return  TRUE;
 190 }
 191
 192 static void
 193 bilateral_filter (GeglBuffer          *src,
 194                   const GeglRectangle *src_rect,
 195                   GeglBuffer          *dst,
 196                   const GeglRectangle *dst_rect,
 197                   gdouble              radius,
 198                   gdouble              preserve,
 199                   const Babl          *format)
 200 {
 201   gfloat *gauss;
 202   gint x,y;
 203   gint offset;
 204   gfloat *src_buf;
 205   gfloat *dst_buf;
 206   gint width = (gint) radius * 2 + 1;
 207   gint iradius = radius;
 208   gint src_width = src_rect->width;
 209   gint src_height = src_rect->height;
 210
 211   gauss = g_newa (gfloat, width * width);
 212   src_buf = g_new0 (gfloat, src_rect->width * src_rect->height * 4);
 213   dst_buf = g_new0 (gfloat, dst_rect->width * dst_rect->height * 4);
 214
 215   gegl_buffer_get (src, src_rect, 1.0, format, src_buf, GEGL_AUTO_ROWSTRIDE,
 216                    GEGL_ABYSS_NONE);
 217
 218   offset = 0;
 219
 220 #define POW2(a) ((a)*(a))
 221   for (y=-iradius;y<=iradius;y++)
 222     for (x=-iradius;x<=iradius;x++)
 223       {
 224         gauss[x+(int)radius + (y+(int)radius)*width] = exp(- 0.5*(POW2(x)+POW2(y))/radius   );
 225       }
 226
 227   for (y=0; y<dst_rect->height; y++)
 228     for (x=0; x<dst_rect->width; x++)
 229       {
 230         gint u,v;
 231         gfloat *center_pix = src_buf + ((x+iradius)+((y+iradius) * src_width)) * 4;
 232         gfloat  accumulated[4]={0,0,0,0};
 233         gfloat  count=0.0;
 234
 235         for (v=-iradius;v<=iradius;v++)
 236           for (u=-iradius;u<=iradius;u++)
 237             {
 238               gint i,j;
 239               i = x + radius + u;
 240               j = y + radius + v;
 241               if (i >= 0 && i < src_width &&
 242                   j >= 0 && j < src_height)
 243                 {
 244                   gint c;
 245
 246                   gfloat *src_pix = src_buf + (i + j * src_width) * 4;
 247
 248                   gfloat diff_map   = exp (- (POW2(center_pix[0] - src_pix[0])+
 249                                               POW2(center_pix[1] - src_pix[1])+
 250                                               POW2(center_pix[2] - src_pix[2])) * preserve
 251                                           );
 252                   gfloat gaussian_weight;
 253                   gfloat weight;
 254
 255                   gaussian_weight = gauss[u+(int)radius+(v+(int)radius)*width];
 256
 257                   weight = diff_map * gaussian_weight;
 258
 259                   for (c=0;c<4;c++)
 260                     {
 261                       accumulated[c] += src_pix[c] * weight;
 262                     }
 263                   count += weight;
 264                 }
 265             }
 266
 267         for (u=0; u<4;u++)
 268           dst_buf[offset*4+u] = accumulated[u]/count;
 269         offset++;
 270       }
 271   gegl_buffer_set (dst, dst_rect, 0, format, dst_buf,
 272                    GEGL_AUTO_ROWSTRIDE);
 273   g_free (src_buf);
 274   g_free (dst_buf);
 275 }
 276
 277
 278 static void
 279 gegl_op_class_init (GeglOpClass *klass)
 280 {
 281   GeglOperationClass       *operation_class;
 282   GeglOperationFilterClass *filter_class;
 283
 284   operation_class  = GEGL_OPERATION_CLASS (klass);
 285   filter_class     = GEGL_OPERATION_FILTER_CLASS (klass);
 286
 287   filter_class->process   = process;
 288   operation_class->prepare = prepare;
 289
 290   operation_class->opencl_support = TRUE;
 291
 292   gegl_operation_class_set_keys (operation_class,
 293            "name", "gegl:bilateral-filter",
 294            "title", _("Bilateral Filter"),
 295            "categories", "enhance:noise-reduction",
 296            "reference-hash", "e19155c5d467143d11f259cc25fdec80",
 297            "description",
 298            _("Like a gaussian blur; but where the contribution for each neighborhood "
 299           "pixel is also weighted by the color difference with the original center pixel."),
 300            NULL);
 301 }
 302
 303 #endif