domain-transform.c 20.2 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 18 19 20 21 22 23
 *
 * Copyright 2018 Felipe Einsfeld Kersting <fekersting@inf.ufrgs.br>
 */

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

#ifdef GEGL_PROPERTIES

24 25 26
property_int (n_iterations, _("Quality"), 3)
  description(_("Number of filtering iterations. "
                "A value between 2 and 4 is usually enough."))
27 28 29
  value_range (1, 5)

property_double (spatial_factor, _("Blur radius"),  30.0)
30 31
  description  (_("Spatial standard deviation of the blur kernel, "
                  "measured in pixels."))
32 33 34
  value_range (0, 1000.0)

property_double (edge_preservation, _("Edge preservation"),  0.8)
35 36 37
  description  (_("Amount of edge preservation. This quantity is inversely "
                  "proportional to the range standard deviation of the blur "
                  "kernel."))
38 39 40 41 42 43 44 45 46 47 48
  value_range (0, 1.0)

#else

#define GEGL_OP_FILTER
#define GEGL_OP_NAME     domain_transform
#define GEGL_OP_C_SOURCE domain-transform.c

#include "gegl-op.h"
#include <math.h>

49 50 51 52 53 54
/* This should be 768, since we have 2^8 possible options for each channel.
 * Since domain transform is given by:
 * 1 + (s_s / s_r) * (diff_channel_R + diff_channel_G + diff_channel_B)
 * We will have 3 * 2^8 different possibilities for the transform of each pixel.
 * 3 * 2^8 = 768
 */
55 56 57 58
#define RF_TABLE_SIZE 768
#define SQRT3 1.7320508075f
#define SQRT2 1.4142135623f
#define BLOCK_STRIDE 1
59
#define REPORT_PROGRESS_TIME 0.5  /* time to report gegl_operation_progress */
60 61 62 63 64 65 66 67

static gint16
absolute (gint16 x)
{
  return (x < 0) ? -x : x;
}

static void
68 69 70
report_progress (GeglOperation  *operation,
                 gdouble         progress,
                 GTimer         *timer)
71 72 73 74 75 76
{
  static gboolean reported = FALSE;

  if (progress == 0.0)
    reported = FALSE;
    
77
  if (g_timer_elapsed (timer, NULL) >= REPORT_PROGRESS_TIME && !reported)
78 79
    {
      reported = TRUE;
80
      gegl_operation_progress (operation, 0.0, "");
81 82 83 84 85 86 87
    }

  if (reported)
    gegl_operation_progress (operation, progress, "");
}

static gint
88 89 90 91 92 93 94 95 96
domain_transform (GeglOperation  *operation,
                  gint            width,
                  gint            height,
                  gint            n_chan,
                  gfloat          spatial_factor,
                  gfloat          range_factor,
                  gint            n_iterations,
                  GeglBuffer     *input,
                  GeglBuffer     *output)
