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 2020 Thomas Manni <thomas.manni@free.fr>
17  *
18  * expected input buffers:
19  *   - input : current selection
20  *   - aux   : color image
21  *   - aux2  : scribbles
22  */
23 
24 #include "config.h"
25 #include <glib/gi18n-lib.h>
26 
27 #ifdef GEGL_PROPERTIES
28 
29 enum_start (gegl_paint_select_mode_type)
30   enum_value (GEGL_PAINT_SELECT_MODE_ADD,      "add",      N_("Add"))
31   enum_value (GEGL_PAINT_SELECT_MODE_SUBTRACT, "subtract", N_("Subtract"))
32 enum_end (GeglPaintSelectModeType)
33 
34 property_enum (mode, _("Mode"),
35                GeglPaintSelectModeType, gegl_paint_select_mode_type,
36                GEGL_PAINT_SELECT_MODE_ADD)
37     description (_("Either to add to or subtract from the mask"))
38 
39 property_boolean (use_local_region, _("Use local region"), FALSE)
40     description (_("Perform graphcut in a local region"))
41 
42 property_int (region_x, _("region-x"), 0)
43     value_range (0, G_MAXINT)
44 
45 property_int (region_y, _("region-y"), 0)
46     value_range (0, G_MAXINT)
47 
48 property_int (region_width, _("region-width"), 0)
49     value_range (0, G_MAXINT)
50 
51 property_int (region_height, _("region-height"), 0)
52     value_range (0, G_MAXINT)
53 
54 #else
55 
56 #define GEGL_OP_COMPOSER3
57 #define GEGL_OP_NAME      paint_select
58 #define GEGL_OP_C_SOURCE  paint-select.cc
59 
60 #include <maxflow/graph.h>
61 #include <math.h>
62 #include "gegl-op.h"
63 
64 #define SELECTION_FORMAT "Y float"
65 #define SCRIBBLES_FORMAT "Y float"
66 #define COLORS_FORMAT    "R'G'B' float"
67 
68 #define POW2(x) ((x)*(x))
69 #define BIG_CAPACITY 100.f
70 #define EPSILON 0.05f
71 #define N_GLOBAL_SAMPLES 1200
72 #define N_BINS 64
73 #define LOCAL_REGION_DILATE 40
74 
75 #define FG_MASK      1.f
76 #define BG_MASK      0.f
77 #define FG_SCRIBBLE  1.f
78 #define BG_SCRIBBLE  0.f
79 
80 typedef maxflow::Graph<gfloat,gfloat,gfloat> GraphType;
81 typedef int node_id;
82 
83 typedef struct
84 {
85   gfloat  rgb[3];
86   gint    x;
87   gint    y;
88 } ColorsSample;
89 
90 typedef struct
91 {
92   GArray  *samples;
93   gfloat   bins[N_BINS][N_BINS][N_BINS];
94 } ColorsModel;
95 
96 typedef enum
97 {
98   SEED_NONE,
99   SEED_SOURCE,
100   SEED_SINK,
101 } SeedType;
102 
103 typedef struct
104 {
105   GeglPaintSelectModeType  mode;
106   GeglRectangle            roi;
107   GeglRectangle            extent;
108 
109   gfloat      *mask;        /* selection mask */
110   gfloat      *colors;      /* input image    */
111   gfloat      *scribbles;   /* user scribbles */
112 } PaintSelect;
113 
114 typedef struct
115 {
116   ColorsModel  *fg_colors;
117   ColorsModel  *bg_colors;
118 } PaintSelectPrivate;
119 
120 typedef struct
121 {
122   gfloat        mask_value_seed; /* Value of the mask where seeds are needed. */
123   SeedType      mask_seed_type;  /* Corresponding seed type.                  */
124   ColorsModel  *source_model;    /* Colors model used to compute terminal     */
125   ColorsModel  *sink_model;      /*  links weights.                           */
126   ColorsModel  *local_colors;
127 
128   gboolean      boundary_top;    /* Seed type added to the local region       */
129   gboolean      boundary_left;   /*  boundaries to avoid the selection to     */
130   gboolean      boundary_right;  /*  align with them.                         */
131   gboolean      boundary_bottom;
132   SeedType      boundary_seed_type;
133 } PaintSelectContext;
134 
135 
136 /* colors model */
137 
138 static ColorsModel *
139 colors_model_new_global (gfloat  *pixels,
140                          gfloat  *mask,
141                          gint     width,
142                          gint     height,
143                          gfloat   mask_value)
144 {
145   ColorsModel *model = g_slice_new0 (ColorsModel);
146   GRand  *gr = g_rand_new_with_seed (0);
147   gint    i = 0;
148 
149   model->samples = g_array_sized_new (FALSE, FALSE, sizeof (ColorsSample), N_GLOBAL_SAMPLES);
150 
151   while (i < N_GLOBAL_SAMPLES)
152     {
153       ColorsSample  sample;
154       gint          m_offset;
155 
156       sample.x = g_rand_int_range (gr, 0, width);
157       sample.y = g_rand_int_range (gr, 0, height);
158       m_offset = sample.x + sample.y * width;
159 
160       if (mask[m_offset] == mask_value)
161         {
162           gint p_offset = m_offset * 3;
163           gint b1, b2, b3;
164 
165           sample.rgb[0] = CLAMP (pixels[p_offset], 0.f, 1.f);
166           sample.rgb[1] = CLAMP (pixels[p_offset+1], 0.f, 1.f);
167           sample.rgb[2] = CLAMP (pixels[p_offset+2], 0.f, 1.f);
168 
169           b1 = (gint) (sample.rgb[0] * (N_BINS - 1));
170           b2 = (gint) (sample.rgb[1] * (N_BINS - 1));
171           b3 = (gint) (sample.rgb[2] * (N_BINS - 1));
172 
173           model->bins[b1][b2][b3]++;
174           g_array_append_val (model->samples, sample);
175           i++;
176         }
177     }
178 
179   g_rand_free (gr);
180   return model;
181 }
182 
183 static void
184 colors_model_update_global (ColorsModel  *model,
185                             gfloat       *pixels,
186                             gfloat       *mask,
187                             gint          width,
188                             gint          height,
189                             gfloat        mask_value)
190 {
191   GRand  *gr = g_rand_new_with_seed (0);
192   guint   i;
193 
194   for (i = 0; i < model->samples->len; i++)
195     {
196       ColorsSample  *sample = &g_array_index (model->samples, ColorsSample, i);
197       gint           m_offset = sample->x + sample->y * width;
198 
199       if (mask[m_offset] != mask_value)
200         {
201           gboolean changed = FALSE;
202           gint b1, b2, b3;
203 
204           b1 = (gint) (sample->rgb[0] * (N_BINS - 1));
205           b2 = (gint) (sample->rgb[1] * (N_BINS - 1));
206           b3 = (gint) (sample->rgb[2] * (N_BINS - 1));
207 
208           model->bins[b1][b2][b3]--;
209 
210           while (! changed)
211             {
212               sample->x = g_rand_int_range (gr, 0, width);
213               sample->y = g_rand_int_range (gr, 0, height);
214               m_offset = sample->x + sample->y * width;
215 
216               if (mask[m_offset] == mask_value)
217                 {
218                   gint p_offset = m_offset * 3;
219 
220                   sample->rgb[0] = pixels[p_offset];
221                   sample->rgb[1] = pixels[p_offset+1];
222                   sample->rgb[2] = pixels[p_offset+2];
223 
224                   b1 = (gint) (CLAMP(sample->rgb[0], 0.f, 1.f) * (N_BINS - 1));
225                   b2 = (gint) (CLAMP(sample->rgb[1], 0.f, 1.f) * (N_BINS - 1));
226                   b3 = (gint) (CLAMP(sample->rgb[2], 0.f, 1.f) * (N_BINS - 1));
227 
228                   model->bins[b1][b2][b3]++;
229                   changed = TRUE;
230                 }
231             }
232         }
233     }
234 
235   g_rand_free (gr);
236 }
237 
238 static ColorsModel *
239 colors_model_new_local (gfloat         *pixels,
240                         gfloat         *mask,
241                         gfloat         *scribbles,
242                         gint            width,
243                         gint            height,
244                         GeglRectangle  *region,
245                         gfloat          mask_value,
246                         gfloat          scribble_value)
247 {
248   gint  x, y;
249   ColorsModel *model = g_slice_new0 (ColorsModel);
250   model->samples = g_array_sized_new (FALSE, FALSE, sizeof (ColorsSample), N_GLOBAL_SAMPLES);
251 
252   for (y = region->y; y < region->y + region->height; y++)
253     {
254       for (x = region->x; x < region->x + region->width; x++)
255         {
256           gint offset = x + y * width;
257 
258           if (scribbles[offset] == scribble_value ||
259               mask[offset] == mask_value)
260             {
261               ColorsSample  sample;
262               gint b1, b2, b3;
263               gint p_offset = offset * 3;
264 
265               sample.x = x;
266               sample.y = y;
267               sample.rgb[0] = CLAMP (pixels[p_offset], 0.f, 1.f);
268               sample.rgb[1] = CLAMP (pixels[p_offset+1], 0.f, 1.f);
269               sample.rgb[2] = CLAMP (pixels[p_offset+2], 0.f, 1.f);
270 
271               b1 = (gint) (sample.rgb[0] * (N_BINS - 1));
272               b2 = (gint) (sample.rgb[1] * (N_BINS - 1));
273               b3 = (gint) (sample.rgb[2] * (N_BINS - 1));
274 
275               model->bins[b1][b2][b3]++;
276               g_array_append_val (model->samples, sample);
277             }
278         }
279     }
280 
281   return model;
282 }
283 
284 static void
285 colors_model_free (ColorsModel  *model)
286 {
287   if (model->samples)
288     g_array_free (model->samples, TRUE);
289 
290   g_slice_free (ColorsModel, model);
291 }
292 
293 static gfloat
294 colors_model_get_likelyhood (ColorsModel  *model,
295                              gfloat       *color)
296 {
297   gint b1 = (gint) (CLAMP(color[0], 0.f, 1.f) * (N_BINS - 1));
298   gint b2 = (gint) (CLAMP(color[1], 0.f, 1.f) * (N_BINS - 1));
299   gint b3 = (gint) (CLAMP(color[2], 0.f, 1.f) * (N_BINS - 1));
300 
301   return model->bins[b1][b2][b3] / (gfloat) model->samples->len;
302 }
303 
304 /* fluctuations removal */
305 
306 static void
307 push_segment (GQueue *segment_queue,
308               gint    y,
309               gint    old_y,
310               gint    start,
311               gint    end,
312               gint    new_y,
313               gint    new_start,
314               gint    new_end)
315 {
316   /* To avoid excessive memory allocation (y, old_y, start, end) tuples are
317    * stored in interleaved format:
318    *
319    * [y1] [old_y1] [start1] [end1] [y2] [old_y2] [start2] [end2]
320    */
321 
322   if (new_y != old_y)
323     {
324       /* If the new segment's y-coordinate is different than the old (source)
325        * segment's y-coordinate, push the entire segment.
326        */
327       g_queue_push_tail (segment_queue, GINT_TO_POINTER (new_y));
328       g_queue_push_tail (segment_queue, GINT_TO_POINTER (y));
329       g_queue_push_tail (segment_queue, GINT_TO_POINTER (new_start));
330       g_queue_push_tail (segment_queue, GINT_TO_POINTER (new_end));
331     }
332   else
333     {
334       /* Otherwise, only push the set-difference between the new segment and
335        * the source segment (since we've already scanned the source segment.)
336        */
337       if (new_start < start)
338         {
339           g_queue_push_tail (segment_queue, GINT_TO_POINTER (new_y));
340           g_queue_push_tail (segment_queue, GINT_TO_POINTER (y));
341           g_queue_push_tail (segment_queue, GINT_TO_POINTER (new_start));
342           g_queue_push_tail (segment_queue, GINT_TO_POINTER (start));
343         }
344 
345       if (new_end > end)
346         {
347           g_queue_push_tail (segment_queue, GINT_TO_POINTER (new_y));
348           g_queue_push_tail (segment_queue, GINT_TO_POINTER (y));
349           g_queue_push_tail (segment_queue, GINT_TO_POINTER (end));
350           g_queue_push_tail (segment_queue, GINT_TO_POINTER (new_end));
351         }
352     }
353 }
354 
355 static void
356 pop_segment (GQueue *segment_queue,
357              gint   *y,
358              gint   *old_y,
359              gint   *start,
360              gint   *end)
361 {
362   *y     = GPOINTER_TO_INT (g_queue_pop_head (segment_queue));
363   *old_y = GPOINTER_TO_INT (g_queue_pop_head (segment_queue));
364   *start = GPOINTER_TO_INT (g_queue_pop_head (segment_queue));
365   *end   = GPOINTER_TO_INT (g_queue_pop_head (segment_queue));
366 }
367 
368 static gboolean
369 find_contiguous_segment (gfloat  *mask,
370                          gfloat  *diff,
371                          gint     width,
372                          gint     initial_x,
373                          gint     initial_y,
374                          gint    *start,
375                          gint    *end)
376 {
377   gint  offset = initial_x + initial_y * width;
378 
379   /* check the starting pixel */
380   if (diff[offset] == 0.f)
381     return FALSE;
382 
383   mask[offset] = 1.f;
384 
385   *start = initial_x - 1;
386 
387   while (*start >= 0)
388     {
389       offset = *start + initial_y * width;
390 
391       if (diff[offset] == 0.f)
392         break;
393 
394       mask[offset] = 1.f;
395 
396       (*start)--;
397     }
398 
399   *end = initial_x + 1;
400 
401   while (*end < width)
402     {
403       offset = *end + initial_y * width;
404 
405       if (diff[offset] == 0.f)
406         break;
407 
408       mask[offset] = 1.f;
409 
410       (*end)++;
411     }
412 
413   return TRUE;
414 }
415 
416 static void
417 paint_select_remove_fluctuations (gfloat  *mask,
418                                   gfloat  *diff,
419                                   gint     width,
420                                   gint     height,
421                                   gint     x,
422                                   gint     y)
423 {
424   gint                 old_y;
425   gint                 start, end;
426   gint                 new_start, new_end;
427   GQueue              *segment_queue;
428 
429   /* mask buffer will hold the result and need to be clean first */
430   memset (mask, 0.f, sizeof (gfloat) * (width * height));
431 
432   segment_queue = g_queue_new ();
433 
434   push_segment (segment_queue,
435                 y, /* dummy values: */ -1, 0, 0,
436                 y, x - 1, x + 1);
437 
438   do
439     {
440       pop_segment (segment_queue,
441                    &y, &old_y, &start, &end);
442 
443       for (x = start + 1; x < end; x++)
444         {
445           gfloat val = mask[x + y * width];
446 
447           if (val != 0.f)
448             {
449               /* If the current pixel is selected, then we've already visited
450                * the next pixel.  (Note that we assume that the maximal image
451                * width is sufficiently low that `x` won't overflow.)
452                */
453               x++;
454               continue;
455             }
456 
457           if (! find_contiguous_segment (mask, diff, width,
458                                          x, y,
459                                          &new_start, &new_end))
460             continue;
461 
462           /* We can skip directly to `new_end + 1` on the next iteration, since
463            * we've just selected all pixels in the range `[x, new_end)`, and
464            * the pixel at `new_end` is above threshold.  (Note that we assume
465            * that the maximal image width is sufficiently low that `x` won't
466            * overflow.)
467            */
468           x = new_end;
469 
470           if (y + 1 < height)
471             {
472               push_segment (segment_queue,
473                             y, old_y, start, end,
474                             y + 1, new_start, new_end);
475             }
476 
477           if (y - 1 >= 0)
478             {
479               push_segment (segment_queue,
480                             y, old_y, start, end,
481                             y - 1, new_start, new_end);
482             }
483 
484         }
485     }
486   while (! g_queue_is_empty (segment_queue));
487 
488   g_queue_free (segment_queue);
489 }
490 
491 /* graphcut
492  *
493  * SOURCE terminal always represents selected region */
494 
495 static inline gfloat
496 pixels_distance (gfloat  *p1,
497                  gfloat  *p2)
498 {
499   return sqrtf (POW2(p1[0] - p2[0]) +
500                 POW2(p1[1] - p2[1]) +
501                 POW2(p1[2] - p2[2]));
502 }
503 
504 static void
505 paint_select_compute_adjacent_costs (gfloat  *pixels,
506                                      gint     width,
507                                      gint     height,
508                                      gfloat  *h_costs,
509                                      gfloat  *v_costs,
510                                      gfloat  *mean)
511 {
512   gint  x, y;
513   gint  n_h_costs = (width - 1) * height;
514   gint  n_v_costs = width * (height - 1);
515   gfloat sum      = 0.f;
516 
517   /* horizontal */
518 
519   for (y = 0; y < height; y++)
520     {
521       for (x = 0; x < width - 1; x++)
522         {
523           gint cost_offset = x + y * (width - 1);
524           gint p1_offset   = (x + y * width) * 3;
525           gint p2_offset   = p1_offset + 3;
526 
527           gfloat d = pixels_distance (pixels + p1_offset, pixels + p2_offset);
528           h_costs[cost_offset] = d;
529           sum += d;
530         }
531     }
532 
533   /* vertical */
534 
535   for (x = 0; x < width; x++)
536     {
537       for (y = 0; y < height - 1; y++)
538         {
539           gint cost_offset = x + y * width;
540           gint p1_offset    = (x + y * width) * 3;
541           gint p2_offset    = p1_offset + width * 3;
542 
543           gfloat d = pixels_distance (pixels + p1_offset, pixels + p2_offset);
544           v_costs[cost_offset] = d;
545           sum += d;
546         }
547     }
548 
549   /* compute mean costs */
550 
551   *mean = sum / (gfloat) (n_h_costs + n_v_costs);
552 }
553 
554 static guint8 *
555 paint_select_compute_seeds_map (gfloat              *mask,
556                                 gfloat              *scribbles,
557                                 gint                 width,
558                                 gint                 height,
559                                 PaintSelectContext  *context)
560 {
561   guint8  *seeds;
562   gint     x, y;
563 
564   seeds = g_new (guint8, width * height);
565 
566   for (y = 0; y < height; y++)
567     {
568       for (x = 0; x < width; x++)
569         {
570           gint offset = x + y * width;
571 
572           gfloat *p_mask      = mask + offset;
573           gfloat *p_scribbles = scribbles + offset;
574           guint8 *p_seeds     = seeds + offset;
575 
576           *p_seeds = SEED_NONE;
577 
578           if (*p_mask == context->mask_value_seed)
579             {
580               *p_seeds = context->mask_seed_type;
581             }
582           else if (*p_scribbles == FG_SCRIBBLE)
583             {
584               *p_seeds = SEED_SOURCE;
585             }
586 
587           else if (*p_scribbles == BG_SCRIBBLE)
588             {
589               *p_seeds = SEED_SINK;
590             }
591         }
592     }
593 
594   /* put boundary seeds if needed */
595 
596   if (context->boundary_top)
597     {
598       y = 0;
599       for (x = 0; x < width; x++)
600         {
601           gint offset = x + y * width;
602           if (seeds[offset] == SEED_NONE)
603             seeds[offset] = context->boundary_seed_type;
604         }
605     }
606 
607   if (context->boundary_left)
608     {
609       x = 0;
610       for (y = 0; y < height; y++)
611         {
612           gint offset = x + y * width;
613           if (seeds[offset] == SEED_NONE)
614             seeds[offset] = context->boundary_seed_type;
615         }
616     }
617 
618   if (context->boundary_right)
619     {
620       x = width - 1;
621       for (y = 0; y < height; y++)
622         {
623           gint offset = x + y * width;
624           if (seeds[offset] == SEED_NONE)
625             seeds[offset] = context->boundary_seed_type;
626         }
627     }
628 
629   if (context->boundary_bottom)
630     {
631       y = height - 1;
632       for (x = 0; x < width; x++)
633         {
634           gint offset = x + y * width;
635           if (seeds[offset] == SEED_NONE)
636             seeds[offset] = context->boundary_seed_type;
637         }
638     }
639 
640   return seeds;
641 }
642 
643 static inline gboolean
644 paint_select_seed_is_boundary (guint8  *seeds,
645                                gint     width,
646                                gint     height,
647                                gint     x,
648                                gint     y)
649 {
650   gint offset = x + y * width;
651   gint neighbor_offset;
652   gint x2, y2;
653 
654   x2 = x - 1;
655   if (x2 >= 0)
656     {
657       neighbor_offset = offset - 1;
658       if (seeds[neighbor_offset] != seeds[offset])
659         return TRUE;
660     }
661 
662   x2 = x + 1;
663   if (x2 < width)
664     {
665       neighbor_offset = offset + 1;
666       if (seeds[neighbor_offset] != seeds[offset])
667         return TRUE;
668     }
669 
670   y2 = y - 1;
671   if (y2 >= 0)
672     {
673       neighbor_offset = offset - width;
674       if (seeds[neighbor_offset] != seeds[offset])
675         return TRUE;
676     }
677 
678   y2 = y + 1;
679   if (y2 < height)
680     {
681       neighbor_offset = offset + width;
682       if (seeds[neighbor_offset] != seeds[offset])
683         return TRUE;
684     }
685 
686   return FALSE;
687 }
688 
689 static void
690 paint_select_graph_init_nodes_and_tlinks (GraphType           *graph,
691                                           gfloat              *pixels,
692                                           guint8              *seeds,
693                                           gint                *nodes,
694                                           gint                 width,
695                                           gint                 height,
696                                           PaintSelectContext  *context)
697 {
698   gint  x, y;
699 
700   for (y = 0; y < height; y++)
701     {
702       for (x = 0; x < width; x++)
703         {
704           node_id  id      = -1;
705           gint     offset  = x + y * width;
706           guint8  *p_seeds = seeds + offset;
707           gint    *p_nodes = nodes + offset;
708 
709           if (*p_seeds == SEED_NONE)
710             {
711               gfloat  source_weight;
712               gfloat  sink_weight;
713               gfloat *color = pixels + offset * 3;
714 
715               id = graph->add_node();
716 
717               sink_weight   = - logf (colors_model_get_likelyhood (context->sink_model, color) + 0.0001f);
718               source_weight = - logf (colors_model_get_likelyhood (context->source_model, color) + 0.0001f);
719 
720               graph->add_tweights (id, source_weight, sink_weight);
721             }
722           else if (paint_select_seed_is_boundary (seeds, width, height, x, y))
723             {
724               id = graph->add_node();
725 
726               if (*p_seeds == SEED_SOURCE)
727                 {
728                   graph->add_tweights (id, BIG_CAPACITY, 0.f);
729                 }
730               else
731                 {
732                   graph->add_tweights (id, 0.f, BIG_CAPACITY);
733                 }
734             }
735 
736           *p_nodes = id;
737         }
738     }
739 }
740 
741 static void
742 paint_select_graph_init_nlinks (GraphType  *graph,
743                                 gint       *nodes,
744                                 gfloat     *h_costs,
745                                 gfloat     *v_costs,
746                                 gfloat      mean_costs,
747                                 gint        width,
748                                 gint        height)
749 {
750   gint x, y;
751 
752   /* horizontal */
753 
754   for (y = 0; y < height; y++)
755     {
756       for (x = 0; x < width - 1; x++)
757         {
758           node_id id1 = nodes[x + y * width];
759           node_id id2 = nodes[x + 1 + y * width];
760 
761           if (id1 != -1 && id2 != -1)
762             {
763               gint costs_off = x + y * (width - 1);
764               gfloat weight  = 60.f * mean_costs / (h_costs[costs_off] + EPSILON);
765               g_assert (weight >= 0.f);
766               graph->add_edge (id1, id2, weight, weight);
767             }
768         }
769     }
770 
771   /* vertical */
772 
773   for (x = 0; x < width; x++)
774     {
775       for (y = 0; y < height - 1; y++)
776         {
777           node_id id1 = nodes[x + y * width];
778           node_id id2 = nodes[x + (y+1) * width];
779 
780           if (id1 != -1 && id2 != -1)
781             {
782               gint costs_off = x + y * width;
783               gfloat weight  = 60.f * mean_costs / (v_costs[costs_off] + EPSILON);
784               g_assert (weight >= 0.f);
785               graph->add_edge (id1, id2, weight, weight);
786             }
787         }
788     }
789 }
790 
791 static gfloat *
792 paint_select_graph_get_segmentation (GraphType  *graph,
793                                      gint       *nodes,
794                                      guint8     *seeds,
795                                      gint        width,
796                                      gint        height)
797 {
798   gint     size = width * height;
799   gfloat  *segmentation;
800   gint     i;
801 
802   segmentation = g_new (gfloat, size);
803 
804   for (i = 0; i < size; i++)
805     {
806       node_id id = nodes[i];
807 
808       if (id != -1)
809         {
810           if (graph->what_segment(id) == GraphType::SOURCE)
811             segmentation[i] = 1.f;
812           else
813             segmentation[i] = 0.f;
814         }
815       else if (seeds[i] == SEED_SOURCE)
816         {
817           segmentation[i] = 1.f;
818         }
819       else if (seeds[i] == SEED_SINK)
820         {
821           segmentation[i] = 0.f;
822         }
823     }
824 
825   return segmentation;
826 }
827 
828 static gfloat *
829 paint_select_graphcut (gfloat              *pixels,
830                        guint8              *seeds,
831                        gint                 width,
832                        gint                 height,
833                        PaintSelectContext  *context)
834 {
835   GraphType  *graph;
836   gint       *nodes;
837   gfloat     *result;
838   gfloat     *h_costs;
839   gfloat     *v_costs;
840   gfloat      mean_costs;
841   gint        n_nodes;
842   gint        n_edges;
843   gint        n_h_costs;
844   gint        n_v_costs;
845 
846   n_h_costs = (width - 1) * height;
847   n_v_costs = (height - 1) * width;
848 
849   n_nodes = width * height;
850   n_edges = n_h_costs + n_v_costs;
851 
852   h_costs = g_new (gfloat, n_h_costs);
853   v_costs = g_new (gfloat, n_v_costs);
854 
855   paint_select_compute_adjacent_costs (pixels, width, height,
856                                        h_costs, v_costs, &mean_costs);
857 
858   graph = new GraphType (n_nodes, n_edges);
859   nodes = g_new (gint, n_nodes);
860 
861   paint_select_graph_init_nodes_and_tlinks (graph, pixels, seeds, nodes,
862                                             width, height, context);
863 
864   paint_select_graph_init_nlinks (graph, nodes, h_costs, v_costs, mean_costs,
865                                   width, height);
866 
867   graph->maxflow();
868 
869   g_free (h_costs);
870   g_free (v_costs);
871 
872   result = paint_select_graph_get_segmentation (graph, nodes, seeds, width, height);
873 
874   g_free (nodes);
875   delete graph;
876 
877   return result;
878 }
879 
880 /* high level functions */
881 
882 static PaintSelectContext *
883 paint_select_context_new (GeglProperties      *o,
884                           gfloat              *pixels,
885                           gfloat              *mask,
886                           gfloat              *scribbles,
887                           GeglRectangle       *roi,
888                           GeglRectangle       *extent,
889                           GeglRectangle       *region)
890 {
891   PaintSelectPrivate *priv = (PaintSelectPrivate *) o->user_data;
892   PaintSelectContext *context = g_slice_new (PaintSelectContext);
893 
894   if (! priv)
895     {
896       priv = g_slice_new (PaintSelectPrivate);
897       priv->fg_colors = NULL;
898       priv->bg_colors = NULL;
899       o->user_data    = priv;
900     }
901 
902   if (o->use_local_region)
903     {
904       context->boundary_top    = (roi->y > 0);
905       context->boundary_left   = (roi->x > 0);
906       context->boundary_right  = (roi->x + roi->width != extent->width);
907       context->boundary_bottom = (roi->y + roi->height != extent->height);
908     }
909 
910   if (o->mode == GEGL_PAINT_SELECT_MODE_ADD)
911     {
912       if (! priv->bg_colors)
913         {
914           priv->bg_colors = colors_model_new_global (pixels, mask,
915                                                      roi->width,
916                                                      roi->height,
917                                                      BG_MASK);
918         }
919       else
920         {
921           colors_model_update_global (priv->bg_colors, pixels, mask,
922                                       roi->width, roi->height, BG_MASK);
923         }
924 
925       context->local_colors = colors_model_new_local (pixels, mask, scribbles,
926                                                       roi->width, roi->height, region,
927                                                       FG_MASK, FG_SCRIBBLE);
928       context->mask_value_seed = FG_MASK;
929       context->mask_seed_type  = SEED_SOURCE;
930       context->source_model    = priv->bg_colors;
931       context->sink_model      = context->local_colors;
932       context->boundary_seed_type = SEED_SINK;
933     }
934   else
935     {
936       if (! priv->fg_colors)
937         {
938           priv->fg_colors = colors_model_new_global (pixels, mask,
939                                                      roi->width, roi->height, FG_MASK);
940         }
941       else
942         {
943           colors_model_update_global (priv->fg_colors, pixels, mask,
944                                       roi->width, roi->height, FG_MASK);
945         }
946 
947       context->local_colors = colors_model_new_local (pixels, mask, scribbles,
948                                                       roi->width, roi->height, region,
949                                                       BG_MASK, BG_SCRIBBLE);
950       context->mask_value_seed = BG_MASK;
951       context->mask_seed_type  = SEED_SINK;
952       context->source_model    = context->local_colors;
953       context->sink_model      = priv->fg_colors;
954       context->boundary_seed_type = SEED_SOURCE;
955     }
956 
957   return context;
958 }
959 
960 static void
961 paint_select_context_free (PaintSelectContext  *context)
962 {
963   colors_model_free (context->local_colors);
964   g_slice_free (PaintSelectContext, context);
965 }
966 
967 static void
968 paint_select_init_buffers (PaintSelect  *ps,
969                            GeglBuffer   *mask,
970                            GeglBuffer   *colors,
971                            GeglBuffer   *scribbles,
972                            GeglProperties  *o)
973 {
974   gint n_pixels;
975 
976   ps->extent = *gegl_buffer_get_extent (mask);
977 
978   if (o->use_local_region)
979     {
980       ps->roi.x      = o->region_x;
981       ps->roi.y      = o->region_y;
982       ps->roi.width  = o->region_width;
983       ps->roi.height = o->region_height;
984     }
985   else
986     {
987       gegl_rectangle_copy (&ps->roi, &ps->extent);
988     }
989 
990   ps->mode   = o->mode;
991   n_pixels   = ps->roi.width * ps->roi.height;
992 
993   ps->mask      = (gfloat *)  gegl_malloc (sizeof (gfloat) * n_pixels);
994   ps->colors    = (gfloat *)  gegl_malloc (sizeof (gfloat) * n_pixels * 3);
995   ps->scribbles = (gfloat *)  gegl_malloc (sizeof (gfloat) * n_pixels);
996 
997   gegl_buffer_get (mask, &ps->roi, 1.0, babl_format (SELECTION_FORMAT),
998                    ps->mask, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
999 
1000   gegl_buffer_get (colors, &ps->roi, 1.0, babl_format (COLORS_FORMAT),
1001                    ps->colors, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
1002 
1003   gegl_buffer_get (scribbles, &ps->roi, 1.0, babl_format (SCRIBBLES_FORMAT),
1004                    ps->scribbles, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE);
1005 }
1006 
1007 static void
1008 paint_select_free_buffers (PaintSelect  *ps)
1009 {
1010   gegl_free (ps->mask);
1011   gegl_free (ps->scribbles);
1012   gegl_free (ps->colors);
1013 }
1014 
1015 static gboolean
1016 paint_select_get_scribble_region (gfloat                  *mask,
1017                                   gfloat                  *scribbles,
1018                                   gint                     width,
1019                                   gint                     height,
1020                                   GeglRectangle           *region,
1021                                   gint                    *pix_x,
1022                                   gint                    *pix_y,
1023                                   GeglPaintSelectModeType  mode)
1024 {
1025   GeglRectangle extent = {0, 0, width, height};
1026   gfloat scribble_val;
1027   gfloat mask_val;
1028   gint minx, miny;
1029   gint maxx, maxy;
1030   gint x, y;
1031 
1032   minx = width;
1033   miny = height;
1034   maxx = maxy = 0;
1035 
1036   if (mode == GEGL_PAINT_SELECT_MODE_ADD)
1037     {
1038       scribble_val = FG_SCRIBBLE;
1039       mask_val     = BG_MASK;
1040     }
1041   else
1042     {
1043       scribble_val = BG_SCRIBBLE;
1044       mask_val     = FG_MASK;
1045     }
1046 
1047   for (y = 0; y < height; y++)
1048     {
1049       for (x = 0; x < width; x++)
1050         {
1051           gint off = x + y * width;
1052           if (scribbles[off] == scribble_val &&
1053               mask[off] == mask_val)
1054             {
1055               /* keep track of one pixel position located in
1056                * local region ; it will be used later as a seed
1057                * point for fluctuations removal.
1058                */
1059 
1060               *pix_x = x;
1061               *pix_y = y;
1062 
1063               if (x < minx)
1064                 minx = x;
1065               else if (x > maxx)
1066                 maxx = x;
1067 
1068               if (y < miny)
1069                 miny = y;
1070               else if (y > maxy)
1071                 maxy = y;
1072             }
1073         }
1074     }
1075 
1076   region->x = minx;
1077   region->y = miny;
1078   region->width = maxx - minx + 1;
1079   region->height = maxy - miny + 1;
1080 
1081   if (gegl_rectangle_is_empty (region) ||
1082       gegl_rectangle_is_infinite_plane (region))
1083     return FALSE;
1084 
1085   region->x -= LOCAL_REGION_DILATE;
1086   region->y -= LOCAL_REGION_DILATE;
1087   region->width += 2 * LOCAL_REGION_DILATE;
1088   region->height += 2 * LOCAL_REGION_DILATE;
1089 
1090   gegl_rectangle_intersect (region, region, &extent);
1091   return TRUE;
1092 }
1093 
1094 static void
1095 paint_select_compute_diff_mask (gfloat  *mask,
1096                                 gfloat  *result,
1097                                 gint     width,
1098                                 gint     height)
1099 {
1100   gint  i;
1101 
1102   for (i = 0; i < width * height; i++)
1103     {
1104       result[i] = result[i] != mask[i] ? 1.f : 0.f;
1105     }
1106 }
1107 
1108 
1109 /* GEGL operation */
1110 
1111 static void
1112 prepare (GeglOperation *operation)
1113 {
1114   const Babl *space = gegl_operation_get_source_space (operation, "aux");
1115   const Babl *selection  = babl_format (SELECTION_FORMAT);
1116   const Babl *scribbles  = babl_format (SCRIBBLES_FORMAT);
1117   const Babl *colors     = babl_format_with_space (COLORS_FORMAT, space);
1118 
1119   gegl_operation_set_format (operation, "input",  selection);
1120   gegl_operation_set_format (operation, "aux",    colors);
1121   gegl_operation_set_format (operation, "aux2",   scribbles);
1122   gegl_operation_set_format (operation, "output", selection);
1123 }
1124 
1125 static GeglRectangle
1126 get_bounding_box (GeglOperation  *operation)
1127 {
1128   GeglProperties  *o = GEGL_PROPERTIES (operation);
1129   GeglRectangle    result = {0, 0, 0, 0};
1130   GeglRectangle   *in_rect;
1131 
1132   in_rect = gegl_operation_source_get_bounding_box (operation, "input");
1133 
1134   if (o->use_local_region)
1135     {
1136       result.x = o->region_x;
1137       result.y = o->region_y;
1138       result.width  = o->region_width;
1139       result.height = o->region_height;
1140     }
1141   else if (in_rect)
1142     {
1143       result = *in_rect;
1144     }
1145 
1146   return result;
1147 }
1148 
1149 static GeglRectangle
1150 get_cached_region (GeglOperation        *operation,
1151                    const GeglRectangle  *roi)
1152 {
1153   GeglRectangle  empty_r  = {0, };
1154   GeglRectangle *in_r = gegl_operation_source_get_bounding_box (operation,
1155                                                                 "input");
1156   if (in_r)
1157     return *in_r;
1158 
1159   return empty_r;
1160 }
1161 
1162 static gboolean
1163 process (GeglOperation       *operation,
1164          GeglBuffer          *input,
1165          GeglBuffer          *aux,
1166          GeglBuffer          *aux2,
1167          GeglBuffer          *output,
1168          const GeglRectangle *result,
1169          gint                 level)
1170 {
1171   GeglProperties  *o = GEGL_PROPERTIES (operation);
1172   PaintSelect   ps;
1173   GeglRectangle region;
1174   gint          x = 0;
1175   gint          y = 0;
1176 
1177   if (! aux || ! aux2)
1178     {
1179       gegl_buffer_copy (input, NULL, GEGL_ABYSS_NONE, output, NULL);
1180       return TRUE;
1181     }
1182 
1183   /* memory allocations, pixels fetch */
1184 
1185   paint_select_init_buffers (&ps, input, aux, aux2, o);
1186 
1187   /* find the overlap region where scribble value doesn't match mask value */
1188 
1189   if (paint_select_get_scribble_region (ps.mask, ps.scribbles,
1190                                         ps.roi.width, ps.roi.height,
1191                                         &region, &x, &y, ps.mode))
1192     {
1193       PaintSelectContext *context;
1194       guint8             *seeds;
1195       gfloat             *result;
1196 
1197       g_printerr ("scribble region: (%d,%d) %d x %d\n",
1198                   region.x, region.y, region.width, region.height);
1199 
1200       context = paint_select_context_new (o,
1201                                           ps.colors,
1202                                           ps.mask,
1203                                           ps.scribbles,
1204                                           &ps.roi,
1205                                           &ps.extent,
1206                                           &region);
1207 
1208       seeds = paint_select_compute_seeds_map (ps.mask,
1209                                               ps.scribbles,
1210                                               ps.roi.width,
1211                                               ps.roi.height,
1212                                               context);
1213 
1214       result = paint_select_graphcut (ps.colors, seeds,
1215                                       ps.roi.width, ps.roi.height,
1216                                       context);
1217 
1218       /* compute difference between original mask and graphcut result.
1219        * then remove fluctuations */
1220 
1221       paint_select_compute_diff_mask (ps.mask, result, ps.roi.width, ps.roi.height);
1222       paint_select_remove_fluctuations (ps.mask, result, ps.roi.width, ps.roi.height, x, y);
1223 
1224       if (o->use_local_region)
1225         {
1226           gegl_buffer_set (output, &ps.roi, 0, babl_format (SELECTION_FORMAT),
1227                        ps.mask, GEGL_AUTO_ROWSTRIDE);
1228         }
1229       else
1230         {
1231           gegl_buffer_set (output, NULL, 0, babl_format (SELECTION_FORMAT),
1232                        ps.mask, GEGL_AUTO_ROWSTRIDE);
1233         }
1234 
1235       g_free (seeds);
1236       g_free (result);
1237       paint_select_context_free (context);
1238     }
1239 
1240   paint_select_free_buffers (&ps);
1241 
1242   return TRUE;
1243 }
1244 
1245 static void
1246 finalize (GObject *object)
1247 {
1248   GeglProperties *o = GEGL_PROPERTIES (object);
1249 
1250   if (o->user_data)
1251     {
1252       PaintSelectPrivate *priv = (PaintSelectPrivate *) o->user_data;
1253 
1254       if (priv->fg_colors)
1255         colors_model_free (priv->fg_colors);
1256 
1257       if (priv->bg_colors)
1258         colors_model_free (priv->bg_colors);
1259 
1260       g_slice_free (PaintSelectPrivate, o->user_data);
1261       o->user_data = NULL;
1262     }
1263 
1264   G_OBJECT_CLASS (gegl_op_parent_class)->finalize (object);
1265 }
1266 
1267 static void
1268 gegl_op_class_init (GeglOpClass *klass)
1269 {
1270   GObjectClass                *object_class;
1271   GeglOperationClass          *operation_class;
1272   GeglOperationComposer3Class *composer_class;
1273 
1274   object_class    = G_OBJECT_CLASS (klass);
1275   operation_class = GEGL_OPERATION_CLASS (klass);
1276   composer_class  = GEGL_OPERATION_COMPOSER3_CLASS (klass);
1277 
1278   object_class->finalize                     = finalize;
1279   operation_class->get_bounding_box          = get_bounding_box;
1280   operation_class->get_cached_region         = get_cached_region;
1281   operation_class->prepare                   = prepare;
1282   operation_class->threaded                  = FALSE;
1283   composer_class->process                    = process;
1284 
1285   gegl_operation_class_set_keys (operation_class,
1286   "name",         "gegl:paint-select",
1287   "title",        _("Paint Select"),
1288   NULL);
1289 }
1290 
1291 #endif
1292