Code covered by the BSD License  

Highlights from
distance2curve

5.0

5.0 | 7 ratings Rate this file 72 Downloads (last 30 days) File Size: 16.3 KB File ID: #34869
image thumbnail

distance2curve

by

 

01 Feb 2012 (Updated )

Find the closest point on a (n-dimensional) curve to any given point or set of points

| Watch this File

File Information
Description

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.

Acknowledgements

This file inspired Xy2sn.

Required Products MATLAB
MATLAB release MATLAB 7.12 (R2011a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (16)
27 Aug 2014 John D'Errico

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.

27 Aug 2014 Sarah

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?

11 Aug 2014 Genevieve  
01 Feb 2014 John D'Errico

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

31 Jan 2014 John D'Errico

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

31 Jan 2014 Friedrich Levi

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.

29 Dec 2013 Cui

Thanks for you kind reply very much. I should think more about my problem. A very useful tool, Thanks very much.

28 Dec 2013 John D'Errico

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.

28 Dec 2013 Cui

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.

21 Jun 2013 Shaun VanWeelden

This file just saved me SO much time!! Thanks a million!

26 Nov 2012 sara

thaaanx alot :-)

07 Sep 2012 Jean-Yves Tinevez

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

06 Sep 2012 John D'Errico

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.

06 Sep 2012 Jean-Yves Tinevez

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

02 Feb 2012 John D'Errico

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.

02 Feb 2012 Christophe Lauwerys

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, ... ?

Updates
27 Feb 2013

Bug fix: if nargout >= 2

Contact us