97
{
98 99 100
  const Babl *space    = gegl_operation_get_source_space (operation, "input");
  const Babl *formatu8 = babl_format_with_space ("R'G'B' u8", space);
  const Babl *format   = babl_format_with_space ("R'G'B'A float", space);
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
  gfloat  **rf_table;
  guint16  *transforms_buffer;
  gfloat   *buffer;

  gint16    last[3];
  gint16    current[3];
  gfloat    lastf[4];
  gfloat    currentf[4];
  gfloat    sum_channels_difference, a, w, sdt_dev;
  gint      i, j, k, n;
  gint      real_stride, d_x_position, d_y_position, biggest_dimension;
  guint16   d;
  GeglRectangle rect;
  GTimer  *timer;

  timer = g_timer_new ();
117 118

  /* PRE-ALLOC MEMORY */
119 120 121
  biggest_dimension = (width > height) ? width : height;
  buffer            = g_new (gfloat, BLOCK_STRIDE * biggest_dimension * n_chan);
  transforms_buffer = g_new (guint16, BLOCK_STRIDE * biggest_dimension);
122

123 124 125 126
  rf_table = g_new (gfloat *, n_iterations);

  for (i = 0; i < n_iterations; ++i)
    rf_table[i] = g_new (gfloat, RF_TABLE_SIZE);
127 128 129

  /* ************************ */

130
  report_progress (operation, 0.0, timer);
131

132
  /* Pre-calculate RF table */
133
  for (i = 0; i < n_iterations; ++i)
134
    {
135 136 137 138
      /* calculating RF feedback coefficient 'a' from desired variance
       * 'a' will change each iteration while the domain transform will
       * remain constant
       */
139 140 141 142 143 144
      sdt_dev = spatial_factor * SQRT3 *
                          (powf (2.0f, (gfloat)(n_iterations - (i + 1))) /
                             sqrtf (powf (4.0f, (gfloat) n_iterations) - 1));

      a = expf (-SQRT2 / sdt_dev);

145
      for (j = 0; j < RF_TABLE_SIZE; ++j)
146 147 148 149
        {
          rf_table[i][j] = powf (a, 1.0f + (spatial_factor / range_factor) *
                             ((gfloat)j / 255.0f));
        }
150 151
    }

152
  /* Filter Iterations */
153
  for (n = 0; n < n_iterations; ++n)
154
    {
155
      /* Horizontal Pass */
156
      for (i = 0; i < height; i += BLOCK_STRIDE)
157
        {
158
          real_stride = (i + BLOCK_STRIDE > height) ? height - i : BLOCK_STRIDE;
159
    
160 161 162 163
          rect.x = 0;
          rect.y = i;
          rect.width = width;
          rect.height = real_stride;
164
    
165 166
          gegl_buffer_get (input, &rect, 1.0, formatu8, buffer,
                           GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_CLAMP);
167
    
168
          /* Domain Transform */
169
          for (j = 0; j < rect.height; ++j)
170
            {
171 172 173
              last[0] = ((guint8*) buffer)[j * rect.width * 3];
              last[1] = ((guint8*) buffer)[j * rect.width * 3 + 1];
              last[2] = ((guint8*) buffer)[j * rect.width * 3 + 2];
174
      
175
              for (k = 0; k < rect.width; ++k)
176
                {
177 178 179
                  current[0] = ((guint8*) buffer)[j * rect.width * 3 + k * 3];
                  current[1] = ((guint8*) buffer)[j * rect.width * 3 + k * 3 + 1];
                  current[2] = ((guint8*) buffer)[j * rect.width * 3 + k * 3 + 2];
180
        
181 182 183
                  sum_channels_difference = absolute (current[0] - last[0]) +
                                            absolute (current[1] - last[1]) +
                                            absolute (current[2] - last[2]);
184
        
185 186 187 188 189
                  /* @NOTE: 'd' should be 1.0f + s_s / s_r * sum_diff
                   * However, we will store just sum_diff.
                   * 1.0f + s_s / s_r will be calculated later when calculating
                   * the RF table. This is done this way because the sum_diff is
                   * perfect to be used as the index of the RF table.
190 191
                   * d = 1.0f + (vdt_information->spatial_factor /
                   *   vdt_information->range_factor) * sum_channels_difference;
192
                   */
193
                  transforms_buffer[j * rect.width + k] = sum_channels_difference;
194
        
195 196 197
                  last[0] = current[0];
                  last[1] = current[1];
                  last[2] = current[2];
198 199 200 201
                }
            }
    
          if (n == 0)
202 203
            gegl_buffer_get (input, &rect, 1.0, format, buffer,
                             GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_CLAMP);
204
          else
205 206
            gegl_buffer_get (output, &rect, 1.0, format, buffer,
                             GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_CLAMP);
207
    
208
          /* Horizontal Filter (Left-Right) */
209
          for (j = 0; j < rect.height; ++j)
210
            {
211 212 213 214
              lastf[0] = buffer[j * rect.width * n_chan];
              lastf[1] = buffer[j * rect.width * n_chan + 1];
              lastf[2] = buffer[j * rect.width * n_chan + 2];
              lastf[3] = buffer[j * rect.width * n_chan + 3];
215
      
216
              for (k = 0; k < rect.width; ++k)
217
                {
218
                  d = transforms_buffer[j * rect.width + k];
219 220
                  w = rf_table[n][d];
        
221 222 223 224
                  currentf[0] = buffer[j * rect.width * n_chan + k * n_chan];
                  currentf[1] = buffer[j * rect.width * n_chan + k * n_chan + 1];
                  currentf[2] = buffer[j * rect.width * n_chan + k * n_chan + 2];
                  currentf[3] = buffer[j * rect.width * n_chan + k * n_chan + 3];
225
        
226 227 228 229
                  lastf[0] = ((1 - w) * currentf[0] + w * lastf[0]);
                  lastf[1] = ((1 - w) * currentf[1] + w * lastf[1]);
                  lastf[2] = ((1 - w) * currentf[2] + w * lastf[2]);
                  lastf[3] = ((1 - w) * currentf[3] + w * lastf[3]);
230
        
231 232 233 234
                  buffer[j * rect.width * n_chan + k * n_chan]     = lastf[0];
                  buffer[j * rect.width * n_chan + k * n_chan + 1] = lastf[1];
                  buffer[j * rect.width * n_chan + k * n_chan + 2] = lastf[2];
                  buffer[j * rect.width * n_chan + k * n_chan + 3] = lastf[3];
235 236 237
                }
            }
    
238
          /* Horizontal Filter (Right-Left) */
239
          for (j = 0; j < rect.height; ++j)
240
            {
241 242 243 244
              lastf[0] = buffer[j * rect.width * n_chan + (rect.width - 1) * n_chan];
              lastf[1] = buffer[j * rect.width * n_chan + (rect.width - 1) * n_chan + 1];
              lastf[2] = buffer[j * rect.width * n_chan + (rect.width - 1) * n_chan + 2];
              lastf[3] = buffer[j * rect.width * n_chan + (rect.width - 1) * n_chan + 3];
245
              
246
              for (k = rect.width - 1; k >= 0; --k)
247
                {
248 249
                  d_x_position = (k < rect.width - 1) ? k + 1 : k;
                  d = transforms_buffer[j * rect.width + d_x_position];
250 251
                  w = rf_table[n][d];
        
252 253 254 255
                  currentf[0] = buffer[j * rect.width * n_chan + k * n_chan];
                  currentf[1] = buffer[j * rect.width * n_chan + k * n_chan + 1];
                  currentf[2] = buffer[j * rect.width * n_chan + k * n_chan + 2];
                  currentf[3] = buffer[j * rect.width * n_chan + k * n_chan + 3];
256
        
257 258 259 260
                  lastf[0] = ((1 - w) * currentf[0] + w * lastf[0]);
                  lastf[1] = ((1 - w) * currentf[1] + w * lastf[1]);
                  lastf[2] = ((1 - w) * currentf[2] + w * lastf[2]);
                  lastf[3] = ((1 - w) * currentf[3] + w * lastf[3]);
261
        
262 263 264 265
                  buffer[j * rect.width * n_chan + k * n_chan] = lastf[0];
                  buffer[j * rect.width * n_chan + k * n_chan + 1] = lastf[1];
                  buffer[j * rect.width * n_chan + k * n_chan + 2] = lastf[2];
                  buffer[j * rect.width * n_chan + k * n_chan + 3] = lastf[3];
266 267 268
                }
            }
    
269 270
          gegl_buffer_set (output, &rect, 0, format, buffer,
                           GEGL_AUTO_ROWSTRIDE);
271 272
        }
  
273
      report_progress (operation, (2.0 * n + 1.0) / (2.0 * n_iterations), timer);
274
  
275
      /* Vertical Pass */
276
      for (i = 0; i < width; i += BLOCK_STRIDE)
277
        {
278
          real_stride = (i + BLOCK_STRIDE > width) ? width - i : BLOCK_STRIDE;
279
    
280 281 282 283
          rect.x = i;
          rect.y = 0;
          rect.width = real_stride;
          rect.height = height;
284
    
285 286
          gegl_buffer_get (input, &rect, 1.0, formatu8, buffer,
                           GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_CLAMP);
287
    
288
          /* Domain Transform */
289
          for (j = 0; j < rect.width; ++j)
290
            {
291 292 293
              last[0] = ((guint8*) buffer)[j * 3];
              last[1] = ((guint8*) buffer)[j * 3 + 1];
              last[2] = ((guint8*) buffer)[j * 3 + 2];
294
      
295
              for (k = 0; k < rect.height; ++k)
296
                {
297 298 299
                  current[0] = ((guint8*) buffer)[k * rect.width * 3 + j * 3];
                  current[1] = ((guint8*) buffer)[k * rect.width * 3 + j * 3 + 1];
                  current[2] = ((guint8*) buffer)[k * rect.width * 3 + j * 3 + 2];
300
        
301 302 303
                  sum_channels_difference = absolute (current[0] - last[0]) +
                                            absolute (current[1] - last[1]) +
                                            absolute (current[2] - last[2]);
304
        
305 306 307 308 309
                  /* @NOTE: 'd' should be 1.0f + s_s / s_r * sum_diff
                   * However, we will store just sum_diff.
                   * 1.0f + s_s / s_r will be calculated later when calculating
                   * the RF table. This is done this way because the sum_diff is
                   * perfect to be used as the index of the RF table.
310 311
                   * d = 1.0f + (vdt_information->spatial_factor /
                   *   vdt_information->range_factor) * sum_channels_difference;
312
                   */
313
                  transforms_buffer[k * rect.width + j] = sum_channels_difference;
314
        
315 316 317
                  last[0] = current[0];
                  last[1] = current[1];
                  last[2] = current[2];
318 319 320
                }
            }
    
321 322
          gegl_buffer_get (output, &rect, 1.0, format, buffer,
                           GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_CLAMP);
323
    
324
          /* Vertical Filter (Top-Down) */
325
          for (j = 0; j < rect.width; ++j)
326
            {
327 328 329 330
              lastf[0] = buffer[j * n_chan];
              lastf[1] = buffer[j * n_chan + 1];
              lastf[2] = buffer[j * n_chan + 2];
              lastf[3] = buffer[j * n_chan + 3];
331
      
332
              for (k = 0; k < rect.height; ++k)
333
                {
334
                  d = transforms_buffer[k * rect.width + j];
335 336
                  w = rf_table[n][d];
                  
337 338 339 340
                  currentf[0] = buffer[k * rect.width * n_chan + j * n_chan];
                  currentf[1] = buffer[k * rect.width * n_chan + j * n_chan + 1];
                  currentf[2] = buffer[k * rect.width * n_chan + j * n_chan + 2];
                  currentf[3] = buffer[k * rect.width * n_chan + j * n_chan + 3];
341
        
342 343 344 345
                  lastf[0] = ((1 - w) * currentf[0] + w * lastf[0]);
                  lastf[1] = ((1 - w) * currentf[1] + w * lastf[1]);
                  lastf[2] = ((1 - w) * currentf[2] + w * lastf[2]);
                  lastf[3] = ((1 - w) * currentf[3] + w * lastf[3]);
346
        
347 348 349 350
                  buffer[k * rect.width * n_chan + j * n_chan]     = lastf[0];
                  buffer[k * rect.width * n_chan + j * n_chan + 1] = lastf[1];
                  buffer[k * rect.width * n_chan + j * n_chan + 2] = lastf[2];
                  buffer[k * rect.width * n_chan + j * n_chan + 3] = lastf[3];
351 352 353
                }
            }
    
354
          /* Vertical Filter (Bottom-Up) */
355
          for (j = 0; j < rect.width; ++j)
356
            {
357 358 359 360
              lastf[0] = buffer[(rect.height - 1) * rect.width * n_chan + j * n_chan];
              lastf[1] = buffer[(rect.height - 1) * rect.width * n_chan + j * n_chan + 1];
              lastf[2] = buffer[(rect.height - 1) * rect.width * n_chan + j * n_chan + 2];
              lastf[3] = buffer[(rect.height - 1) * rect.width * n_chan + j * n_chan + 3];
361
              
362
              for (k = rect.height - 1; k >= 0; --k)
363
                {
364 365
                  d_y_position = (k < rect.height - 1) ? k + 1 : k;
                  d = transforms_buffer[d_y_position * rect.width + j];
366 367
                  w = rf_table[n][d];
        
368 369 370 371
                  currentf[0] = buffer[k * rect.width * n_chan + j * n_chan];
                  currentf[1] = buffer[k * rect.width * n_chan + j * n_chan + 1];
                  currentf[2] = buffer[k * rect.width * n_chan + j * n_chan + 2];
                  currentf[3] = buffer[k * rect.width * n_chan + j * n_chan + 3];
372
        
373 374 375 376
                  lastf[0] = ((1 - w) * currentf[0] + w * lastf[0]);
                  lastf[1] = ((1 - w) * currentf[1] + w * lastf[1]);
                  lastf[2] = ((1 - w) * currentf[2] + w * lastf[2]);
                  lastf[3] = ((1 - w) * currentf[3] + w * lastf[3]);
377
        
378 379 380 381
                  buffer[k * rect.width * n_chan + j * n_chan]     = lastf[0];
                  buffer[k * rect.width * n_chan + j * n_chan + 1] = lastf[1];
                  buffer[k * rect.width * n_chan + j * n_chan + 2] = lastf[2];
                  buffer[k * rect.width * n_chan + j * n_chan + 3] = lastf[3];
382 383 384
                }
            }
    
385 386
          gegl_buffer_set (output, &rect, 0, format, buffer,
                           GEGL_AUTO_ROWSTRIDE);
387 388
        }
  
389
      report_progress (operation, (2.0 * n + 2.0) / (2.0 * n_iterations), timer);
390 391
    }

392 393 394 395 396
  g_free (transforms_buffer);
  g_free (buffer);

  for (i = 0; i < n_iterations; ++i)
    g_free (rf_table[i]);
397

398 399
  g_free (rf_table);
  g_timer_destroy (timer);
400 401 402 403 404 405 406

  return 0;
}

