Euclidean Distance Transform
sbossuyt
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