operations/common/mean-curvature-blur.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 2018 Thomas Manni <thomas.manni@free.fr>
  17  */
  18
  19 #include "config.h"
  20 #include <glib/gi18n-lib.h>
  21
  22 #ifdef GEGL_PROPERTIES
  23
  24 property_int  (iterations, _("Iterations"), 20)
  25   description (_("Controls the number of iterations"))
  26   value_range (0, 500)
  27   ui_range    (0, 60)
  28
  29 #else
  30
  31 #define GEGL_OP_AREA_FILTER
  32 #define GEGL_OP_NAME         mean_curvature_blur
  33 #define GEGL_OP_C_SOURCE     mean-curvature-blur.c
  34
  35 #define POW2(x) ((x)*(x))
  36
  37 #include "gegl-op.h"
  38
  39 static void
  40 prepare (GeglOperation *operation)
  41 {
  42   const Babl *space = gegl_operation_get_source_space (operation, "input");
  43   GeglOperationAreaFilter *area   = GEGL_OPERATION_AREA_FILTER (operation);
  44   GeglProperties          *o      = GEGL_PROPERTIES (operation);
  45   const Babl              *format = babl_format_with_space ("R'G'B'A float", space);
  46
  47   area->left = area->right = area->top = area->bottom = o->iterations;
  48
  49   gegl_operation_set_format (operation, "input",  format);
  50   gegl_operation_set_format (operation, "output", format);
  51 }
  52
  53 static GeglRectangle
  54 get_bounding_box (GeglOperation *operation)
  55 {
  56   GeglRectangle  result = { 0, 0, 0, 0 };
  57   GeglRectangle *in_rect;
  58
  59   in_rect = gegl_operation_source_get_bounding_box (operation, "input");
  60   if (in_rect)
  61     {
  62       result = *in_rect;
  63     }
  64
  65   return result;
  66 }
  67
  68 static void
  69 mean_curvature_flow (gfloat *src_buf,
  70                      gint    src_stride,
  71                      gfloat *dst_buf,
  72                      gint    dst_width,
  73                      gint    dst_height,
  74                      gint    dst_stride)
  75 {
  76   gint c;
  77   gint x, y;
  78   gfloat *center_pix;
  79
  80 #define O(u,v) (((u)+((v) * src_stride)) * 4)
  81   gint   offsets[8] = { O( -1, -1), O(0, -1), O(1, -1),
  82                         O( -1,  0),           O(1,  0),
  83                         O( -1,  1), O(0, 1),  O(1,  1)};
  84 #undef O
  85
  86 #define LEFT(c) ((center_pix + offsets[3])[c])
  87 #define RIGHT(c) ((center_pix + offsets[4])[c])
  88 #define TOP(c) ((center_pix + offsets[1])[c])
  89 #define BOTTOM(c) ((center_pix + offsets[6])[c])
  90 #define TOPLEFT(c) ((center_pix + offsets[0])[c])
  91 #define TOPRIGHT(c) ((center_pix + offsets[2])[c])
  92 #define BOTTOMLEFT(c) ((center_pix + offsets[5])[c])
  93 #define BOTTOMRIGHT(c) ((center_pix + offsets[7])[c])
  94
  95   for (y = 0; y < dst_height; y++)
  96     {
  97       gint dst_offset = dst_stride * y;
  98       center_pix = src_buf + ((y+1) * src_stride + 1) * 4;
  99
 100       for (x = 0; x < dst_width; x++)
 101         {
 102           for (c = 0; c < 3; c++) /* do each color component individually */
 103             {
 104               gdouble dx  = RIGHT(c) - LEFT(c);
 105               gdouble dy  = BOTTOM(c) - TOP(c);
 106               gdouble magnitude = sqrt (POW2(dx) + POW2(dy));
 107
 108               dst_buf[dst_offset * 4 + c] = center_pix[c];
 109
 110               if (magnitude)
 111                 {
 112                   gdouble dx2 = POW2(dx);
 113                   gdouble dy2 = POW2(dy);
 114
 115                   gdouble dxx = RIGHT(c) + LEFT(c) - 2. * center_pix[c];
 116                   gdouble dyy = BOTTOM(c) + TOP(c) - 2. * center_pix[c];
 117                   gdouble dxy = 0.25 * (BOTTOMRIGHT(c) - TOPRIGHT(c) - BOTTOMLEFT(c) + TOPLEFT(c));
 118
 119                   gdouble n = dx2 * dyy + dy2 * dxx - 2. * dx * dy * dxy;
 120                   gdouble d = sqrt (pow (dx2 + dy2, 3.));
 121                   gdouble mean_curvature = n / d;
 122
 123                   dst_buf[dst_offset * 4 + c] += (0.25 * magnitude * mean_curvature);
 124                 }
 125             }
 126
 127           dst_buf[dst_offset * 4 + 3] = center_pix[3];
 128
 129           dst_offset++;
 130           center_pix += 4;
 131         }
 132     }
 133
 134 #undef LEFT
 135 #undef RIGHT
 136 #undef TOP
 137 #undef BOTTOM
 138 #undef TOPLEFT
 139 #undef TOPRIGHT
 140 #undef BOTTOMLEFT
 141 #undef BOTTOMRIGHT
 142
 143 }
 144
 145 static gboolean
 146 process (GeglOperation       *operation,
 147          GeglBuffer          *input,
 148          GeglBuffer          *output,
 149          const GeglRectangle *roi,
 150          gint                 level)
 151 {
 152   GeglProperties  *o = GEGL_PROPERTIES (operation);
 153   const Babl      *format = gegl_operation_get_format (operation, "output");
 154
 155   gint iteration;
 156   gint stride;
 157   gfloat *src_buf;
 158   gfloat *dst_buf;
 159   GeglRectangle rect;
 160
 161   rect = *roi;
 162
 163   rect.x      -= o->iterations;
 164   rect.y      -= o->iterations;
 165   rect.width  += o->iterations*2;
 166   rect.height += o->iterations*2;
 167
 168   stride = roi->width + o->iterations * 2;
 169
 170   src_buf = g_new (gfloat,
 171          (stride) * (roi->height + o->iterations * 2) * 4);
 172
 173   dst_buf = g_new0 (gfloat,
 174          (stride) * (roi->height + o->iterations * 2) * 4);
 175
 176   gegl_buffer_get (input, &rect, 1.0, format, src_buf, stride * 4 * 4,
 177                    GEGL_ABYSS_CLAMP);
 178
 179   for (iteration = 0; iteration < o->iterations; iteration++)
 180     {
 181       mean_curvature_flow (src_buf, stride,
 182                            dst_buf,
 183                            roi->width  + (o->iterations - 1 - iteration) * 2,
 184                            roi->height + (o->iterations - 1 - iteration) * 2,
 185                            stride);
 186
 187       { /* swap buffers */
 188         gfloat *tmp = src_buf;
 189         src_buf = dst_buf;
 190         dst_buf = tmp;
 191       }
 192     }
 193
 194   gegl_buffer_set (output, roi, 0, format, src_buf, stride * 4 * 4);
 195
 196   g_free (src_buf);
 197   g_free (dst_buf);
 198
 199   return TRUE;
 200 }
 201
 202 /* Pass-through when iterations parameter is set to zero */
 203
 204 static gboolean
 205 operation_process (GeglOperation        *operation,
 206                    GeglOperationContext *context,
 207                    const gchar          *output_prop,
 208                    const GeglRectangle  *result,
 209                    gint                  level)
 210 {
 211   GeglOperationClass  *operation_class;
 212   GeglProperties      *o = GEGL_PROPERTIES (operation);
 213
 214   operation_class = GEGL_OPERATION_CLASS (gegl_op_parent_class);
 215
 216   if (! o->iterations)
 217     {
 218       gpointer in = gegl_operation_context_get_object (context, "input");
 219       gegl_operation_context_take_object (context, "output",
 220                                           g_object_ref (G_OBJECT (in)));
 221       return TRUE;
 222     }
 223
 224   return operation_class->process (operation, context, output_prop, result,
 225                                    gegl_operation_context_get_level (context));
 226 }
 227
 228 static void
 229 gegl_op_class_init (GeglOpClass *klass)
 230 {
 231   GeglOperationClass       *operation_class;
 232   GeglOperationFilterClass *filter_class;
 233
 234   operation_class = GEGL_OPERATION_CLASS (klass);
 235   filter_class    = GEGL_OPERATION_FILTER_CLASS (klass);
 236
 237   filter_class->process             = process;
 238   operation_class->process          = operation_process;
 239   operation_class->prepare          = prepare;
 240   operation_class->get_bounding_box = get_bounding_box;
 241   operation_class->opencl_support   = FALSE;
 242
 243   gegl_operation_class_set_keys (operation_class,
 244     "name",           "gegl:mean-curvature-blur",
 245     "title",          _("Mean Curvature Blur"),
 246     "categories",     "blur",
 247     "reference-hash", "8856d371c39a439e501dc2f2a74d6417",
 248     "description",    _("Regularize geometry at a speed proportional to the local mean curvature value"),
 249     NULL);
 250 }
 251
 252 #endif