static void
prepare (GeglOperation *operation)
{
407 408 409 410 411
  const Babl *space  = gegl_operation_get_source_space (operation, "input");
  const Babl *format = babl_format_with_space ("R'G'B'A float", space);

  gegl_operation_set_format (operation, "input", format);
  gegl_operation_set_format (operation, "output", format);
412 413 414 415 416 417 418
}

static GeglRectangle
get_required_for_output (GeglOperation       *operation,
                         const gchar         *input_pad,
                         const GeglRectangle *roi)
{
419 420
  GeglRectangle *result =
      gegl_operation_source_get_bounding_box (operation, "input");
421 422

  /* Don't request an infinite plane */
423
  if (result && gegl_rectangle_is_infinite_plane (result))
424 425
    return *roi;

426
  return *result;
427 428 429 430 431 432
}

static GeglRectangle
get_cached_region (GeglOperation       *operation,
                   const GeglRectangle *roi)
{
433 434
  GeglRectangle *result =
       gegl_operation_source_get_bounding_box (operation, "input");
435

436
  if (result && gegl_rectangle_is_infinite_plane (result))
437 438
    return *roi;

439
  return *result;
440 441 442 443 444 445 446 447 448
}

static gboolean
process (GeglOperation       *operation,
         GeglBuffer          *input,
         GeglBuffer          *output,
         const GeglRectangle *result,
         gint                 level)
{
449 450
  GeglProperties  *o = GEGL_PROPERTIES (operation);
  gfloat range_factor;
451

452 453
  if (o->edge_preservation != 0.0)
    range_factor = (1.0 / o->edge_preservation) - 1.0f;
454
  else
455 456 457 458 459 460 461 462 463 464 465 466
    range_factor = G_MAXFLOAT;

  /* Buffer is ready for domain transform */
  domain_transform (operation,
                    result->width,
                    result->height,
                    4,
                    o->spatial_factor,
                    range_factor,
                    o->n_iterations,
                    input,
                    output);
467 468 469 470 471

  return TRUE;
}

