Leptonica  1.83.1
Image processing and image analysis suite
watershed.c
Go to the documentation of this file.
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 
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 
129 {
130  l_int32 x;
131  l_int32 y;
132 };
133 typedef struct L_NewPixel L_NEWPIXEL;
134 
136 struct L_WSPixel
137 {
138  l_float32 val;
139  l_int32 x;
140  l_int32 y;
141  l_int32 index;
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  *-----------------------------------------------------------------------*/
206 L_WSHED *
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 
251 void
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 
305 l_ok
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  *-----------------------------------------------------------------------*/
552 static void
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 
594 static l_int32
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 
705 static l_int32
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 
767 static l_int32
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  *-----------------------------------------------------------------------*/
1017 l_ok
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 
1039 PIX *
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 
1073 PIX *
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 }
#define SET_DATA_BIT(pdata, n)
Definition: arrayaccess.h:127
#define SET_DATA_FOUR_BYTES(pdata, n, val)
Definition: arrayaccess.h:235
#define GET_DATA_BYTE(pdata, n)
Definition: arrayaccess.h:188
#define GET_DATA_FOUR_BYTES(pdata, n)
Definition: arrayaccess.h:231
#define GET_DATA_BIT(pdata, n)
Definition: arrayaccess.h:123
BOX * boxCreate(l_int32 x, l_int32 y, l_int32 w, l_int32 h)
boxCreate()
Definition: boxbasic.c:171
void lheapDestroy(L_HEAP **plh, l_int32 freeflag)
lheapDestroy()
Definition: heap.c:152
L_HEAP * lheapCreate(l_int32 n, l_int32 direction)
lheapCreate()
Definition: heap.c:112
l_int32 lheapGetCount(L_HEAP *lh)
lheapGetCount()
Definition: heap.c:273
l_ok lheapAdd(L_HEAP *lh, void *item)
lheapAdd()
Definition: heap.c:189
void * lheapRemove(L_HEAP *lh)
lheapRemove()
Definition: heap.c:243
l_ok numaAddNumber(NUMA *na, l_float32 val)
numaAddNumber()
Definition: numabasic.c:460
NUMA * numaCreate(l_int32 n)
numaCreate()
Definition: numabasic.c:193
NUMA * numaClone(NUMA *na)
numaClone()
Definition: numabasic.c:414
void numaDestroy(NUMA **pna)
numaDestroy()
Definition: numabasic.c:357
l_ok numaSetValue(NUMA *na, l_int32 index, l_float32 val)
numaSetValue()
Definition: numabasic.c:750
l_int32 numaGetCount(NUMA *na)
numaGetCount()
Definition: numabasic.c:630
l_ok numaGetIValue(NUMA *na, l_int32 index, l_int32 *pival)
numaGetIValue()
Definition: numabasic.c:720
l_int32 * numaGetIArray(NUMA *na)
numaGetIArray()
Definition: numabasic.c:807
l_ok numaJoin(NUMA *nad, NUMA *nas, l_int32 istart, l_int32 iend)
numaJoin()
Definition: numafunc1.c:3521
NUMA * numaMakeConstant(l_float32 val, l_int32 size)
numaMakeConstant()
Definition: numafunc1.c:820
NUMA * numaMakeSequence(l_float32 startval, l_float32 increment, l_int32 size)
numaMakeSequence()
Definition: numafunc1.c:792
void pixDestroy(PIX **ppix)
pixDestroy()
Definition: pix1.c:608
l_ok pixGetDimensions(const PIX *pix, l_int32 *pw, l_int32 *ph, l_int32 *pd)
pixGetDimensions()
Definition: pix1.c:1074
PIX * pixCopy(PIX *pixd, const PIX *pixs)
pixCopy()
Definition: pix1.c:689
PIX * pixCreate(l_int32 width, l_int32 height, l_int32 depth)
pixCreate()
Definition: pix1.c:315
PIX * pixClone(PIX *pixs)
pixClone()
Definition: pix1.c:582
void ** pixGetLinePtrs(PIX *pix, l_int32 *psize)
pixGetLinePtrs()
Definition: pix1.c:1844
l_ok pixSetPixel(PIX *pix, l_int32 x, l_int32 y, l_uint32 val)
pixSetPixel()
Definition: pix2.c:263
l_ok pixGetPixel(PIX *pix, l_int32 x, l_int32 y, l_uint32 *pval)
pixGetPixel()
Definition: pix2.c:192
l_ok pixSetAllArbitrary(PIX *pix, l_uint32 val)
pixSetAllArbitrary()
Definition: pix2.c:929
l_ok pixPaintThroughMask(PIX *pixd, PIX *pixm, l_int32 x, l_int32 y, l_uint32 val)
pixPaintThroughMask()
Definition: pix3.c:618
l_ok pixCombineMasked(PIX *pixd, PIX *pixs, PIX *pixm)
pixCombineMasked()
Definition: pix3.c:378
PIX * pixClipRectangle(PIX *pixs, BOX *box, BOX **pboxc)
pixClipRectangle()
Definition: pix5.c:994
#define PIX_DST
Definition: pix.h:445
@ L_CLONE
Definition: pix.h:506
@ L_INSERT
Definition: pix.h:504
#define PIX_SRC
Definition: pix.h:444
@ L_SORT_INCREASING
Definition: pix.h:522
l_ok pixaAddPix(PIXA *pixa, PIX *pix, l_int32 copyflag)
pixaAddPix()
Definition: pixabasic.c:493
void pixaDestroy(PIXA **ppixa)
pixaDestroy()
Definition: pixabasic.c:404
PIXA * pixaCreate(l_int32 n)
pixaCreate()
Definition: pixabasic.c:167
PIXA * pixaCopy(PIXA *pixa, l_int32 copyflag)
pixaCopy()
Definition: pixabasic.c:442
l_int32 pixaGetCount(PIXA *pixa)
pixaGetCount()
Definition: pixabasic.c:629
l_ok pixaAddBox(PIXA *pixa, BOX *box, l_int32 copyflag)
pixaAddBox()
Definition: pixabasic.c:540
l_ok pixaGetBoxGeometry(PIXA *pixa, l_int32 index, l_int32 *px, l_int32 *py, l_int32 *pw, l_int32 *ph)
pixaGetBoxGeometry()
Definition: pixabasic.c:800
PIX * pixaGetPix(PIXA *pixa, l_int32 index, l_int32 accesstype)
pixaGetPix()
Definition: pixabasic.c:647
PIX * pixaDisplay(PIXA *pixa, l_int32 w, l_int32 h)
pixaDisplay()
Definition: pixafunc2.c:191
PIX * pixaDisplayRandomCmap(PIXA *pixa, l_int32 w, l_int32 h)
pixaDisplayRandomCmap()
Definition: pixafunc2.c:267
PIX * pixConvertTo32(PIX *pixs)
pixConvertTo32()
Definition: pixconv.c:3246
l_ok ptaGetIPt(PTA *pta, l_int32 index, l_int32 *px, l_int32 *py)
ptaGetIPt()
Definition: ptabasic.c:527
l_int32 ptaGetCount(PTA *pta)
ptaGetCount()
Definition: ptabasic.c:480
void ptaDestroy(PTA **ppta)
ptaDestroy()
Definition: ptabasic.c:191
PIX * pixGenerateFromPta(PTA *pta, l_int32 w, l_int32 h)
pixGenerateFromPta()
Definition: ptafunc1.c:1962
l_int32 lqueueGetCount(L_QUEUE *lq)
lqueueGetCount()
Definition: queue.c:276
void lqueueDestroy(L_QUEUE **plq, l_int32 freeflag)
lqueueDestroy()
Definition: queue.c:132
void * lqueueRemove(L_QUEUE *lq)
lqueueRemove()
Definition: queue.c:249
l_ok lqueueAdd(L_QUEUE *lq, void *item)
lqueueAdd()
Definition: queue.c:184
L_QUEUE * lqueueCreate(l_int32 nalloc)
lqueueCreate()
Definition: queue.c:93
l_ok pixRasterop(PIX *pixd, l_int32 dx, l_int32 dy, l_int32 dw, l_int32 dh, l_int32 op, PIX *pixs, l_int32 sx, l_int32 sy)
pixRasterop()
Definition: rop.c:204
l_ok pixSelectMinInConnComp(PIX *pixs, PIX *pixm, PTA **ppta, NUMA **pnav)
pixSelectMinInConnComp()
Definition: seedfill.c:3266
l_ok pixLocalExtrema(PIX *pixs, l_int32 maxmin, l_int32 minmax, PIX **ppixmin, PIX **ppixmax)
pixLocalExtrema()
Definition: seedfill.c:2975
PIX * pixRemoveSeededComponents(PIX *pixd, PIX *pixs, PIX *pixm, l_int32 connectivity, l_int32 bordersize)
pixRemoveSeededComponents()
Definition: seedfill.c:3377
l_int32 lstackGetCount(L_STACK *lstack)
lstackGetCount()
Definition: stack.c:242
void lstackDestroy(L_STACK **plstack, l_int32 freeflag)
lstackDestroy()
Definition: stack.c:122
l_ok lstackAdd(L_STACK *lstack, void *item)
lstackAdd()
Definition: stack.c:166
void * lstackRemove(L_STACK *lstack)
lstackRemove()
Definition: stack.c:196
L_STACK * lstackCreate(l_int32 n)
lstackCreate()
Definition: stack.c:82
Definition: heap.h:78
l_int32 x
Definition: watershed.c:130
l_int32 y
Definition: watershed.c:131
Definition: queue.h:65
struct L_Stack * stack
Definition: queue.h:71
Definition: stack.h:60
l_int32 y
Definition: watershed.c:140
l_int32 x
Definition: watershed.c:139
l_float32 val
Definition: watershed.c:138
l_int32 index
Definition: watershed.c:141
struct Pta * ptas
Definition: watershed.h:50
l_int32 arraysize
Definition: watershed.h:59
struct Pixa * pixad
Definition: watershed.h:49
void ** linet1
Definition: watershed.h:48
void ** linelab32
Definition: watershed.h:47
struct Numa * nash
Definition: watershed.h:52
struct Pix * pixm
Definition: watershed.h:41
struct Numa * nasi
Definition: watershed.h:51
l_int32 mindepth
Definition: watershed.h:42
l_int32 debug
Definition: watershed.h:60
struct Numa * namh
Definition: watershed.h:53
struct Pix * pixs
Definition: watershed.h:40
void ** lines8
Definition: watershed.h:45
struct Pix * pixt
Definition: watershed.h:44
l_int32 * lut
Definition: watershed.h:57
struct Numa * nalevels
Definition: watershed.h:54
l_int32 nseeds
Definition: watershed.h:55
void ** linem1
Definition: watershed.h:46
struct Numa ** links
Definition: watershed.h:58
struct Pix * pixlab
Definition: watershed.h:43
l_int32 nother
Definition: watershed.h:56
void lept_stderr(const char *fmt,...)
lept_stderr()
Definition: utils1.c:306
PIX * wshedRenderFill(L_WSHED *wshed)
wshedRenderFill()
Definition: watershed.c:1040
static l_int32 wshedGetHeight(L_WSHED *wshed, l_int32 val, l_int32 label, l_int32 *pheight)
wshedGetHeight()
Definition: watershed.c:768
static l_int32 identifyWatershedBasin(L_WSHED *wshed, l_int32 index, l_int32 level, BOX **pbox, PIX **ppixd)
identifyWatershedBasin()
Definition: watershed.c:595
L_WSHED * wshedCreate(PIX *pixs, PIX *pixm, l_int32 mindepth, l_int32 debugflag)
wshedCreate()
Definition: watershed.c:207
void wshedDestroy(L_WSHED **pwshed)
wshedDestroy()
Definition: watershed.c:252
static l_int32 mergeLookup(L_WSHED *wshed, l_int32 sindex, l_int32 dindex)
mergeLookup()
Definition: watershed.c:706
l_ok wshedBasins(L_WSHED *wshed, PIXA **ppixa, NUMA **pnalevels)
wshedBasins()
Definition: watershed.c:1018
static void wshedSaveBasin(L_WSHED *wshed, l_int32 index, l_int32 level)
wshedSaveBasin()
Definition: watershed.c:553
l_ok wshedApply(L_WSHED *wshed)
wshedApply()
Definition: watershed.c:306
PIX * wshedRenderColors(L_WSHED *wshed)
wshedRenderColors()
Definition: watershed.c:1074