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