version 2.0.0.0 (6 KB) by
Michael Yoshpe

Find minimum distances from points to a polyline or to a closed polygon

**Editor's Note:** This file was selected as MATLAB Central Pick of the Week

The following files are included:

p_poly_dist.m - Compute the distances from each one of a set of np points p(1), p(2),... p(np) on a 2D plane to a polyline or a closed polygon. Polyline is defined as a set of nv-1 segments connecting nv ordered vertices v(1), v(2), ..., v(nv). The polyline can optionally be treated as a closed polygon.

Distance from point j to a segment k is defined as a distance from this point to a straight line passing through vertices v(k) and v(k+1), when the projection of point j on this line falls INSIDE of segment k; and to the closest of v(k) or v(k+1) vertices, when the projection falls OUTSIDE of segment k. Distance from point j to a polyline is defined as a minimum of this point's distances to all segments.

In a case when the projected point fall OUTSIDE of ALL polyline segments, the distance to a closest vertex of the polyline is returned

test_p_poly_dist.m - a simple unit test for p_poly_dist. Plots the results of a call to p_poly_dist function (see help for usage example)

p_poly_dist1.m – previous version (ver. 1.0) of p_poly_dist.m

Notes on version 2.0:

Oct 2, 2015 - version 2.0 (Michael Yoshpe). The original ver. 1.0 function was completely redesigned. The main changes introduced in Ver. 2.0 are:

1. Distances to polyline (instead of a closed polygon in Ver. 1.0) are returned. The polyline can optionally be treated as a closed polygon.

2. Distances from multiple points to the same polyline can be found.

3. The algorithm for finding the closest points is now based on coordinate system transformation. The new algorithm avoids numerical problems that Ver. 1.0 algorithm could experience in "ill-conditioned" cases.

4. Many optional output variables were added. In particular, the closest points on polyline can be returned.

5. Added input validity checks

