[mmcthin] [Up] [mmskiz] Thinning And Thickening

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] Python