noise-reduction.c 15.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
/* 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
14
 * License along with GEGL; if not, see <https://www.gnu.org/licenses/>.
15 16 17
 *
 * Ali Alsam, Hans Jakob Rivertz, Øyvind Kolås (c) 2011
 */
18 19
#include "config.h"
#include <glib/gi18n-lib.h>
20

21 22 23
#include "config.h"
#include <glib/gi18n-lib.h>

24
#ifdef GEGL_PROPERTIES
25

26 27
property_int  (iterations, _("Strength"), 4)
  description (_("Controls the number of iterations; lower values give less plastic results"))
28 29
  value_range (0, 32)
  ui_range    (0, 8)
30 31 32

#else

33
#define GEGL_OP_AREA_FILTER
34
#define GEGL_OP_NAME     noise_reduction
35
#define GEGL_OP_C_SOURCE noise-reduction.c
36

37
#include "gegl-op.h"
38 39
#include <math.h>

40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
/* The core noise_reduction function, which is implemented as
 * portable C - this is the function where most cpu time goes
 */
static void
noise_reduction (float *src_buf,     /* source buffer, one pixel to the left
                                        and up from the starting pixel */
                 int    src_stride,  /* stridewidth of buffer in pixels */
                 float *dst_buf,     /* destination buffer */
                 int    dst_width,   /* width to render */
                 int    dst_height,  /* height to render */
                 int    dst_stride)  /* stride of target buffer */
{
  int c;
  int x,y;
  int dst_offset;

#define NEIGHBOURS 8
#define AXES       (NEIGHBOURS/2)
58

59
#define POW2(a) ((a)*(a))
60 61
/* core code/formulas to be tweaked for the tuning the implementation */
#define GEN_METRIC(before, center, after) \
62
                   POW2((center) * 2 - (before) - (after))
63 64

/* Condition used to bail diffusion from a direction */
65
#define BAIL_CONDITION(new,original) ((new) > (original))
66

67
#define SYMMETRY(a)  (NEIGHBOURS - (a) - 1) /* point-symmetric neighbour pixel */
68

69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
#define O(u,v) (((u)+((v) * src_stride)) * 4)
  int   offsets[NEIGHBOURS] = {  /* array of the relative distance i float
                                  * pointers to each of neighbours
                                  * in source buffer, allows quick referencing.
                                  */
              O( -1, -1), O(0, -1), O(1, -1),
              O( -1,  0),           O(1,  0),
              O( -1,  1), O(0, 1),  O(1,  1)};
#undef O

  dst_offset = 0;
  for (y=0; y<dst_height; y++)
    {
      float *center_pix = src_buf + ((y+1) * src_stride + 1) * 4;
      dst_offset = dst_stride * y;
      for (x=0; x<dst_width; x++)
        {
          for (c=0; c<3; c++) /* do each color component individually */
            {
              float  metric_reference[AXES];
              int    axis;
              int    direction;
              float  sum;
              int    count;

              for (axis = 0; axis < AXES; axis++)
                { /* initialize original metrics for the horizontal, vertical
                     and 2 diagonal metrics */
                  float *before_pix  = center_pix + offsets[axis];
                  float *after_pix   = center_pix + offsets[SYMMETRY(axis)];

                  metric_reference[axis] =
                    GEN_METRIC (before_pix[c], center_pix[c], after_pix[c]);
                }

              sum   = center_pix[c];
              count = 1;

              /* try smearing in data from all neighbours */
              for (direction = 0; direction < NEIGHBOURS; direction++)
                {
                  float *pix   = center_pix + offsets[direction];
                  float  value = pix[c] * 0.5 + center_pix[c] * 0.5;
                  int    axis;
                  int    valid;

                  /* check if the non-smoothing operating check is true if
                   * smearing from this direction for any of the axes */
                  valid = 1; /* assume it will be valid */
                  for (axis = 0; axis < AXES; axis++)
                    {
                      float *before_pix = center_pix + offsets[axis];
                      float *after_pix  = center_pix + offsets[SYMMETRY(axis)];
                      float  metric_new =
                             GEN_METRIC (before_pix[c], value, after_pix[c]);

                      if (BAIL_CONDITION(metric_new, metric_reference[axis]))
                        {
                          valid = 0; /* mark as not a valid smoothing, and .. */
                          break;     /* .. break out of loop */
                        }
                    }
                  if (valid) /* we were still smooth in all axes */
                    {        /* add up contribution to final result  */
                      sum += value;
                      count ++;
                    }
                }
              dst_buf[dst_offset*4+c] = sum / count;
            }
          dst_buf[dst_offset*4+3] = center_pix[3]; /* copy alpha unmodified */
          dst_offset++;
          center_pix += 4;
        }
    }
}
145 146 147

