轮廓跟踪及检测

效果图如上所示，主要完成轮廓的跟踪及内外轮廓检测，轮廓跟踪代码如下所示：

function pt =BoundaryTrack(bw, pstart)
% 单个边缘跟踪，pstart为联通区右下方的像素点 当边缘位于边界处，无法跟踪，需将
% 图像人为扩大2*2像素并置0，并将最后结果做边界处理
[row, col]= size(bw);
k=0;
flag = 1;
ptcurrent = pstart;
off=[1 0; 1 1; 0 1; -1 1; -1 0; -1 -1; 0 -1; 1 -1];

pt=[pstart.x pstart.y];
while flag
flag = 0;   
for i =1:8
    k=bitand(k,7);
    x = ptcurrent.x + off(k+1,1);
    y = ptcurrent.y + off(k+1,2);
    if x>=1 && x <=col && y>=1 && y <=row
        if bw(y,x)
            for j = 1:2:8
                ptemp.x = x + off(j,1);
                ptemp.y = y + off(j,2);
                if ptemp.x >=1 && ptemp.x <=col && ptemp.y >=1 && ptemp.y <=row
                    if bw(ptemp.y, ptemp.x) == 0
                        flag = 1;
                        ptcurrent.x = x;
                        ptcurrent.y = y;
                        pt=[pt;x y];
                        break;
                    end
                end
            end
        end
    end
    if flag
        if pstart.x == ptcurrent.x && pstart.y == ptcurrent.y
            flag = 0;
        end
        break;
    end
    k = k +1;
end
k = k + 6;
end

为了完成多个区域的轮廓检测及跟踪，利用形态学算法，对图像取反，便能获取内轮廓代码如下：

function boundaries = MorBoundaries(bwt)
[row,col] = size(bwt);
%% calculate the outside bandaries
outbw = bwt;
outbw=[zeros(1,col);outbw;zeros(1,col)];
outbw=[zeros(row+2,1) outbw zeros(row+2,1)];

[L,num1] = bwlabel(outbw,8);
for i = 1: num1
    [pty, ptx] = find(L == i);
    npt=[ptx pty];
    dex = find(npt(:,2) == max(pty));
    pstart.x = max(ptx(dex));
    pstart.y = max(pty);
    pt = BoundaryTrack(outbw, pstart) - 1;
    xx = pt(:,1); yy = pt(:,2);
    xx(xx > col) = col;
    yy(yy > row) = row;

    bout{i} = [xx yy];
end

%% calculate the inside boundaries
inbw = ~bwt;
inbw=[zeros(1,col);inbw;zeros(1,col)];
inbw=[zeros(row+2,1) inbw zeros(row+2,1)];

[L,num2] = bwlabel(inbw,8);
for i = 1: num2
    [pty, ptx] = find(L == i);
    npt=[ptx pty];
    dex = find(npt(:,2) == max(pty));
    pstart.x = max(ptx(dex));
    pstart.y = max(pty);
    pt = BoundaryTrack(inbw, pstart) - 1;
    xx = pt(:,1); yy = pt(:,2);
    xx(xx > col) = col;
    yy(yy > row) = row;

    bin{i} = [xx yy];
end

%% 分解为内外轮廓
for i = 1: num1
    boundaries(i).out = bout{i};
    boundaries(i).in={};
    for j = 1:length(bin)
        pt = bin{j};
        if IsPointInPolygon(bout{i},pt(1,:))
            boundaries(i).in=[boundaries(i).in, pt];
        end
    end
end

其中，用到了判断点时候在多边形内部算法，这个实现由很多种算法，包括面积和、角度和、射线法，此处，利用射线法判断点时候在多边形内部，过点p做水平线，判断单边交点个数，个数为奇数，则在多边形内部，为偶数，则在多边形外部，代码如下：

function isIn = IsPointInPolygon(pt, point) 
% 求解通过该点的水平线与多边形个边的交点 单边交点个数为奇数，成立
nCross = 0;
for i = 1 : size(pt,1) - 1
p1 = pt(i,:);
p2 = pt(i+1,:);
if p1(2) == p2(2)
    continue;
end
if point(2) < min(p1(2),p2(2))
    continue;
end
if point(2) >= max(p1(2),p2(2))
    continue;
end
x = (point(2)- p1(2))*(p2(1)-p1(1))/(p2(2)-p1(2))+p1(1);
if x > point(1)
    nCross = nCross + 1;
end
end
isIn = mod(nCross,2) == 1;