matlab - 减少多边形的边数

Question

假设我们有一个简单图形的图像，并且我们知道它是一个多边形，略微扭曲。有没有一种图像处理方式来逼近图形对象的原始参数？

下面的矩阵由代码创建，然后缩小以显示第五个感兴趣的区域：

EnumeratedMask=bwlabel(Zdata<-.06);

0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   5   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   5   5   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   5   5   5   5   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   5   5   5   5   5   0   0
0   0   0   5   5   5   5   5   5   5   5   5   5   5   5   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   5   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0
0   0   0   0   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   5   5   0   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   5   5   5   5   5   0   0   0   0   0   0   0   0   0   0   0
0   0   0   5   5   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0

对于下一步，我需要ABC/ABCD坐标来跨由这些点进一步定义的线获取 z 轮廓。

Answer 1

在 Patrick 发布他的教育答案之前，我已经开始工作，并且我遇到了 Ramer-Douglas-Peucker 算法的一个“问题”：根据定义，它保留了第一点和最后一点。convhull和函数都boundary从某个地方开始，并不总是在角落里，这会触发误报。第三和第四步解决了这个问题——另一点更有可能是真正的角落。

算法：

1. 检测凸包/外点。
2. 应用更紧密拟合的 Ramer-Douglas-Peucker 算法 (RDP)
3. 倒回循环
4. 应用 RDP 算法去除拟合较松的假阳性角

代码：

eps1=2;
eps2=4;
img=[...
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 5 5 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 0 0
0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 5 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0
0 0 0 5 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

[X,Y]=find(matr==5); % // find image coordinates of points of interest
is=boundary(X,Y,1);  % // find points on the boundary

xy=[X(is),Y(is)];    % // boundary points (step 1)

% // apply RDP algorithm (step 2)
xy=RDP(xy,eps1);

% // rewind the loop (step 3)
xy=xy([2:end-1 1:2],:);

% // apply the DRP algorithm (step 4)
xy=RDP(xy,eps2);

function[hull]=RDP(hull,eps)
    sp=hull(1,:);
    ep=hull(end,:);
    ip=hull(2:end-1,:);
    % // calculate distances of inner points from the first-last line
    dst=PerpDist(sp,ep,ip);
    % // find the point furthest from the f-l line
    [mx,mi]=max(dst);
    if mx>eps % // furthest point does not fit in - split and apply RDP recursively
        lp=[sp;ip(1:mi,:)];
        if size(lp,1)>2 % // there are points left to assess
            lp=RDP(lp,eps);
        end
        rp=[ip(mi:end,:);ep];
        if size(rp,1)>2 % // there are points left to assess
            rp=RDP(rp,eps);
        end
        hull=[lp;rp(2:end,:)]; % // concatenate the branches
    else % // the furthest poit fits in the limit, drop all inner points
        hull=[sp;ep];
    end
end

function[D]=PerpDist(A,B,C)
    D=nan(size(C,1),1);
    if A==B % // edge is defined by one point, use euclidean distance between points
        for PDi=1:size(C,1)
            D(PDi)=norm(C(PDi,:)-A);
        end
    else % // edge is a line, use eucleidean distance from a line
        for PDi=1:size(C,1)
            D(PDi)=abs(A(1)*(B(2)-C(PDi,2))+B(1)*(C(PDi,2)-A(2))+C(PDi,1)*(A(2)-B(2)))/norm(B-A);
        end
    end
end

点：形状内的点。
红色方块：boundary函数返回的第一个和最后一个点。
绿线：第一次 RDP 简化的结果。
洋红色线：最终三角形（原始形状为三角形）。

编辑：

1. 留下找到两个最远点来分割循环的想法，因为它甚至对于三角形也会触发误报。从任意点开始并使用 RDP 两次就可以了。

学分：

Steve Eddins 为max_feret 直径
Cris Luengo 为Ramer-Douglas-Peucker 算法评论

Answer 2

这是上面 Cris Luengo 评论中建议的 Ramer-Douglas-Peucker 算法的实现。

这是答案的第一个版本的完整编辑，用于edge查找对象的边界。正如 Cris Luengo 在评论中指出的那样，bwboundaries它是二值图像的更好选择。bwboundaries返回排序后的点，这大大简化了代码。

以下代码执行以下操作：

1) 使用 bwboundaries 找到对象的边缘。它们已经排序。

2) 使用 Ramer-Douglas-Pecker 算法简化点列表。

由于我需要一些视觉提示来进行调试，因此代码会打开一个显示正在发生的事情的图形。

请注意，代码远未得到正确测试。

img = [...
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 5 5 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 5 5 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 5 5 5 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
    0 0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 5 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 5 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

watch = true;

if watch
    f = figure;
    subplot(1,2,1);
    imagesc(img);
end


% first, find the edges

b = bwboundaries(img);

% note that the first and the last point are the same.
% However they are already sorted.
x = b{1}(:,1);
y = b{1}(:,2);

edges = zeros(size(img));
edges(sub2ind(size(e), x,y)) = 1;

if watch
    ax = subplot(1,2,2);

    img_h = imagesc(edges);
    hold on;

end 

title('Performing Douglas-Peucker algorithm');

% Omit the last point for the algorithm.

points = l_DouglasPeucker( [x(1:end-1), y(1:end-1)], 1, watch);

title('Final result');
plot([points(:,2); points(1,2)], [points(:,1); points(1,1)]);



function res = l_DouglasPeucker( points, ep, watch )

    if nargin < 3
        watch = false;
    end

    if watch
        subplot(1,2,2);
        hold on;
        hp = plot(points(:,2), points(:,1), 'ko-');
        hp2 = plot([points(1,2) points(end,2)], [points(1,1) points(end,1)], 'r-');
    end
    distances = zeros(size(points,1),1);
    for i = 1:size(points,1)
        distances(i) = l_distance_to_line(points(1,:), points(end,:), points(i,:));
    end

    idx = find(distances == max(distances),1);


    if watch
        hp4 = plot(points(idx,2), points(idx,1), 'mo', 'MarkerFaceColor', [1,1,1]);
        pause(.5);
        delete(hp);
        delete(hp2);
        delete(hp4);
    end

    if max(distances) > ep
        res = [l_DouglasPeucker(points(1:idx,:), ep, watch); l_DouglasPeucker(points(idx:end,:), ep, watch)];
    else
        res = [points(1,:); points(end,:)];
    end


end

function d = l_distance_to_line(p1,p2,p)
% computes the distance of p to the line through by p1 and p2
% There might be much better implementations of this...

% we need 3-dimensional data for this

p1 = [p1(1), p1(2), 0];
p2 = [p2(1), p2(2), 0];
p = [p(1,1) p(1,2) 0];

a = p2 - p1;
b = p - p1;

d = norm(cross(a,b)) ./ norm(a);
end