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_int (radius, _("Radius"), 8)
27     description(_("Radius of square pixel region, (width and height will be radius*2+1)"))
28     value_range (0, 100)
29     ui_range    (0, 40)
30     ui_gamma    (1.5)
31     ui_meta     ("unit", "pixel-distance")
32 
33 property_int (pairs, _("Pairs"), 2)
34   description(_("Number of pairs; higher number preserves more acute features"))
35   value_range (1, 2)
36 
37 #else
38 
39 #define GEGL_OP_AREA_FILTER
40 #define GEGL_OP_NAME     snn_mean
41 #define GEGL_OP_C_SOURCE snn-mean.c
42 
43 #include "gegl-op.h"
44 
45 static void
46 snn_mean (GeglBuffer          *src,
47           const GeglRectangle *src_rect,
48           GeglBuffer          *dst,
49           const GeglRectangle *dst_rect,
50           gdouble              radius,
51           gint                 pairs,
52           gint                 level,
53           const Babl          *space);
54 
55 
56 static void prepare (GeglOperation *operation)
57 {
58   const Babl *space = gegl_operation_get_source_space (operation, "input");
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->radius);
63   gegl_operation_set_format (operation, "input",
64                              babl_format_with_space ("RGBA float", space));
65   gegl_operation_set_format (operation, "output",
66                              babl_format_with_space ("RGBA float", space));
67 }
68 
69 static gboolean
70 cl_process (GeglOperation       *operation,
71             GeglBuffer          *input,
72             GeglBuffer          *output,
73             const GeglRectangle *result);
74 
75 static gboolean
76 process (GeglOperation       *operation,
77          GeglBuffer          *input,
78          GeglBuffer          *output,
79          const GeglRectangle *result,
80          gint                 level)
81 {
82   GeglProperties      *o = GEGL_PROPERTIES (operation);
83   GeglRectangle        compute;
84 
85   if (gegl_operation_use_opencl (operation))
86     if (cl_process (operation, input, output, result))
87       return TRUE;
88 
89   compute = gegl_operation_get_required_for_output (operation, "input", result);
90 
91   if (o->radius < 1.0)
92     {
93       output = g_object_ref (input);
94     }
95   else
96     {
97       snn_mean (input, &compute, output, result, o->radius, o->pairs, level, gegl_operation_get_format (operation, "output"));
98     }
99 
100   return  TRUE;
101 }
102 
103 #define POW2(a)((a)*(a))
104 
105 static inline gfloat colordiff (gfloat *pixA,
106                                 gfloat *pixB)
107 {
108   return POW2(pixA[0]-pixB[0])+
109          POW2(pixA[1]-pixB[1])+
110          POW2(pixA[2]-pixB[2]);
111 }
112 
113 
114 static void
115 snn_mean (GeglBuffer          *src,
116           const GeglRectangle *src_rect,
117           GeglBuffer          *dst,
118           const GeglRectangle *dst_rect,
119           gdouble              dradius,
120           gint                 pairs,
121           gint                 level,
122           const Babl          *space)
123 {
124   gint x,y;
125   gint offset;
126   gfloat *src_buf;
127   gfloat *dst_buf;
128   gint radius = dradius;
129   GeglRectangle src_rect_scaled, dst_rect_scaled;
130   if (level)
131   {
132     src_rect_scaled = *src_rect;
133     dst_rect_scaled = *dst_rect;
134     src_rect_scaled.x >>= level;
135     src_rect_scaled.y >>= level;
136     src_rect_scaled.width >>= level;
137     src_rect_scaled.height >>= level;
138     dst_rect_scaled.x >>= level;
139     dst_rect_scaled.y >>= level;
140     dst_rect_scaled.width >>= level;
141     dst_rect_scaled.height >>= level;
142     src_rect = &src_rect_scaled;
143     dst_rect = &dst_rect_scaled;
144     dradius /= (1<<level);
145   }
146   radius = dradius;
147 
148 
149   src_buf = g_new0 (gfloat, src_rect->width * src_rect->height * 4);
150   dst_buf = g_new0 (gfloat, dst_rect->width * dst_rect->height * 4);
151 
152   gegl_buffer_get (src, src_rect, 1.0/(1<<level), babl_format_with_space ("RGBA float", space), src_buf,
153                    GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
154 
155   offset = 0;
156 
157   for (y=0; y<dst_rect->height; y++)
158     {
159       gfloat *center_pix;
160 
161       center_pix = src_buf + ((radius) + (y+radius)* src_rect->width)*4;
162 
163       for (x=0; x<dst_rect->width; x++)
164         {
165           gint u,v;
166 
167           gfloat  accumulated[4]={0.0f, 0.0f, 0.0f, 0.0f};
168           gint    count=0;
169 
170           /* iterate through the upper left quarter of pixels */
171           for (v=-radius;v<=0;v++)
172             for (u=-radius;u<= (pairs==1?radius:0);u++)
173               {
174                 gfloat *selected_pix = center_pix;
175                 gfloat  best_diff = 1000.0;
176                 gint    i;
177 
178                 /* skip computations for the center pixel */
179                 if (u != 0 &&
180                     v != 0)
181                   {
182                     /* compute the coordinates of the symmetric pairs for
183                      * this locaiton in the quadrant
184                      */
185                     gint xs[4], ys[4];
186 
187                     xs[0] = x+u+radius;
188                     xs[1] = x-u+radius;
189                     xs[2] = x-u+radius;
190                     xs[3] = x+u+radius;
191                     ys[0] = y+v+radius;
192                     ys[1] = y-v+radius;
193                     ys[2] = y+v+radius;
194                     ys[3] = y-v+radius;
195 
196                     /* check which member of the symmetric quadruple to use */
197                     for (i=0;i<pairs*2;i++)
198                       {
199                         if (xs[i] >= 0 && xs[i] < src_rect->width &&
200                             ys[i] >= 0 && ys[i] < src_rect->height)
201                           {
202                             gfloat *tpix = src_buf + (xs[i]+ys[i]* src_rect->width)*4;
203                             gfloat diff = colordiff (tpix, center_pix);
204                             if (diff < best_diff)
205                               {
206                                 best_diff = diff;
207                                 selected_pix = tpix;
208                               }
209                           }
210                       }
211                   }
212 
213                 /* accumulate the components of the best sample from
214                  * the symmetric quadruple
215                  */
216                 for (i=0;i<4;i++)
217                   {
218                     accumulated[i] += selected_pix[i];
219                   }
220                 count++;
221 
222                 if (u==0 && v==0)
223                   break; /* to avoid doubly processing when using only 1 pair */
224               }
225           for (u=0; u<4; u++)
226             dst_buf[offset*4+u] = accumulated[u]/count;
227           offset++;
228 
229           center_pix += 4;
230         }
231     }
232   gegl_buffer_set (dst, dst_rect, level, babl_format_with_space ("RGBA float", space), dst_buf,
233                    GEGL_AUTO_ROWSTRIDE);
234   g_free (src_buf);
235   g_free (dst_buf);
236 }
237 
238 
239 #include "opencl/gegl-cl.h"
240 #include "gegl-buffer-cl-iterator.h"
241 
242 #include "opencl/snn-mean.cl.h"
243 
244 static GeglClRunData *cl_data = NULL;
245 
246 static gboolean
247 cl_snn_mean (cl_mem                in_tex,
248              cl_mem                out_tex,
249              const GeglRectangle  *src_rect,
250              const GeglRectangle  *roi,
251              gint                  radius,
252              gint                  pairs)
253 {
254   cl_int cl_err = 0;
255   size_t global_ws[2];
256 
257   if (!cl_data)
258     {
259       const char *kernel_name[] = {"snn_mean", NULL};
260       cl_data = gegl_cl_compile_and_build (snn_mean_cl_source, kernel_name);
261     }
262   if (!cl_data) return TRUE;
263 
264 
265   global_ws[0] = roi->width;
266   global_ws[1] = roi->height;
267 
268   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 0, sizeof(cl_mem),   (void*)&in_tex);
269   CL_CHECK;
270   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 1, sizeof(cl_int),   (void*)&src_rect->width);
271   CL_CHECK;
272   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 2, sizeof(cl_int),   (void*)&src_rect->height);
273   CL_CHECK;
274   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 3, sizeof(cl_mem),   (void*)&out_tex);
275   CL_CHECK;
276   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 4, sizeof(cl_int),   (void*)&radius);
277   CL_CHECK;
278   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 5, sizeof(cl_int),   (void*)&pairs);
279   CL_CHECK;
280 
281   cl_err = gegl_clEnqueueNDRangeKernel(gegl_cl_get_command_queue (),
282                                         cl_data->kernel[0], 2,
283                                         NULL, global_ws, NULL,
284                                         0, NULL, NULL);
285   CL_CHECK;
286 
287   return FALSE;
288 
289 error:
290   return TRUE;
291 }
292 
293 static gboolean
294 cl_process (GeglOperation       *operation,
295             GeglBuffer          *input,
296             GeglBuffer          *output,
297             const GeglRectangle *result)
298 {
299   const Babl *in_format  = gegl_operation_get_format (operation, "input");
300   const Babl *out_format = gegl_operation_get_format (operation, "output");
301   gint err;
302 
303   GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
304   GeglProperties *o = GEGL_PROPERTIES (operation);
305 
306   GeglBufferClIterator *i = gegl_buffer_cl_iterator_new (output,
307                                                          result,
308                                                          out_format,
309                                                          GEGL_CL_BUFFER_WRITE);
310 
311   gint read = gegl_buffer_cl_iterator_add_2 (i,
312                                              input,
313                                              result,
314                                              in_format,
315                                              GEGL_CL_BUFFER_READ,
316                                              op_area->left,
317                                              op_area->right,
318                                              op_area->top,
319                                              op_area->bottom,
320                                              GEGL_ABYSS_NONE);
321 
322   while (gegl_buffer_cl_iterator_next (i, &err))
323     {
324       if (err) return FALSE;
325 
326       err = cl_snn_mean(i->tex[read],
327                         i->tex[0],
328                         &i->roi[read],
329                         &i->roi[0],
330                         ceil(o->radius),
331                         o->pairs);
332 
333       if (err) return FALSE;
334     }
335 
336   return TRUE;
337 }
338 
339 /* Pass-through when radius parameter is set to zero */
340 
341 static gboolean
342 operation_process (GeglOperation        *operation,
343                    GeglOperationContext *context,
344                    const gchar          *output_prop,
345                    const GeglRectangle  *result,
346                    gint                  level)
347 {
348   GeglOperationClass  *operation_class;
349   GeglProperties      *o = GEGL_PROPERTIES (operation);
350 
351   operation_class = GEGL_OPERATION_CLASS (gegl_op_parent_class);
352 
353   if (! o->radius)
354     {
355       gpointer in = gegl_operation_context_get_object (context, "input");
356       gegl_operation_context_take_object (context, "output",
357                                           g_object_ref (G_OBJECT (in)));
358       return TRUE;
359     }
360 
361   return operation_class->process (operation, context, output_prop, result,
362                                    gegl_operation_context_get_level (context));
363 }
364 
365 static void
366 gegl_op_class_init (GeglOpClass *klass)
367 {
368   GeglOperationClass       *operation_class;
369   GeglOperationFilterClass *filter_class;
370 
371   operation_class = GEGL_OPERATION_CLASS (klass);
372   filter_class    = GEGL_OPERATION_FILTER_CLASS (klass);
373 
374   filter_class->process    = process;
375   operation_class->prepare = prepare;
376   operation_class->process = operation_process;
377 
378   operation_class->opencl_support = TRUE;
379 
380   gegl_operation_class_set_keys (operation_class,
381     "name"       , "gegl:snn-mean",
382     "categories" , "enhance:noise-reduction",
383     "title",       _("Symmetric Nearest Neighbour"),
384     "reference-hash", "dc342e01d1016f92f21d1b2196dc1a7f",
385     "reference-hashB", "1f5c30085011311cf743ddc91a44f1f0",
386     "description",
387         _("Noise reducing edge preserving blur filter based "
388           "on Symmetric Nearest Neighbours"),
389         NULL);
390 }
391 
392 #endif
393