/*====================================================================*
 -  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;
}