noise-reduction.c 14 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
/* 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/>.
 *
 * 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 28 29
property_int  (iterations, _("Strength"), 4)
  description (_("Controls the number of iterations; lower values give less plastic results"))
  value_range (1, 32)
  ui_range    (1, 8)
30 31 32

#else

33 34
#define GEGL_OP_AREA_FILTER
#define GEGL_OP_C_FILE       "noise-reduction.c"
35

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

39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
/* 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)
57

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

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

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

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
#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;
        }
    }
}
144 145 146 147

static void prepare (GeglOperation *operation)
{
  GeglOperationAreaFilter *area = GEGL_OPERATION_AREA_FILTER (operation);
148
  GeglProperties              *o = GEGL_PROPERTIES (operation);
149 150 151 152 153 154

  area->left = area->right = area->top = area->bottom = o->iterations;
  gegl_operation_set_format (operation, "input",  babl_format ("R'G'B'A float"));
  gegl_operation_set_format (operation, "output", babl_format ("R'G'B'A float"));
}

155 156 157
#include "opencl/gegl-cl.h"
#include "buffer/gegl-buffer-cl-iterator.h"

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

160
static GeglClRunData *cl_data = NULL;
161

Victor Oliveira's avatar
Victor Oliveira committed
162
static gboolean
163 164 165 166 167 168 169 170 171 172 173 174 175 176
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;

177 178
  cl_mem temp_tex = NULL;
  cl_mem tmptex   = NULL;
179 180 181 182 183 184

  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
185
      cl_data = gegl_cl_compile_and_build(noise_reduction_cl_source, kernel_name);
186
    }
Victor Oliveira's avatar
Victor Oliveira committed
187
  if (!cl_data)  return TRUE;
188 189 190 191 192

  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
193
  CL_CHECK;
194 195 196 197 198

  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
199
  CL_CHECK;
200 201 202 203 204 205 206 207 208 209 210 211 212

  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
213 214 215 216 217 218 219 220
      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;
221 222 223 224

      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
225
      CL_CHECK;
226 227 228 229 230
    }

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

Victor Oliveira's avatar
Victor Oliveira committed
231 232 233 234 235 236
  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;
237 238 239 240

  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
241
  CL_CHECK;
242 243

  cl_err = gegl_clFinish(gegl_cl_get_command_queue());
Victor Oliveira's avatar
Victor Oliveira committed
244
  CL_CHECK;
245

246 247 248 249 250
  if (tmptex)
    {
      cl_err = gegl_clReleaseMemObject (tmptex);
      CL_CHECK_ONLY (cl_err);
    }
251

Victor Oliveira's avatar
Victor Oliveira committed
252 253 254
  return FALSE;

error:
255 256 257
  if (tmptex)
    gegl_clReleaseMemObject (tmptex);

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

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

  GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
272
  GeglProperties *o = GEGL_PROPERTIES (operation);
273

274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
  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);

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

298
  while (gegl_buffer_cl_iterator_next (i, &err) && !err)
299
    {
300 301 302 303 304 305 306 307 308 309 310 311 312
      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;
        }
313 314
    }

315
  return !err;
316 317
}

318 319
#define INPLACE 1

320 321 322 323
static gboolean
process (GeglOperation       *operation,
         GeglBuffer          *input,
         GeglBuffer          *output,
324 325
         const GeglRectangle *result,
         gint                 level)
326
{
327
  GeglProperties   *o = GEGL_PROPERTIES (operation);
328

329
  int iteration;
330
  int stride;
331
  float *src_buf;
332
#ifndef INPLACE
333
  float *dst_buf;
334
#endif
335
  GeglRectangle rect;
336

337
  if (gegl_operation_use_opencl (operation))
338 339 340
    if(cl_process(operation, input, output, result))
      return TRUE;

341
  rect = *result;
342

343
  stride = result->width + o->iterations * 2;
344

345 346
  src_buf = g_new0 (float,
         (stride) * (result->height + o->iterations * 2) * 4);
347
#ifndef INPLACE
348 349
  dst_buf = g_new0 (float,
         (stride) * (result->height + o->iterations * 2) * 4);
350
#endif
351

352 353 354 355 356
  {
    rect.x      -= o->iterations;
    rect.y      -= o->iterations;
    rect.width  += o->iterations*2;
    rect.height += o->iterations*2;
357
    gegl_buffer_get (input, &rect, 1.0, babl_format ("R'G'B'A float"),
358
                     src_buf, stride * 4 * 4, GEGL_ABYSS_NONE);
359
  }
360

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

381
  gegl_buffer_set (output, result, 0, babl_format ("R'G'B'A float"),
382
#ifndef INPLACE
383
                   src_buf,
384 385 386
#else
                   src_buf + ((stride +1) * 4) * o->iterations,
#endif
387 388
                   stride * 4 * 4);

389
  g_free (src_buf);
390
#ifndef INPLACE
391
  g_free (dst_buf);
392
#endif
393 394

  return  TRUE;
395 396
}

397

398 399 400 401 402 403 404 405 406 407 408
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;
}

409
static void
410
gegl_op_class_init (GeglOpClass *klass)
411 412 413 414 415 416 417 418 419
{
  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;
420 421
  operation_class->opencl_support = TRUE;

422
  operation_class->get_bounding_box = get_bounding_box;
423

424
  gegl_operation_class_set_keys (operation_class,
425
    "title",       _("Noise Reduction"),
426
    "name"       , "gegl:noise-reduction",
427
    "categories" , "enhance:noise-reduction",
428
    "description", _("Anisotropic smoothing operation"),
429
    NULL);
430 431 432
}

#endif