Distance Between Points and a Line Segment in 3D

Hi
I have an mx3 array of point cloud data (X,Y,Z) and a 3D line. The line is defined by a start point (1x3) and an end point (1x3).
I'd like to find the orthogonal distance between each point and the line, but the line must only extend between the start and end point (i.e. not between +/- infinity).
Any solutions to this problem would be greatly appreciated.
Many thanks!
Jack

 Réponse acceptée

Rik
Rik le 12 Sep 2017
This problem can be tackled by separating the points in 3 parts. Let's define some distances for clarity: dsp is the distance between the start and the current point, dep is the distance between the end point and the current point, and dse is the distance between start and end.
  • type 1 is all points where sqrt(dse^2+dep^2)>=dsp
  • type 2 is all points where sqrt(dse^2+dsp^2)>=dep
  • type 3 is all other points
For type 1, the distance to the line segment is simply dep. For type 2, the distance to the line segment is simply dsp. (to see this, draw the line as one side of a right-sided triangle and apply Pythagoras)
Type 3 is between start and end, so there we need the orthogonal distance. We can slightly modify a staff-provided answer:
function d = point_to_line(pt, v1, v2)
% pt should be nx3
% v1 and v2 are vertices on the line (each 1x3)
% d is a nx1 vector with the orthogonal distances
v1 = repmat(v1,size(pt,1),1);
v2 = repmat(v2,size(pt,1),1);
a = v1 - v2;
b = pt - v2;
d = sqrt(sum(cross(a,b,2).^2,2)) ./ sqrt(sum(a.^2,2));
You can easily separate these three cases with logical indexing.

4 commentaires

Songqiu
Songqiu le 4 Déc 2019
Shouldn't be?
  • type 1 is all points where sqrt(dse^2+dep^2)<=dsp
  • type 2 is all points where sqrt(dse^2+dsp^2)<=dep
Rik
Rik le 4 Déc 2019
Modifié(e) : Rik le 4 Déc 2019
Yes, you are correct, I mixed up the direction of the inequality in my initial description. If we sketch the situation, it becomes more clear. If you consider the example below, the distance to the line segment is equal to dep if the point P is either at this location, or to the right. Currently, sqrt(dse^2+dep^2)=dsp. Moving P to the right will increase dsp (while dse stays the same and the new P can be chosen such that dep stays the same as well), without changing the fact that the output should still be dep. That means the inequalities should indeed be as you describe.
sketch.png
As far as the implementation: I don't think it is smart to merge the two concepts, because that would make the FEX submission less clear. Although they are related, I think they are better suited to be separate functions. If you are unable to implent it yourself, post a new question and post a link here.
Below you find the implementation for this situation. I barely tested this, so there may still be some issues. If you find more bugs, please let me know.
function distance=point_to_line_segment_distance(pt, v1, v2)
%assume nx3,1x3,1x3 inputs
distance=zeros(size(pt,1),1);
%type 1 are all points closer to the end-point than to the rest of the line
%type 2 is closer to the start-point
%type 3 is the rest
%for type 3 we can use the point_to_line_distance FEX submission
%https://www.mathworks.com/matlabcentral/fileexchange/64396-point-to-line-distance
%dsp is the distance between the start and points
%dep is the distance between the end and points
%dse is the distance between start and end
%dpl is the distance between P and the projection on the entire line
dsp=sqrt( (pt(:,1)-v1(1)).^2 + (pt(:,2)-v1(2)).^2 + (pt(:,3)-v1(3)).^2 );
dep=sqrt( (pt(:,1)-v2(1)).^2 + (pt(:,2)-v2(2)).^2 + (pt(:,3)-v2(3)).^2 );
dse=sqrt( sum( (v1-v2).^2 ) );
dpl=point_to_line_distance(pt, v1, v2);
type1= sqrt(dse^2+dep.^2)<=dsp ;
type2= sqrt(dse^2+dsp.^2)<=dep ;
type3= ~type1 & ~type2;
distance(type1)=dep(type1);
distance(type2)=dsp(type2);
distance(type3)=dpl(type3);
end
Great, thanks so much! But there seems to be a mistake as I noticed in my tests that the distance gets lower when the segment length gets bigger.
.. Gravitation effect? ;)
Rik
Rik le 20 Fév 2020
Do you have an example I could test later?

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by