现有的两个答案有很大的缺点:
Durkee's method 仅在后续点之间的距离完全相同时才有效。点必须具有可以完美表示为浮点值的坐标,以便可以发现一条线上后续点之间的距离是相同的。如果点不等距,则该方法不执行任何操作。此外,多边形的起点和终点不会一起检查,因此如果在多边形的起点/终点形成一条直线,则会留下太多的点。
ShadowMan's method 更好,因为距离不需要相同,并且正确处理了穿过多边形起点/终点的线。但是,它也使用浮点相等比较,这通常不起作用。只有使用整数坐标,此方法才能正常工作。此外,它使用vecnorm(计算平方根)和除法,两者都是相对昂贵的操作(与此处显示的方法相比)。
要查看三个点是否形成一条直线,可以使用一个简单的算术规则。假设我们有积分p0、p1 和p2。从p0 到p1 的向量和从p0 到p2 的向量构成平行四边形的基础,其面积可以通过the cross product of the two vectors 计算(在二维中,叉积被理解为使用@987654335 @,结果向量具有x=0 和y=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