static void prepare (GeglOperation *operation)
{
148
  const Babl *space = gegl_operation_get_source_space (operation, "input");
149
  GeglOperationAreaFilter *area = GEGL_OPERATION_AREA_FILTER (operation);
150
  GeglProperties              *o = GEGL_PROPERTIES (operation);
151 152

  area->left = area->right = area->top = area->bottom = o->iterations;
153 154
  gegl_operation_set_format (operation, "input",  babl_format_with_space ("R~G~B~A float", space));
  gegl_operation_set_format (operation, "output", babl_format_with_space ("R~G~B~A float", space));
155 156
}

157
#include "opencl/gegl-cl.h"
158
#include "gegl-buffer-cl-iterator.h"
159

Victor Oliveira's avatar
Victor Oliveira committed
160
#include "opencl/noise-reduction.cl.h"
161

162
static GeglClRunData *cl_data = NULL;
163

Victor Oliveira's avatar
Victor Oliveira committed
164
static gboolean
165 166 167 168 169 170 171 172 173 174 175 176 177 178
cl_noise_reduction (cl_mem                in_tex,
                    cl_mem                aux_tex,
                    cl_mem                out_tex,
                    size_t                global_worksize,
                    const GeglRectangle  *src_roi,
                    const GeglRectangle  *roi,
                    const int             iterations)
{
  int i = 0;
  size_t gbl_size_tmp[2];

  cl_int n_src_stride  = roi->width + iterations * 2;
  cl_int cl_err = 0;

179 180
  cl_mem temp_tex = NULL;
  cl_mem tmptex   = NULL;
181 182 183 184 185 186

  gint stride = 16; /*R'G'B'A float*/

  if (!cl_data)
    {
      const char *kernel_name[] ={"noise_reduction_cl","transfer", NULL};
Victor Oliveira's avatar
Victor Oliveira committed
187
      cl_data = gegl_cl_compile_and_build(noise_reduction_cl_source, kernel_name);
188
    }
Victor Oliveira's avatar
Victor Oliveira committed
189
  if (!cl_data)  return TRUE;
190 191 192 193 194

  temp_tex = gegl_clCreateBuffer (gegl_cl_get_context(),
                                  CL_MEM_READ_WRITE,
                                  src_roi->width * src_roi->height * stride,
                                  NULL, &cl_err);
Victor Oliveira's avatar
Victor Oliveira committed
195
  CL_CHECK;
196 197 198 199 200

  cl_err = gegl_clEnqueueCopyBuffer(gegl_cl_get_command_queue(),
                                    in_tex , temp_tex , 0 , 0 ,
                                    src_roi->width * src_roi->height * stride,
                                    0, NULL, NULL);
Victor Oliveira's avatar
Victor Oliveira committed
201
  CL_CHECK;
202 203 204 205 206 207 208 209 210 211 212 213 214

  tmptex = temp_tex;
  for (i = 0;i<iterations;i++)
    {
      if (i > 0)
      {
        cl_mem temp = aux_tex;
        aux_tex     = temp_tex;
        temp_tex    = temp;
      }
      gbl_size_tmp[0] = roi->width  + 2 * (iterations - 1 -i);
      gbl_size_tmp[1] = roi->height + 2 * (iterations - 1 -i);

Victor Oliveira's avatar
Victor Oliveira committed
215 216 217 218 219 220 221 222
      cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 0, sizeof(cl_mem), (void*)&temp_tex);
      CL_CHECK;
      cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 1, sizeof(cl_int), (void*)&n_src_stride);
      CL_CHECK;
      cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 2, sizeof(cl_mem), (void*)&aux_tex);
      CL_CHECK;
      cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 3, sizeof(cl_int), (void*)&n_src_stride);
      CL_CHECK;
