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_double (blur_radius, _("Blur radius"), 4.0)
27   description(_("Radius of square pixel region, (width and height will be radius*2+1)."))
28   value_range   (0.0, 1000.0)
29   ui_range      (0.0, 100.0)
30   ui_gamma      (1.5)
31 
32 property_double (edge_preservation, _("Edge preservation"), 8.0)
33   description   (_("Amount of edge preservation"))
34   value_range   (0.0, 100.0)
35 
36 #else
37 
38 #define GEGL_OP_AREA_FILTER
39 #define GEGL_OP_NAME     bilateral_filter
40 #define GEGL_OP_C_SOURCE bilateral-filter.c
41 
42 #include "gegl-op.h"
43 
44 static void
45 bilateral_filter (GeglBuffer          *src,
46                   const GeglRectangle *src_rect,
47                   GeglBuffer          *dst,
48                   const GeglRectangle *dst_rect,
49                   gdouble              radius,
50                   gdouble              preserve,
51                   const Babl          *format);
52 
53 #include <stdio.h>
54 
55 static void prepare (GeglOperation *operation)
56 {
57   const Babl *space = gegl_operation_get_source_space (operation, "input");
58   const Babl *format = babl_format_with_space ("RGBA float", space);
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->blur_radius);
63   gegl_operation_set_format (operation, "input", format);
64   gegl_operation_set_format (operation, "output", format);
65 }
66 
67 #include "opencl/gegl-cl.h"
68 #include "gegl-buffer-cl-iterator.h"
69 
70 #include "opencl/bilateral-filter.cl.h"
71 
72 static GeglClRunData *cl_data = NULL;
73 
74 static gboolean
75 cl_bilateral_filter (cl_mem                in_tex,
76                      cl_mem                out_tex,
77                      size_t                global_worksize,
78                      const GeglRectangle  *roi,
79                      gfloat                radius,
80                      gfloat                preserve)
81 {
82   cl_int cl_err = 0;
83   size_t global_ws[2];
84 
85   if (!cl_data)
86     {
87       const char *kernel_name[] = {"bilateral_filter", NULL};
88       cl_data = gegl_cl_compile_and_build (bilateral_filter_cl_source, kernel_name);
89     }
90   if (!cl_data) return TRUE;
91 
92   global_ws[0] = roi->width;
93   global_ws[1] = roi->height;
94 
95   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 0, sizeof(cl_mem),   (void*)&in_tex);
96   CL_CHECK;
97   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 1, sizeof(cl_mem),   (void*)&out_tex);
98   CL_CHECK;
99   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 2, sizeof(cl_float), (void*)&radius);
100   CL_CHECK;
101   cl_err = gegl_clSetKernelArg(cl_data->kernel[0], 3, sizeof(cl_float), (void*)&preserve);
102   CL_CHECK;
103 
104   cl_err = gegl_clEnqueueNDRangeKernel(gegl_cl_get_command_queue (),
105                                        cl_data->kernel[0], 2,
106                                        NULL, global_ws, NULL,
107                                        0, NULL, NULL);
108   CL_CHECK;
109 
110   return FALSE;
111 
112 error:
113   return TRUE;
114 }
115 
116 static gboolean
117 cl_process (GeglOperation       *operation,
118             GeglBuffer          *input,
119             GeglBuffer          *output,
120             const GeglRectangle *result)
121 {
122   const Babl *in_format  = gegl_operation_get_format (operation, "input");
123   const Babl *out_format = gegl_operation_get_format (operation, "output");
124   gint err;
125 
126   GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
127   GeglProperties *o = GEGL_PROPERTIES (operation);
128 
129   GeglBufferClIterator *i = gegl_buffer_cl_iterator_new (output,
130                                                          result,
131                                                          out_format,
132                                                          GEGL_CL_BUFFER_WRITE);
133 
134   gint read = gegl_buffer_cl_iterator_add_2 (i,
135                                              input,
136                                              result,
137                                              in_format,
138                                              GEGL_CL_BUFFER_READ,
139                                              op_area->left,
140                                              op_area->right,
141                                              op_area->top,
142                                              op_area->bottom,
143                                              GEGL_ABYSS_NONE);
144 
145   while (gegl_buffer_cl_iterator_next (i, &err))
146     {
147       if (err) return FALSE;
148 
149       err = cl_bilateral_filter(i->tex[read],
150                                 i->tex[0],
151                                 i->size[0],
152                                 &i->roi[0],
153                                 ceil(o->blur_radius),
154                                 o->edge_preservation);
155 
156       if (err) return FALSE;
157     }
158 
159   return TRUE;
160 }
161 
162 static gboolean
163 process (GeglOperation       *operation,
164          GeglBuffer          *input,
165          GeglBuffer          *output,
166          const GeglRectangle *result,
167          gint                 level)
168 {
169   GeglProperties *o = GEGL_PROPERTIES (operation);
170   GeglRectangle compute;
171   const Babl *format = gegl_operation_get_format (operation, "output");
172 
173   if (o->blur_radius >= 1.0 && gegl_operation_use_opencl (operation))
174     if (cl_process (operation, input, output, result))
175       return TRUE;
176 
177   compute = gegl_operation_get_required_for_output (operation, "input",result);
178 
179   if (o->blur_radius < 1.0)
180     {
181       gegl_buffer_copy (input, result, GEGL_ABYSS_NONE,
182                         output, result);
183     }
184   else
185     {
186       bilateral_filter (input, &compute, output, result, o->blur_radius, o->edge_preservation, format);
187     }
188 
189   return  TRUE;
190 }
191 
192 static void
193 bilateral_filter (GeglBuffer          *src,
194                   const GeglRectangle *src_rect,
195                   GeglBuffer          *dst,
196                   const GeglRectangle *dst_rect,
197                   gdouble              radius,
198                   gdouble              preserve,
199                   const Babl          *format)
200 {
201   gfloat *gauss;
202   gint x,y;
203   gint offset;
204   gfloat *src_buf;
205   gfloat *dst_buf;
206   gint width = (gint) radius * 2 + 1;
207   gint iradius = radius;
208   gint src_width = src_rect->width;
209   gint src_height = src_rect->height;
210 
211   gauss = g_newa (gfloat, width * width);
212   src_buf = g_new0 (gfloat, src_rect->width * src_rect->height * 4);
213   dst_buf = g_new0 (gfloat, dst_rect->width * dst_rect->height * 4);
214 
215   gegl_buffer_get (src, src_rect, 1.0, format, src_buf, GEGL_AUTO_ROWSTRIDE,
216                    GEGL_ABYSS_NONE);
217 
218   offset = 0;
219 
220 #define POW2(a) ((a)*(a))
221   for (y=-iradius;y<=iradius;y++)
222     for (x=-iradius;x<=iradius;x++)
223       {
224         gauss[x+(int)radius + (y+(int)radius)*width] = exp(- 0.5*(POW2(x)+POW2(y))/radius   );
225       }
226 
227   for (y=0; y<dst_rect->height; y++)
228     for (x=0; x<dst_rect->width; x++)
229       {
230         gint u,v;
231         gfloat *center_pix = src_buf + ((x+iradius)+((y+iradius) * src_width)) * 4;
232         gfloat  accumulated[4]={0,0,0,0};
233         gfloat  count=0.0;
234 
235         for (v=-iradius;v<=iradius;v++)
236           for (u=-iradius;u<=iradius;u++)
237             {
238               gint i,j;
239               i = x + radius + u;
240               j = y + radius + v;
241               if (i >= 0 && i < src_width &&
242                   j >= 0 && j < src_height)
243                 {
244                   gint c;
245 
246                   gfloat *src_pix = src_buf + (i + j * src_width) * 4;
247 
248                   gfloat diff_map   = exp (- (POW2(center_pix[0] - src_pix[0])+
249                                               POW2(center_pix[1] - src_pix[1])+
250                                               POW2(center_pix[2] - src_pix[2])) * preserve
251                                           );
252                   gfloat gaussian_weight;
253                   gfloat weight;
254 
255                   gaussian_weight = gauss[u+(int)radius+(v+(int)radius)*width];
256 
257                   weight = diff_map * gaussian_weight;
258 
259                   for (c=0;c<4;c++)
260                     {
261                       accumulated[c] += src_pix[c] * weight;
262                     }
263                   count += weight;
264                 }
265             }
266 
267         for (u=0; u<4;u++)
268           dst_buf[offset*4+u] = accumulated[u]/count;
269         offset++;
270       }
271   gegl_buffer_set (dst, dst_rect, 0, format, dst_buf,
272                    GEGL_AUTO_ROWSTRIDE);
273   g_free (src_buf);
274   g_free (dst_buf);
275 }
276 
277 
278 static void
279 gegl_op_class_init (GeglOpClass *klass)
280 {
281   GeglOperationClass       *operation_class;
282   GeglOperationFilterClass *filter_class;
283 
284   operation_class  = GEGL_OPERATION_CLASS (klass);
285   filter_class     = GEGL_OPERATION_FILTER_CLASS (klass);
286 
287   filter_class->process   = process;
288   operation_class->prepare = prepare;
289 
290   operation_class->opencl_support = TRUE;
291 
292   gegl_operation_class_set_keys (operation_class,
293            "name", "gegl:bilateral-filter",
294            "title", _("Bilateral Filter"),
295            "categories", "enhance:noise-reduction",
296            "reference-hash", "ae2daa632761f829c4a59225f17bf211",
297            "reference-hashB", "5cfcdea9b2f5917f48c54a8972374d8a",
298            "description",
299            _("Like a gaussian blur; but where the contribution for each neighbourhood "
300           "pixel is also weighted by the color difference with the original center pixel."),
301            NULL);
302 }
303 
304 #endif
305