Mercurial > hgrepos > Python2 > PyMuPDF

comparison mupdf-source/thirdparty/leptonica/src/watershed.c @ 2:b50eed0cc0ef upstream

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