| /*====================================================================* |
| - Copyright (C) 2001 Leptonica. All rights reserved. |
| - This software is distributed in the hope that it will be |
| - useful, but with NO WARRANTY OF ANY KIND. |
| - No author or distributor accepts responsibility to anyone for the |
| - consequences of using this software, or for whether it serves any |
| - particular purpose or works at all, unless he or she says so in |
| - writing. Everyone is granted permission to copy, modify and |
| - redistribute this source code, for commercial or non-commercial |
| - purposes, with the following restrictions: (1) the origin of this |
| - source code must not be misrepresented; (2) modified versions must |
| - be plainly marked as such; and (3) this notice may not be removed |
| - or altered from any source or modified source distribution. |
| *====================================================================*/ |
| |
| /* |
| * watershed.c |
| * |
| * Top-level |
| * L_WSHED *wshedCreate() |
| * void wshedDestroy() |
| * l_int32 wshedApply() |
| * |
| * Helpers |
| * static l_int32 identifyWatershedBasin() |
| * static l_int32 mergeLookup() |
| * static l_int32 wshedGetHeight() |
| * static void pushNewPixel() |
| * static void popNewPixel() |
| * static void pushWSPixel() |
| * static void popWSPixel() |
| * static void debugPrintLUT() |
| * static void debugWshedMerge() |
| * |
| * Output |
| * l_int32 wshedBasins() |
| * PIX *wshedRenderFill() |
| * PIX *wshedRenderColors() |
| * |
| * The watershed function identifies the "catch basins" of the input |
| * 8 bpp image, with respect to the specified seeds or "markers". |
| * The use is in segmentation, but the selection of the markers is |
| * critical to getting meaningful results. |
| * |
| * How are the markers selected? You can't simply use the local |
| * minima, because a typical image has sufficient noise so that |
| * a useful catch basin can easily have multiple local minima. However |
| * they are selected, the question for the watershed function is |
| * how to handle local minima that are not markers. The reason |
| * this is important is because of the algorithm used to find the |
| * watersheds, which is roughly like this: |
| * |
| * (1) Identify the markers and the local minima, and enter them |
| * into a priority queue based on the pixel value. Each marker |
| * is shrunk to a single pixel, if necessary, before the |
| * operation starts. |
| * (2) Feed the priority queue with neighbors of pixels that are |
| * popped off the queue. Each of these queue pixels is labelled |
| * with the index value of its parent. |
| * (3) Each pixel is also labelled, in a 32-bit image, with the marker |
| * or local minimum index, from which it was originally derived. |
| * (4) There are actually 3 classes of labels: seeds, minima, and |
| * fillers. The fillers are labels of regions that have already |
| * been identified as watersheds and are continuing to fill, for |
| * the purpose of finding higher watersheds. |
| * (5) When a pixel is popped that has already been labelled in the |
| * 32-bit image and that label differs from the label of its |
| * parent (stored in the queue pixel), a boundary has been crossed. |
| * There are several cases: |
| * (a) Both parents are derived from markers but at least one |
| * is not deep enough to become a watershed. Absorb the |
| * shallower basin into the deeper one, fixing the LUT to |
| * redirect the shallower index to the deeper one. |
| * (b) Both parents are derived from markers and both are deep |
| * enough. Identify and save the watershed for each marker. |
| * (c) One parent was derived from a marker and the other from |
| * a minima: absorb the minima basin into the marker basin. |
| * (d) One parent was derived from a marker and the other is |
| * a filler: identify and save the watershed for the marker. |
| * (e) Both parents are derived from minima: merge them. |
| * (f) One parent is a filler and the other is derived from a |
| * minima: merge the minima into the filler. |
| * (6) The output of the watershed operation consists of: |
| * - a pixa of the basins |
| * - a pta of the markers |
| * - a numa of the watershed levels |
| * |
| * Typical usage: |
| * L_WShed *wshed = wshedCreate(pixs, pixseed, mindepth, 0); |
| * wshedApply(wshed); |
| * |
| * wshedBasins(wshed, &pixa, &nalevels); |
| * ... do something with pixa, nalevels ... |
| * pixaDestroy(&pixa); |
| * numaDestroy(&nalevels); |
| * |
| * Pix *pixd = wshedRenderFill(wshed); |
| * |
| * wshedDestroy(&wshed); |
| */ |
| |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include "allheaders.h" |
| |
| #ifndef NO_CONSOLE_IO |
| #define DEBUG_WATERSHED 0 |
| #endif /* ~NO_CONSOLE_IO */ |
| |
| static const l_uint32 MAX_LABEL_VALUE = 0x7fffffff; /* largest l_int32 */ |
| |
| struct L_NewPixel |
| { |
| l_int32 x; |
| l_int32 y; |
| }; |
| typedef struct L_NewPixel L_NEWPIXEL; |
| |
| struct L_WSPixel |
| { |
| l_float32 val; /* pixel value */ |
| l_int32 x; |
| l_int32 y; |
| l_int32 index; /* label for set to which pixel belongs */ |
| }; |
| typedef struct L_WSPixel L_WSPIXEL; |
| |
| |
| /* Static functions for obtaining bitmap of watersheds */ |
| static void wshedSaveBasin(L_WSHED *wshed, l_int32 index, l_int32 level); |
| |
| static l_int32 identifyWatershedBasin(L_WSHED *wshed, |
| l_int32 index, l_int32 level, |
| BOX **pbox, PIX **ppixd); |
| |
| /* Static function for merging lut and backlink arrays */ |
| static l_int32 mergeLookup(L_WSHED *wshed, l_int32 sindex, l_int32 dindex); |
| |
| /* Static function for finding the height of the current pixel |
| above its seed or minima in the watershed. */ |
| static l_int32 wshedGetHeight(L_WSHED *wshed, l_int32 val, l_int32 label, |
| l_int32 *pheight); |
| |
| /* Static accessors for NewPixel on a queue */ |
| static void pushNewPixel(L_QUEUE *lq, l_int32 x, l_int32 y, |
| l_int32 *pminx, l_int32 *pmaxx, |
| l_int32 *pminy, l_int32 *pmaxy); |
| static void popNewPixel(L_QUEUE *lq, l_int32 *px, l_int32 *py); |
| |
| /* Static accessors for WSPixel on a heap */ |
| static void pushWSPixel(L_HEAP *lh, L_STACK *stack, l_int32 val, |
| l_int32 x, l_int32 y, l_int32 index); |
| static void popWSPixel(L_HEAP *lh, L_STACK *stack, l_int32 *pval, |
| l_int32 *px, l_int32 *py, l_int32 *pindex); |
| |
| /* Static debug print output */ |
| static void debugPrintLUT(l_int32 *lut, l_int32 size, l_int32 debug); |
| |
| static void debugWshedMerge(L_WSHED *wshed, char *descr, l_int32 x, |
| l_int32 y, l_int32 label, l_int32 index); |
| |
| |
| /*-----------------------------------------------------------------------* |
| * Top-level watershed * |
| *-----------------------------------------------------------------------*/ |
| /*! |
| * wshedCreate() |
| * |
| * Input: pixs (8 bpp source) |
| * pixm (1 bpp 'marker' seed) |
| * mindepth (minimum depth; anything less is not saved) |
| * debugflag (1 for debug output) |
| * Return: WShed, or null on error |
| * |
| * Notes: |
| * (1) It is not necessary for the fg pixels in the seed image |
| * be at minima, or that they be isolated. We extract a |
| * single pixel from each connected component, and a seed |
| * anywhere in a watershed will eventually label the watershed |
| * when the filling level reaches it. |
| * (2) Set mindepth to some value to ignore noise in pixs that |
| * can create small local minima. Any watershed shallower |
| * than mindepth, even if it has a seed, will not be saved; |
| * It will either be incorporated in another watershed or |
| * eliminated. |
| */ |
| L_WSHED * |
| wshedCreate(PIX *pixs, |
| PIX *pixm, |
| l_int32 mindepth, |
| l_int32 debugflag) |
| { |
| l_int32 w, h; |
| L_WSHED *wshed; |
| |
| PROCNAME("wshedCreate"); |
| |
| if (!pixs) |
| return (L_WSHED *)ERROR_PTR("pixs is not defined", procName, NULL); |
| if (pixGetDepth(pixs) != 8) |
| return (L_WSHED *)ERROR_PTR("pixs is not 8 bpp", procName, NULL); |
| if (!pixm) |
| return (L_WSHED *)ERROR_PTR("pixm is not defined", procName, NULL); |
| if (pixGetDepth(pixm) != 1) |
| return (L_WSHED *)ERROR_PTR("pixm is not 1 bpp", procName, NULL); |
| pixGetDimensions(pixs, &w, &h, NULL); |
| if (pixGetWidth(pixm) != w || pixGetHeight(pixm) != h) |
| return (L_WSHED *)ERROR_PTR("pixs/m sizes are unequal", procName, NULL); |
| |
| if ((wshed = (L_WSHED *)CALLOC(1, sizeof(L_WSHED))) == NULL) |
| return (L_WSHED *)ERROR_PTR("wshed not made", procName, NULL); |
| |
| wshed->pixs = pixClone(pixs); |
| wshed->pixm = pixClone(pixm); |
| wshed->mindepth = L_MAX(1, mindepth); |
| wshed->pixlab = pixCreate(w, h, 32); |
| pixSetAllArbitrary(wshed->pixlab, MAX_LABEL_VALUE); |
| wshed->pixt = pixCreate(w, h, 1); |
| wshed->lines8 = pixGetLinePtrs(pixs, NULL); |
| wshed->linem1 = pixGetLinePtrs(pixm, NULL); |
| wshed->linelab32 = pixGetLinePtrs(wshed->pixlab, NULL); |
| wshed->linet1 = pixGetLinePtrs(wshed->pixt, NULL); |
| wshed->debug = debugflag; |
| return wshed; |
| } |
| |
| |
| /*! |
| * wshedDestroy() |
| * |
| * Input: &wshed (<will be set to null before returning>) |
| * Return: void |
| */ |
| void |
| wshedDestroy(L_WSHED **pwshed) |
| { |
| l_int32 i; |
| L_WSHED *wshed; |
| |
| PROCNAME("wshedDestroy"); |
| |
| if (pwshed == NULL) { |
| L_WARNING("ptr address is null!", procName); |
| return; |
| } |
| |
| if ((wshed = *pwshed) == NULL) |
| return; |
| |
| pixDestroy(&wshed->pixs); |
| pixDestroy(&wshed->pixm); |
| pixDestroy(&wshed->pixlab); |
| pixDestroy(&wshed->pixt); |
| if (wshed->lines8) FREE(wshed->lines8); |
| if (wshed->linem1) FREE(wshed->linem1); |
| if (wshed->linelab32) FREE(wshed->linelab32); |
| if (wshed->linet1) FREE(wshed->linet1); |
| pixaDestroy(&wshed->pixad); |
| ptaDestroy(&wshed->ptas); |
| numaDestroy(&wshed->nash); |
| numaDestroy(&wshed->nasi); |
| numaDestroy(&wshed->namh); |
| numaDestroy(&wshed->nalevels); |
| if (wshed->lut) |
| FREE(wshed->lut); |
| if (wshed->links) { |
| for (i = 0; i < wshed->arraysize; i++) |
| numaDestroy(&wshed->links[i]); |
| FREE(wshed->links); |
| } |
| FREE(wshed); |
| *pwshed = NULL; |
| return; |
| } |
| |
| |
| /*! |
| * wshedApply() |
| * |
| * Input: wshed (generated from wshedCreate()) |
| * Return: 0 if OK, 1 on error |
| * |
| * Iportant note: |
| * (1) This is buggy. It seems to locate watersheds that are |
| * duplicates. The watershed extraction after complete fill |
| * grabs some regions belonging to existing watersheds. |
| * See prog/watershedtest.c for testing. |
| */ |
| l_int32 |
| wshedApply(L_WSHED *wshed) |
| { |
| char two_new_watersheds[] = "Two new watersheds"; |
| char seed_absorbed_into_seeded_basin[] = "Seed absorbed into seeded basin"; |
| char one_new_watershed_label[] = "One new watershed (label)"; |
| char one_new_watershed_index[] = "One new watershed (index)"; |
| char minima_absorbed_into_seeded_basin[] = |
| "Minima absorbed into seeded basin"; |
| char minima_absorbed_by_filler_or_another[] = |
| "Minima absorbed by filler or another"; |
| l_int32 nseeds, nother, nboth, arraysize; |
| l_int32 i, j, val, x, y, w, h, index, mindepth; |
| l_int32 imin, imax, jmin, jmax, cindex, clabel, nindex; |
| l_int32 hindex, hlabel, hmin, hmax, minhindex, maxhindex; |
| l_int32 *lut; |
| l_uint32 ulabel, uval; |
| void **lines8, **linelab32; |
| NUMA *nalut, *nalevels, *nash, *namh, *nasi; |
| NUMA **links; |
| L_HEAP *lh; |
| PIX *pixmin, *pixsd; |
| PIXA *pixad; |
| L_STACK *rstack; |
| PTA *ptas, *ptao; |
| |
| PROCNAME("wshedApply"); |
| |
| if (!wshed) |
| return ERROR_INT("wshed not defined", procName, 1); |
| |
| /* ------------------------------------------------------------ * |
| * Initialize priority queue and pixlab with seeds and minima * |
| * ------------------------------------------------------------ */ |
| |
| lh = lheapCreate(0, L_SORT_INCREASING); /* remove lowest values first */ |
| rstack = lstackCreate(0); /* for reusing the WSPixels */ |
| pixGetDimensions(wshed->pixs, &w, &h, NULL); |
| lines8 = wshed->lines8; /* wshed owns this */ |
| linelab32 = wshed->linelab32; /* ditto */ |
| |
| /* Identify seed (marker) pixels, 1 for each c.c. in pixm */ |
| ptas = pixSelectMinInConnComp(wshed->pixs, wshed->pixm, &nash); |
| pixsd = pixGenerateFromPta(ptas, w, h); |
| nseeds = ptaGetCount(ptas); |
| for (i = 0; i < nseeds; i++) { |
| ptaGetIPt(ptas, i, &x, &y); |
| uval = GET_DATA_BYTE(lines8[y], x); |
| pushWSPixel(lh, rstack, (l_int32)uval, x, y, i); |
| } |
| wshed->ptas = ptas; |
| nasi = numaMakeConstant(1, nseeds); /* indicator array */ |
| wshed->nasi = nasi; |
| wshed->nash = nash; |
| wshed->nseeds = nseeds; |
| |
| /* Identify minima that are not seeds. Use these 4 steps: |
| * (1) Get the local minima, which can have components |
| * of arbitrary size. This will be a clipping mask. |
| * (2) Get the image of the actual seeds (pixsd) |
| * (3) Remove all elements of the clipping mask that have a seed. |
| * (4) Shrink each of the remaining elements of the minima mask |
| * to a single pixel. */ |
| pixLocalExtrema(wshed->pixs, 200, 0, &pixmin, NULL); |
| pixRemoveSeededComponents(pixmin, pixsd, pixmin, 8, 2); |
| ptao = pixSelectMinInConnComp(wshed->pixs, pixmin, &namh); |
| nother = ptaGetCount(ptao); |
| for (i = 0; i < nother; i++) { |
| ptaGetIPt(ptao, i, &x, &y); |
| uval = GET_DATA_BYTE(lines8[y], x); |
| pushWSPixel(lh, rstack, (l_int32)uval, x, y, nseeds + i); |
| } |
| wshed->namh = namh; |
| |
| /* ------------------------------------------------------------ * |
| * Initialize merging lookup tables * |
| * ------------------------------------------------------------ */ |
| |
| /* nalut should always give the current after-merging index. |
| * links are effectively backpointers: they are numas associated with |
| * a dest index of all indices in nalut that point to that index. */ |
| mindepth = wshed->mindepth; |
| nboth = nseeds + nother; |
| arraysize = 2 * nboth; |
| wshed->arraysize = arraysize; |
| nalut = numaMakeSequence(0, 1, arraysize); |
| lut = numaGetIArray(nalut); |
| wshed->lut = lut; /* wshed owns this */ |
| links = (NUMA **)CALLOC(arraysize, sizeof(NUMA *)); |
| wshed->links = links; /* wshed owns this */ |
| nindex = nseeds + nother; /* the next unused index value */ |
| |
| /* ------------------------------------------------------------ * |
| * Fill the basins, using the priority queue * |
| * ------------------------------------------------------------ */ |
| |
| pixad = pixaCreate(nseeds); |
| wshed->pixad = pixad; /* wshed owns this */ |
| nalevels = numaCreate(nseeds); |
| wshed->nalevels = nalevels; /* wshed owns this */ |
| L_INFO_INT2("nseeds = %d, nother = %d\n", procName, nseeds, nother); |
| while (lheapGetCount(lh) > 0) { |
| popWSPixel(lh, rstack, &val, &x, &y, &index); |
| /* fprintf(stderr, "x = %d, y = %d, index = %d\n", x, y, index); */ |
| ulabel = GET_DATA_FOUR_BYTES(linelab32[y], x); |
| if (ulabel == MAX_LABEL_VALUE) |
| clabel = ulabel; |
| else |
| clabel = lut[ulabel]; |
| cindex = lut[index]; |
| if (clabel == cindex) continue; /* have already seen this one */ |
| if (clabel == MAX_LABEL_VALUE) { /* new one; assign index and try to |
| * propagate to all neighbors */ |
| SET_DATA_FOUR_BYTES(linelab32[y], x, cindex); |
| imin = L_MAX(0, y - 1); |
| imax = L_MIN(h - 1, y + 1); |
| jmin = L_MAX(0, x - 1); |
| jmax = L_MIN(w - 1, x + 1); |
| for (i = imin; i <= imax; i++) { |
| for (j = jmin; j <= jmax; j++) { |
| if (i == y && j == x) continue; |
| uval = GET_DATA_BYTE(lines8[i], j); |
| pushWSPixel(lh, rstack, (l_int32)uval, j, i, cindex); |
| } |
| } |
| } |
| else { /* this pixel is already labeled (differently); must resolve */ |
| |
| /* If both indices are seeds, check if the min height is |
| * greater than mindepth. If so, we have two new watersheds; |
| * locate them and assign to both regions a new index |
| * for further waterfill. If not, absorb the shallower |
| * watershed into the deeper one and continue filling it. */ |
| pixGetPixel(pixsd, x, y, &uval); |
| if (clabel < nseeds && cindex < nseeds) { |
| wshedGetHeight(wshed, val, clabel, &hlabel); |
| wshedGetHeight(wshed, val, cindex, &hindex); |
| hmin = L_MIN(hlabel, hindex); |
| hmax = L_MAX(hlabel, hindex); |
| if (hmin == hmax) { |
| hmin = hlabel; |
| hmax = hindex; |
| } |
| if (wshed->debug) { |
| fprintf(stderr, "clabel,hlabel = %d,%d\n", clabel, hlabel); |
| fprintf(stderr, "hmin = %d, hmax = %d\n", hmin, hmax); |
| fprintf(stderr, "cindex,hindex = %d,%d\n", cindex, hindex); |
| if (hmin < mindepth) |
| fprintf(stderr, "Too shallow!\n"); |
| } |
| |
| if (hmin >= mindepth) { |
| debugWshedMerge(wshed, two_new_watersheds, |
| x, y, clabel, cindex); |
| wshedSaveBasin(wshed, cindex, val - 1); |
| wshedSaveBasin(wshed, clabel, val - 1); |
| numaSetValue(nasi, cindex, 0); |
| numaSetValue(nasi, clabel, 0); |
| |
| if (wshed->debug) fprintf(stderr, "nindex = %d\n", nindex); |
| debugPrintLUT(lut, nindex, wshed->debug); |
| mergeLookup(wshed, clabel, nindex); |
| debugPrintLUT(lut, nindex, wshed->debug); |
| mergeLookup(wshed, cindex, nindex); |
| debugPrintLUT(lut, nindex, wshed->debug); |
| nindex++; |
| } |
| else /* extraneous seed within seeded basin; absorb */ |
| debugWshedMerge(wshed, seed_absorbed_into_seeded_basin, |
| x, y, clabel, cindex); |
| maxhindex = clabel; |
| minhindex = cindex; |
| if (hindex > hlabel) { |
| maxhindex = cindex; |
| minhindex = clabel; |
| } |
| mergeLookup(wshed, minhindex, maxhindex); |
| } |
| |
| /* If one index is a seed and the other is a merge of |
| * 2 watersheds, generate a single watershed. */ |
| else if (clabel < nseeds && cindex >= nboth) { |
| debugWshedMerge(wshed, one_new_watershed_label, |
| x, y, clabel, cindex); |
| wshedSaveBasin(wshed, clabel, val - 1); |
| numaSetValue(nasi, clabel, 0); |
| mergeLookup(wshed, clabel, cindex); |
| } |
| else if (cindex < nseeds && clabel >= nboth) { |
| debugWshedMerge(wshed, one_new_watershed_index, |
| x, y, clabel, cindex); |
| wshedSaveBasin(wshed, cindex, val - 1); |
| numaSetValue(nasi, cindex, 0); |
| mergeLookup(wshed, cindex, clabel); |
| } |
| /* If one index is a seed and the other is from a minimum, |
| * merge the minimum wshed into the seed wshed. */ |
| else if (clabel < nseeds) { /* cindex from minima; absorb */ |
| debugWshedMerge(wshed, minima_absorbed_into_seeded_basin, |
| x, y, clabel, cindex); |
| mergeLookup(wshed, cindex, clabel); |
| } |
| else if (cindex < nseeds) { /* clabel from minima; absorb */ |
| debugWshedMerge(wshed, minima_absorbed_into_seeded_basin, |
| x, y, clabel, cindex); |
| mergeLookup(wshed, clabel, cindex); |
| } |
| /* If neither index is a seed, just merge */ |
| else { |
| debugWshedMerge(wshed, minima_absorbed_by_filler_or_another, |
| x, y, clabel, cindex); |
| mergeLookup(wshed, clabel, cindex); |
| } |
| } |
| } |
| |
| #if 0 |
| /* Use the indicator array to save any watersheds that fill |
| * to the maximum value. This seems to screw things up! */ |
| for (i = 0; i < nseeds; i++) { |
| numaGetIValue(nasi, i, &ival); |
| if (ival == 1) { |
| wshedSaveBasin(wshed, lut[i], val - 1); |
| numaSetValue(nasi, i, 0); |
| } |
| } |
| #endif |
| |
| numaDestroy(&nalut); |
| pixDestroy(&pixmin); |
| pixDestroy(&pixsd); |
| ptaDestroy(&ptao); |
| lheapDestroy(&lh, TRUE); |
| lstackDestroy(&rstack, TRUE); |
| return 0; |
| } |
| |
| |
| /*-----------------------------------------------------------------------* |
| * Helpers * |
| *-----------------------------------------------------------------------*/ |
| /*! |
| * wshedSaveBasin() |
| * |
| * Input: wshed |
| * index (index of basin to be located) |
| * level (filling level reached at the time this function |
| * is called) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) This identifies a single watershed. It does not change |
| * the LUT, which must be done subsequently. |
| * (2) The fill level of a basin is taken to be @level - 1. |
| */ |
| static void |
| wshedSaveBasin(L_WSHED *wshed, |
| l_int32 index, |
| l_int32 level) |
| { |
| BOX *box; |
| PIX *pix; |
| |
| PROCNAME("wshedSaveBasin"); |
| |
| if (!wshed) |
| return ERROR_VOID("wshed not defined", procName); |
| |
| if (identifyWatershedBasin(wshed, index, level, &box, &pix) == 0) { |
| pixaAddPix(wshed->pixad, pix, L_INSERT); |
| pixaAddBox(wshed->pixad, box, L_INSERT); |
| numaAddNumber(wshed->nalevels, level - 1); |
| } |
| return; |
| } |
| |
| |
| /*! |
| * identifyWatershedBasin() |
| * |
| * Input: wshed |
| * index (index of basin to be located) |
| * level (of basin at point at which the two basins met) |
| * &box (<return> bounding box of basin) |
| * &pixd (<return> pix of basin, cropped to its bounding box) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) This is a static function, so we assume pixlab, pixs and pixt |
| * exist and are the same size. |
| * (2) It selects all pixels that have the label @index in pixlab |
| * and that have a value in pixs that is less than @level. |
| * (3) It is used whenever two seeded basins meet (typically at a saddle), |
| * or when one seeded basin meets a 'filler'. All identified |
| * basins are saved as a watershed. |
| */ |
| static l_int32 |
| identifyWatershedBasin(L_WSHED *wshed, |
| l_int32 index, |
| l_int32 level, |
| BOX **pbox, |
| PIX **ppixd) |
| { |
| l_int32 imin, imax, jmin, jmax, minx, miny, maxx, maxy; |
| l_int32 bw, bh, i, j, w, h, x, y; |
| l_int32 *lut; |
| l_uint32 label, bval, lval; |
| void **lines8, **linelab32, **linet1; |
| BOX *box; |
| PIX *pixs, *pixlab, *pixt, *pixd; |
| L_QUEUE *lq; |
| |
| PROCNAME("identifyWatershedBasin"); |
| |
| if (!pbox) |
| return ERROR_INT("&box not defined", procName, 1); |
| *pbox = NULL; |
| if (!ppixd) |
| return ERROR_INT("&pixd not defined", procName, 1); |
| *ppixd = NULL; |
| if (!wshed) |
| return ERROR_INT("wshed not defined", procName, 1); |
| |
| /* Make a queue and an auxiliary stack */ |
| lq = lqueueCreate(0); |
| lq->stack = lstackCreate(0); |
| |
| pixs = wshed->pixs; |
| pixlab = wshed->pixlab; |
| pixt = wshed->pixt; |
| lines8 = wshed->lines8; |
| linelab32 = wshed->linelab32; |
| linet1 = wshed->linet1; |
| lut = wshed->lut; |
| pixGetDimensions(pixs, &w, &h, NULL); |
| |
| /* Prime the queue with the seed pixel for this watershed. */ |
| minx = miny = 1000000; |
| maxx = maxy = 0; |
| ptaGetIPt(wshed->ptas, index, &x, &y); |
| pixSetPixel(pixt, x, y, 1); |
| pushNewPixel(lq, x, y, &minx, &maxx, &miny, &maxy); |
| if (wshed->debug) fprintf(stderr, "prime: (x,y) = (%d, %d)\n", x, y); |
| |
| /* Each pixel in a spreading breadth-first search is inspected. |
| * It is accepted as part of this watershed, and pushed on |
| * the search queue, if: |
| * (1) It has a label value equal to @index |
| * (2) The pixel value is less than @level, the overflow |
| * height at which the two basins join. |
| * (3) It has not yet been seen in this search. */ |
| while (lqueueGetCount(lq) > 0) { |
| popNewPixel(lq, &x, &y); |
| imin = L_MAX(0, y - 1); |
| imax = L_MIN(h - 1, y + 1); |
| jmin = L_MAX(0, x - 1); |
| jmax = L_MIN(w - 1, x + 1); |
| for (i = imin; i <= imax; i++) { |
| for (j = jmin; j <= jmax; j++) { |
| if (j == x && i == y) continue; /* parent */ |
| label = GET_DATA_FOUR_BYTES(linelab32[i], j); |
| if (label == MAX_LABEL_VALUE || lut[label] != index) continue; |
| bval = GET_DATA_BIT(linet1[i], j); |
| if (bval == 1) continue; /* already seen */ |
| lval = GET_DATA_BYTE(lines8[i], j); |
| if (lval >= level) continue; /* too high */ |
| SET_DATA_BIT(linet1[i], j); |
| pushNewPixel(lq, j, i, &minx, &maxx, &miny, &maxy); |
| } |
| } |
| } |
| |
| /* Extract the box and pix, and clear pixt */ |
| bw = maxx - minx + 1; |
| bh = maxy - miny + 1; |
| box = boxCreate(minx, miny, bw, bh); |
| pixd = pixClipRectangle(pixt, box, NULL); |
| pixRasterop(pixt, minx, miny, bw, bh, PIX_SRC ^ PIX_DST, pixd, 0, 0); |
| *pbox = box; |
| *ppixd = pixd; |
| |
| lqueueDestroy(&lq, 1); |
| return 0; |
| } |
| |
| |
| /*! |
| * mergeLookup() |
| * |
| * Input: wshed |
| * sindex (primary index being changed in the merge) |
| * dindex (index that @sindex will point to after the merge) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) The links are a sparse array of Numas showing current back-links. |
| * The lut gives the current index (of the seed or the minima |
| * for the wshed in which it is located. |
| * (2) Think of each entry in the lut. There are two types: |
| * owner: lut[index] = index |
| * redirect: lut[index] != index |
| * (3) This is called each time a merge occurs. It puts the lut |
| * and backlinks in a canonical form after the merge, where |
| * all entries in the lut point to the current "owner", which |
| * has all backlinks. That is, every "redirect" in the lut |
| * points to an "owner". The lut always gives the index of |
| * the current owner. |
| */ |
| static l_int32 |
| mergeLookup(L_WSHED *wshed, |
| l_int32 sindex, |
| l_int32 dindex) |
| { |
| l_int32 i, n, size, index; |
| l_int32 *lut; |
| NUMA *na; |
| NUMA **links; |
| |
| PROCNAME("mergeLookup"); |
| |
| if (!wshed) |
| return ERROR_INT("wshed not defined", procName, 1); |
| size = wshed->arraysize; |
| if (sindex < 0 || sindex >= size) |
| return ERROR_INT("invalid sindex", procName, 1); |
| if (dindex < 0 || dindex >= size) |
| return ERROR_INT("invalid dindex", procName, 1); |
| |
| /* Redirect links in the lut */ |
| n = 0; |
| links = wshed->links; |
| lut = wshed->lut; |
| if ((na = links[sindex]) != NULL) { |
| n = numaGetCount(na); |
| for (i = 0; i < n; i++) { |
| numaGetIValue(na, i, &index); |
| lut[index] = dindex; |
| } |
| } |
| lut[sindex] = dindex; |
| |
| /* Shift the backlink arrays from sindex to dindex. |
| * sindex should have no backlinks because all entries in the |
| * lut that were previously pointing to it have been redirected |
| * to dindex. */ |
| if (!links[dindex]) |
| links[dindex] = numaCreate(n); |
| numaJoin(links[dindex], links[sindex], 0, 0); |
| numaAddNumber(links[dindex], sindex); |
| numaDestroy(&links[sindex]); |
| |
| return 0; |
| } |
| |
| |
| /*! |
| * wshedGetHeight() |
| * |
| * Input: wshed (array of current indices) |
| * val (value of current pixel popped off queue) |
| * label (of pixel or 32 bpp label image) |
| * &height (<return> height of current value from seed |
| * or minimum of watershed) |
| * Return: 0 if OK, 1 on error |
| * |
| * Notes: |
| * (1) It is only necessary to find the height for a watershed |
| * that is indexed by a seed or a minima. This function should |
| * not be called on a finished watershed (that continues to fill). |
| */ |
| static l_int32 |
| wshedGetHeight(L_WSHED *wshed, |
| l_int32 val, |
| l_int32 label, |
| l_int32 *pheight) |
| { |
| l_int32 minval; |
| |
| PROCNAME("wshedGetHeight"); |
| |
| if (!pheight) |
| return ERROR_INT("&height not defined", procName, 1); |
| *pheight = 0; |
| if (!wshed) |
| return ERROR_INT("wshed not defined", procName, 1); |
| |
| if (label < wshed->nseeds) |
| numaGetIValue(wshed->nash, label, &minval); |
| else if (label < wshed->nseeds + wshed->nother) |
| numaGetIValue(wshed->namh, label, &minval); |
| else |
| return ERROR_INT("finished watershed; should not call", procName, 1); |
| |
| *pheight = val - minval; |
| return 0; |
| } |
| |
| |
| /* |
| * pushNewPixel() |
| * |
| * Input: lqueue |
| * x, y (pixel coordinates) |
| * &minx, &maxx, &miny, &maxy (<return> bounding box update) |
| * Return: void |
| * |
| * Notes: |
| * (1) This is a wrapper for adding a NewPixel to a queue, which |
| * updates the bounding box for all pixels on that queue and |
| * uses the storage stack to retrieve a NewPixel. |
| */ |
| static void |
| pushNewPixel(L_QUEUE *lq, |
| l_int32 x, |
| l_int32 y, |
| l_int32 *pminx, |
| l_int32 *pmaxx, |
| l_int32 *pminy, |
| l_int32 *pmaxy) |
| { |
| L_NEWPIXEL *np; |
| |
| PROCNAME("pushNewPixel"); |
| |
| if (!lq) |
| return ERROR_VOID(procName, "queue not defined"); |
| |
| /* Adjust bounding box */ |
| *pminx = L_MIN(*pminx, x); |
| *pmaxx = L_MAX(*pmaxx, x); |
| *pminy = L_MIN(*pminy, y); |
| *pmaxy = L_MAX(*pmaxy, y); |
| |
| /* Get a newpixel to use */ |
| if (lstackGetCount(lq->stack) > 0) |
| np = (L_NEWPIXEL *)lstackRemove(lq->stack); |
| else |
| np = (L_NEWPIXEL *)CALLOC(1, sizeof(L_NEWPIXEL)); |
| |
| np->x = x; |
| np->y = y; |
| lqueueAdd(lq, np); |
| return; |
| } |
| |
| |
| /* |
| * popNewPixel() |
| * |
| * Input: lqueue |
| * &x, &y (<return> pixel coordinates) |
| * Return: void |
| * |
| * Notes: |
| * (1) This is a wrapper for removing a NewPixel from a queue, |
| * which returns the pixel coordinates and saves the NewPixel |
| * on the storage stack. |
| */ |
| static void |
| popNewPixel(L_QUEUE *lq, |
| l_int32 *px, |
| l_int32 *py) |
| { |
| L_NEWPIXEL *np; |
| |
| PROCNAME("popNewPixel"); |
| |
| if (!lq) |
| return ERROR_VOID(procName, "lqueue not defined"); |
| |
| if ((np = (L_NEWPIXEL *)lqueueRemove(lq)) == NULL) |
| return; |
| *px = np->x; |
| *py = np->y; |
| lstackAdd(lq->stack, np); /* save for re-use */ |
| return; |
| } |
| |
| |
| /* |
| * pushWSPixel() |
| * |
| * Input: lh (priority queue) |
| * stack (of reusable WSPixels) |
| * val (pixel value: used for ordering the heap) |
| * x, y (pixel coordinates) |
| * index (label for set to which pixel belongs) |
| * Return: void |
| * |
| * Notes: |
| * (1) This is a wrapper for adding a WSPixel to a heap. It |
| * uses the storage stack to retrieve a WSPixel. |
| */ |
| static void |
| pushWSPixel(L_HEAP *lh, |
| L_STACK *stack, |
| l_int32 val, |
| l_int32 x, |
| l_int32 y, |
| l_int32 index) |
| { |
| L_WSPIXEL *wsp; |
| |
| PROCNAME("pushWSPixel"); |
| |
| if (!lh) |
| return ERROR_VOID(procName, "heap not defined"); |
| if (!stack) |
| return ERROR_VOID(procName, "stack not defined"); |
| |
| /* Get a wspixel to use */ |
| if (lstackGetCount(stack) > 0) |
| wsp = (L_WSPIXEL *)lstackRemove(stack); |
| else |
| wsp = (L_WSPIXEL *)CALLOC(1, sizeof(L_WSPIXEL)); |
| |
| wsp->val = (l_float32)val; |
| wsp->x = x; |
| wsp->y = y; |
| wsp->index = index; |
| lheapAdd(lh, wsp); |
| return; |
| } |
| |
| |
| /* |
| * popWSPixel() |
| * |
| * Input: lh (priority queue) |
| * stack (of reusable WSPixels) |
| * &val (<return> pixel value) |
| * &x, &y (<return> pixel coordinates) |
| * &index (<return> label for set to which pixel belongs) |
| * Return: void |
| * |
| * Notes: |
| * (1) This is a wrapper for removing a WSPixel from a heap, |
| * which returns the WSPixel data and saves the WSPixel |
| * on the storage stack. |
| */ |
| static void |
| popWSPixel(L_HEAP *lh, |
| L_STACK *stack, |
| l_int32 *pval, |
| l_int32 *px, |
| l_int32 *py, |
| l_int32 *pindex) |
| { |
| L_WSPIXEL *wsp; |
| |
| PROCNAME("popWSPixel"); |
| |
| if (!lh) |
| return ERROR_VOID(procName, "lheap not defined"); |
| if (!stack) |
| return ERROR_VOID(procName, "stack not defined"); |
| if (!pval || !px || !py || !pindex) |
| return ERROR_VOID(procName, "data can't be returned"); |
| |
| if ((wsp = (L_WSPIXEL *)lheapRemove(lh)) == NULL) |
| return; |
| *pval = (l_int32)wsp->val; |
| *px = wsp->x; |
| *py = wsp->y; |
| *pindex = wsp->index; |
| lstackAdd(stack, wsp); /* save for re-use */ |
| return; |
| } |
| |
| |
| static void |
| debugPrintLUT(l_int32 *lut, |
| l_int32 size, |
| l_int32 debug) |
| { |
| l_int32 i; |
| |
| if (!debug) return; |
| fprintf(stderr, "lut: "); |
| for (i = 0; i < size; i++) |
| fprintf(stderr, "%d ", lut[i]); |
| fprintf(stderr, "\n"); |
| return; |
| } |
| |
| |
| static void |
| debugWshedMerge(L_WSHED *wshed, |
| char *descr, |
| l_int32 x, |
| l_int32 y, |
| l_int32 label, |
| l_int32 index) |
| { |
| if (!wshed || (wshed->debug == 0)) |
| return; |
| fprintf(stderr, "%s:\n", descr); |
| fprintf(stderr, " (x, y) = (%d, %d)\n", x, y); |
| fprintf(stderr, " clabel = %d, cindex = %d\n", label, index); |
| return; |
| } |
| |
| |
| /*-----------------------------------------------------------------------* |
| * Output * |
| *-----------------------------------------------------------------------*/ |
| /*! |
| * wshedBasins() |
| * |
| * Input: wshed |
| * &pixa (<optional return> mask of watershed basins) |
| * &nalevels (<optional return> watershed levels) |
| * Return: 0 if OK, 1 on error |
| */ |
| l_int32 |
| wshedBasins(L_WSHED *wshed, |
| PIXA **ppixa, |
| NUMA **pnalevels) |
| { |
| PROCNAME("wshedBasins"); |
| |
| if (!wshed) |
| return ERROR_INT("wshed not defined", procName, 1); |
| |
| if (ppixa) |
| *ppixa = pixaCopy(wshed->pixad, L_CLONE); |
| if (pnalevels) |
| *pnalevels = numaClone(wshed->nalevels); |
| return 0; |
| } |
| |
| |
| /*! |
| * wshedRenderFill() |
| * |
| * Input: wshed |
| * Return: pixd (initial image with all basins filled), or null on error |
| */ |
| PIX * |
| wshedRenderFill(L_WSHED *wshed) |
| { |
| l_int32 i, n, level, bx, by; |
| NUMA *na; |
| PIX *pix, *pixd; |
| PIXA *pixa; |
| |
| PROCNAME("wshedRenderFill"); |
| |
| if (!wshed) |
| return (PIX *)ERROR_PTR("wshed not defined", procName, NULL); |
| |
| wshedBasins(wshed, &pixa, &na); |
| pixd = pixCopy(NULL, wshed->pixs); |
| n = pixaGetCount(pixa); |
| for (i = 0; i < n; i++) { |
| pix = pixaGetPix(pixa, i, L_CLONE); |
| pixaGetBoxGeometry(pixa, i, &bx, &by, NULL, NULL); |
| numaGetIValue(na, i, &level); |
| pixPaintThroughMask(pixd, pix, bx, by, level); |
| pixDestroy(&pix); |
| } |
| |
| pixaDestroy(&pixa); |
| numaDestroy(&na); |
| return pixd; |
| } |
| |
| |
| /*! |
| * wshedRenderColors() |
| * |
| * Input: wshed |
| * Return: pixd (initial image with all basins filled), or null on error |
| */ |
| PIX * |
| wshedRenderColors(L_WSHED *wshed) |
| { |
| l_int32 w, h; |
| PIX *pixg, *pixt, *pixc, *pixm, *pixd; |
| PIXA *pixa; |
| |
| PROCNAME("wshedRenderColors"); |
| |
| if (!wshed) |
| return (PIX *)ERROR_PTR("wshed not defined", procName, NULL); |
| |
| wshedBasins(wshed, &pixa, NULL); |
| pixg = pixCopy(NULL, wshed->pixs); |
| pixGetDimensions(wshed->pixs, &w, &h, NULL); |
| pixd = pixConvertTo32(pixg); |
| pixt = pixaDisplayRandomCmap(pixa, w, h); |
| pixc = pixConvertTo32(pixt); |
| pixm = pixaDisplay(pixa, w, h); |
| pixCombineMasked(pixd, pixc, pixm); |
| |
| pixDestroy(&pixg); |
| pixDestroy(&pixt); |
| pixDestroy(&pixc); |
| pixDestroy(&pixm); |
| pixaDestroy(&pixa); |
| return pixd; |
| } |
| |
