Euclidean Distance Transform

Average rating
(1 vote)

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.

Back to top