欧氏距离变换以后再添加。
下面分别给出计算欧式距离（EDT）的matlab和C代码。
参考网址
不清楚上面的网址是什么情况，但是有很多的matlab的mex文件的源码。
matlab：
1 % 下面的为欧式距离变换的代码 2 function Dr = mybwdist(binaryimg) 3 [height, width]=size(binaryimg); 4 D=[]; 5 for k = 1: width % 由于matlab按列存储，所以内层遍历行，外层遍历列。 6 D=[D calcFirstPassD1(binaryimg(:,k), height)]; 7 end 8 9 Dr=processDimN([height, width], D); 10 11 end12 13 function D=calcFirstPassD1(I, vectorLength)14 D1=zeros(vectorLength,1);15 16 % // Initialize D1 - Feature points get set to zero, -1 otherwise17 % D1(i) = (I(col*nrows+i) == 1) ? 0.0f : -1.0f;18 D1(I==0)= (-1);19 20 % // Process column 21 D=voronoiEDT(D1); 22 23 end24 25 function D=voronoiEDT(D)26 27 % // Note: g and h are working vectors allocated in calling function28 [ns, g, h] = constructPartialVoronoi(D);29 if ns == -130 return;31 end32 33 D=queryPartialVoronoi(g, h, ns, D);34 35 end 36 37 function [el, g, h]=constructPartialVoronoi(D)38 % //39 % // Construct partial voronoi diagram (see Maurer et al., 2003, Fig. 5, lines 1-14)40 % // Note: short variable names are used to mimic the notation of the paper41 % //42 g=zeros(1,length(D));43 h=zeros(1,length(D));44 el = -1;45 for k = 0:length(D)-1 % 在matlab中对应数组的高度46 fk = D(k+1);47 if fk ~= -148 if (el < 1)49 el=el+1;50 g(el+1) = fk;51 h(el+1) = k;52 else 53 while ( (el >= 1) && RemoveEDT(g(el), g(el+1), fk, h(el), h(el+1), k))54 el=el-1;55 end 56 el=el+1;57 g(el+1) = fk;58 h(el+1) = k;59 end60 end61 end62 63 end64 65 function res = RemoveEDT(du, dv, dw, u, v, w)66 a = v - u;67 b = w - v;68 c = w - u;69 70 % // See Eq. 2 from Maurer et al., 200371 res=(c * dv - b * du - a * dw) > (a * b * c);72 end73 74 function D=queryPartialVoronoi(g, h, ns, D)75 76 % //77 % // Query partial Voronoi diagram (see Maurer et al., 2003, Fig. 5, lines 18-24)78 % //79 el = 0;80 for k = 0:length(D)-181 while ( (el < ns) && ((g(el+1) + ((h(el+1) - k)*(h(el+1) - k))) > (g(el+2) + ((h(el+2) - k)*(h(el+2) - k)))) ) 82 el=el+1;83 end84 D(k+1) = g(el+1) + (h(el+1) - k)*(h(el+1) - k);85 end86 end87 88 function D3=processDimN(dims, D)89 % // Create temporary vectors for processing along dimension N90 D3=[];91 for k = 0 : dims(1)-1 % // 若为二维数组，则得到第一维元素总数，即行数。注意，matlab按列存储，C按行存储。92 % // Generate indices to points along vector path93 D3=[D3;voronoiEDT(D(k+1,:))];94 end95 end
分别使用matlab自带的bwdist、上面的代码进行测试，测试结果一样。
C++代码（进行了简化，将不需要的代码合并了，同时，有些代码并非完全对应，但是结果正确。注意，matlab调用时需要取反，下面的C代码不再需要取反，如107行）：
1 // Removal Function - FVs that are the Voronoi sites of Voronoi cells that do not intersect Rd are removed 2 inline bool RemoveEDT(const double du, const double dv, const double dw, const double u, const double v, const double w) 3 { 4 double a = v - u; 5 double b = w - v; 6 double c = w - u; 7 8 // See Eq. 2 from Maurer et al., 2003 9 return ((c * dv - b * du - a * dw) > (a * b * c)); 10 } 11 12 inline int constructPartialVoronoi(double* D, double* g, int* h, int length) 13 { 14 // 15 // Construct partial voronoi diagram (see Maurer et al., 2003, Fig. 5, lines 1-14) 16 // Note: short variable names are used to mimic the notation of the paper 17 // 18 int el = -1; 19 20 for (int k = 0; k < length; k++) 21 { 22 double fk = D[k]; 23 if (fk != -1.0f) 24 { 25 if (el < 1) 26 { 27 el++; 28 g[el] = fk; 29 h[el] = k; 30 } 31 else 32 { 33 while ((el >= 1) && RemoveEDT(g[el - 1], g[el], fk, h[el - 1], h[el], k)) 34 { 35 el--; 36 } 37 el++; 38 g[el] = fk; 39 h[el] = k; 40 } 41 } 42 } 43 44 return(el); 45 } 46 47 inline void queryPartialVoronoi(const double* g, const int* h, const int ns, double* D, int length) 48 { 49 // 50 // Query partial Voronoi diagram (see Maurer et al., 2003, Fig. 5, lines 18-24) 51 // 52 int el = 0; 53 54 for (int k = 0; k < length; k++) 55 { 56 while ((el < ns) && ((g[el] + ((h[el] - k)*(h[el] - k))) >(g[el + 1] + ((h[el + 1] - k)*(h[el + 1] - k))))) 57 { 58 el++; 59 } 60 D[k] = g[el] + (h[el] - k)*(h[el] - k); 61 } 62 } 63 64 void voronoiEDT(double* D, double* g, int* h, int length) 65 { 66 // Note: g and h are working vectors allocated in calling function 67 int ns = constructPartialVoronoi(D, g, h, length); 68 if (ns == -1) 69 return; 70 71 queryPartialVoronoi(g, h, ns, D, length); 72 73 return; 74 } 75 76 inline void processDimN(int width, int height/*const mwSize *dims, const mwSize ndims*/, double *D) 77 { 78 // Create temporary vectors for processing along dimension N 79 80 double* g = new double[width]; 81 int* h = new int[width]; 82 83 // 若为二维数组，则得到第一维元素总数，即行数。注意，matlab按列存储，C按行存储。 84 for (int k = 0; k < height; k++) 85 { 86 // Process vector 87 voronoiEDT(&D[k*width], g, h, width); 88 } 89 90 delete[] g; 91 delete[] h; 92 } 93 94 // mxI为输入，unsigned char*类型，mxD为输出，double* 类型 95 // 注意，matlab中mxI是逻辑类型，需要1-mxI才是真正的输入；此处省略了这一步。 96 // 输入为二值化图像，大于阈值的不为0，小于阈值的为0。 97 void computeEDT(double *mxD, const unsigned char *mxI, int width, int height) 98 { 99 double* D1 = new double[width];100 double* g = new double[width];101 int* h = new int[width];102 103 for (int k = 0; k < width; k++)104 {105 for (int i = 0; i < height; i++)106 {107 D1[i] = (mxI[i*width + k] == 0) ? 0.0f : -1.0f;108 }109 110 voronoiEDT(D1, g, h, height);111 112 for (int i = 0; i < height; i++)113 {114 mxD[i*width + k] = D1[i];115 }116 }117 118 delete[] D1;119 delete[] g;120 delete[] h;121 122 processDimN(width, height, mxD);123 } matlab代码对320*240的图像处理，耗时200ms左右（记不清了）；当然，matlab使用的自带的代码耗时2ms左右。C中debug时耗时9ms，release耗时2ms。