Euclidean Distance Transform

Thu, 04/18/2013 - 06:55 pm

Euclidean Distance Transform calculated using the "Fast euclidean distance transformation in two scans using a 3x3 neighborhood" algorithm from F.Y. Shih and Y.-T. Wu, Comput. Vis. Image Underst. 93 (2004) 195-205.
This function takes an image mask M_fg and returns a matrix M_DistanceTransform of the same size, containing the Euclidean distance from each foreground pixel to its nearest background pixel. The matrix M_DistanceTransform is created in the current folder, overwriting any existing wave with that name (much like what happens with many forms of the ImageTransform operation). The image mask may be a binary image, or a grayscale image with an optional threshold parameter to specify the background level. The other optional parameters specify the use of periodic boundary conditions or a skeletonization of the EDT, which returns only the values of the distance for those pixels where moving away from the nearest edge pixel changes which is the nearest edge pixel. Typically, the skeletonization gives a centerline of the object with side branches at sharp corners of the object. You can pass a negative number for the skeletonization parameter to attempt to detect only the centerline.

function/WAVE make_FastEDT(M_fg, [bg, wrap,skel])
    WAVE M_fg       // foreground mask indicating pixels to calculate Euclidean distance from
    Variable bg     // background level, below which pixels are considered background
    Variable wrap   // flag to use periodic boundary conditions
    Variable skel   // flag to do skeletonization (0 to skip, >0 for skeletonization, <0 for centerline only)
 
    bg= ParamIsDefault(bg)? 0 : bg
    wrap= ParamIsDefault(wrap)? 0 : wrap
    skel= ParamIsDefault(skel)? 0 : skel
 
    DFREF saveDFR= getDataFolderDFR()
    SetDataFolder GetWavesDataFolderDFR(M_fg)
 
    String dfName= UniqueName("EDT_",11,0)
    NewDataFolder/O/S $dfName
 
    Variable dX= DimDelta(M_fg,0)
    Variable nX= DimSize(M_fg,0)
    Variable dY= DimDelta(M_fg,1)
    Variable nY= DimSize(M_fg,1)
    Variable bigNum= (nX*dX)^2 + (nY*dY)^2
    Make/O/N=(nX+2,nY+2) F=bigNum, Rx=0 , Ry=0, C = NaN
    F[1,nX][1,nY]= (M_fg[p-1][q-1] <= bg )? 0 : bigNum
    if (wrap)
        F[1,nX][0]= F[p][nY]
        F[1,nX][nY+1]= F[p][1]
        F[0][]= F[nX][q]
        F[nX+1][]= F[1][q]
    endif
    Variable iX,iY // indices of current pixel
    Variable pF,qF // current value and test value of squared distance
    Variable pRx,pRy // current value of relative location to nearest bg pixel
    for ( iY=1 ; iY<=nY ; iY+=1 )
        for ( iX=1 ; iX<=nX ; iX+=1 )
            pF= F[iX][iY]
            // immediately skip masked pixel
            if ( pF<=bg )
                Rx[iX][iY]= 0
                Ry[iX][iY]= 0
                C[iX][iY]= 0
                Continue
            endif
            // upper neighbour
            qF= F[iX][iY-1] - 2*Ry[iX][iY-1] + dY
            if ( pF > qF )
                pF= qF
                pRx= Rx[iX][iY-1]
                pRy= Ry[iX][iY-1]-dY
            endif
            // left neighbour
            qF= F[iX-1][iY] - 2*Rx[iX-1][iY] + dX
            if ( pF > qF )
                pF= qF
                pRx= Rx[iX-1][iY]-dX
                pRy= Ry[iX-1][iY]
            endif
            // upper-left neighbour
            qF= F[iX-1][iY-1] - 2*Rx[iX-1][iY-1] + dX - 2*Ry[iX-1][iY-1] + dY
            if ( pF > qF )
                pF= qF
                pRx= Rx[iX-1][iY-1]-dX
                pRy= Ry[iX-1][iY-1]-dY
            endif
            // upper-right neighbour
            qF= F[iX+1][iY-1] + 2*Rx[iX+1][iY-1] + dX - 2*Ry[iX+1][iY-1] + dY
            if ( pF > qF )
                pF= qF
                pRx= Rx[iX+1][iY-1]+dX
                pRy= Ry[iX+1][iY-1]-dY
            endif
            F[iX][iY]= pF
            Rx[iX][iY]= pRx
            Ry[iX][iY]= pRy
        endfor
    endfor
    for ( iY=nY ; iY>=1 ; iY-=1 )
        for ( iX=nX ; iX>=1 ; iX-=1 )
            pF= F[iX][iY]
            // immediately skip masked pixel
            if ( pF==0 )
                Continue
            endif
            pRx= Rx[iX][iY]
            pRy= Ry[iX][iY]
            // lower neighbour
            qF= F[iX][iY+1] + 2*Ry[iX][iY+1] + dY
            if ( pF >= qF )
                pF= qF
                pRx= Rx[iX][iY+1]
                pRy= Ry[iX][iY+1]+dY
            endif
            // right neighbour
            qF= F[iX+1][iY] + 2*Rx[iX+1][iY] + dX
            if ( pF > qF )
                pF= qF
                pRx= Rx[iX+1][iY]+dX
                pRy= Ry[iX+1][iY]
            endif
            // lower-right neighbour
            qF= F[iX+1][iY+1] + 2*Rx[iX+1][iY+1] + dX + 2*Ry[iX+1][iY+1] + dY
            if ( pF > qF )
                pF= qF
                pRx= Rx[iX+1][iY+1]+dX
                pRy= Ry[iX+1][iY+1]+dY
            endif
            // lower-left neighbour
            qF= F[iX-1][iY+1] - 2*Rx[iX-1][iY+1] + dX + 2*Ry[iX-1][iY+1] + dY
            if ( pF > qF )
                pF= qF
                pRx= Rx[iX-1][iY+1]-dX
                pRy= Ry[iX-1][iY+1]+dY
            endif
            F[iX][iY]= pF
            Rx[iX][iY]= pRx
            Ry[iX][iY]= pRy
        endfor
    endfor

    Variable sX,sY // location of pixel to check for skeletonization
    if ( skel ) // check whether the next pixel away from nearest  edge pixel uses different edge
        for ( iY=1 ; iY<=nY ; iY+=1 )
            for ( iX=1 ; iX<=nX ; iX+=1 )
                pF= F[iX][iY]
                if ( pF == 0 )
                    Continue
                endif
                pRx= Rx[iX][iY]
                pRy= Ry[iX][iY]
                if ( abs(pRx) >= abs(pRy) )
                    sX = iX - sign(pRx)
                    sY = iY
                else
                    sX = iX
                    sY = iY - sign(pRy)
                endif
                if ( skel > 0  && (pRx-Rx[sX][sY])^2+(pRy-Ry[sX][sY])^2 > skel^2)
                    // keep values where nearest edge pixels is are not near to each other
                    C[iX][iY] = F
                elseif ( skel < 0 &&(pRx)^2+(pRy)^2 > skel^2 && (pRx+Rx[sX][sY])^2+(pRy+Ry[sX][sY])^2 < skel^2 *pF)
                    // only keep values where nearest edge pixel is (near) diametrically opposed
                    C[iX][iY] = F
                endif
            endfor
        endfor

        MatrixOP/O M_DistanceTransform= sqrt(C)
    else
        MatrixOP/O M_DistanceTransform= sqrt(F)
    endif

    Duplicate/O/R=[1,nX][1,nY] M_DistanceTransform ::M_DistanceTransform
 
    WAVE wResult= ::M_DistanceTransform
    CopyScales M_fg wResult
    KillDataFolder :
 
    SetDataFolder saveDFR
    Return wResult
end


This algorithm assumes that pixelated Voronoi neighborhoods are 8-connected. This is not always true, e.g., for a pixel whose nearest background pixel is at offset (-4,1) with other background pixels at (-3,3) and (-5,0). However, the discrepancy with the true EDT is never more than a small fraction of a pixel.

Forum

Support

Gallery

Igor Pro 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More