223 224 225 226

      cl_err = gegl_clEnqueueNDRangeKernel(gegl_cl_get_command_queue(), cl_data->kernel[0],
                                           2, NULL, gbl_size_tmp, NULL,
                                           0, NULL, NULL);
Victor Oliveira's avatar
Victor Oliveira committed
227
      CL_CHECK;
228 229 230 231 232
    }

  gbl_size_tmp[0] = roi->width ;
  gbl_size_tmp[1] = roi->height;

Victor Oliveira's avatar
Victor Oliveira committed
233 234 235 236 237 238
  cl_err = gegl_clSetKernelArg(cl_data->kernel[1], 0, sizeof(cl_mem), (void*)&aux_tex);
  CL_CHECK;
  cl_err = gegl_clSetKernelArg(cl_data->kernel[1], 1, sizeof(cl_int), (void*)&n_src_stride);
  CL_CHECK;
  cl_err = gegl_clSetKernelArg(cl_data->kernel[1], 2, sizeof(cl_mem), (void*)&out_tex);
  CL_CHECK;
239 240 241 242

  cl_err = gegl_clEnqueueNDRangeKernel(gegl_cl_get_command_queue(), cl_data->kernel[1],
                                       2, NULL, gbl_size_tmp, NULL,
                                       0, NULL, NULL);
Victor Oliveira's avatar
Victor Oliveira committed
243
  CL_CHECK;
244 245

  cl_err = gegl_clFinish(gegl_cl_get_command_queue());
Victor Oliveira's avatar
Victor Oliveira committed
246
  CL_CHECK;
247

248 249 250 251 252
  if (tmptex)
    {
      cl_err = gegl_clReleaseMemObject (tmptex);
      CL_CHECK_ONLY (cl_err);
    }
253

Victor Oliveira's avatar
Victor Oliveira committed
254 255 256
  return FALSE;

error:
257 258 259
  if (tmptex)
    gegl_clReleaseMemObject (tmptex);

Victor Oliveira's avatar
Victor Oliveira committed
260
  return TRUE;
261 262 263 264 265 266 267 268 269 270
}

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");
271
  gint err = 0;
272 273

  GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
274
  GeglProperties *o = GEGL_PROPERTIES (operation);
275

276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291
  GeglBufferClIterator *i = gegl_buffer_cl_iterator_new (output,
                                                         result,
                                                         out_format,
                                                         GEGL_CL_BUFFER_WRITE);

  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);

292 293 294 295 296 297 298
  gint aux = gegl_buffer_cl_iterator_add_aux (i,
                                              result,
                                              in_format,
                                              op_area->left,
                                              op_area->right,
                                              op_area->top,
                                              op_area->bottom);
299

300
  while (gegl_buffer_cl_iterator_next (i, &err) && !err)
301
    {
302 303 304 305 306 307 308 309 310 311 312 313 314
      err = cl_noise_reduction (i->tex[read],
                                i->tex[aux],
                                i->tex[0],
                                i->size[0],
                                &i->roi[read],
                                &i->roi[0],
                                o->iterations);

      if (err)
        {
          gegl_buffer_cl_iterator_stop (i);
          break;
        }
315 316
    }

317
  return !err;
318 319
}

320 321
#define INPLACE 1

322 323 324 325
static gboolean
process (GeglOperation       *operation,
         GeglBuffer          *input,
         GeglBuffer          *output,
326 327
         const GeglRectangle *result,
         gint                 level)