Michael Yoshpe (2021). Distance from points to polyline or polygon (https://www.mathworks.com/matlabcentral/fileexchange/12744-distance-from-points-to-polyline-or-polygon), MATLAB Central File Exchange. Retrieved .

Created with
R12

Compatible with any release

**Inspired:**
Minimum distance between two polylines, Distance from a point to polygon

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

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Michael YoshpeHello, Rebecca, thank you very much for your rating!

I haven't tested the function with self-intersected polygons, but I think it should work, because the algorithm computes the distance to each segment independently, and then finds the minimum. It won't work for multiply connected polygon though, simply because there is currently no code inside the function that would separate the polygons, like in Matlab inpolygon function.

Regards

Michael

Rebecca JiangThis function is great -- it is efficient, accurate, and very useful.

One question: does the polygon have to be simply connected for this function to compute the distance accurately?

Thank you!

Gerald KoudijsGreat function, it save me a lot of time. As I needed did for a lot of data points, I needed to smartly select the points and crop the polyline to make to improve the computation time. I think the function itself could be optimized a bit as well to improve the computation time. Functionally however perfect.

Michael YoshpeHello, Wai Kei

I am glad you've found my function useful, but I don't think I fully understand what do you mean by saying "why need input the same dimensions of points and polygon?". As stated in the description, the function was specifically designed to find the distances of MULTIPLE point to a SINGLE polyline or polygon. Please see the description of the input variables:

% xp - vector of points X coordinates (1 X np)

% yp - vector of points Y coordinates (1 X np)

% xv - vector of polygon vertices' X coordinates (1 X nv)

% yv - vector of polygon vertices' Y coordinates (1 X nv)

The np and nv DO NOT have to be of the same value.

I suggest you have a look at test_p_poly_dist function. Please uncomment the lines after usage example and run the code, i think it might help you.

Best regards

Michael

wai keiImportant code, and help me work a lot.

but one question, why need input the same dimensions of points and polygon? If I want to lots of points, but only one polygon, do I need to split the points to the same dimensions of polygon,then use the command?

thanks a lot!

Juan MarrugoThis code was very useful to me. I needed the distance from a fault line to a set of points. At the beginning I thought of measuring the distance to the vertices in the fault line, but then I realized I needed to account for the closest distance to be in the segments as well. This code saved me a lot of time, thanks a lot. JC

SongqiuCould you please improve your codes for 3D points and polylines?

Muhammad UsmanDear Micheal,

The code is very useful. But I want your comment on another code written for the same purpose.

The code name is "dpoly". It also find the distance between point(s) and a polygon. You can find it here

https://people.sc.fsu.edu/~jburkardt/m_src/dist_plot/dpoly.m

Can you please comment on this?

I would greatly appreciate it if you kindly give me some feedback on this.

Regards

Muhammad Usman

Muhammad UsmanSamuel LondnerVery useful function, fast and exact.

fobooDear Michael,

Thanks so much for this great function!

I have one question/feature request.

Based on the current output of the function, I can derive what I call the "closest line segment function", i.e., the line segment on which the closest point has been located.

Let me denote this function as

f(t) = a + t*(b-a)

where "t" is between 0 and 1, and where the points "a" and "b" can already be deduced from the current output of the function.

However, I am also looking for the value of "t" which corresponds to the smallest distance (as computed).

Currently, I am trying to compute value of "t" via

t = ([x*, y*] - a )/ (b-a),

where x*,y* are the closest points on the polyline (as returned by the function).

This is cumbersome, as it gives me two possible solutions for "t".

Would be possible that the function itself also returns the value of "t" corresponding to the smallest distance?

fobooAddendum to my previous comment:

I did not mean

t = (distance - a )/ (b-a).

Instead, I meant

t = ([x*, y*] - a )/ (b-a),

where x*,y* are the closest points on the polyline (as returned by the function).

snoop58Hi Michael,

thanks to your reply, I started digging and found that my polyline generator does not always make consecutive connected points. Moreover, the function trisurf() that I have used to display the polyline, plots a confusing edge-only graph, even when the points are not consecutive!

Anyhow, after passing a convhull() filter over my polylines, I was able to apply your function like a charm. Thank you!

Michael YoshpeHello, snoop58

Thank you very much for your comment.

Unfortunately, I could not reproduce the problem you mentioned entirely, since I could only load the polylines without the test code (the code itself is in pdf file instead of .m file).

Still, I have a few remarks.

First, setting the optional find_in_out flag to 'true' treats any OPEN polyline as A CLOSED one. I quote from the function's help section:

% find_in_out - logical flag. When true, the polyline is treated as a

% closed polygon, and each input point is checked wether it is inside or

% outside of this polygon. In such case, an open polyline is automatically

% closed by adding a last point identical to a first one.

So, using find_in_out = 'true', effectively adds a segment to the polyline which connects last point to the first one, and the polylines in your figures are plotted as OPEN without this closing segment.

Second, I could not reproduce the polylines themselves as they are plotted in your document. For example, I used the following command to plot the A polyline for the wrong case:

polylines = load('polylines');

S2 = polylines.S2;

plot(S2{1}.vertices(:,1), S2{1}.vertices(:,2), 'm-')

and the line I got is very different from what you plotted (try this command yourself). Note that xv,yv input vectors are treated as the coordinates of consecutive connected points.

I am sorry I could not be of more help to you. I suggest you read the help section again and check the input data. If you still have problems with the results, please give both the data AND the code as .m function, so I could try to debug it further.

Best regards

Michael

snoop58Hello Michael. Great work, thank you.

The function is very nice and does exactly what I want; however, it fails in some cases. I attach a PDF and mat file, which contain a couple of screenshots of correct and incorrect distance calculation, a mat file with polyline triangulation, and a code to run the data.

https://www.dropbox.com/s/dg77yepl0btsykg/distance_test.pdf?dl=0

https://www.dropbox.com/s/nt17douftlyhj8u/polylines.mat?dl=0

In short: there are 2 sets of 4 open polylines and a single point in fixed location.

In the first case, the point is outside all 4 polylines, and the shortest distances to each polyline are calculated correctly.

In the second case, the point is outside top & bottom polylines (distance is correct) but inside the left & right polylines: the min distance is wrong, as well as projection coordinates.

I have tried various polyline combinations, and I always get wrong distance calculation when the left/right polylines cover the point of interest. I have also tried to keep the polylines open or closed (find_in_out = true/false) with identical results.

Can you please help me trace and fix the error? Thank you!

Michael YoshpeHello, Philippe

Thank you very much for your comments

I ran the unit test test_p_poly_dist (included with my submission) on Matlab R2015a, and it was OK.

Please let me know which version of Matlab caused the problems you observed, and what were the inputs, specifically the coordinates of vertices (xv, yv), and points (xp, yp).

Best regards

Michael

Philippe Lebelalso, on line 291 (or so...)

d_min(is_vertex) = dpv_min(is_vertex);

Philippe LebelHi, I found out about some bugs that might have been introduced by new matlab releases.

When trying to compute multiple points, the function did not work at all.

here are the fixes i had to implement in order to make it work:

at line 219:

if length(xp)<= 1

Ppr1 = Ppr1 - a1';

Ppr2 = Ppr2 - a2';

else

Ppr1 = Ppr1 - a1;

Ppr2 = Ppr2 - a2;

end

then at line 316:

I_cr_min_temp = I_cr_min(~is_vertex);

if ~isempty(I_cr_min_temp)

I_cr_min_temp = I_cr_min_temp + [0:length(I_cr_min_temp)-1]'*(nv - 1);

xc_temp = xc';

yc_temp = yc';

x_d_min = xc_temp(I_cr_min_temp);

y_d_min = yc_temp(I_cr_min_temp);

end

thank you for your work!

Turgut S.David SHi Michael, thanks for this very useful function. Here is an example of a group of polygons that form "donut" holes. Interestingly, the code correctly parses the right donut hole, but at the left one, it adds a line connecting the outer polygon to the inner one. Is there a way to avoid this? https://www.dropbox.com/s/9yxilcfx635bvvi/p_poly_dist_example.m?dl=0

Michael YoshpeCharles, you might want to have a look on a p_spoly_dist function I just added (http://www.mathworks.com/matlabcentral/fileexchange/52734-p-spoly-dist).

I think it's pretty close to what you want, EXCEPT it DOES NOT for now return an in/out indication for a point in a closed spherical polygon. I simply could not find a reliable algorithm for this task.

Pablo NiloVery very useful, perfect for geographical data filtering and screening. Thanks a lot!!

Just a question from an Earth Observartion investigator: Could it be faster? We need to deal with satellite tracks (120.000 records each orbit).

charlesvery helpful, thank you. however, is 2D. is there a comparable function for finding the distance on surface of a sphere, from point on the sphere's surface to polygon on the sphere's surface? point and vertices given in lat-long?

IlyaTo Eric Schmitz's comment:

why we need "clear id"? By the way, it consumes quite a lot of time.

TreySorry, my previous comment is in regards to Bastian's update, not the original by Michael Yoshpe

TreyHi, thanks for the code. It's fast, but the 23 July 2013 code has a bug.

There's an 'X' in

uY = [Y(2:end, :); X(1, :)] - Y;

that should be 'Y':

uY = [Y(2:end, :); Y(1, :)] - Y;

BastianI've changed the code a little bit. So it can be used with the same syntax as the MATLAB-Funktion inpolygon.

function d = distancepolygon(x, y, X, Y)

m = size(X, 1) * size(X, 2);

n = size(x, 1) * size(x, 2);

X = X(:) * ones(1, n);

Y = Y(:) * ones(1, n);

x = ones(m, 1) * x(:)';

y = ones(m, 1) * y(:)';

% projection

uX = [X(2:end, :); X(1, :)] - X;

uY = [Y(2:end, :); X(1, :)] - Y;

LAMBDA = ((x - X) .* uX + (y - Y) .* uY) ./ (uX.^2 + uY.^2);

PX = X + LAMBDA .* uX;

PY = Y + LAMBDA .* uY;

% distance to projections points which is inside a segment

PD = NaN(size(x));

b = ((LAMBDA > 0) & (LAMBDA < 1));

PD(b) = sqrt((PX(b) - x(b)).^2 + (PY(b) - y(b)).^2);

% distance to polygon point

D = sqrt((X - x).^2 + (Y - y).^2);

% distance

d = min([D; PD]);

% change sign if the point is inside the polygon

b = inpolygon(x(1, :), y(1, :), X(:, 1), Y(:, 1));

d(b) = -d(b);

Xavier NievesHi, I have been trying to create a function similar to this one... if I were to have a random polygon but always have as an input the x,y coordinates of the vertices's and and the x,y coordinate of a point. can I use this function to determine the minimum distance of the point to the polygon, and if I can should express the inputs in what way... ?? nice function

Cristian DThanks for this function! it was very usefull!

Alejandro WeinsteinI added the suggested test from Eric Schmitz, and also return the point of the polygon closest to the point. It is available here:

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=19398&objectType=file

Eric SchmitzFor the case where a polygon rib is either horizontal or vertical, the test for checking if the projected point is inside a segment may not work all the time because of numerical accuracy. I added a test to force either the x or y coordinate of projected point to be set equal to the corresponding rib x or y coordinate (before the section computing idx_x and idx_y):

id = find(diff(xv)==0);

xp(id)=xv(id);

clear id

id = find(diff(yv)==0);

yp(id)=yv(id);

Jesse Collinsthis function would benefit greatly from also returning the point on the given polygon closest to the point from which the distance is to be found.

Kardi Teknomo