（原）欧式距离变换

欧氏距离变换以后再添加。

下面分别给出计算欧式距离（EDT）的matlab和C代码。

参考网址

https://lost-contact.mit.edu/afs/cs.stanford.edu/package/matlab-r2009b/matlab/r2009b/toolbox/images/images/private/mexsrc/bwdistComputeEDT/

不清楚上面的网址是什么情况，但是有很多的matlab的mex文件的源码。

matlab：

% 下面的为欧式距离变换的代码
function Dr = mybwdist(binaryimg)  
[height, width]=size(binaryimg); 
D=[];
for  k = 1: width   % 由于matlab按列存储，所以内层遍历行，外层遍历列。
    D=[D calcFirstPassD1(binaryimg(:,k), height)];
end 
 
Dr=processDimN([height, width], D); 

end

function D=calcFirstPassD1(I, vectorLength)
D1=zeros(vectorLength,1);

%     // Initialize D1 - Feature points get set to zero, -1 otherwise
%     D1(i) = (I(col*nrows+i) == 1) ? 0.0f : -1.0f;
D1(I==0)= (-1);

%     // Process column 
D=voronoiEDT(D1);  
 
end

function D=voronoiEDT(D)
 
%     // Note: g and h are working vectors allocated in calling function
[ns, g, h] = constructPartialVoronoi(D);
if ns == -1
    return;
end

D=queryPartialVoronoi(g, h, ns, D);

end 

function [el, g, h]=constructPartialVoronoi(D)
%     //
%     // Construct partial voronoi diagram (see Maurer et al., 2003, Fig. 5, lines 1-14)
%     // Note: short variable names are used to mimic the notation of the paper
%     //
g=zeros(1,length(D));
h=zeros(1,length(D));
el = -1;
for k = 0:length(D)-1   % 在matlab中对应数组的高度
    fk = D(k+1);
    if fk ~= -1
        if (el < 1)
            el=el+1;
            g(el+1) = fk;
            h(el+1) = k;
        else 
            while ( (el >= 1) && RemoveEDT(g(el), g(el+1), fk, h(el), h(el+1), k))
                el=el-1;
            end  
            el=el+1;
            g(el+1) = fk;
            h(el+1) = k;
        end
    end
end

end

function res = RemoveEDT(du, dv, dw, u, v, w)
a = v - u;
b = w - v;
c = w - u;

%     // See Eq. 2 from Maurer et al., 2003
res=(c * dv - b * du - a * dw) > (a * b * c);
end

function D=queryPartialVoronoi(g, h, ns, D)

%     //
%     // Query partial Voronoi diagram (see Maurer et al., 2003, Fig. 5, lines 18-24)
%     //
el = 0;
for k = 0:length(D)-1
    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)))) ) 
        el=el+1;
    end
    D(k+1) = g(el+1) + (h(el+1) - k)*(h(el+1) - k);
end
end

function D3=processDimN(dims, D)
%     // Create temporary vectors for processing along dimension N
D3=[];
for k = 0 : dims(1)-1         % // 若为二维数组，则得到第一维元素总数，即行数。注意，matlab按列存储，C按行存储。
% // Generate indices to points along vector path
    D3=[D3;voronoiEDT(D(k+1,:))];
end
end

分别使用matlab自带的bwdist、上面的代码进行测试，测试结果一样。

C++代码（进行了简化，将不需要的代码合并了，同时，有些代码并非完全对应，但是结果正确。注意，matlab调用时需要取反，下面的C代码不再需要取反，如107行）：

// Removal Function - FVs that are the Voronoi sites of Voronoi cells that do not intersect Rd are removed
inline bool RemoveEDT(const double du, const double dv, const double dw, const double u, const double v, const double w)
{
    double a = v - u;
    double b = w - v;
    double c = w - u;

    // See Eq. 2 from Maurer et al., 2003
    return ((c * dv - b * du - a * dw) > (a * b * c));
}

inline int constructPartialVoronoi(double* D, double* g, int* h, int length)
{
    //
    // Construct partial voronoi diagram (see Maurer et al., 2003, Fig. 5, lines 1-14)
    // Note: short variable names are used to mimic the notation of the paper
    //
    int el = -1;

    for (int k = 0; k < length; k++)
    {
        double fk = D[k];
        if (fk != -1.0f)
        {
            if (el < 1)
            {
                el++;
                g[el] = fk;
                h[el] = k;
            }
            else
            {
                while ((el >= 1) && RemoveEDT(g[el - 1], g[el], fk, h[el - 1], h[el], k))
                {
                    el--;
                }
                el++;
                g[el] = fk;
                h[el] = k;
            }
        }
    }

    return(el);
}

inline void queryPartialVoronoi(const double* g, const int* h, const int ns, double* D, int length)
{
    //
    // Query partial Voronoi diagram (see Maurer et al., 2003, Fig. 5, lines 18-24)
    //
    int el = 0;

    for (int k = 0; k < length; k++)
    {
        while ((el < ns) && ((g[el] + ((h[el] - k)*(h[el] - k))) >(g[el + 1] + ((h[el + 1] - k)*(h[el + 1] - k))))) 
        {
            el++;
        }
        D[k] = g[el] + (h[el] - k)*(h[el] - k);
    }
}

void voronoiEDT(double* D, double* g, int* h, int length)
{
    // Note: g and h are working vectors allocated in calling function
    int ns = constructPartialVoronoi(D, g, h, length);
    if (ns == -1)
        return;

    queryPartialVoronoi(g, h, ns, D, length);

    return;
}

inline void processDimN(int width, int height/*const mwSize *dims, const mwSize ndims*/, double *D)
{
    // Create temporary vectors for processing along dimension N

    double* g = new double[width];
    int* h = new int[width];

    // 若为二维数组，则得到第一维元素总数，即行数。注意，matlab按列存储，C按行存储。
    for (int k = 0; k < height; k++)
    {
        // Process vector
        voronoiEDT(&D[k*width], g, h, width);
    }

    delete[] g;
    delete[] h;
}

// mxI为输入，unsigned char*类型，mxD为输出，double* 类型
// 注意，matlab中mxI是逻辑类型，需要1-mxI才是真正的输入；此处省略了这一步。
// 输入为二值化图像，大于阈值的不为0，小于阈值的为0。
void computeEDT(double *mxD, const unsigned char *mxI, int width, int height)
{
    double* D1 = new double[width];
    double* g = new double[width];
    int* h = new int[width];

    for (int k = 0; k < width; k++)
    {
        for (int i = 0; i < height; i++)
        {
            D1[i] = (mxI[i*width + k] == 0) ? 0.0f : -1.0f;
        }

        voronoiEDT(D1, g, h, height);

        for (int i = 0; i < height; i++)
        {
            mxD[i*width + k] = D1[i];
        }
    }

    delete[] D1;
    delete[] g;
    delete[] h;

    processDimN(width, height, mxD);
}

matlab代码对320*240的图像处理，耗时200ms左右（记不清了）；当然，matlab使用的自带的代码耗时2ms左右。C中debug时耗时9ms，release耗时2ms。