Find the closest point on a (n-dimensional) curve to any given point or set of points
I've seen many people ask for a way to find the closest point on a curve from some given point in space. If the curve is a piecewise linear one, this is not too difficult, since this reduces to finding the closest point on a line segment, and then testing each line segment. (With care, you need not test EVERY line segment.)
For a cubic spline however, this becomes more difficult, but still doable in a mathematical sense without an explicit optimization. Not too much more than roots is necessary.
Distance2curve allows you to specify a set of general points in an n-dimensional space as a connected space curve. A spline (or pchip) is then fit in terms of arc length along the curve, and the closest point identified.
1.1 | Bug fix: if nargout >= 2 |
Joshua Diamond (view profile)
Hey John,
Thanks for all your work on this. Unfortunately, I am having some problems. The result is always literally a point of curvexy, rather than a point on the surface expressed by curvexy. Does that make sense? Here is my code. Thank you!
curvexy=[0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.299972441470565,0.0100000000000000;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.283719107177895,0.107400892521526;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.236720405266962,0.194246887007779;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.164069374349813,0.261126872401726;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.0736388809232166,0.300793364536036;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.0247715278763768,0.308947883498916;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.120497557260707,0.284706760630339;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.203165806603595,0.230696897732891;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.263817888358687,0.152771101500882;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.295881208353804,0.0593738410993592;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.295881208353804,-0.0393738410993591;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.263817888358687,-0.132771101500882;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.203165806603595,-0.210696897732891;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.120497557260707,-0.264706760630339;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,-0.0247715278763769,-0.288947883498916;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.0736388809232165,-0.280793364536036;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.164069374349813,-0.241126872401726;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.236720405266962,-0.174246887007778;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.283719107177895,-0.0874008925215257;0.0350000000000000,0,0.0100000000000000;0.0525582744890453,0.299972441470565,0.00999999999999993]
mapxy=[.042 .3 .3]
xy = distance2curve(curvexy,mapxy,'linear');
plot3(curvexy(:,1),curvexy(:,2),curvexy(:,3),'k-o')
hold on
plot3(mapxy(:,1),mapxy(:,2),mapxy(:,3),'r*')
line([mapxy(:,1),xy(:,1)]',[mapxy(:,2),xy(:,2)]',[mapxy(:,3),xy(:,3)]','color',[0 0 1])
Pei-Ann Lin (view profile)
thank you for a great script. easy to understand and adapt!
Tillmann (view profile)
Jonathan Ahlstedt (view profile)
Bosheng (view profile)
Mithun Nallana (view profile)
Setafania Vaga (view profile)
Perfect! Exactly what I needed! Thank you so much!!!!!
Kevin Hughes (view profile)
Thanks John another great program. As usual it was just what I needed.
Diogo (view profile)
David (view profile)
TYo (view profile)
Great tool!
I just have an inquiry about your explanation about the interpolated curve , in response to Jean.
When I try your approach to evaluate a spline defined by xx and yy in a 1x30, I get:
Error using distance2curve (line 149)
curvexy and mapxy do not appear to live in the same dimension spaces
Any ideas what goes wrong?
Newsam Niu (view profile)
John D'Errico (view profile)
I'm sorry, but I simply refuse to do consulting in the comments of my submissions. This question has nothing directly to do with distance2curve anyway, except that you are trying to compute some kind of distance that involves curves.
As for your problem, it is basic calculus, once you generate the respective spline fits. And no, I cannot show you how to solve it here.
Finally, please do not contact me by mail with this question, as I also do not answer direct questions. Were I to do so, I would quickly become overwhelmed with the volume of questions.
Ajabofu Augoye (view profile)
Hi John,
Thanks for your contribution so far. I have two discrete curves and I need to find the normal distance from each discrete point of one of the curves to the points where this normal intersects the other curve. Both curves are completely different functions and intersects at some points.
Thank you.
Kobi
John D'Errico (view profile)
If x is the SECOND index of the array and y the first, so that when we plot it, the plot looks like we expect with x on the x axis. I'll also assume that y=0 is at the bottom of the resulting image array when plotted.
curvexy = [30 0;50 15;60 30;70 15;90 0];
[x,y] = meshgrid(1:201);
xy = [x(:),y(:)];
[xyc,distance] = distance2curve(curvexy,xy,'linear');
k = (distance < 3);
sum(k)
ans =
493
201*201
ans =
40401
A = flipud(accumarray(xy(k,[2 1]),1,[201,201]));
image(A)
Markus Rühlmann (view profile)
Hi John,
Thank you for the nice programm.
My name is Markus and I've a matlab-problem. I think I could solve it with your program.
Following:
I have a 201x201 double field.
The individual points of this field have either the value 1 or 0.
My goal now is to programm a parabel or a curve of ones in this field.
And my idea was to construct a curve (with serveral lines) through this field and all points with a certain distance to the curve should assume the value 1.
My code for that propose is:
curvexy = [30 0;50 15;60 30;70 15;90 0];
for i = 1:1:201
for n = 1:1:201
mapxy = [i n];
[xy,distance,t] = distance2curve(curvexy,mapxy,'linear');
if distance < 3
Field(mapxy)=1;
end;
end;
end;
But it doesn't work. So I try to explain my Code with words:
In curvexy I've built the curve with serval lines.
In mapxy are all points of my 201x201 double field.
If the distance of curvexy to the points of my field is less than 3, the points of my field should get the value 1.
I hope you understand, what is my problem.
Please help me. :)
Kind regards
Markus
John D'Errico (view profile)
Hi Diego,
distance2curve returns the closest point on a curve to a given point. There are two possibilities (assuming that one uses this function with a spline interpolant):
1. If that closest point is at some internal point on the curve, then the line connecting those points IS perpendicular to the curve.
2. Only if the closest point is one of the end points of the curve, then that line need not be perpendicular to the curve.
Of course, if you specify a linear interpolant, then it makes no sense to require the line of minimum distance to be always normal to a non-smooth curve.
So given the spline or pchip interpolating methods that distance2curve offers, you WILL get a normal line as long as that closest point is not an endpoint of the curve. Were that closest point of approach not perpendicular, then there would be another point on the curve that was closer to your given point.
If you are worried about the endpoints of the curve sometimes being returned as the point of minimum distance, yet that does not yield a perpendicular, what exactly would you have me return in that case? The point of minimum distance is what this tool is designed to return.
For example, suppose I were to pose a curve that is a semi-circle, half of a circle of radius 1, where the center of the circle is at the origin, (0,0). I'll set the end points of the semi-circle to be (1,0) and (-1,0), a circular arc that goes from 0 degrees to 180 degrees, that goes through the point (0,1). This semi-circle lives in quadrants 1 and 2 of the (x,y) plane.
Now, choose the point at (0+eps,-5). What is the closest point on that semi-circle? It will be an end-point of the curve at (1,0), since I offset the point by eps in the x direction to yield a unique solution. The distance to the circular arc is essentially:
sqrt(25 + 1)
ans =
5.09901951359278
But the line connecting those two points is not perpendicular to the curve. It cannot be so.
There is only one point on the curve where the connecting line is indeed perpendicular to the curve, and that point is the MAXIMUM distance to the curve. So the line segment that connects the two points (0+eps,-5) and (0,1) has length 6. There is no point on the curve that is farther away, yet that point of MAXIMUM distance is the only point with a perpendicular line that connects the point and the curve.
So, what you are asking for is what you already get in context of a curve that has at least C1 differentiability, WHENEVER that is possible. (Either spline or pchip as the interpolant.)
You may perhaps be confused IF you are drawing plots without the axes being equal. In that case, it can appear as if the line segment is not indeed perpendicular to the curve. The fix for this is to use
axis equal
after you use plot. This sets the aspect ratio for the plot to unity, so that a tick mark along each axis will represent the same distance. Otherwise, perpendicular lines will not look perpendicular.
For example, try this:
plot([0 1],[0 2],'r-',[0 6],[1 -2],'g-')
the two line segments are perpendicular, but they will not appear to look as if they are. The first line segment has slope 2. The second line segment has slope -0.5, so with slopes that are negative reciprocals, they must be perpendicular. Unless you also use
axis equal
the two line segments will not appear perpendicular.
Diego (view profile)
Thanks for the useful code!
Is there an easy way to calculate the normal (perpendicular) distance between mapxy and curvexy instead of the minimum distance?
I would like to use this code in an image analysis application, where the normal distance between two lines in needed (to calculate a thickness distribution).
Diego (view profile)
Stephanos Kythreotis (view profile)
Very useful program. Thanks!
Michael (view profile)
It would be awfully convenient if the code also returned the total arc length (for px,py pair). That way fractional arc length could be trivially converted to actual pixel units. This length will depend on interpolation method of course. And yes, it is often a one-liner to compute separately hence the comment about sheer convenience.
John D'Errico (view profile)
I suppose it would have been nice to add that as a capability, to pass in a spline itself to this code, but it accepts only a list of points. There are several different forms of splines, more than I wanted to worry about. So distance2curve just builds its own spline from scratch.
If you have a spline function of the form
y = f(x)
then it is simple. f(x) is merely a piecewise (often cubic) polynomial. So in each segment to find the closest point to the point (x0,y0), compute the polynomial
(x-x0)^2 + (f(x) - y0)^2
differentiate wrt x, and set to zero, using roots to solve. Then pick the real roots that lie in the interval in question, and choose the point on the curve with the smallest distance. Compare the distance found to the distance to the point at the endpoints of the interval. Then just loop over the intervals.
Other spline forms will be similar. For example, distance2curve builds a spline of the form (x(t),y(t)), where t is a cumulative linear arc length parameter. This allows the curve to be a completely general one that need not be a single valued function, and it works in any number of dimensions. The general idea is the same as above though.
Sarah (view profile)
Hi John,i learned a lot from your program.Thank you very much.But i have a question,how can i find the closest point on a B-Spline?
Genevieve (view profile)
John D'Errico (view profile)
To completely resolve the question, perhaps we might look for a symbolic solution.
syms x
D = (x - 4)^2 + (sin(x) - 0)^2;
Differentiating, and setting to zero
dx = diff(D,x)
dx =
2*x + 2*cos(x)*sin(x) - 8
Which we might recognize reduces to:
2*x + sin(2*x) - 8 = 0
I don't see any obvious solution to the mixed trig equation, but it is easy enough to throw the symbolic toolbox at it.
solve('2*x + sin(2*x) - 8 = 0')
ans =
3.6019704355966679057142560063871
See that this agrees with the distance2curve prediction (3.6019704316906) out to about 9 significant digits.
Or, we could use fzero, or for the numerical gear heads like me, write this as a fixed point iteration.
fp = @(x) (8 - sin(2*x))/2;
x = [inf,1];
iter = 0;
while abs(diff(x)) > (10*eps)
iter = iter + 1;
x = [x(2),fp(x(2))];
end
x(2)
ans =
3.60197043559667
iter
iter =
65
John D'Errico (view profile)
Hi Friedrich,
Ok, I tried your test. You did not tell me what set of points on the curve you tried, so I was arbitrary. It will be close to your result. A different spacing, i.e., (0:.1:6) gave me essentially the same results as I got here.
x = 0:.01:6;
y = sin(x);
[xy dd]=distance2curve([x(:) y(:)], [4, 0])
xy =
3.6023 -0.44458
dd =
0.5965
plot(x,y,'g-',[xy(1),4],[xy(2),0],'r-o')
axis equal
Sorry, but I have a very difficult time believing that it did NOT find the closest point on that curve. Perhaps you did not use "axis equal" if you plotted it. My distance2curve code is not perfect here, because it uses an interpolant.
Note that the line connecting the point it found and your chosen point is orthogonal to the curve. I could even compute the tangent line at that point to convince you that it is an accurate estimate.
That is, the slope at the point distance2curve found was
cos(xy(1))
ans =
-0.89574
This is the slope of the tangent line. The normal vector to that point would have a slope that is the negative reciprocal.
-1./cos(xy(1))
ans =
1.1164
And if I now compute the slope of the line that connects the two points, we get...
(xy(2) - 0)./(xy(1) - 4)
ans =
1.1179
Which seems close enough for government work, yielding a number that is surprisingly close to that which was predicted. Finally, note that distance2curve with NO third argument uses a piecewise linear interpolant, which is fast, but not quite as accurate. The spline interpolant will be considerably more accurate.
format long g
[xy dd]=distance2curve([x(:) y(:)], [4, 0],'spline')
xy =
3.6019704316906 -0.444286584307468
dd =
0.59650490881816
(xy(2) - 0)./(xy(1) - 4)
ans =
1.11621502441275
-1./cos(xy(1))
ans =
1.11621504194962
See that the slopes now agree to within 8 significant digits, about as accurate as I could hope for.
Oh, and just in case you want more affirmation, try this further test.
fun = @(x) norm([x,sin(x)] - [4 0]);
[xmin,dval] = fminbnd(fun,0,6)
xmin =
3.6019676945136
dval =
0.596504908814437
So I do hate to say it, but I think you are mistaken. Distance2curve clearly nailed it. I could not have done better with a hammer.
John
Friedrich Levi (view profile)
Dear John, I tested this function with y=sin(x), and then
[xy dd]=distance2curve([x(:) y(:)], [4, 0])
The problem is that the result is
xy=[3.6297 -0.4561].
While xy is in the curve y(x), it is clearly not the closest point in y(x) to [4 0]. Could you please let me know what's going on? Thanks.
Cui (view profile)
Thanks for you kind reply very much. I should think more about my problem. A very useful tool, Thanks very much.
John D'Errico (view profile)
Cui - I'm sorry, but it cannot do magic, or even just things it was not programmed to do. It does only interpolate a set of points, then efficiently finds the closest point on that interpolated curve. This is done efficiently to find a globally closest point since the interpolated curve is a piecewise polynomial.
If you insist on not allowing the curve to be interpolated, then you must write your own code to find a minimal distance. Beware then that you do not find a locally closest point as opposed to the globally closest point, since the optimization problem for a general curve may also have local minima.
Cui (view profile)
hi, if I have the equation of the curve, how can I get the closet point.
If I use your program, then I must create many points of the curve first, because you input is points of the curve rather than the equation.
Shaun VanWeelden (view profile)
This file just saved me SO much time!! Thanks a million!
sara (view profile)
thaaanx alot :-)
Jean-Yves Tinevez (view profile)
Thanks John. I tested with interparc and managed to produce correct curves. I also understood that the discrepancies were due to me using the cumulative interchordal length instead of the proper arc length along the target interpolated curve.
Thanks again
jy
John D'Errico (view profile)
The curves themselves are parametric splines, so if we define the cumulative (piecewise linear) chordal arc length as t, then the splines are defined as (px(t),py(t)). So essentially we have a pair of parametric splines (or more in higher dimensions.) As such, its not very useful if you really want the function, but you can interpolate points along that curve using interparc. I've tried to make the two tools consistent in this respect.
You can create those explicit splines themselves, but be careful, as distance2curve and interparc both interface to you in terms of the actual distance along that spline curve, as opposed to that cumulative chordal arclength.
Anyway, test it by picking some random curve in the plane.
px = rand(10,1);
py = rand(10,1);
Now find the distance to the origin, and the closest point on that curve to the origin.
[xy,distance,t_a] = distance2curve([px,py],[0 0],'spline')
xy =
0.16771 0.29285
distance =
0.33747
t_a =
0.95696
Note that distance2curve actually gives you the interpolated value along the curve at that closest point.
Now, can we interpolate this point on the curve using interparc? I'd hope that interparc can generate the same point for the value t_a.
interparc(t_a,px,py,'spline')
ans =
0.16771 0.29285
Actually, interparc is different down in the 7th decimal digit, which without thinking about it deeply, probably reflects some numerical algorithmic issues. For example, interparc uses an odesolver (ode45), triggering at zero crossings of a variable. Instead, distance2curve uses roots as a root finder on a polynomial that defines the derivative of the square of the distance.
The distinct algorithms, coupled with the tolerances employed there are probably sufficient to explain that tiny variation.
If this does not answer your question, send me an e-mail and we can carry the discussion further offline.
Jean-Yves Tinevez (view profile)
Hi John,
Again a great submission in a flawless catalogue.
I just have a side question: How would I do to get the interpolated curve we measure the distance to in your contribution? My attempts to approach the implicit curve we compute the distance to here have failed (see fex #22948).
Best
jy
John D'Errico (view profile)
Is it easy? I suppose the definition of easy depends on the person doing the writing.
I have written a tool that finds the closest point on a general manifold in 2 or 3 dimensions, or for low dimensional manifolds in higher dimensions. So if one wanted to find the closest point on a triangulated surface in 3-d, this would be no problem. It is even reasonably efficient. That tool could also solve for the closest point on a triangulated 2-manifold in a 4 or higher dimensional space. For higher dimensions than tessellated 2 manifolds, this is still eminently doable, but I had to stop somewhere. This tool is part of my simplicialcomplex toolbox, something that never seems to be ready to release to the general public.
However, that code is also limited to piecewise linear manifolds, since I needed to use roots to solve for the cubic case in distance2curve. Of course, roots is not available for more than one variable, so it would get a bit more nasty in higher dimensions.
The other question is how to find the minimum distance between two general curves, or even from a curve to a line. I am planning to write a tool for the minimal distance from a curve to a straight line. The minimum distance between two general space curves (1-manifolds) is somewhat more difficult, because it becomes essentially a nonlinear optimization in two variables.
Christophe Lauwerys (view profile)
As always, great submission, thanks.
Could one easily extend this algorithm to finding the closest point on a surface from given a point in space? Or from another given curve, surface, ... ?