# mmcwatershed
Detection of watershed from markers.

## Synopsis

Y = mmcwatershed( f, g, Bc = None, LINEREG = "LINES" )

Implemented in Python.

### Input

f | Image Gray-scale (uint8 or uint16) image.
g | Image Gray-scale (uint8 or uint16) or binary image. marker image: binary or labeled.
Bc | Structuring Element (watershed connectivity) Default: `None` (3x3 elementary cross)
LINEREG | String 'LINES' or ' REGIONS'. Default: `"LINES"`

### Output

Y | Image Gray-scale (uint8 or uint16) or binary image.

## Description

mmcwatershed creates the image `y` by detecting the domain of the catchment basins of `f` indicated by the marker image `g`, according to the connectivity defined by `Bc`. According to the flag LINEREG `y` will be a labeled image of the catchment basins domain or just a binary image that presents the watershed lines. To know more about watershed and watershed from markers, see BeucMeye:93. The implementation of this function is based on LotuFalc:00. WARNING: There is a common mistake related to the marker image `g`. If this image contains only zeros and ones, but it is not a binary image, the result will be an image with all ones. If the marker image is binary, you have to set this explicitly using the `logical` function.

## Examples

    >>> a = uint8([\
        [10,   10,   10,   10,   10,   10,   10],\
        [10,    9,    6,   18,    6,    5,   10],\
        [10,    9,    6,   18,    6,    8,   10],\
        [10,    9,    9,   15,    9,    9,   10],\
        [10,    9,    9,   15,   12,   10,   10],\
        [10,   10,   10,   10,   10,   10,   10]])

    >>> b = mmcmp(a,'==',uint8(6))

    >>> print mmcwatershed(a,b)

    [[0 0 0 1 0 0 0]
     [0 0 0 1 0 0 0]
     [0 0 0 1 0 0 0]
     [0 0 0 1 0 0 0]
     [0 0 0 1 0 0 0]
     [0 0 0 0 1 0 0]]

    >>> print mmcwatershed(a,b,mmsecross(),'REGIONS')

    [[1 1 1 1 2 2 2]
     [1 1 1 1 2 2 2]
     [1 1 1 1 2 2 2]
     [1 1 1 1 2 2 2]
     [1 1 1 1 2 2 2]
     [1 1 1 1 2 2 2]]

    >>> f=mmreadgray('astablet.tif')

    >>> grad=mmgradm(f)

    >>> mark=mmregmin(mmhmin(grad,17))

    >>> w=mmcwatershed(grad,mark)

    >>> mmshow(grad)

    >>> mmshow(mark)

    >>> mmshow(w)

grad | mark |

w |

## Equation

Equation (watershed regions):
Minimum length of a point to a set:
Length of a point to a set:
Where:
is the path between the point x and subset X under connectivity Bc

## Source Code

    def mmcwatershed(f, g, Bc=None, LINEREG="LINES"):
        from Numeric import ones, zeros, nonzero, array, put, take, argmin, transpose, compress, concatenate
        if Bc is None: Bc = mmsecross()
        return g
        print 'starting'
        withline = (LINEREG == 'LINES')
        if mmis(g,'binary'):
            g = mmlabel(g,Bc)
        print 'before 1. mmpad4n'
        status = mmpad4n(uint8(zeros(f.shape)),Bc, 3)
        f = mmpad4n( f,Bc,0)                 #pad input image
        print 'before 2. mmpad4n'
        y = mmpad4n( g,Bc,0)                  # pad marker image with 0
        if withline:
            y1 = mmintersec(mmbinary(y), 0)
        costM = mmlimits(f)[1] * ones(f.shape)  # cummulative cost function image
        mi = nonzero(mmgradm(y,mmsebox(0),Bc).flat)  # 1D index of internal contour of marker
        print 'before put costM'
        put(costM.flat,mi, 0)
        HQueue=transpose([mi, take(costM.flat, mi)])       # init hierarquical queue: index,value
        print 'before mmse2list0'
        Bi=mmse2list0(f,Bc)                # get 1D displacement neighborhood pixels
        x,v = mmmat2set(Bc)
        while HQueue:
            print 'Hq=',HQueue
            i = argmin(HQueue[:,1])           # i is the index of minimum value
            print 'imin=',i
            pi = HQueue[i,0]
            print 'pi=',pi
            ii = ones(HQueue.shape[0])
            ii[i] = 0
            print 'ii=',ii
            HQueue = transpose(array([compress(ii,HQueue[:,0]),
                                      compress(ii,HQueue[:,1])])) # remove this pixel from queue
            print 'H=',HQueue
            put(status.flat, pi, 1)          # make it a permanent label
            for qi in pi+Bi :                # for each neighbor of pi
                if (status.flat[qi] != 3):          # not image border
                    if (status.flat[qi] != 1):        # if not permanent
                        cost_M = max(costM.flat[pi], f.flat[qi])
                        if cost_M < costM.flat[qi]:
                            print 'qi=',qi
                            costM.flat[qi] = cost_M
                            y.flat[qi] = y.flat[pi]                  # propagate the label
                            aux = zeros(array(HQueue.shape) + [1,0])
                            aux[:-1,:] = HQueue
                            aux[-1,:]=[qi, cost_M]
                            HQueue = aux # insert pixel in the queue
                            print 'insert H=',HQueue
                    elif (withline        and
                         (y.flat[qi] != y.flat[pi]) and
                         (y1.flat[pi] == 0)    and
                         (y1.flat[qi] == 0)     ):
                        y1.flat[pi] = 1
        if withline:
            Y = y1
        else:
            Y = y
        return Y

## See also

mmfreedom | Control automatic data type conversion.
mmwatershed | Watershed detection.
mmsebox | Create a box structuring element.
mmsecross | Diamond structuring element and elementary 3x3 cross.
mmswatershed | Detection of similarity-based watershed from markers.
mmdcalc | Extract the keys of a calculator.

[`mmcthin`] [Up] [`mmskiz`] |
Copyright (c) 2003, Roberto A. Lotufo, UNICAMP-University of Campinas; Rubens C. Machado, CenPRA-Renato Archer Research Center.