/* Pass-through when trying to perform a reduction on an infinite plane
472
 * or when edge preservation is one.
473 474 475 476 477 478 479 480
 */
static gboolean
operation_process (GeglOperation        *operation,
                   GeglOperationContext *context,
                   const gchar          *output_prop,
                   const GeglRectangle  *result,
                   gint                  level)
{
481
  GeglProperties      *o = GEGL_PROPERTIES (operation);
482 483 484 485 486 487 488
  GeglOperationClass  *operation_class;

  const GeglRectangle *in_rect =
    gegl_operation_source_get_bounding_box (operation, "input");

  operation_class = GEGL_OPERATION_CLASS (gegl_op_parent_class);

489 490
  if ((in_rect && gegl_rectangle_is_infinite_plane (in_rect)) ||
       o->edge_preservation == 1.0)
491 492 493 494 495
    {
      gpointer in = gegl_operation_context_get_object (context, "input");
      gegl_operation_context_take_object (context, "output",
                                          g_object_ref (G_OBJECT (in)));
      return TRUE;
496 497 498
    }

  /* chain up, which will create the needed buffers for our actual
499 500
   * process function
   */
501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518
  return operation_class->process (operation, context, output_prop, result,
    gegl_operation_context_get_level (context));
}

/* This is called at the end of the gobject class_init function.
 *
 * Here we override the standard passthrough options for the rect
 * computations.
 */
static void
gegl_op_class_init (GeglOpClass *klass)
{
  GeglOperationClass       *operation_class;
  GeglOperationFilterClass *filter_class;

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

519 520 521 522
  filter_class->process                    = process;
  operation_class->prepare                 = prepare;
  operation_class->threaded                = FALSE;
  operation_class->process                 = operation_process;
523
  operation_class->get_required_for_output = get_required_for_output;
524 525
  operation_class->get_cached_region       = get_cached_region;
  operation_class->opencl_support          = FALSE;
526 527

  gegl_operation_class_set_keys (operation_class,
528 529 530 531 532 533 534
    "name",        "gegl:domain-transform",
    "title",       _("Smooth by Domain Transform"),
    "categories" , "enhance:noise-reduction",
    "description", _("An edge-preserving smoothing filter implemented with the "
                     "Domain Transform recursive technique. Similar to a "
                     "bilateral filter, but faster to compute."),
    NULL);
535 536 537
}

#endif