# Simple distance transform function using Chamfer metrics
# DT(x,y) with integer or fractional metric (Borgefors/Butt and Maragos) and 5x5 template, 101x101 cells
# Programd yields a symmetric (hexadecagon) surface about the sampled point (50,50)
# M J de Smith - 5/12/08 for Python (as a template prog)

# Libraries required:
# Numpy - Numerical Python - see http://numpy.scipy.org/
# Matplotlib - MatLab-like plotting lib - see http://matplotlib.sourceforge.net/index.html

import numeric as np
import matplotlib.pyplot as plt

# define arrays and dimensions
base = np.array((999.0))
DT=np.resize(base,(101,101))
DT[50,50]=0.0 # target point in centre of array
# Local Distance Metric (LDM) and mask arrays
LDM=np.zeros((9));DX=np.zeros((9));DY=np.zeros((9))

# define chamfer values
a1,a2,a3 = 2.2062, 1.4141, 0.9866

# forward scan
LDM[0]=a1; LDM[1]=a1; LDM[2]=a1; LDM[3]=a2; LDM[4]=a3; LDM[5]=a2; LDM[6]=a1; LDM[7]=a3; LDM[8]=0
DX[0]=-2; DX[1]=-2; DX[2]=-1; DX[3]=-1; DX[4]=-1; DX[5]=-1; DX[6]=-1; DX[7]=0;  DX[8]=0
DY[0]=-1; DY[1]=1;  DY[2]=-2; DY[3]=-1; DY[4]=0;  DY[5]=1;  DY[6]=2;  DY[7]=-1; DY[8]=0

for i in range(2,99):
    for j in range(2,99):
        d0=DT[i,j]
        for k in range(0,9):
                r=int(i+DX[k])
                c=int(j+DY[k])
                d=DT[r,c]
                d1=d+LDM[k]
                d0=min(d1,d0)
            # k
        DT[i,j]=d0
    # j
# i

# backwards scan
LDM[0]=0; LDM[1]=a3; LDM[2]=a1; LDM[3]=a2; LDM[4]=a3; LDM[5]=a2; LDM[6]=a1; LDM[7]=a1; LDM[8]=a1
DX[0]=0; DX[1]=0; DX[2]=1;  DX[3]=1;  DX[4]=1; DX[5]=1; DX[6]=1; DX[7]=2;  DX[8]=2
DY[0]=0; DY[1]=1; DY[2]=-2; DY[3]=-1; DY[4]=0; DY[5]=1; DY[6]=2; DY[7]=-1; DY[8]=1

for i in range(98,1,-1):
    for j in range(98,1,-1):
        d0=DT[i,j]
        for k in range(0,9):
                r=int(i+DX[k])
                c=int(j+DY[k])
                d=DT[r,c]+LDM[k]
                d0=min(d,d0)
             # k
        DT[i,j]=d0
    # j
# i
# trim off edges
DT=DT[2:99,2:99]

# plot result
plt.contourf(DT)
cbar = plt.colorbar()
plt.show()