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 2015 Thomas Manni <thomas.manni@free.fr>
17  *
18  */
19 
20 #include "config.h"
21 #include <glib/gi18n-lib.h>
22 #include <stdlib.h>
23 
24 #ifdef GEGL_PROPERTIES
25 
26 enum_start (gegl_median_blur_neighborhood)
27   enum_value (GEGL_MEDIAN_BLUR_NEIGHBORHOOD_SQUARE,  "square",  N_("Square"))
28   enum_value (GEGL_MEDIAN_BLUR_NEIGHBORHOOD_CIRCLE,  "circle",  N_("Circle"))
29   enum_value (GEGL_MEDIAN_BLUR_NEIGHBORHOOD_DIAMOND, "diamond", N_("Diamond"))
30 enum_end (GeglMedianBlurNeighborhood)
31 
32 enum_start (gegl_median_blur_abyss_policy)
33    enum_value (GEGL_MEDIAN_BLUR_ABYSS_NONE,  "none",  N_("None"))
34    enum_value (GEGL_MEDIAN_BLUR_ABYSS_CLAMP, "clamp", N_("Clamp"))
35 enum_end (GeglMedianBlurAbyssPolicy)
36 
37 property_enum (neighborhood, _("Neighborhood"),
38                GeglMedianBlurNeighborhood, gegl_median_blur_neighborhood,
39                GEGL_MEDIAN_BLUR_NEIGHBORHOOD_CIRCLE)
40   description (_("Neighborhood type"))
41 
42 property_int  (radius, _("Radius"), 3)
43   value_range (-400, 400)
44   ui_range    (0, 100)
45   ui_meta     ("unit", "pixel-distance")
46   description (_("Neighborhood radius, a negative value will calculate with inverted percentiles"))
47 
48 property_double  (percentile, _("Percentile"), 50)
49   value_range (0, 100)
50   description (_("Neighborhood color percentile"))
51 
52 property_double  (alpha_percentile, _("Alpha percentile"), 50)
53   value_range (0, 100)
54   description (_("Neighborhood alpha percentile"))
55 
56 property_enum (abyss_policy, _("Abyss policy"), GeglMedianBlurAbyssPolicy,
57                gegl_median_blur_abyss_policy, GEGL_MEDIAN_BLUR_ABYSS_CLAMP)
58   description (_("How image edges are handled"))
59 
60 property_boolean (high_precision, _("High precision"), FALSE)
61   description (_("Avoid clipping and quantization (slower)"))
62 
63 #else
64 
65 #define GEGL_OP_AREA_FILTER
66 #define GEGL_OP_NAME         median_blur
67 #define GEGL_OP_C_SOURCE     median-blur.c
68 
69 #include "gegl-op.h"
70 
71 #define DEFAULT_N_BINS   256
72 #define MAX_CHUNK_WIDTH  128
73 #define MAX_CHUNK_HEIGHT 128
74 
75 #define SAFE_CLAMP(x, min, max) ((x) > (min) ? (x) < (max) ? (x) : (max) : (min))
76 
77 static gfloat        default_bin_values[DEFAULT_N_BINS];
78 static gint          default_alpha_values[DEFAULT_N_BINS];
79 static volatile gint default_values_initialized = FALSE;
80 
81 typedef struct
82 {
83   gboolean  quantize;
84   gint     *neighborhood_outline;
85 } UserData;
86 
87 typedef struct
88 {
89   gfloat value;
90   gint   index;
91 } InputValue;
92 
93 typedef struct
94 {
95   gint   *bins;
96   gfloat *bin_values;
97   gint    last_median;
98   gint    last_median_sum;
99 } HistogramComponent;
100 
101 typedef struct
102 {
103   HistogramComponent  components[4];
104   gint               *alpha_values;
105   gint                count;
106   gint                size;
107   gint                n_components;
108   gint                n_color_components;
109 } Histogram;
110 
111 typedef enum
112 {
113   LEFT_TO_RIGHT,
114   RIGHT_TO_LEFT,
115   TOP_TO_BOTTOM
116 } Direction;
117 
118 static inline gfloat
119 histogram_get_median (Histogram *hist,
120                       gint       component,
121                       gdouble    percentile)
122 {
123   gint                count = hist->count;
124   HistogramComponent *comp  = &hist->components[component];
125   gint                i     = comp->last_median;
126   gint                sum   = comp->last_median_sum;
127 
128   if (component == hist->n_color_components)
129     count = hist->size;
130 
131   if (count == 0)
132     return 0.0f;
133 
134   count = (gint) ceil (count * percentile);
135   count = MAX (count, 1);
136 
137   if (sum < count)
138     {
139       while ((sum += comp->bins[++i]) < count);
140     }
141   else
142     {
143       while ((sum -= comp->bins[i--]) >= count);
144       sum += comp->bins[++i];
145     }
146 
147   comp->last_median     = i;
148   comp->last_median_sum = sum;
149 
150   return comp->bin_values[i];
151 }
152 
153 static inline void
154 histogram_modify_val (Histogram    *hist,
155                       const gint32 *src,
156                       gint          diff,
157                       gint          n_color_components,
158                       gboolean      has_alpha)
159 {
160   gint alpha = diff;
161   gint c;
162 
163   if (has_alpha)
164     alpha *= hist->alpha_values[src[n_color_components]];
165 
166   for (c = 0; c < n_color_components; c++)
167     {
168       HistogramComponent *comp = &hist->components[c];
169       gint                bin  = src[c];
170 
171       comp->bins[bin] += alpha;
172 
173       /* this is shorthand for:
174        *
175        *   if (bin <= comp->last_median)
176        *     comp->last_median_sum += alpha;
177        *
178        * but with a notable speed boost.
179        */
180       comp->last_median_sum += (bin <= comp->last_median) * alpha;
181     }
182 
183   if (has_alpha)
184     {
185       HistogramComponent *comp = &hist->components[n_color_components];
186       gint                bin  = src[n_color_components];
187 
188       comp->bins[bin] += diff;
189 
190       comp->last_median_sum += (bin <= comp->last_median) * diff;
191     }
192 
193   hist->count += alpha;
194 }
195 
196 static inline void
197 histogram_modify_vals (Histogram    *hist,
198                        const gint32 *src,
199                        gint          stride,
200                        gint          xmin,
201                        gint          ymin,
202                        gint          xmax,
203                        gint          ymax,
204                        gint          diff)
205 {
206   gint     n_components       = hist->n_components;
207   gint     n_color_components = hist->n_color_components;
208   gboolean has_alpha          = n_color_components < n_components;
209   gint x;
210   gint y;
211 
212   if (xmin > xmax || ymin > ymax)
213     return;
214 
215   src += ymin * stride + xmin * n_components;
216 
217   if (n_color_components == 3)
218     {
219       if (has_alpha)
220         {
221           for (y = ymin; y <= ymax; y++, src += stride)
222             {
223               const gint32 *pixel = src;
224 
225               for (x = xmin; x <= xmax; x++, pixel += n_components)
226                 {
227                   histogram_modify_val (hist, pixel, diff, 3, TRUE);
228                 }
229             }
230         }
231       else
232         {
233           for (y = ymin; y <= ymax; y++, src += stride)
234             {
235               const gint32 *pixel = src;
236 
237               for (x = xmin; x <= xmax; x++, pixel += n_components)
238                 {
239                   histogram_modify_val (hist, pixel, diff, 3, FALSE);
240                 }
241             }
242         }
243     }
244   else
245     {
246       if (has_alpha)
247         {
248           for (y = ymin; y <= ymax; y++, src += stride)
249             {
250               const gint32 *pixel = src;
251 
252               for (x = xmin; x <= xmax; x++, pixel += n_components)
253                 {
254                   histogram_modify_val (hist, pixel, diff, 1, TRUE);
255                 }
256             }
257         }
258       else
259         {
260           for (y = ymin; y <= ymax; y++, src += stride)
261             {
262               const gint32 *pixel = src;
263 
264               for (x = xmin; x <= xmax; x++, pixel += n_components)
265                 {
266                   histogram_modify_val (hist, pixel, diff, 1, FALSE);
267                 }
268             }
269         }
270     }
271 }
272 
273 static inline void
274 histogram_update (Histogram                  *hist,
275                   const gint32               *src,
276                   gint                        stride,
277                   GeglMedianBlurNeighborhood  neighborhood,
278                   gint                        radius,
279                   const gint                 *neighborhood_outline,
280                   Direction                   dir)
281 {
282   gint i;
283 
284   switch (neighborhood)
285     {
286     case GEGL_MEDIAN_BLUR_NEIGHBORHOOD_SQUARE:
287       switch (dir)
288         {
289           case LEFT_TO_RIGHT:
290             histogram_modify_vals (hist, src, stride,
291                                    -radius - 1, -radius,
292                                    -radius - 1, +radius,
293                                    -1);
294 
295             histogram_modify_vals (hist, src, stride,
296                                    +radius, -radius,
297                                    +radius, +radius,
298                                    +1);
299             break;
300 
301           case RIGHT_TO_LEFT:
302             histogram_modify_vals (hist, src, stride,
303                                    +radius + 1, -radius,
304                                    +radius + 1, +radius,
305                                    -1);
306 
307             histogram_modify_vals (hist, src, stride,
308                                    -radius, -radius,
309                                    -radius, +radius,
310                                    +1);
311             break;
312 
313           case TOP_TO_BOTTOM:
314             histogram_modify_vals (hist, src, stride,
315                                    -radius, -radius - 1,
316                                    +radius, -radius - 1,
317                                    -1);
318 
319             histogram_modify_vals (hist, src, stride,
320                                    -radius, +radius,
321                                    +radius, +radius,
322                                    +1);
323             break;
324         }
325       break;
326 
327     default:
328       switch (dir)
329         {
330           case LEFT_TO_RIGHT:
331             for (i = 0; i < radius; i++)
332               {
333                 histogram_modify_vals (hist, src, stride,
334                                        -i - 1, -neighborhood_outline[i],
335                                        -i - 1, -neighborhood_outline[i + 1] - 1,
336                                        -1);
337                 histogram_modify_vals (hist, src, stride,
338                                        -i - 1, +neighborhood_outline[i + 1] + 1,
339                                        -i - 1, +neighborhood_outline[i],
340                                        -1);
341 
342                 histogram_modify_vals (hist, src, stride,
343                                        +i, -neighborhood_outline[i],
344                                        +i, -neighborhood_outline[i + 1] - 1,
345                                        +1);
346                 histogram_modify_vals (hist, src, stride,
347                                        +i, +neighborhood_outline[i + 1] + 1,
348                                        +i, +neighborhood_outline[i],
349                                        +1);
350               }
351 
352             histogram_modify_vals (hist, src, stride,
353                                    -i - 1, -neighborhood_outline[i],
354                                    -i - 1, +neighborhood_outline[i],
355                                    -1);
356 
357             histogram_modify_vals (hist, src, stride,
358                                    +i, -neighborhood_outline[i],
359                                    +i, +neighborhood_outline[i],
360                                    +1);
361 
362             break;
363 
364           case RIGHT_TO_LEFT:
365             for (i = 0; i < radius; i++)
366               {
367                 histogram_modify_vals (hist, src, stride,
368                                        +i + 1, -neighborhood_outline[i],
369                                        +i + 1, -neighborhood_outline[i + 1] - 1,
370                                        -1);
371                 histogram_modify_vals (hist, src, stride,
372                                        +i + 1, +neighborhood_outline[i + 1] + 1,
373                                        +i + 1, +neighborhood_outline[i],
374                                        -1);
375 
376                 histogram_modify_vals (hist, src, stride,
377                                        -i, -neighborhood_outline[i],
378                                        -i, -neighborhood_outline[i + 1] - 1,
379                                        +1);
380                 histogram_modify_vals (hist, src, stride,
381                                        -i, +neighborhood_outline[i + 1] + 1,
382                                        -i, +neighborhood_outline[i],
383                                        +1);
384               }
385 
386             histogram_modify_vals (hist, src, stride,
387                                    +i + 1, -neighborhood_outline[i],
388                                    +i + 1, +neighborhood_outline[i],
389                                    -1);
390 
391             histogram_modify_vals (hist, src, stride,
392                                    -i, -neighborhood_outline[i],
393                                    -i, +neighborhood_outline[i],
394                                    +1);
395 
396             break;
397 
398           case TOP_TO_BOTTOM:
399             for (i = 0; i < radius; i++)
400               {
401                 histogram_modify_vals (hist, src, stride,
402                                        -neighborhood_outline[i],         -i - 1,
403                                        -neighborhood_outline[i + 1] - 1, -i - 1,
404                                        -1);
405                 histogram_modify_vals (hist, src, stride,
406                                        +neighborhood_outline[i + 1] + 1, -i - 1,
407                                        +neighborhood_outline[i],         -i - 1,
408                                        -1);
409 
410                 histogram_modify_vals (hist, src, stride,
411                                        -neighborhood_outline[i],         +i,
412                                        -neighborhood_outline[i + 1] - 1, +i,
413                                        +1);
414                 histogram_modify_vals (hist, src, stride,
415                                        +neighborhood_outline[i + 1] + 1, +i,
416                                        +neighborhood_outline[i],         +i,
417                                        +1);
418               }
419 
420             histogram_modify_vals (hist, src, stride,
421                                    -neighborhood_outline[i], -i - 1,
422                                    +neighborhood_outline[i], -i - 1,
423                                    -1);
424 
425             histogram_modify_vals (hist, src, stride,
426                                    -neighborhood_outline[i], +i,
427                                    +neighborhood_outline[i], +i,
428                                    +1);
429 
430             break;
431         }
432       break;
433     }
434 }
435 
436 static void
437 init_neighborhood_outline (GeglMedianBlurNeighborhood  neighborhood,
438                            gint                        radius,
439                            gint                       *neighborhood_outline)
440 {
441   gint i;
442 
443   for (i = 0; i <= radius; i++)
444     {
445       switch (neighborhood)
446         {
447         case GEGL_MEDIAN_BLUR_NEIGHBORHOOD_SQUARE:
448           neighborhood_outline[i] = radius;
449           break;
450 
451         case GEGL_MEDIAN_BLUR_NEIGHBORHOOD_CIRCLE:
452           neighborhood_outline[i] =
453             (gint) sqrt ((radius + .5) * (radius + .5) - i * i);
454           break;
455 
456         case GEGL_MEDIAN_BLUR_NEIGHBORHOOD_DIAMOND:
457           neighborhood_outline[i] = radius - i;
458           break;
459         }
460     }
461 }
462 
463 static void
464 sort_input_values (InputValue **values,
465                    InputValue **scratch,
466                    gint         n_values)
467 {
468   const InputValue *in     = *values;
469   InputValue       *out    = *scratch;
470   gint              size;
471 
472   for (size = 1; size < n_values; size *= 2)
473     {
474       InputValue *temp;
475       gint        i;
476 
477       for (i = 0; i < n_values; )
478         {
479           gint l     = i;
480           gint l_end = MIN (l + size, n_values);
481           gint r     = l_end;
482           gint r_end = MIN (r + size, n_values);
483 
484           while (l < l_end && r < r_end)
485             {
486               if (in[r].value < in[l].value)
487                 out[i] = in[r++];
488               else
489                 out[i] = in[l++];
490 
491               i++;
492             }
493 
494           memcpy (&out[i], &in[l], (l_end - l) * sizeof (InputValue));
495           i += l_end - l;
496 
497           memcpy (&out[i], &in[r], (r_end - r) * sizeof (InputValue));
498           i += r_end - r;
499         }
500 
501       temp = (InputValue *) in;
502       in   = out;
503       out  = temp;
504     }
505 
506   *values  = (InputValue *) in;
507   *scratch = out;
508 }
509 
510 static void
511 convert_values_to_bins (Histogram *hist,
512                         gint32    *src,
513                         gint       n_pixels,
514                         gboolean   quantize)
515 {
516   gint     n_components       = hist->n_components;
517   gint     n_color_components = hist->n_color_components;
518   gboolean has_alpha          = n_color_components < n_components;
519   gint     i;
520   gint     c;
521 
522   if (quantize)
523     {
524       for (c = 0; c < n_components; c++)
525         {
526           hist->components[c].bins       = g_new0 (gint, DEFAULT_N_BINS);
527           hist->components[c].bin_values = default_bin_values;
528         }
529 
530       while (n_pixels--)
531         {
532           for (c = 0; c < n_components; c++)
533             {
534               gfloat value = ((gfloat *) src)[c];
535               gint   bin;
536 
537               bin = floorf (SAFE_CLAMP (value, 0.0f, 1.0f) * (DEFAULT_N_BINS - 1) + 0.5f);
538 
539               src[c] = bin;
540             }
541 
542           src += n_components;
543         }
544 
545       hist->alpha_values = default_alpha_values;
546     }
547   else
548     {
549       InputValue *values       = g_new (InputValue, n_pixels);
550       InputValue *scratch      = g_new (InputValue, n_pixels);
551       gint       *alpha_values = NULL;
552 
553       if (has_alpha)
554         {
555           alpha_values = g_new (gint, n_pixels);
556 
557           hist->alpha_values = alpha_values;
558         }
559 
560       for (i = 0; i < n_pixels; i++)
561         {
562           values[i].value = ((gfloat *) src)[i * n_components];
563           values[i].index = i;
564         }
565 
566       for (c = 0; c < n_components; c++)
567         {
568           gint    bin = 0;
569           gfloat  prev_value;
570           gfloat *bin_values;
571 
572           sort_input_values (&values, &scratch, n_pixels);
573 
574           prev_value = values[0].value;
575 
576           bin_values    = g_new (gfloat, n_pixels);
577           bin_values[0] = prev_value;
578           if (c == n_color_components)
579             {
580               alpha_values[0] = floorf (SAFE_CLAMP (prev_value, 0.0f, 1.0f) *
581                                         (gfloat) (1 << 10) + 0.5f);
582             }
583 
584           for (i = 0; i < n_pixels; i++)
585             {
586               gint32 *p = &src[values[i].index * n_components + c];
587 
588               if (values[i].value != prev_value)
589                 {
590                   bin++;
591                   prev_value      = values[i].value;
592                   bin_values[bin] = prev_value;
593                   if (c == n_color_components)
594                     {
595                       alpha_values[bin] = floorf (SAFE_CLAMP (prev_value, 0.0f, 1.0f) *
596                                                   (gfloat) (1 << 10) + 0.5f);
597                     }
598                 }
599 
600               *p = bin;
601               if (c < n_components - 1)
602                 values[i].value = ((gfloat *) p)[1];
603             }
604 
605           hist->components[c].bins       = g_new0 (gint, bin + 1);
606           hist->components[c].bin_values = bin_values;
607         }
608 
609       g_free (scratch);
610       g_free (values);
611     }
612 }
613 
614 static void
615 prepare (GeglOperation *operation)
616 {
617   GeglOperationAreaFilter *area      = GEGL_OPERATION_AREA_FILTER (operation);
618   GeglProperties          *o         = GEGL_PROPERTIES (operation);
619   const Babl              *in_format = gegl_operation_get_source_format (operation, "input");
620   const Babl              *format    = NULL;
621   UserData                *data;
622   gint                     radius    = abs (o->radius);
623 
624   area->left   =
625   area->right  =
626   area->top    =
627   area->bottom = radius;
628 
629   if (! o->user_data)
630     o->user_data = g_slice_new0 (UserData);
631 
632   data                       = o->user_data;
633   data->quantize             = ! o->high_precision;
634   data->neighborhood_outline = g_renew (gint, data->neighborhood_outline,
635                                         radius + 1);
636   init_neighborhood_outline (o->neighborhood, radius,
637                              data->neighborhood_outline);
638 
639   if (in_format)
640     {
641       const Babl *model = babl_format_get_model (in_format);
642 
643       if (o->high_precision)
644         {
645           if (babl_model_is (model, "Y"))
646             format = babl_format_with_space ("Y float", in_format);
647           else if (babl_model_is (model, "Y'"))
648             format = babl_format_with_space ("Y' float", in_format);
649           else if (babl_model_is (model, "YA") || babl_model_is (model, "YaA"))
650             format = babl_format_with_space ("YA float", in_format);
651           else if (babl_model_is (model, "Y'A") || babl_model_is (model, "Y'aA"))
652             format = babl_format_with_space ("Y'A float", in_format);
653           else if (babl_model_is (model, "RGB"))
654             format = babl_format_with_space ("RGB float", in_format);
655           else if (babl_model_is (model, "R'G'B'"))
656             format = babl_format_with_space ("R'G'B' float", in_format);
657           else if (babl_model_is (model, "RGBA") || babl_model_is (model, "RaGaBaA"))
658             format = babl_format_with_space ("RGBA float", in_format);
659           else if (babl_model_is (model, "R'G'B'A") || babl_model_is (model, "R'aG'aB'aA"))
660             format = babl_format_with_space ("R'G'B'A float", in_format);
661 
662           if (format)
663             {
664               gint n_components = babl_format_get_n_components (in_format);
665               gint i;
666 
667               data->quantize = TRUE;
668 
669               for (i = 0; i < n_components; i++)
670                 {
671                   if (babl_format_get_type (in_format, i) != babl_type ("u8"))
672                     {
673                       data->quantize = FALSE;
674                       break;
675                     }
676                 }
677             }
678         }
679       else
680         {
681           if (babl_model_is (model, "Y") || babl_model_is (model, "Y'"))
682             format = babl_format_with_space ("Y' float", in_format);
683           else if (babl_model_is (model, "YA")  || babl_model_is (model, "YaA") ||
684                    babl_model_is (model, "Y'A") || babl_model_is (model, "Y'aA"))
685             format = babl_format_with_space ("Y'A float", in_format);
686           else if (babl_model_is (model, "RGB") || babl_model_is (model, "R'G'B'"))
687             format = babl_format_with_space ("R'G'B' float", in_format);
688           else if (babl_model_is (model, "RGBA")    || babl_model_is (model, "RaGaBaA") ||
689                    babl_model_is (model, "R'G'B'A") || babl_model_is (model, "R'aG'aB'aA"))
690             format = babl_format_with_space ("R'G'B'A float", in_format);
691         }
692 
693       if (format == NULL)
694         {
695           if (babl_format_has_alpha (in_format))
696             format = babl_format_with_space ("R'G'B'A float", in_format);
697           else
698             format = babl_format_with_space ("R'G'B' float", in_format);
699         }
700     }
701   else
702     {
703       if (o->high_precision)
704         format = babl_format_with_space ("RGBA float", in_format);
705       else
706         format = babl_format_with_space ("R'G'B'A float", in_format);
707     }
708 
709   if (data->quantize && ! g_atomic_int_get (&default_values_initialized))
710     {
711       gint i;
712 
713       for (i = 0; i < DEFAULT_N_BINS; i++)
714         {
715           default_bin_values[i]   = (gfloat) i / (gfloat) (DEFAULT_N_BINS - 1);
716           default_alpha_values[i] = i;
717         }
718 
719       g_atomic_int_set (&default_values_initialized, TRUE);
720     }
721 
722   gegl_operation_set_format (operation, "input", format);
723   gegl_operation_set_format (operation, "output", format);
724 }
725 
726 static GeglRectangle
727 get_bounding_box (GeglOperation *operation)
728 {
729   GeglProperties *o = GEGL_PROPERTIES (operation);
730 
731   if (o->abyss_policy != GEGL_MEDIAN_BLUR_ABYSS_NONE)
732     {
733       GeglRectangle  result = { 0, 0, 0, 0 };
734       GeglRectangle *in_rect;
735 
736       in_rect = gegl_operation_source_get_bounding_box (operation, "input");
737 
738       if (in_rect)
739         result = *in_rect;
740 
741       return result;
742     }
743 
744   return GEGL_OPERATION_CLASS (gegl_op_parent_class)->get_bounding_box (operation);
745 }
746 
747 static GeglAbyssPolicy
748 get_abyss_policy (GeglOperation *operation,
749                   const gchar   *input_pad)
750 {
751   GeglProperties *o = GEGL_PROPERTIES (operation);
752 
753   switch (o->abyss_policy)
754     {
755     case GEGL_MEDIAN_BLUR_ABYSS_NONE:  return GEGL_ABYSS_NONE;
756     case GEGL_MEDIAN_BLUR_ABYSS_CLAMP: return GEGL_ABYSS_CLAMP;
757     }
758 
759   g_return_val_if_reached (GEGL_ABYSS_NONE);
760 }
761 
762 static gboolean
763 process (GeglOperation       *operation,
764          GeglBuffer          *input,
765          GeglBuffer          *output,
766          const GeglRectangle *roi,
767          gint                 level)
768 {
769   GeglProperties *o    = GEGL_PROPERTIES (operation);
770   UserData       *data = o->user_data;
771 
772   gint            radius               = abs (o->radius);
773   gdouble         percentile           = o->percentile       / 100.0;
774   gdouble         alpha_percentile     = o->alpha_percentile / 100.0;
775   const gint     *neighborhood_outline = data->neighborhood_outline;
776 
777   const Babl     *format               = gegl_operation_get_format (operation, "input");
778   gint            n_components         = babl_format_get_n_components (format);
779   gint            n_color_components   = n_components;
780   gboolean        has_alpha            = babl_format_has_alpha (format);
781 
782   G_STATIC_ASSERT (sizeof (gint32) == sizeof (gfloat));
783   gint32         *src_buf;
784   gfloat         *dst_buf;
785   GeglRectangle   src_rect;
786   gint            src_stride;
787   gint            dst_stride;
788   gint            n_src_pixels;
789   gint            n_dst_pixels;
790 
791   Histogram      *hist;
792 
793   const gint32   *src;
794   gfloat         *dst;
795   gint            dst_x, dst_y;
796   Direction       dir;
797 
798   gint            i;
799   gint            c;
800 
801   if (o->radius < 0)
802     {
803       percentile       = 1.0 - percentile;
804       alpha_percentile = 1.0 - alpha_percentile;
805     }
806 
807   if (! data->quantize &&
808       (roi->width > MAX_CHUNK_WIDTH || roi->height > MAX_CHUNK_HEIGHT))
809     {
810       gint n_x = (roi->width  + MAX_CHUNK_WIDTH  - 1) / MAX_CHUNK_WIDTH;
811       gint n_y = (roi->height + MAX_CHUNK_HEIGHT - 1) / MAX_CHUNK_HEIGHT;
812       gint x;
813       gint y;
814 
815       for (y = 0; y < n_y; y++)
816         {
817           for (x = 0; x < n_x; x++)
818             {
819               GeglRectangle chunk;
820 
821               chunk.x      = roi->x + roi->width  * x       / n_x;
822               chunk.y      = roi->y + roi->height * y       / n_y;
823               chunk.width  = roi->x + roi->width  * (x + 1) / n_x - chunk.x;
824               chunk.height = roi->y + roi->height * (y + 1) / n_y - chunk.y;
825 
826               if (! process (operation, input, output, &chunk, level))
827                 return FALSE;
828             }
829         }
830 
831       return TRUE;
832     }
833 
834   if (has_alpha)
835     n_color_components--;
836 
837   g_return_val_if_fail (n_color_components == 1 || n_color_components == 3, FALSE);
838 
839   hist = g_slice_new0 (Histogram);
840 
841   hist->n_components       = n_components;
842   hist->n_color_components = n_color_components;
843 
844   src_rect     = gegl_operation_get_required_for_output (operation, "input", roi);
845   src_stride   = src_rect.width * n_components;
846   dst_stride   = roi->width * n_components;
847   n_src_pixels = src_rect.width * src_rect.height;
848   n_dst_pixels = roi->width * roi->height;
849   src_buf = g_new (gint32, n_src_pixels * n_components);
850   dst_buf = g_new (gfloat, n_dst_pixels * n_components);
851 
852   gegl_buffer_get (input, &src_rect, 1.0, format, src_buf,
853                    GEGL_AUTO_ROWSTRIDE, get_abyss_policy (operation, "input"));
854   convert_values_to_bins (hist, src_buf, n_src_pixels, data->quantize);
855 
856   src = src_buf + radius * (src_rect.width + 1) * n_components;
857   dst = dst_buf;
858 
859   /* compute the first window */
860 
861   for (i = -radius; i <= radius; i++)
862     {
863       histogram_modify_vals (hist, src, src_stride,
864                              i, -neighborhood_outline[abs (i)],
865                              i, +neighborhood_outline[abs (i)],
866                              +1);
867 
868       hist->size += 2 * neighborhood_outline[abs (i)] + 1;
869     }
870 
871   for (c = 0; c < n_color_components; c++)
872     dst[c] = histogram_get_median (hist, c, percentile);
873   if (has_alpha)
874     dst[c] = histogram_get_median (hist, c, alpha_percentile);
875 
876   dst_x = 0;
877   dst_y = 0;
878 
879   n_dst_pixels--;
880   dir = LEFT_TO_RIGHT;
881 
882   while (n_dst_pixels--)
883     {
884       /* move the src coords based on current direction and positions */
885       if (dir == LEFT_TO_RIGHT)
886         {
887           if (dst_x != roi->width - 1)
888             {
889               dst_x++;
890               src += n_components;
891               dst += n_components;
892             }
893           else
894             {
895               dst_y++;
896               src += src_stride;
897               dst += dst_stride;
898               dir = TOP_TO_BOTTOM;
899             }
900         }
901       else if (dir == TOP_TO_BOTTOM)
902         {
903           if (dst_x == 0)
904             {
905               dst_x++;
906               src += n_components;
907               dst += n_components;
908               dir = LEFT_TO_RIGHT;
909             }
910           else
911             {
912               dst_x--;
913               src -= n_components;
914               dst -= n_components;
915               dir = RIGHT_TO_LEFT;
916             }
917         }
918       else if (dir == RIGHT_TO_LEFT)
919         {
920           if (dst_x != 0)
921             {
922               dst_x--;
923               src -= n_components;
924               dst -= n_components;
925             }
926           else
927             {
928               dst_y++;
929               src += src_stride;
930               dst += dst_stride;
931               dir = TOP_TO_BOTTOM;
932             }
933         }
934 
935       histogram_update (hist, src, src_stride,
936                         o->neighborhood, radius, neighborhood_outline,
937                         dir);
938 
939       for (c = 0; c < n_color_components; c++)
940         dst[c] = histogram_get_median (hist, c, percentile);
941       if (has_alpha)
942         dst[c] = histogram_get_median (hist, c, alpha_percentile);
943     }
944 
945   gegl_buffer_set (output, roi, 0, format, dst_buf, GEGL_AUTO_ROWSTRIDE);
946 
947   for (c = 0; c < n_components; c++)
948     {
949       g_free (hist->components[c].bins);
950 
951       if (! data->quantize)
952         g_free (hist->components[c].bin_values);
953     }
954 
955   if (! data->quantize && has_alpha)
956     g_free (hist->alpha_values);
957 
958   g_slice_free (Histogram, hist);
959   g_free (dst_buf);
960   g_free (src_buf);
961 
962   return TRUE;
963 }
964 
965 static void
966 finalize (GObject *object)
967 {
968   GeglOperation  *op = (void*) object;
969   GeglProperties *o  = GEGL_PROPERTIES (op);
970 
971   if (o->user_data)
972     {
973       UserData *data = o->user_data;
974 
975       g_free (data->neighborhood_outline);
976       g_slice_free (UserData, data);
977     }
978 
979   G_OBJECT_CLASS (gegl_op_parent_class)->finalize (object);
980 }
981 
982 static void
983 gegl_op_class_init (GeglOpClass *klass)
984 {
985   GObjectClass                 *object_class;
986   GeglOperationClass           *operation_class;
987   GeglOperationFilterClass     *filter_class;
988   GeglOperationAreaFilterClass *area_class;
989 
990   object_class    = G_OBJECT_CLASS (klass);
991   operation_class = GEGL_OPERATION_CLASS (klass);
992   filter_class    = GEGL_OPERATION_FILTER_CLASS (klass);
993   area_class      = GEGL_OPERATION_AREA_FILTER_CLASS (klass);
994 
995   object_class->finalize            = finalize;
996   filter_class->process             = process;
997   operation_class->prepare          = prepare;
998   operation_class->get_bounding_box = get_bounding_box;
999   area_class->get_abyss_policy      = get_abyss_policy;
1000 
1001   gegl_operation_class_set_keys (operation_class,
1002     "name",           "gegl:median-blur",
1003     "title",          _("Median Blur"),
1004     "categories",     "blur",
1005     "reference-hash", "1865918d2f3b95690359534bbd58b513",
1006     "description",    _("Blur resulting from computing the median "
1007                         "color in the neighborhood of each pixel."),
1008     NULL);
1009 }
1010 
1011 #endif
1012