• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*====================================================================*
2  -  Copyright (C) 2001 Leptonica.  All rights reserved.
3  -  This software is distributed in the hope that it will be
4  -  useful, but with NO WARRANTY OF ANY KIND.
5  -  No author or distributor accepts responsibility to anyone for the
6  -  consequences of using this software, or for whether it serves any
7  -  particular purpose or works at all, unless he or she says so in
8  -  writing.  Everyone is granted permission to copy, modify and
9  -  redistribute this source code, for commercial or non-commercial
10  -  purposes, with the following restrictions: (1) the origin of this
11  -  source code must not be misrepresented; (2) modified versions must
12  -  be plainly marked as such; and (3) this notice may not be removed
13  -  or altered from any source or modified source distribution.
14  *====================================================================*/
15 
16 /*
17  *  watershed.c
18  *
19  *      Top-level
20  *            L_WSHED         *wshedCreate()
21  *            void             wshedDestroy()
22  *            l_int32          wshedApply()
23  *
24  *      Helpers
25  *            static l_int32   identifyWatershedBasin()
26  *            static l_int32   mergeLookup()
27  *            static l_int32   wshedGetHeight()
28  *            static void      pushNewPixel()
29  *            static void      popNewPixel()
30  *            static void      pushWSPixel()
31  *            static void      popWSPixel()
32  *            static void      debugPrintLUT()
33  *            static void      debugWshedMerge()
34  *
35  *      Output
36  *            l_int32          wshedBasins()
37  *            PIX             *wshedRenderFill()
38  *            PIX             *wshedRenderColors()
39  *
40  *  The watershed function identifies the "catch basins" of the input
41  *  8 bpp image, with respect to the specified seeds or "markers".
42  *  The use is in segmentation, but the selection of the markers is
43  *  critical to getting meaningful results.
44  *
45  *  How are the markers selected?  You can't simply use the local
46  *  minima, because a typical image has sufficient noise so that
47  *  a useful catch basin can easily have multiple local minima.  However
48  *  they are selected, the question for the watershed function is
49  *  how to handle local minima that are not markers.  The reason
50  *  this is important is because of the algorithm used to find the
51  *  watersheds, which is roughly like this:
52  *
53  *    (1) Identify the markers and the local minima, and enter them
54  *        into a priority queue based on the pixel value.  Each marker
55  *        is shrunk to a single pixel, if necessary, before the
56  *        operation starts.
57  *    (2) Feed the priority queue with neighbors of pixels that are
58  *        popped off the queue.  Each of these queue pixels is labelled
59  *        with the index value of its parent.
60  *    (3) Each pixel is also labelled, in a 32-bit image, with the marker
61  *        or local minimum index, from which it was originally derived.
62  *    (4) There are actually 3 classes of labels: seeds, minima, and
63  *        fillers.  The fillers are labels of regions that have already
64  *        been identified as watersheds and are continuing to fill, for
65  *        the purpose of finding higher watersheds.
66  *    (5) When a pixel is popped that has already been labelled in the
67  *        32-bit image and that label differs from the label of its
68  *        parent (stored in the queue pixel), a boundary has been crossed.
69  *        There are several cases:
70  *         (a) Both parents are derived from markers but at least one
71  *             is not deep enough to become a watershed.  Absorb the
72  *             shallower basin into the deeper one, fixing the LUT to
73  *             redirect the shallower index to the deeper one.
74  *         (b) Both parents are derived from markers and both are deep
75  *             enough.  Identify and save the watershed for each marker.
76  *         (c) One parent was derived from a marker and the other from
77  *             a minima: absorb the minima basin into the marker basin.
78  *         (d) One parent was derived from a marker and the other is
79  *             a filler: identify and save the watershed for the marker.
80  *         (e) Both parents are derived from minima: merge them.
81  *         (f) One parent is a filler and the other is derived from a
82  *             minima: merge the minima into the filler.
83  *    (6) The output of the watershed operation consists of:
84  *         - a pixa of the basins
85  *         - a pta of the markers
86  *         - a numa of the watershed levels
87  *
88  *  Typical usage:
89  *      L_WShed *wshed = wshedCreate(pixs, pixseed, mindepth, 0);
90  *      wshedApply(wshed);
91 *
92  *      wshedBasins(wshed, &pixa, &nalevels);
93  *        ... do something with pixa, nalevels ...
94  *      pixaDestroy(&pixa);
95  *      numaDestroy(&nalevels);
96  *
97  *      Pix *pixd = wshedRenderFill(wshed);
98  *
99  *      wshedDestroy(&wshed);
100  */
101 
102 #include <stdio.h>
103 #include <stdlib.h>
104 #include "allheaders.h"
105 
106 #ifndef  NO_CONSOLE_IO
107 #define   DEBUG_WATERSHED     0
108 #endif  /* ~NO_CONSOLE_IO */
109 
110 static const l_uint32  MAX_LABEL_VALUE = 0x7fffffff;  /* largest l_int32 */
111 
112 struct L_NewPixel
113 {
114     l_int32    x;
115     l_int32    y;
116 };
117 typedef struct L_NewPixel  L_NEWPIXEL;
118 
119 struct L_WSPixel
120 {
121     l_float32  val;    /* pixel value */
122     l_int32    x;
123     l_int32    y;
124     l_int32    index;  /* label for set to which pixel belongs */
125 };
126 typedef struct L_WSPixel  L_WSPIXEL;
127 
128 
129     /* Static functions for obtaining bitmap of watersheds  */
130 static void wshedSaveBasin(L_WSHED *wshed, l_int32 index, l_int32 level);
131 
132 static l_int32 identifyWatershedBasin(L_WSHED *wshed,
133                                       l_int32 index, l_int32 level,
134                                       BOX **pbox, PIX **ppixd);
135 
136     /* Static function for merging lut and backlink arrays */
137 static l_int32 mergeLookup(L_WSHED *wshed, l_int32 sindex, l_int32 dindex);
138 
139     /* Static function for finding the height of the current pixel
140        above its seed or minima in the watershed.  */
141 static l_int32 wshedGetHeight(L_WSHED *wshed, l_int32 val, l_int32 label,
142                               l_int32 *pheight);
143 
144     /* Static accessors for NewPixel on a queue */
145 static void pushNewPixel(L_QUEUE *lq, l_int32 x, l_int32 y,
146                          l_int32 *pminx, l_int32 *pmaxx,
147                          l_int32 *pminy, l_int32 *pmaxy);
148 static void popNewPixel(L_QUEUE *lq, l_int32 *px, l_int32 *py);
149 
150     /* Static accessors for WSPixel on a heap */
151 static void pushWSPixel(L_HEAP *lh, L_STACK *stack, l_int32 val,
152                         l_int32 x, l_int32 y, l_int32 index);
153 static void popWSPixel(L_HEAP *lh, L_STACK *stack, l_int32 *pval,
154                        l_int32 *px, l_int32 *py, l_int32 *pindex);
155 
156     /* Static debug print output */
157 static void debugPrintLUT(l_int32 *lut, l_int32 size, l_int32 debug);
158 
159 static void debugWshedMerge(L_WSHED *wshed, char *descr, l_int32 x,
160                             l_int32 y, l_int32 label, l_int32 index);
161 
162 
163 /*-----------------------------------------------------------------------*
164  *                        Top-level watershed                            *
165  *-----------------------------------------------------------------------*/
166 /*!
167  *  wshedCreate()
168  *
169  *      Input:  pixs  (8 bpp source)
170  *              pixm  (1 bpp 'marker' seed)
171  *              mindepth (minimum depth; anything less is not saved)
172  *              debugflag (1 for debug output)
173  *      Return: WShed, or null on error
174  *
175  *  Notes:
176  *      (1) It is not necessary for the fg pixels in the seed image
177  *          be at minima, or that they be isolated.  We extract a
178  *          single pixel from each connected component, and a seed
179  *          anywhere in a watershed will eventually label the watershed
180  *          when the filling level reaches it.
181  *      (2) Set mindepth to some value to ignore noise in pixs that
182  *          can create small local minima.  Any watershed shallower
183  *          than mindepth, even if it has a seed, will not be saved;
184  *          It will either be incorporated in another watershed or
185  *          eliminated.
186  */
187 L_WSHED *
wshedCreate(PIX * pixs,PIX * pixm,l_int32 mindepth,l_int32 debugflag)188 wshedCreate(PIX     *pixs,
189             PIX     *pixm,
190             l_int32  mindepth,
191             l_int32  debugflag)
192 {
193 l_int32   w, h;
194 L_WSHED  *wshed;
195 
196     PROCNAME("wshedCreate");
197 
198     if (!pixs)
199         return (L_WSHED *)ERROR_PTR("pixs is not defined", procName, NULL);
200     if (pixGetDepth(pixs) != 8)
201         return (L_WSHED *)ERROR_PTR("pixs is not 8 bpp", procName, NULL);
202     if (!pixm)
203         return (L_WSHED *)ERROR_PTR("pixm is not defined", procName, NULL);
204     if (pixGetDepth(pixm) != 1)
205         return (L_WSHED *)ERROR_PTR("pixm is not 1 bpp", procName, NULL);
206     pixGetDimensions(pixs, &w, &h, NULL);
207     if (pixGetWidth(pixm) != w || pixGetHeight(pixm) != h)
208         return (L_WSHED *)ERROR_PTR("pixs/m sizes are unequal", procName, NULL);
209 
210     if ((wshed = (L_WSHED *)CALLOC(1, sizeof(L_WSHED))) == NULL)
211         return (L_WSHED *)ERROR_PTR("wshed not made", procName, NULL);
212 
213     wshed->pixs = pixClone(pixs);
214     wshed->pixm = pixClone(pixm);
215     wshed->mindepth = L_MAX(1, mindepth);
216     wshed->pixlab = pixCreate(w, h, 32);
217     pixSetAllArbitrary(wshed->pixlab, MAX_LABEL_VALUE);
218     wshed->pixt = pixCreate(w, h, 1);
219     wshed->lines8 = pixGetLinePtrs(pixs, NULL);
220     wshed->linem1 = pixGetLinePtrs(pixm, NULL);
221     wshed->linelab32 = pixGetLinePtrs(wshed->pixlab, NULL);
222     wshed->linet1 = pixGetLinePtrs(wshed->pixt, NULL);
223     wshed->debug = debugflag;
224     return wshed;
225 }
226 
227 
228 /*!
229  *  wshedDestroy()
230  *
231  *      Input:  &wshed (<will be set to null before returning>)
232  *      Return: void
233  */
234 void
wshedDestroy(L_WSHED ** pwshed)235 wshedDestroy(L_WSHED  **pwshed)
236 {
237 l_int32   i;
238 L_WSHED  *wshed;
239 
240     PROCNAME("wshedDestroy");
241 
242     if (pwshed == NULL) {
243         L_WARNING("ptr address is null!", procName);
244         return;
245     }
246 
247     if ((wshed = *pwshed) == NULL)
248         return;
249 
250     pixDestroy(&wshed->pixs);
251     pixDestroy(&wshed->pixm);
252     pixDestroy(&wshed->pixlab);
253     pixDestroy(&wshed->pixt);
254     if (wshed->lines8) FREE(wshed->lines8);
255     if (wshed->linem1) FREE(wshed->linem1);
256     if (wshed->linelab32) FREE(wshed->linelab32);
257     if (wshed->linet1) FREE(wshed->linet1);
258     pixaDestroy(&wshed->pixad);
259     ptaDestroy(&wshed->ptas);
260     numaDestroy(&wshed->nash);
261     numaDestroy(&wshed->nasi);
262     numaDestroy(&wshed->namh);
263     numaDestroy(&wshed->nalevels);
264     if (wshed->lut)
265          FREE(wshed->lut);
266     if (wshed->links) {
267         for (i = 0; i < wshed->arraysize; i++)
268             numaDestroy(&wshed->links[i]);
269         FREE(wshed->links);
270     }
271     FREE(wshed);
272     *pwshed = NULL;
273     return;
274 }
275 
276 
277 /*!
278  *  wshedApply()
279  *
280  *      Input:  wshed (generated from wshedCreate())
281  *      Return: 0 if OK, 1 on error
282  *
283  *  Iportant note:
284  *      (1) This is buggy.  It seems to locate watersheds that are
285  *          duplicates.  The watershed extraction after complete fill
286  *          grabs some regions belonging to existing watersheds.
287  *          See prog/watershedtest.c for testing.
288  */
289 l_int32
wshedApply(L_WSHED * wshed)290 wshedApply(L_WSHED  *wshed)
291 {
292 char      two_new_watersheds[] = "Two new watersheds";
293 char      seed_absorbed_into_seeded_basin[] = "Seed absorbed into seeded basin";
294 char      one_new_watershed_label[] = "One new watershed (label)";
295 char      one_new_watershed_index[] = "One new watershed (index)";
296 char      minima_absorbed_into_seeded_basin[] =
297                  "Minima absorbed into seeded basin";
298 char      minima_absorbed_by_filler_or_another[] =
299                  "Minima absorbed by filler or another";
300 l_int32   nseeds, nother, nboth, arraysize;
301 l_int32   i, j, val, x, y, w, h, index, mindepth;
302 l_int32   imin, imax, jmin, jmax, cindex, clabel, nindex;
303 l_int32   hindex, hlabel, hmin, hmax, minhindex, maxhindex;
304 l_int32  *lut;
305 l_uint32  ulabel, uval;
306 void    **lines8, **linelab32;
307 NUMA     *nalut, *nalevels, *nash, *namh, *nasi;
308 NUMA    **links;
309 L_HEAP   *lh;
310 PIX      *pixmin, *pixsd;
311 PIXA     *pixad;
312 L_STACK  *rstack;
313 PTA      *ptas, *ptao;
314 
315     PROCNAME("wshedApply");
316 
317     if (!wshed)
318         return ERROR_INT("wshed not defined", procName, 1);
319 
320     /* ------------------------------------------------------------ *
321      *  Initialize priority queue and pixlab with seeds and minima  *
322      * ------------------------------------------------------------ */
323 
324     lh = lheapCreate(0, L_SORT_INCREASING);  /* remove lowest values first */
325     rstack = lstackCreate(0);  /* for reusing the WSPixels */
326     pixGetDimensions(wshed->pixs, &w, &h, NULL);
327     lines8 = wshed->lines8;  /* wshed owns this */
328     linelab32 = wshed->linelab32;  /* ditto */
329 
330         /* Identify seed (marker) pixels, 1 for each c.c. in pixm */
331     ptas = pixSelectMinInConnComp(wshed->pixs, wshed->pixm, &nash);
332     pixsd = pixGenerateFromPta(ptas, w, h);
333     nseeds = ptaGetCount(ptas);
334     for (i = 0; i < nseeds; i++) {
335         ptaGetIPt(ptas, i, &x, &y);
336         uval = GET_DATA_BYTE(lines8[y], x);
337         pushWSPixel(lh, rstack, (l_int32)uval, x, y, i);
338     }
339     wshed->ptas = ptas;
340     nasi = numaMakeConstant(1, nseeds);  /* indicator array */
341     wshed->nasi = nasi;
342     wshed->nash = nash;
343     wshed->nseeds = nseeds;
344 
345         /* Identify minima that are not seeds.  Use these 4 steps:
346          *  (1) Get the local minima, which can have components
347          *      of arbitrary size.  This will be a clipping mask.
348          *  (2) Get the image of the actual seeds (pixsd)
349          *  (3) Remove all elements of the clipping mask that have a seed.
350          *  (4) Shrink each of the remaining elements of the minima mask
351          *      to a single pixel.  */
352     pixLocalExtrema(wshed->pixs, 200, 0, &pixmin, NULL);
353     pixRemoveSeededComponents(pixmin, pixsd, pixmin, 8, 2);
354     ptao = pixSelectMinInConnComp(wshed->pixs, pixmin, &namh);
355     nother = ptaGetCount(ptao);
356     for (i = 0; i < nother; i++) {
357         ptaGetIPt(ptao, i, &x, &y);
358         uval = GET_DATA_BYTE(lines8[y], x);
359         pushWSPixel(lh, rstack, (l_int32)uval, x, y, nseeds + i);
360     }
361     wshed->namh = namh;
362 
363     /* ------------------------------------------------------------ *
364      *                Initialize merging lookup tables              *
365      * ------------------------------------------------------------ */
366 
367         /* nalut should always give the current after-merging index.
368          * links are effectively backpointers: they are numas associated with
369          * a dest index of all indices in nalut that point to that index. */
370     mindepth = wshed->mindepth;
371     nboth = nseeds + nother;
372     arraysize = 2 * nboth;
373     wshed->arraysize = arraysize;
374     nalut = numaMakeSequence(0, 1, arraysize);
375     lut = numaGetIArray(nalut);
376     wshed->lut = lut;  /* wshed owns this */
377     links = (NUMA **)CALLOC(arraysize, sizeof(NUMA *));
378     wshed->links = links;  /* wshed owns this */
379     nindex = nseeds + nother;  /* the next unused index value */
380 
381     /* ------------------------------------------------------------ *
382      *              Fill the basins, using the priority queue       *
383      * ------------------------------------------------------------ */
384 
385     pixad = pixaCreate(nseeds);
386     wshed->pixad = pixad;  /* wshed owns this */
387     nalevels = numaCreate(nseeds);
388     wshed->nalevels = nalevels;  /* wshed owns this */
389     L_INFO_INT2("nseeds = %d, nother = %d\n", procName, nseeds, nother);
390     while (lheapGetCount(lh) > 0) {
391         popWSPixel(lh, rstack, &val, &x, &y, &index);
392 /*        fprintf(stderr, "x = %d, y = %d, index = %d\n", x, y, index); */
393         ulabel = GET_DATA_FOUR_BYTES(linelab32[y], x);
394         if (ulabel == MAX_LABEL_VALUE)
395             clabel = ulabel;
396         else
397             clabel = lut[ulabel];
398         cindex = lut[index];
399         if (clabel == cindex) continue;  /* have already seen this one */
400         if (clabel == MAX_LABEL_VALUE) {  /* new one; assign index and try to
401                                            * propagate to all neighbors */
402             SET_DATA_FOUR_BYTES(linelab32[y], x, cindex);
403             imin = L_MAX(0, y - 1);
404             imax = L_MIN(h - 1, y + 1);
405             jmin = L_MAX(0, x - 1);
406             jmax = L_MIN(w - 1, x + 1);
407             for (i = imin; i <= imax; i++) {
408                 for (j = jmin; j <= jmax; j++) {
409                     if (i == y && j == x) continue;
410                     uval = GET_DATA_BYTE(lines8[i], j);
411                     pushWSPixel(lh, rstack, (l_int32)uval, j, i, cindex);
412                 }
413             }
414         }
415         else {  /* this pixel is already labeled (differently); must resolve */
416 
417                 /* If both indices are seeds, check if the min height is
418                  * greater than mindepth.  If so, we have two new watersheds;
419                  * locate them and assign to both regions a new index
420                  * for further waterfill.  If not, absorb the shallower
421                  * watershed into the deeper one and continue filling it. */
422             pixGetPixel(pixsd, x, y, &uval);
423             if (clabel < nseeds && cindex < nseeds) {
424                 wshedGetHeight(wshed, val, clabel, &hlabel);
425                 wshedGetHeight(wshed, val, cindex, &hindex);
426                 hmin = L_MIN(hlabel, hindex);
427                 hmax = L_MAX(hlabel, hindex);
428                 if (hmin == hmax) {
429                     hmin = hlabel;
430                     hmax = hindex;
431                 }
432                 if (wshed->debug) {
433                     fprintf(stderr, "clabel,hlabel = %d,%d\n", clabel, hlabel);
434                     fprintf(stderr, "hmin = %d, hmax = %d\n", hmin, hmax);
435                     fprintf(stderr, "cindex,hindex = %d,%d\n", cindex, hindex);
436                     if (hmin < mindepth)
437                         fprintf(stderr, "Too shallow!\n");
438                 }
439 
440                 if (hmin >= mindepth) {
441                     debugWshedMerge(wshed, two_new_watersheds,
442                                     x, y, clabel, cindex);
443                     wshedSaveBasin(wshed, cindex, val - 1);
444                     wshedSaveBasin(wshed, clabel, val - 1);
445                     numaSetValue(nasi, cindex, 0);
446                     numaSetValue(nasi, clabel, 0);
447 
448                     if (wshed->debug) fprintf(stderr, "nindex = %d\n", nindex);
449                     debugPrintLUT(lut, nindex, wshed->debug);
450                     mergeLookup(wshed, clabel, nindex);
451                     debugPrintLUT(lut, nindex, wshed->debug);
452                     mergeLookup(wshed, cindex, nindex);
453                     debugPrintLUT(lut, nindex, wshed->debug);
454                     nindex++;
455                 }
456                 else  /* extraneous seed within seeded basin; absorb */
457                     debugWshedMerge(wshed, seed_absorbed_into_seeded_basin,
458                                     x, y, clabel, cindex);
459                     maxhindex = clabel;
460                     minhindex = cindex;
461                     if (hindex > hlabel) {
462                         maxhindex = cindex;
463                         minhindex = clabel;
464                     }
465                     mergeLookup(wshed, minhindex, maxhindex);
466             }
467 
468                 /* If one index is a seed and the other is a merge of
469                  * 2 watersheds, generate a single watershed. */
470             else if (clabel < nseeds && cindex >= nboth) {
471                 debugWshedMerge(wshed, one_new_watershed_label,
472                                 x, y, clabel, cindex);
473                 wshedSaveBasin(wshed, clabel, val - 1);
474                 numaSetValue(nasi, clabel, 0);
475                 mergeLookup(wshed, clabel, cindex);
476             }
477             else if (cindex < nseeds && clabel >= nboth) {
478                 debugWshedMerge(wshed, one_new_watershed_index,
479                                 x, y, clabel, cindex);
480                 wshedSaveBasin(wshed, cindex, val - 1);
481                 numaSetValue(nasi, cindex, 0);
482                 mergeLookup(wshed, cindex, clabel);
483             }
484                 /* If one index is a seed and the other is from a minimum,
485                  * merge the minimum wshed into the seed wshed. */
486             else if (clabel < nseeds) {  /* cindex from minima; absorb */
487                 debugWshedMerge(wshed, minima_absorbed_into_seeded_basin,
488                                 x, y, clabel, cindex);
489                 mergeLookup(wshed, cindex, clabel);
490             }
491             else if (cindex < nseeds) {  /* clabel from minima; absorb */
492                 debugWshedMerge(wshed, minima_absorbed_into_seeded_basin,
493                                 x, y, clabel, cindex);
494                 mergeLookup(wshed, clabel, cindex);
495             }
496                 /* If neither index is a seed, just merge */
497             else {
498                 debugWshedMerge(wshed, minima_absorbed_by_filler_or_another,
499                                 x, y, clabel, cindex);
500                 mergeLookup(wshed, clabel, cindex);
501             }
502         }
503     }
504 
505 #if 0
506         /*  Use the indicator array to save any watersheds that fill
507          *  to the maximum value.  This seems to screw things up!  */
508     for (i = 0; i < nseeds; i++) {
509         numaGetIValue(nasi, i, &ival);
510         if (ival == 1) {
511             wshedSaveBasin(wshed, lut[i], val - 1);
512             numaSetValue(nasi, i, 0);
513         }
514     }
515 #endif
516 
517     numaDestroy(&nalut);
518     pixDestroy(&pixmin);
519     pixDestroy(&pixsd);
520     ptaDestroy(&ptao);
521     lheapDestroy(&lh, TRUE);
522     lstackDestroy(&rstack, TRUE);
523     return 0;
524 }
525 
526 
527 /*-----------------------------------------------------------------------*
528  *                               Helpers                                 *
529  *-----------------------------------------------------------------------*/
530 /*!
531  *  wshedSaveBasin()
532  *
533  *      Input:  wshed
534  *              index (index of basin to be located)
535  *              level (filling level reached at the time this function
536  *                     is called)
537  *      Return: 0 if OK, 1 on error
538  *
539  *  Notes:
540  *      (1) This identifies a single watershed.  It does not change
541  *          the LUT, which must be done subsequently.
542  *      (2) The fill level of a basin is taken to be @level - 1.
543  */
544 static void
wshedSaveBasin(L_WSHED * wshed,l_int32 index,l_int32 level)545 wshedSaveBasin(L_WSHED  *wshed,
546                l_int32   index,
547                l_int32   level)
548 {
549 BOX  *box;
550 PIX  *pix;
551 
552     PROCNAME("wshedSaveBasin");
553 
554     if (!wshed)
555         return ERROR_VOID("wshed not defined", procName);
556 
557     if (identifyWatershedBasin(wshed, index, level, &box, &pix) == 0) {
558         pixaAddPix(wshed->pixad, pix, L_INSERT);
559         pixaAddBox(wshed->pixad, box, L_INSERT);
560         numaAddNumber(wshed->nalevels, level - 1);
561     }
562     return;
563 }
564 
565 
566 /*!
567  *  identifyWatershedBasin()
568  *
569  *      Input:  wshed
570  *              index (index of basin to be located)
571  *              level (of basin at point at which the two basins met)
572  *              &box (<return> bounding box of basin)
573  *              &pixd (<return> pix of basin, cropped to its bounding box)
574  *      Return: 0 if OK, 1 on error
575  *
576  *  Notes:
577  *      (1) This is a static function, so we assume pixlab, pixs and pixt
578  *          exist and are the same size.
579  *      (2) It selects all pixels that have the label @index in pixlab
580  *          and that have a value in pixs that is less than @level.
581  *      (3) It is used whenever two seeded basins meet (typically at a saddle),
582  *          or when one seeded basin meets a 'filler'.  All identified
583  *          basins are saved as a watershed.
584  */
585 static l_int32
identifyWatershedBasin(L_WSHED * wshed,l_int32 index,l_int32 level,BOX ** pbox,PIX ** ppixd)586 identifyWatershedBasin(L_WSHED  *wshed,
587                        l_int32   index,
588                        l_int32   level,
589                        BOX     **pbox,
590                        PIX     **ppixd)
591 {
592 l_int32   imin, imax, jmin, jmax, minx, miny, maxx, maxy;
593 l_int32   bw, bh, i, j, w, h, x, y;
594 l_int32  *lut;
595 l_uint32  label, bval, lval;
596 void    **lines8, **linelab32, **linet1;
597 BOX      *box;
598 PIX      *pixs, *pixlab, *pixt, *pixd;
599 L_QUEUE  *lq;
600 
601     PROCNAME("identifyWatershedBasin");
602 
603     if (!pbox)
604         return ERROR_INT("&box not defined", procName, 1);
605     *pbox = NULL;
606     if (!ppixd)
607         return ERROR_INT("&pixd not defined", procName, 1);
608     *ppixd = NULL;
609     if (!wshed)
610         return ERROR_INT("wshed not defined", procName, 1);
611 
612         /* Make a queue and an auxiliary stack */
613     lq = lqueueCreate(0);
614     lq->stack = lstackCreate(0);
615 
616     pixs = wshed->pixs;
617     pixlab = wshed->pixlab;
618     pixt = wshed->pixt;
619     lines8 = wshed->lines8;
620     linelab32 = wshed->linelab32;
621     linet1 = wshed->linet1;
622     lut = wshed->lut;
623     pixGetDimensions(pixs, &w, &h, NULL);
624 
625         /* Prime the queue with the seed pixel for this watershed. */
626     minx = miny = 1000000;
627     maxx = maxy = 0;
628     ptaGetIPt(wshed->ptas, index, &x, &y);
629     pixSetPixel(pixt, x, y, 1);
630     pushNewPixel(lq, x, y, &minx, &maxx, &miny, &maxy);
631     if (wshed->debug) fprintf(stderr, "prime: (x,y) = (%d, %d)\n", x, y);
632 
633         /* Each pixel in a spreading breadth-first search is inspected.
634          * It is accepted as part of this watershed, and pushed on
635          * the search queue, if:
636          *     (1) It has a label value equal to @index
637          *     (2) The pixel value is less than @level, the overflow
638          *         height at which the two basins join.
639          *     (3) It has not yet been seen in this search.  */
640     while (lqueueGetCount(lq) > 0) {
641         popNewPixel(lq, &x, &y);
642         imin = L_MAX(0, y - 1);
643         imax = L_MIN(h - 1, y + 1);
644         jmin = L_MAX(0, x - 1);
645         jmax = L_MIN(w - 1, x + 1);
646         for (i = imin; i <= imax; i++) {
647             for (j = jmin; j <= jmax; j++) {
648                 if (j == x && i == y) continue;  /* parent */
649                 label = GET_DATA_FOUR_BYTES(linelab32[i], j);
650                 if (label == MAX_LABEL_VALUE || lut[label] != index) continue;
651                 bval = GET_DATA_BIT(linet1[i], j);
652                 if (bval == 1) continue;  /* already seen */
653                 lval = GET_DATA_BYTE(lines8[i], j);
654                 if (lval >= level) continue;  /* too high */
655                 SET_DATA_BIT(linet1[i], j);
656                 pushNewPixel(lq, j, i, &minx, &maxx, &miny, &maxy);
657             }
658         }
659     }
660 
661         /* Extract the box and pix, and clear pixt */
662     bw = maxx - minx + 1;
663     bh = maxy - miny + 1;
664     box = boxCreate(minx, miny, bw, bh);
665     pixd = pixClipRectangle(pixt, box, NULL);
666     pixRasterop(pixt, minx, miny, bw, bh, PIX_SRC ^ PIX_DST, pixd, 0, 0);
667     *pbox = box;
668     *ppixd = pixd;
669 
670     lqueueDestroy(&lq, 1);
671     return 0;
672 }
673 
674 
675 /*!
676  *  mergeLookup()
677  *
678  *      Input:  wshed
679  *              sindex (primary index being changed in the merge)
680  *              dindex (index that @sindex will point to after the merge)
681  *      Return: 0 if OK, 1 on error
682  *
683  *  Notes:
684  *      (1) The links are a sparse array of Numas showing current back-links.
685  *          The lut gives the current index (of the seed or the minima
686  *          for the wshed  in which it is located.
687  *      (2) Think of each entry in the lut.  There are two types:
688  *             owner:     lut[index] = index
689  *             redirect:  lut[index] != index
690  *      (3) This is called each time a merge occurs.  It puts the lut
691  *          and backlinks in a canonical form after the merge, where
692  *          all entries in the lut point to the current "owner", which
693  *          has all backlinks.  That is, every "redirect" in the lut
694  *          points to an "owner".  The lut always gives the index of
695  *          the current owner.
696  */
697 static l_int32
mergeLookup(L_WSHED * wshed,l_int32 sindex,l_int32 dindex)698 mergeLookup(L_WSHED  *wshed,
699             l_int32   sindex,
700             l_int32   dindex)
701 {
702 l_int32   i, n, size, index;
703 l_int32  *lut;
704 NUMA     *na;
705 NUMA    **links;
706 
707     PROCNAME("mergeLookup");
708 
709     if (!wshed)
710         return ERROR_INT("wshed not defined", procName, 1);
711     size = wshed->arraysize;
712     if (sindex < 0 || sindex >= size)
713         return ERROR_INT("invalid sindex", procName, 1);
714     if (dindex < 0 || dindex >= size)
715         return ERROR_INT("invalid dindex", procName, 1);
716 
717         /* Redirect links in the lut */
718     n = 0;
719     links = wshed->links;
720     lut = wshed->lut;
721     if ((na = links[sindex]) != NULL) {
722         n = numaGetCount(na);
723         for (i = 0; i < n; i++) {
724             numaGetIValue(na, i, &index);
725             lut[index] = dindex;
726         }
727     }
728     lut[sindex] = dindex;
729 
730         /* Shift the backlink arrays from sindex to dindex.
731          * sindex should have no backlinks because all entries in the
732          * lut that were previously pointing to it have been redirected
733          * to dindex. */
734     if (!links[dindex])
735         links[dindex] = numaCreate(n);
736     numaJoin(links[dindex], links[sindex], 0, 0);
737     numaAddNumber(links[dindex], sindex);
738     numaDestroy(&links[sindex]);
739 
740     return 0;
741 }
742 
743 
744 /*!
745  *  wshedGetHeight()
746  *
747  *      Input:  wshed (array of current indices)
748  *              val (value of current pixel popped off queue)
749  *              label (of pixel or 32 bpp label image)
750  *              &height (<return> height of current value from seed
751  *                       or minimum of watershed)
752  *      Return: 0 if OK, 1 on error
753  *
754  *  Notes:
755  *      (1) It is only necessary to find the height for a watershed
756  *          that is indexed by a seed or a minima.  This function should
757  *          not be called on a finished watershed (that continues to fill).
758  */
759 static l_int32
wshedGetHeight(L_WSHED * wshed,l_int32 val,l_int32 label,l_int32 * pheight)760 wshedGetHeight(L_WSHED  *wshed,
761                l_int32   val,
762                l_int32   label,
763                l_int32  *pheight)
764 {
765 l_int32  minval;
766 
767     PROCNAME("wshedGetHeight");
768 
769     if (!pheight)
770         return ERROR_INT("&height not defined", procName, 1);
771     *pheight = 0;
772     if (!wshed)
773         return ERROR_INT("wshed not defined", procName, 1);
774 
775     if (label < wshed->nseeds)
776         numaGetIValue(wshed->nash, label, &minval);
777     else if (label < wshed->nseeds + wshed->nother)
778         numaGetIValue(wshed->namh, label, &minval);
779     else
780         return ERROR_INT("finished watershed; should not call", procName, 1);
781 
782     *pheight = val - minval;
783     return 0;
784 }
785 
786 
787 /*
788  *  pushNewPixel()
789  *
790  *      Input:  lqueue
791  *              x, y   (pixel coordinates)
792  *              &minx, &maxx, &miny, &maxy  (<return> bounding box update)
793  *      Return: void
794  *
795  *  Notes:
796  *      (1) This is a wrapper for adding a NewPixel to a queue, which
797  *          updates the bounding box for all pixels on that queue and
798  *          uses the storage stack to retrieve a NewPixel.
799  */
800 static void
pushNewPixel(L_QUEUE * lq,l_int32 x,l_int32 y,l_int32 * pminx,l_int32 * pmaxx,l_int32 * pminy,l_int32 * pmaxy)801 pushNewPixel(L_QUEUE  *lq,
802              l_int32   x,
803              l_int32   y,
804              l_int32  *pminx,
805              l_int32  *pmaxx,
806              l_int32  *pminy,
807              l_int32  *pmaxy)
808 {
809 L_NEWPIXEL  *np;
810 
811     PROCNAME("pushNewPixel");
812 
813     if (!lq)
814         return ERROR_VOID(procName, "queue not defined");
815 
816         /* Adjust bounding box */
817     *pminx = L_MIN(*pminx, x);
818     *pmaxx = L_MAX(*pmaxx, x);
819     *pminy = L_MIN(*pminy, y);
820     *pmaxy = L_MAX(*pmaxy, y);
821 
822         /* Get a newpixel to use */
823     if (lstackGetCount(lq->stack) > 0)
824         np = (L_NEWPIXEL *)lstackRemove(lq->stack);
825     else
826         np = (L_NEWPIXEL *)CALLOC(1, sizeof(L_NEWPIXEL));
827 
828     np->x = x;
829     np->y = y;
830     lqueueAdd(lq, np);
831     return;
832 }
833 
834 
835 /*
836  *  popNewPixel()
837  *
838  *      Input:  lqueue
839  *              &x, &y   (<return> pixel coordinates)
840  *      Return: void
841  *
842  *   Notes:
843  *       (1) This is a wrapper for removing a NewPixel from a queue,
844  *           which returns the pixel coordinates and saves the NewPixel
845  *           on the storage stack.
846  */
847 static void
popNewPixel(L_QUEUE * lq,l_int32 * px,l_int32 * py)848 popNewPixel(L_QUEUE  *lq,
849             l_int32  *px,
850             l_int32  *py)
851 {
852 L_NEWPIXEL  *np;
853 
854     PROCNAME("popNewPixel");
855 
856     if (!lq)
857         return ERROR_VOID(procName, "lqueue not defined");
858 
859     if ((np = (L_NEWPIXEL *)lqueueRemove(lq)) == NULL)
860         return;
861     *px = np->x;
862     *py = np->y;
863     lstackAdd(lq->stack, np);  /* save for re-use */
864     return;
865 }
866 
867 
868 /*
869  *  pushWSPixel()
870  *
871  *      Input:  lh  (priority queue)
872  *              stack  (of reusable WSPixels)
873  *              val  (pixel value: used for ordering the heap)
874  *              x, y  (pixel coordinates)
875  *              index  (label for set to which pixel belongs)
876  *      Return: void
877  *
878  *  Notes:
879  *      (1) This is a wrapper for adding a WSPixel to a heap.  It
880  *          uses the storage stack to retrieve a WSPixel.
881  */
882 static void
pushWSPixel(L_HEAP * lh,L_STACK * stack,l_int32 val,l_int32 x,l_int32 y,l_int32 index)883 pushWSPixel(L_HEAP   *lh,
884             L_STACK  *stack,
885             l_int32   val,
886             l_int32   x,
887             l_int32   y,
888             l_int32   index)
889 {
890 L_WSPIXEL  *wsp;
891 
892     PROCNAME("pushWSPixel");
893 
894     if (!lh)
895         return ERROR_VOID(procName, "heap not defined");
896     if (!stack)
897         return ERROR_VOID(procName, "stack not defined");
898 
899         /* Get a wspixel to use */
900     if (lstackGetCount(stack) > 0)
901         wsp = (L_WSPIXEL *)lstackRemove(stack);
902     else
903         wsp = (L_WSPIXEL *)CALLOC(1, sizeof(L_WSPIXEL));
904 
905     wsp->val = (l_float32)val;
906     wsp->x = x;
907     wsp->y = y;
908     wsp->index = index;
909     lheapAdd(lh, wsp);
910     return;
911 }
912 
913 
914 /*
915  *  popWSPixel()
916  *
917  *      Input:  lh  (priority queue)
918  *              stack  (of reusable WSPixels)
919  *              &val  (<return> pixel value)
920  *              &x, &y  (<return> pixel coordinates)
921  *              &index  (<return> label for set to which pixel belongs)
922  *      Return: void
923  *
924  *   Notes:
925  *       (1) This is a wrapper for removing a WSPixel from a heap,
926  *           which returns the WSPixel data and saves the WSPixel
927  *           on the storage stack.
928  */
929 static void
popWSPixel(L_HEAP * lh,L_STACK * stack,l_int32 * pval,l_int32 * px,l_int32 * py,l_int32 * pindex)930 popWSPixel(L_HEAP   *lh,
931            L_STACK  *stack,
932            l_int32  *pval,
933            l_int32  *px,
934            l_int32  *py,
935            l_int32  *pindex)
936 {
937 L_WSPIXEL  *wsp;
938 
939     PROCNAME("popWSPixel");
940 
941     if (!lh)
942         return ERROR_VOID(procName, "lheap not defined");
943     if (!stack)
944         return ERROR_VOID(procName, "stack not defined");
945     if (!pval || !px || !py || !pindex)
946         return ERROR_VOID(procName, "data can't be returned");
947 
948     if ((wsp = (L_WSPIXEL *)lheapRemove(lh)) == NULL)
949         return;
950     *pval = (l_int32)wsp->val;
951     *px = wsp->x;
952     *py = wsp->y;
953     *pindex = wsp->index;
954     lstackAdd(stack, wsp);  /* save for re-use */
955     return;
956 }
957 
958 
959 static void
debugPrintLUT(l_int32 * lut,l_int32 size,l_int32 debug)960 debugPrintLUT(l_int32  *lut,
961               l_int32   size,
962               l_int32   debug)
963 {
964 l_int32  i;
965 
966     if (!debug) return;
967     fprintf(stderr, "lut: ");
968     for (i = 0; i < size; i++)
969         fprintf(stderr, "%d ", lut[i]);
970     fprintf(stderr, "\n");
971     return;
972 }
973 
974 
975 static void
debugWshedMerge(L_WSHED * wshed,char * descr,l_int32 x,l_int32 y,l_int32 label,l_int32 index)976 debugWshedMerge(L_WSHED *wshed,
977                 char    *descr,
978                 l_int32  x,
979                 l_int32  y,
980                 l_int32  label,
981                 l_int32  index)
982 {
983     if (!wshed || (wshed->debug == 0))
984          return;
985     fprintf(stderr, "%s:\n", descr);
986     fprintf(stderr, "   (x, y) = (%d, %d)\n", x, y);
987     fprintf(stderr, "   clabel = %d, cindex = %d\n", label, index);
988     return;
989 }
990 
991 
992 /*-----------------------------------------------------------------------*
993  *                                 Output                                *
994  *-----------------------------------------------------------------------*/
995 /*!
996  *  wshedBasins()
997  *
998  *      Input:  wshed
999  *              &pixa  (<optional return> mask of watershed basins)
1000  *              &nalevels   (<optional return> watershed levels)
1001  *      Return: 0 if OK, 1 on error
1002  */
1003 l_int32
wshedBasins(L_WSHED * wshed,PIXA ** ppixa,NUMA ** pnalevels)1004 wshedBasins(L_WSHED  *wshed,
1005             PIXA    **ppixa,
1006             NUMA    **pnalevels)
1007 {
1008     PROCNAME("wshedBasins");
1009 
1010     if (!wshed)
1011         return ERROR_INT("wshed not defined", procName, 1);
1012 
1013     if (ppixa)
1014         *ppixa = pixaCopy(wshed->pixad, L_CLONE);
1015     if (pnalevels)
1016         *pnalevels = numaClone(wshed->nalevels);
1017     return 0;
1018 }
1019 
1020 
1021 /*!
1022  *  wshedRenderFill()
1023  *
1024  *      Input:  wshed
1025  *      Return: pixd (initial image with all basins filled), or null on error
1026  */
1027 PIX *
wshedRenderFill(L_WSHED * wshed)1028 wshedRenderFill(L_WSHED  *wshed)
1029 {
1030 l_int32  i, n, level, bx, by;
1031 NUMA    *na;
1032 PIX     *pix, *pixd;
1033 PIXA    *pixa;
1034 
1035     PROCNAME("wshedRenderFill");
1036 
1037     if (!wshed)
1038         return (PIX *)ERROR_PTR("wshed not defined", procName, NULL);
1039 
1040     wshedBasins(wshed, &pixa, &na);
1041     pixd = pixCopy(NULL, wshed->pixs);
1042     n = pixaGetCount(pixa);
1043     for (i = 0; i < n; i++) {
1044         pix = pixaGetPix(pixa, i, L_CLONE);
1045         pixaGetBoxGeometry(pixa, i, &bx, &by, NULL, NULL);
1046         numaGetIValue(na, i, &level);
1047         pixPaintThroughMask(pixd, pix, bx, by, level);
1048         pixDestroy(&pix);
1049     }
1050 
1051     pixaDestroy(&pixa);
1052     numaDestroy(&na);
1053     return pixd;
1054 }
1055 
1056 
1057 /*!
1058  *  wshedRenderColors()
1059  *
1060  *      Input:  wshed
1061  *      Return: pixd (initial image with all basins filled), or null on error
1062  */
1063 PIX *
wshedRenderColors(L_WSHED * wshed)1064 wshedRenderColors(L_WSHED  *wshed)
1065 {
1066 l_int32  w, h;
1067 PIX     *pixg, *pixt, *pixc, *pixm, *pixd;
1068 PIXA    *pixa;
1069 
1070     PROCNAME("wshedRenderColors");
1071 
1072     if (!wshed)
1073         return (PIX *)ERROR_PTR("wshed not defined", procName, NULL);
1074 
1075     wshedBasins(wshed, &pixa, NULL);
1076     pixg = pixCopy(NULL, wshed->pixs);
1077     pixGetDimensions(wshed->pixs, &w, &h, NULL);
1078     pixd = pixConvertTo32(pixg);
1079     pixt = pixaDisplayRandomCmap(pixa, w, h);
1080     pixc = pixConvertTo32(pixt);
1081     pixm = pixaDisplay(pixa, w, h);
1082     pixCombineMasked(pixd, pixc, pixm);
1083 
1084     pixDestroy(&pixg);
1085     pixDestroy(&pixt);
1086     pixDestroy(&pixc);
1087     pixDestroy(&pixm);
1088     pixaDestroy(&pixa);
1089     return pixd;
1090 }
1091 
1092