将多边形的一组点压缩成一组较短的点

Condense a set of points of a polygon into a shorter set of points

我有以下多边形,它只是一组二维点,如下所示:-

poly0=[80    60
    90    60
   100    60
   110    60
   110    50
   120    50
   130    50
   140    50
   150    50
   160    50
   170    50
   180    50
   190    50
   200    50
   210    50
   210    60
   210    70
   210    80
   210    90
   220    90
   220   100
   210   100
   210   110
   200   110
   200   120
   190   120
   180   120
   180   130
   170   130
   160   130
   150   130
   140   130
   130   130
   130   120
   120   120
   110   120
   110   110
   100   110
   100   100
    90   100
    90    90
    90    80
    90    70
    80    70
    80    60];

现在我可以使用它来绘制它了。

>> line(poly0(:,1), poly0(:,2),'Color','k','LineWidth',3,'LineStyle',':');

这清楚地表明一件事,即我原来的多边形点集是高度冗余的。基本上,上面列举了位于同一条直线上的多个点,这是不需要的。我可以开始检查每对点,如果它们在同一条直线上,我可以删除它们。但这将意味着使用许多 for 循环。我想不出一个聪明的矢量化方式。

如何获得一组新的点,这些点的尺寸比以前的要小得多,但仍然表示完全相同的多边形?我应该只拥有与多边形中的顶点一样多的点。那么换句话说如何快速地从上述数据集中找到顶点?

PS:这里的顶角是 90 度,但是如果你给出了一个解决方案,不要试图利用这个事实。我想要一个更笼统的答案。

'vector' 的方式可以做得很优雅。我也尝试了 for 循环方式,你可以用这种方式做同样的事情,但你要求矢量,所以这是我的方式。

我对您的数据所做的唯一更改是在启动此脚本之前删除所有重复项。此外,提供的点应按顺时针或逆时针顺序排列。

    poly0Z = circshift(poly0,1);
    poly0I = circshift(poly0,-1); 
    unitVectIn =(poly0 - poly0I)./vecnorm((poly0 - poly0I),2,2);
    unitVectOut =(poly0Z - poly0)./vecnorm((poly0Z - poly0),2,2)  ;
    cornerIndices = sum(unitVectIn == unitVectOut,2)==0
    corners = poly0(cornerIndices,:)

    line(poly0(:,1), poly0(:,2),'Color','k','LineWidth',2,'LineStyle',':');
    hold on
    scatter(corners(:,1), corners(:,2),'filled')

这个方法的基础是走到每个点,计算进来的单位向量,和出去的单位向量。输入的单位向量与输出的单位向量不匹配的点是角点。

好的,我对此进行了调整以处理非方角。

考虑一个由点识别的三角形

P = [0 0; 1 0; 2 0; 1.5 1; 1 2; .5 1; 0 0];

这是一个 7x2 数组,如果我们然后按照我在评论中提到的问题定义的那样定义 2 个导数向量。

a = logical([1 diff(P(:,1),2)' 1]);
b = logical([1 diff(P(:,2),2)' 1]);

从那里,我们可以将两者结合起来得到一个新的索引变量

idx = or(a,b);

最后,我们可以用它来制作我们的情节

line(P(idx,1), P(idx,2),'Color','k','LineWidth',3,'LineStyle',':');

如果您正在绘制线图,我认为您需要将最后一个变量设置为 false。

idx(end) = false;

现有的两个答案都有很大的不足:

  • 仅在后续点之间的距离完全相同时才有效。点必须具有可以完美表示为浮点值的坐标,以便可以发现一条线上的后续点之间的距离是相同的。如果点不等距,则该方法不执行任何操作。此外,多边形的起点和终点不会一起检查,因此如果在多边形的 start/end 上形成一条直线,将保留一个点太多。

  • 更好,因为距离不需要相同,并且正确处理了穿过多边形 start/end 的线。然而,它也使用浮点数相等比较,这在一般情况下是行不通的。只有使用整数坐标,此方法才能正常工作。此外,它使用 vecnorm(计算平方根)和除法,两者都是相对昂贵的操作(与此处显示的方法相比)。

要看三点是否构成一条直线,可以用一个简单的算术规则。假设我们有点 p0p1p2。从p0p1的向量和从p0p2的向量构成平行四边形的基础,其面积可以通过the cross product of the two vectors计算(在2D ,叉积被理解为使用 z=0,结果向量具有 x=0y=0,只有 z 值有用;因此,我们假设的 2D 叉积产生一个标量值)。可以这样计算:

v1 = p1 - p0;
v2 = p2 - p0;
x = v1(1)*v2(2) - v1(2)*v2(1);
如果两个向量平行,则叉积

x 将为零,这意味着三个点共线。但是等于 0 的测试必须有一定的公差,因为浮点运算是不精确的。我在这里使用 1e-6 作为公差。使用比点之间的距离小几个数量级的值。

给定一组输入点 p,我们可以找到角点:

p1 = p;                                  % point 1
p0 = circshift(p1,1);                    % point 0
v1 = p1 - p0;                            % vector from point 0 to 1
v2 = circshift(p1,-1) - p0;              % vector from point 0 to 2
x = v1(:,1).*v2(:,2) - v1(:,2).*v2(:,1); % cross product
idx = abs(x) > 1e-6;                     % comparison with tolerance
p = p(idx,:);                            % corner points

请注意,如果两个连续的点具有相同的坐标(即其中一个向量的长度为零),则此叉积测试将失败。如果数据可能有重复的点,则需要进行额外的测试。

这是三种方法的结果。我创建了一个具有非平凡坐标且顶点间距不等的多边形。我还将 start/end 间隙放在直边的中间。这些特点是为了展示其他两种方法的缺点。

这是我用来生成图表的代码:

% Make a polygon that will be difficult for the other two methods
p = [0,0 ; 0.5,0 ; 1,0 ; 1,1 ; 0.5,1 ; 0,1];
p = p + rand(size(p))/3;
p(end+1,:) = p(1,:);
q = [];
for ii = 1:size(p,1)-1
   t = p(ii,:) + (p(ii+1,:) - p(ii,:)) .* [0;0.1;1/3;0.45;0.5897545;pi/4;exp(1)/3];
   q = [q;t];
end
q = circshift(q,3,1);

figure
subplot(2,2,1)
plot(q(:,1),q(:,2),'bo-')
axis equal
title('input')

subplot(2,2,2)
res1 = method1(q);
plot(res1(:,1),res1(:,2),'ro-')
axis equal
title('Durkee''s method')

subplot(2,2,3)
res2 = method2(q);
plot(res2(:,1),res2(:,2),'ro-')
axis equal
title('ShadowMan''s method')

subplot(2,2,4)
res3 = method3(q);
plot(res3(:,1),res3(:,2),'go-')
axis equal
title('correct method')

% Durkee's method: 
function P = method1(P)
a = logical([1 diff(P(:,1),2)' 1]);
b = logical([1 diff(P(:,2),2)' 1]);
idx = or(a,b);
P = P(idx,:);
end

% ShadowMan's method: 
function corners = method2(poly0)
poly0Z = circshift(poly0,1);
poly0I = circshift(poly0,-1);
unitVectIn =(poly0 - poly0I)./vecnorm((poly0 - poly0I),2,2);
unitVectOut =(poly0Z - poly0)./vecnorm((poly0Z - poly0),2,2);
cornerIndices = sum(unitVectIn == unitVectOut,2)==0;
corners = poly0(cornerIndices,:);
end
% vecnorm is new to R2017b, I'm still running R2017a.
function p = vecnorm(p,n,d)
% n is always 2
p = sqrt(sum(p.^2,d));
end

function p = method3(p1)
p0 = circshift(p1,1);
v1 = p1 - p0;
v2 = circshift(p1,-1) - p0;
x = v1(:,1).*v2(:,2) - v1(:,2).*v2(:,1);
idx = abs(x) > 1e-6;
p = p1(idx,:);
end