328
{
329
  GeglProperties   *o = GEGL_PROPERTIES (operation);
330 331
  const Babl *in_format  = gegl_operation_get_format (operation, "input");
  const Babl *out_format = gegl_operation_get_format (operation, "output");
332

333
  int iteration;
334
  int stride;
335
  float *src_buf;
336
#ifndef INPLACE
337
  float *dst_buf;
338
#endif
339
  GeglRectangle rect;
340

341
  if (gegl_operation_use_opencl (operation))
342 343 344
    if(cl_process(operation, input, output, result))
      return TRUE;

345
  rect = *result;
346

347
  stride = result->width + o->iterations * 2;
348

349 350
  src_buf = g_new0 (float,
         (stride) * (result->height + o->iterations * 2) * 4);
351
#ifndef INPLACE
352 353
  dst_buf = g_new0 (float,
         (stride) * (result->height + o->iterations * 2) * 4);
354
#endif
355

356 357 358 359 360
  {
    rect.x      -= o->iterations;
    rect.y      -= o->iterations;
    rect.width  += o->iterations*2;
    rect.height += o->iterations*2;
361
    gegl_buffer_get (input, &rect, 1.0, in_format,
362
                     src_buf, stride * 4 * 4, GEGL_ABYSS_NONE);
363
  }
364

365
  for (iteration = 0; iteration < o->iterations; iteration++)
366
    {
367
      noise_reduction (src_buf, stride,
368 369 370
#ifdef INPLACE
                       src_buf + (stride + 1) * 4,
#else
371
                       dst_buf,
372
#endif
373 374 375
                       result->width  + (o->iterations - 1 - iteration) * 2,
                       result->height + (o->iterations - 1 - iteration) * 2,
                       stride);
376
#ifndef INPLACE
377 378 379 380 381
      { /* swap buffers */
        float *tmp = src_buf;
        src_buf = dst_buf;
        dst_buf = tmp;
      }
382
#endif
383 384
    }

385
  gegl_buffer_set (output, result, 0, out_format,
386
#ifndef INPLACE
387
                   src_buf,
388 389 390
#else
                   src_buf + ((stride +1) * 4) * o->iterations,
#endif
391 392
                   stride * 4 * 4);

393
  g_free (src_buf);
394
#ifndef INPLACE
395
  g_free (dst_buf);
396
#endif
397 398

  return  TRUE;
399 400
}

401

402 403 404 405 406 407 408 409 410 411 412
static GeglRectangle
get_bounding_box (GeglOperation *operation)
{
  GeglRectangle  result = {0,0,0,0};
  GeglRectangle *in_rect = gegl_operation_source_get_bounding_box (operation,
                                                                     "input");
  if (!in_rect)
    return result;
  return *in_rect;
}

413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
/* Pass-through when iterations parameter is set to zero */

static gboolean
operation_process (GeglOperation        *operation,
                   GeglOperationContext *context,
                   const gchar          *output_prop,
                   const GeglRectangle  *result,
                   gint                  level)
{
  GeglOperationClass  *operation_class;
  GeglProperties      *o = GEGL_PROPERTIES (operation);

  operation_class = GEGL_OPERATION_CLASS (gegl_op_parent_class);

  if (! o->iterations)
    {
      gpointer in = gegl_operation_context_get_object (context, "input");
      gegl_operation_context_take_object (context, "output",
                                          g_object_ref (G_OBJECT (in)));
      return TRUE;
    }

  return operation_class->process (operation, context, output_prop, result,
                                   gegl_operation_context_get_level (context));
}

439
static void
440
gegl_op_class_init (GeglOpClass *klass)
441 442 443 444 445 446 447
{
  GeglOperationClass       *operation_class;
  GeglOperationFilterClass *filter_class;

  operation_class  = GEGL_OPERATION_CLASS (klass);
  filter_class     = GEGL_OPERATION_FILTER_CLASS (klass);

448 449 450
  filter_class->process           = process;
  operation_class->process        = operation_process;
  operation_class->prepare        = prepare;
451 452
  operation_class->opencl_support = TRUE;

453
  operation_class->get_bounding_box = get_bounding_box;
454

455
  gegl_operation_class_set_keys (operation_class,
456
    "title",       _("Noise Reduction"),
457
    "name"       , "gegl:noise-reduction",
458
    "categories" , "enhance:noise-reduction",
459
    "reference-hash", "7cd18da7f407f4bc6f917c8318e50b59",
460
    "description", _("Anisotropic smoothing operation"),
461
    NULL);
462 463 464
}

#endif