欧式距离变换 EDT_chamfer distance 欧几里得距离

零、建议阅读

文献 & 博客：

1. 李均利,陈爽,王秀英.三维欧氏距离变换快速算法[J].计算机辅助设计与图形学学报,2012,24(12):1559-1567. （适合一开始连概念都不懂的，只看介绍部分就行，算法不用看

2. 欧几里得距离转换（EDT）算法_distance transforms of sampled functions-CSDN博客

3. 图像距离变换_根据图像中的距离计算实际距离的算法或模型-CSDN博客

4. Fabbri R, Costa L D F, Torelli J C, et al. 2D Euclidean distance transform algorithms: A comparative survey[J]. ACM Computing Surveys (CSUR), 2008, 40(1): 1-44. （特别推荐！

5. Meijster A, Roerdink J B T M, Hesselink W H. A general algorithm for computing distance transforms in linear time[J]. Mathematical Morphology and its applications to image and signal processing, 2000: 331-340.

6. Felzenszwalb P F, Huttenlocher D P. Distance transforms of sampled functions[J]. Theory of computing, 2012, 8(1): 415-428.

7. 欧几里得距离转换（EDT）算法记录 - 知乎 (zhihu.com) (看6的文献时可以结合着一起看

代码：

1. euclidean-distance-transform-3d/python/edt.hpp at a38a2a0a434076baba6b43c1844649618f07c790 · seung-lab/euclidean-distance-transform-3d (github.com)

2.distance-transform-cuda/DistanceTransformCUDA.cu at master · 1danielcoelho/distance-transform-cuda (github.com)

3.maplibre-native/src/mbgl/util/tiny_sdf.cpp at 098dca8c5aaefa460e44d0a6643ff2dba1fd74e4 · maplibre/maplibre-native (github.com)

---

以下是个人笔记，没有顺序与逻辑。可能有错误，谨慎观看。

---

一、距离度量（distance metrics）

1. 棋盘距离（D8）
2. 城市街区距离（D4，又叫曼哈顿距离）
3. 欧式距离
4. 倒角距离（chamfer）（不知道这个放在这里对不对

倒角距离使用预先定义的距离模板来近似欧式距离，常见的倒角距离模板有3x3模板和5x5模板（对于2D图像来说，如果是3D图像则是3x3x3和5x5x5）。

倒角方法，具体可查看文献：

1. Shih F Y, Wu Y T. Fast Euclidean distance transformation in two scans using a 3× 3 neighborhood[J]. Computer Vision and Image Understanding, 2004, 93(2): 195-205. （2D的情况）

2. Svensson S, Borgefors G. Digital distance transforms in 3D images using information from neighbourhoods up to 5× 5× 5[J]. Computer Vision and Image Understanding, 2002, 88(1): 24-53. （3D的情况）

二、距离变换（distance transform, DT）算法

有很多方法，也有不同的分类方法。

1. 可以分为：iterative 和 sequential (or recursive)，其中iterative可以并行，sequential不可以

2. 也可以大致分为：(这里的分类是not exclusive，即不是排它的)

1. 有序传播算法（ordered propagation algorithm）
2. 光栅扫描算法（raster scanning algorithm）
3. 独立扫描算法（independent algorithm），也可以叫降维算法（dimensional reduction）

chamfer DT就是属于光栅扫描算法，不可并行。

三、独立扫描算法

如果是2D图像，则包含两部分：1st transformation, 2nd transformation。每个部分负责一个维度。详细内容去看”建议阅读“中的2，3，4.

独立扫描算法也有很多不同的方法，但每种方法的1st tranformation都是相同的。根据2nd tansformation的不同可以分为三种：

1. based on parabola intersection (e.g. Satio's algorithm, Meijster's algorithm, Felzenszwalb‘s algorithm)
2. based on mathmatical morphology
3. based on Voronoi diagram intersections with image line (e.g. Maurer's algorithm)

---

四、尝试写的代码

1. chamfer DT

使用3x3x3的距离模板进行chamfer DT，语言matlab

function F = edt(I,a,b,c)
 
[nx,ny,nz] = size(I);
F = inf(nx,ny,nz);
F(I) = 0;
F = padarray(F,[1 1 1],NaN,'both');
 
% local distance matrix
D = zeros(3,3,3);
D([5,11,13,15,17,23]) = a;
D([2,4,6,8,10,12,16,18,20,22,24,26]) = b;
D([1,3,7,9,19,21,25,27]) = c;
 
% foreward raster scan
for z = 2:nz+1
    for y = 2:ny+1
        for x = 2:nx+1
            F(x,y,z) = min([F(x,y,z), ...
                D(1)+F(x-1,y-1,z-1),D(4)+F(x,y-1,z-1),  D(7)+F(x+1,y-1,z-1), ...
                D(2)+F(x-1,y,z-1),  D(5)+F(x,y,z-1),    D(8)+F(x+1,y,z-1), ...
                D(3)+F(x-1,y+1,z-1),D(6)+F(x,y+1,z-1),  D(9)+F(x+1,y+1,z-1) ...
                D(10)+F(x-1,y-1,z), D(13)+F(x,y-1,z),   D(16)+F(x+1,y-1,z) ...
                D(11)+F(x-1,y,z)]);
        end
    end
end
 
% backward raster scan
for z = nz+1:-1:2
    for y = ny+1:-1:2
        for x = nx+1:-1:2
            F(x,y,z) = min([F(x,y,z), ...
                D(17)+F(x+1,y,z), ...
                D(12)+F(x-1,y+1,z),     D(15)+F(x,y+1,z),   D(18)+F(x+1,y+1,z), ...
                D(19)+F(x-1,y-1,z+1),   D(22)+F(x,y-1,z+1), D(25)+F(x+1,y-1,z+1), ...
                D(20)+F(x-1,y,z+1),     D(23)+F(x,y,z+1),   D(26)+F(x+1,y,z+1), ...
                D(21)+F(x-1,y+1,z+1),   D(24)+F(x,y+1,z+1), D(27)+F(x+1,y+1,z+1)]);
        end
    end
end
 
F = F(2:end-1,2:end-1,2:end-1);
 
end

2. Saito方法

对2D binary图像I进行欧式距离变换，语言matlab

ps1: 参考 github 1danielcoelho/distance-transform-cuda （他说Meijster方法，但其实还是Saito方法

ps2: F是平方距离，I是2D binary image

function F = edt(I)
 
% nx(width), ny(height)
[height,width] = size(I);
% N (infinity)
N = width * height + 1;
 
% first phase
G = zeros(size(I));
for i = 1:width
    % Initialize val to either 0 or 'infinity'
    val = (1 - I(1,i)) * N;
    G(1,i) = val;
 
    % scan 1
    for j = 2:height
        val = (1 - I(j,i))*(1+val);
        G(j,i) = val;
    end
 
    % scan 2
    for j = height-1:-1:1
        if G(j,i) > val
            G(j,i) = 1+val;
        end
        val = G(j,i);
    end
end
 
% second phase
F = G;
x = 1:width;
for j = 1:height
    for i = 1:width
        F(j,i) = min((x-i).*(x-i)+G(j,:).*G(j,:));
    end
end
 
end

3. Felzenszwalb方法

（该方法据说是目前最好的)，对2D binary图像进行欧式距离变换，语言matlab

ps: 参考 github maplibre-native

function G = edt(I)
 
[height,width] = size(I);
G = zeros(size(I));
G(~I) = height*width+1;
 
% col
for x = 1:width
    f = G(:,x);
    d = edt1d(f,height);
    G(:,x) = d;
end
 
% row
for y = 1:height
    f = G(y,:);
    d = edt1d(f,width);
    G(y,:) = d;
end
 
end
 
function d = edt1d(f,n)
d = zeros(1,n);
v = zeros(1,n);
z = zeros(1,n+1);
 
v(1) = 1;
z(1) = -inf;
z(2) = inf;
 
k = 1;
for q = 2:n
    s = ((f(q) + q * q) - (f(v(k)) + v(k) * v(k))) / (2 * q - 2 * v(k));
    while (s <= z(k))
        k = k-1;
        s = ((f(q) + q * q) - (f(v(k)) + v(k) * v(k))) / (2 * q - 2 * v(k));
    end
    k = k+1;
    v(k) = q;
    z(k) = s;
    z(k+1) = inf;
end
 
k = 1;
for q = 1:n
    while (z(k+1)<q), k = k+1; end
    d(q) = (q - v(k))*(q - v(k)) + f(v(k));
end
end

---

2024.6.22