How to get the max orthogonal (perpendicular) distances from a set of coordinates

9 次查看(过去 30 天)
I have arrays of X and Y values that I used to plot on a figure and to find the max distance between the points (see picture). What I want to do now is instead of finding the max distance, I want to find 2 distances that are orthogonal (perpandicular) to each other that is the max possible value for both. Thanks.
  11 个评论
DGM
DGM 2021-6-3
So you're basically looking for the longest inscribed perpendicular to the maximal distance vector. That sounds solvable, but I know my own solutions to geometry problems like this tend to be ... inelegant.
Talha Majeed
Talha Majeed 2021-6-3
Yes that's what I'm looking for I'm completely stuck, any solution you can provide however inelegant would be much appreciated, thanks

请先登录,再进行评论。

采纳的回答

DGM
DGM 2021-6-3
编辑:DGM 2021-6-3
I'm sure this can be done much more efficiently, and I hope someone is motivated to prove that point. I'm cheating by using polyxpoly(), but I hate doing line intersections and dealing with a bunch of conditions. I sure hope you have polyxpoly().
% generate random point ring
a = 0.2; % r variance
b = 0.1; % th variance
r = 100; % maximal radius
th = 0:(2*pi/20):(2*pi*(1-1/20)) + b*rand(1,20);
x = r*(1-a + a*rand(1,20)).*cos(th);
y = r*(1-a + a*rand(1,20)).*sin(th);
% need to close the polygon
x = [x x(1)];
y = [y y(1)];
D = sqrt((x-x.').^2 + (y-y.').^2); % distance to/from all points
D(D<1E-6) = NaN; % remove self-distances
[Dn Nn] = max(D,[],2); % maximize
% Dn is distance to furthest neighbor
% Nn is index of furthest neighbor
% find maximal distance vector
[Dmax,idx] = max(Dn);
mxidx = [idx Nn(idx)];
DVMx = x(mxidx)
DVMy = y(mxidx)
plot(x,y,':o'); hold on;
plot(DVMx,DVMy,'b')
axis equal
% angle of DVM
ang = atan2d(diff(DVMy),diff(DVMx))
% each candidate perpendicular should have at least one endpoint on a vertex
% so consider all vertices V except those belonging to DVM
% find distance from V to opposite edge of polygon
perplen = zeros(numel(x),1);
pix = zeros(numel(x),1);
piy = zeros(numel(x),1);
for v = 1:numel(x)
% skip if V belongs to DVM
if ismember(v,mxidx); continue; end
% project a perpendicular from V
px = x(v) + [1 -1]*Dmax*cosd(ang+90);
py = y(v) + [1 -1]*Dmax*sind(ang+90);
% find intersections between perpendicular and polygon
[xi,yi] = polyxpoly(x,y,px,py);
% find distance to points of intersection from V, maximize
pid = [xi.'-x(v); yi.'-y(v)];
[perplen(v) mpididx] = max(sqrt(pid(1,:).^2 + pid(2,:).^2));
% this is the other end of the candidate perpendicular
pix(v) = xi(mpididx);
piy(v) = yi(mpididx);
plot([x(v) pix(v)],[y(v) piy(v)],'m')
end
% find maximal perpendicular vector
[~,idx] = max(perplen);
DVPx = [x(idx) pix(idx)];
DVPy = [y(idx) piy(idx)];
plot(DVPx,DVPy,'k--','linewidth',2)
This naive method is not robust against cases where the polygon has deep concave features or is self-intersecting. I don't know if point ordering even matters for your usage.
Actually, that's a good point. This code presumes x and y are ordered. If they're not, the polygon drawn between the points will likely be self-intersecting nonsense.
  7 个评论
DGM
DGM 2021-6-4
What are you using instead of polyxpoly()? I imagine that might be part of it. Your example set of points works for me.
Well, it works as well as suggested (note the intersection with the concave feature).
The fact that points are projected at different angles is curious. I imagine that the projected lines are still perpendicular:
plot(px,py,'b:') % add after px,py are calculated in the loop
but there's probably some problem with how the results from the intersection calculation are used to calculate pix, piy, or how those points are associated with the vertex list (x and y). It's like some of the okay lines are behaving as intended with one endpoint on a vertex and the opposite endpoint at whatever intersection is required, while the strange lines all seem to have a vertex at both ends -- and they seem to like sharing one vertex. I'd have to see the changes to know why it's different.

请先登录,再进行评论。

更多回答(1 个)

Matt J
Matt J 2021-6-4
编辑:Matt J 2021-6-4
I created an alternative to polyxpoly here,
which finds the intersection of a single line and a polyshape boundary.
Once you have this, it's not too difficult to write code that will find the perpendicular lines even for polygons that have very deep concavities, though I don't know that that scenario is actually relevant to you. I illustrate this below with the same golf-club shaped region as in my earlier comment.
load('polygon','pgon')
[xl,yl] = boundingbox(pgon); tol=min(diff(xl), diff(yl))/1e6;
[xl,yl] = ndgrid(xl,yl);
warning off
bbox=convhull(polyshape([xl(:),yl(:)]));
warning on
dpgon=subtract(bbox,pgon);
V=pgon.Vertices; nv=size(V,1);
IJ=nchoosek(1:nv,2);
I=IJ(:,1); J=IJ(:,2);
XY1=V(I(:),:).'; XY2=V(J(:),:).';
N=numel(I);
Dists=nan(nv,1);
for i=1:N
xy1=XY1(:,i); xy2=XY2(:,i);
out=linexlines2D(dpgon,xy1,xy2);
out=uniquetol(out.',tol,'ByRows',1).';
if size(out,2)>2
continue;
else
Dists(i)=norm(xy1-xy2);
end
end
[L,imax]=max(Dists);
xy1=XY1(:,imax); xy2=XY2(:,imax); %maximally separated vertices
dvec=[0 -1; 1 0]*(xy2-xy1); %perpendicular direction vector
dvec(end+1)=0;
dvec=dvec/norm(dvec);
T=linspace(0,1,100);
fun=@(t)perpWidth(xy1+t*(xy2-xy1), pgon,dvec);
[~,imax]=min( arrayfun(fun , T ));
tmin=fminsearch(fun, T(imax));
[fval, out]=fun(tmin);
xy3=out(:,1); xy4=out(:,2); %maximally separated perpendicular points.
%% Plot the results
plot(pgon);
hold on
plot([xy1(1),xy2(1)] , [xy1(2),xy2(2)],'-om','MarkerFaceColor','m')
plot([xy3(1),xy4(1)] , [xy3(2),xy4(2)],'-om','MarkerFaceColor','m')
axis equal
hold off
function [fval, out]=perpWidth(xy,pgon,dvec)
E=cross([xy;1],dvec);
out=linexlines2D(pgon,E);
if size(out,2)<2, fval=0; out=[xy,xy]; return; end
[~,idx]=mink( vecnorm(out-xy,2,1) , 2);
out=out(:,idx);
fval=-vecnorm(out(:,1)-out(:,2),2,1);
end

类别

Help CenterFile Exchange 中查找有关 Annotations 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by