距离变换[Rosenfeld and Pfaltz,1966,1968]

坐标为P(px,py)和Q(qx,qy)的两点间距离可以定义为几种形式，主要有三种形式：欧式距离(Euclidean distance)、城市距离(city block distance)D4和棋盘距离(chessboard distance)D8。

    欧式距离被定义为：

    城市距离D4被定义为：

    棋盘距离D8被定义为：

对于距离度量(D4/D8/DE)，Rosenfeld and Pfaltz提出了一种计算距离变换的算法，该算法用一个小的掩模先从左上遍历图像后用另一个掩模反向遍历。第一次使用的掩模AL和第二次使用的掩模BR如下图所示。

按照指定的距离度量，对M*N的图像I中的一个任意子集S(I中的部分像素组成的区域)计算距离变换，同时建立一个M*N大小的二维数组F，并置F在S有像素的位置为0，其他为无穷(Inf)。

第一次从左至右，从上到下用掩模AL遍历I，对于任意一点P，依次计算F(P)的值 （Q为P的邻接像素）：

第二次从右至左，从下到上用掩模BR遍历I，对于任意一点P，依次计算F(P)的值 （Q为P的邻接像素）

最后得到的F即为S的斜切。

主函数代码如下：

I=imread('lena.jpg');
if ndims(I)==3%如果为RGB图像则转为灰度
    I=rgb2gray(I);
end
[m,n]=size(I);
S=I(100:200,100:200);
F=zeros(m,n);
D_para=0;
for i=1:m
    for j=1:n
        if i<100 || i>200 || j<100 || j>200
            F(i,j)=Inf;
        end
    end
end
imshow(F);
figure,
for i=1:m%第一次遍历
    for j=1:n
        for k=1:4%四个AL邻接点
            [qx,qy]=returnpoint(i,j,k);
            if ifavailable(qx,qy,m,n)==1
                F(i,j)=min(F(i,j),D(i,j,qx,qy,D_para)+F(qx,qy));
            end
        end
    end
end
for i=m:-1:1%第二次遍历
    for j=n:-1:1
        for k=5:8
            [qx,qy]=returnpoint(i,j,k);
            if ifavailable(qx,qy,m,n)==1
                F(i,j)=min(F(i,j),D(i,j,qx,qy,D_para)+F(qx,qy));
            end
        end
    end
end
F=uint8(F);
imshow(F);

在主代码中，设子集为图像I中第100~200行，第100~200列的矩形区域，returnpoint(i,j,k)为自定义函数，返回点(i,j)的第k个邻接点，k取值为1~8，当k分别为1~8时，该函数分别返回左上、上、左、左下、右上、右、下、右下的坐标。D_para为选择距离度量的类型，0为欧式距离，1为D4，2为D8。D(i,j,qx,qy,D_para)为自定义函数，根据D_para选择的类型计算(i,j)和(qx,qy)返回两点的距离。ifavailable为自定义函数，判断(qx,qy)点是否超出图像I的区域。自定义函数分别如下：

function [qx,qy]=returnpoint(px,py,k) 
%根据p和k返回q
%k = 1/2/3/4 左上/上/左/左下
%k = 5/6/7/8 右上/右/下/右下
    switch k
        case 1 
            qx=px-1;
            qy=py-1;
        case 2 
            qx=px-1;
            qy=py;
        case 3 
            qx=px;
            qy=py-1;
        case 4 
            qx=px+1;
            qy=py-1;    
        case 5 
            qx=px-1;
            qy=py+1;
        case 6 
            qx=px;
            qy=py+1;
        case 7 
            qx=px+1;
            qy=py;
        case 8 
            qx=px+1;
            qy=py+1;  
    end
end

returnpoint函数直接采用switch-case-end，该关键字和C语言用法有较大出入。详情help switch。

function dis=D(px,py,qx,qy,D_para)
    if D_para==1
        dis=abs(px-qx)+abs(py-qy);
    elseif D_para==2
        dis=max(abs(px-qx),abs(py-qy));
    elseif D_para==0
        dis=sqrt((px-qx)^2+(py-qy)^2);
    end
end

function temp=ifavailable(x,y,m,n)
%这个坐标是否在图内
    if x<1 || x>m || y<1 || y>n
        temp=0;
    else
        temp=1;
    end
end

S(白色区域)：

运行结果：(越黑代表距离越近)

DE:

D4:

D8: