【问题标题】:Condense a set of points of a polygon into a shorter set of points将多边形的一组点压缩为一组较短的点
【发布时间】:2019-08-31 07:21:06
【问题描述】:

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

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 度,但如果你给出了一个解决方案,不要试图利用这个事实。我想要一个更笼统的答案。

【问题讨论】:

  • 这只是一个循环:对于每3个连续点,如果它们在一条直线上,则删除中间的一个。如果您删除一个点,请再次检查相同的第一个点,添加一个新的第三个点。在循环之后,您还需要检查点end-1end1,以及end12
  • @Durkee Ha,工作!但我认为它利用了多边形角度为 90 度的事实。请参阅上面的 PS。
  • @CrisLuengo 我猜你的答案是最通用的答案。如果用一些代码做成答案,我会接受。
  • 积分是否强制连续?

标签: matlab polygon


【解决方案1】:

“矢量”方式可以非常优雅地完成。我也尝试了 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')

这个方法的基础是去每个点,计算进来的单位向量,出去的单位向量。单位向量 in 与单位向量 out 不匹配的点就是角点。

【讨论】:

  • 这个方法效果很好,除了相等比较。在比较归一化向量(abs(unitVectIn-unitVectOut)<1e-6 或类似的)时,您应该允许容差。
【解决方案2】:

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

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

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

这是一个 7x2 数组,如果我们随后定义我在 cmets 中提到的问题所定义的 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;

【讨论】:

  • 这种方法比我的计算速度要快得多,而且编写代码的时间也快:D。
  • 如果第一个点和最后一个点形成一条直线,则您没有删除足够的点。我认为,另一个答案确实可以通过 circshift 做到这一点。
  • 是的,视情况而定,这可能无关紧要。或者,可以使用 if 语句来修复它。
  • 我刚刚意识到这个方法还有一个缺点。您假设顶点是等距的。如果直线上的点的距离不相等,它们将保留在输出中。
【解决方案3】:

现有的两个答案有很大的缺点:

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

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

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

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

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

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

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

% 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: https://stackoverflow.com/a/55603145/7328782
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: https://stackoverflow.com/a/55603040/7328782
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

【讨论】:

  • 克里斯,感谢您的反馈。我试图使用你的紧张点在我的脚本中实现容忍。 (abs(unitVectIn-unitVectOut)<1e-6 它并没有完全做到这一点,因为它给我留下了一个向量。 sum((abs(unitVectIn-unitVectOut)<1e-6) == 2 是否合适,因为现在我不是比较浮点数,而是比较第一个布尔值结果的整数?如果您知道任何好的资源,我想更多地了解您的批评背后的原则。
  • @ShadowMan:是的,我就是这个意思。它是关于用abs(a-b)<tol 替换a==bThis is the work that everyone always references when discussing issues about floating-point representation.
  • 您有什么理由选择手动计算叉积而不是使用cross()
  • @ShadowMan: cross 需要 3D 矢量。我必须添加一列零(对于第三维),然后它会计算我已经知道为零的 x 和 y 值。所以效率会更低。我不知道 MATLAB 是否有计算二维叉积的函数。此外,我从已有的 C++ 代码中修改了代码,我确实努力使用我在您的答案中看到的 circshift 替换循环,但没有看是否可以简化其他任何事情。所以留下这个手动交叉产品是一个简单的解决方案。 :)
  • 感谢这里的提示,即使在“抖动”数据上更正我的脚本中的浮点比较也会产生良好的结果。从这里开始,我将始终使用我的代码来寻找这个问题:D
猜你喜欢
  • 2010-10-24
  • 2021-06-28
  • 1970-01-01
  • 1970-01-01
  • 2013-07-14
  • 1970-01-01
  • 1970-01-01
  • 2016-10-28
  • 1970-01-01
相关资源
最近更新 更多