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.
Comment only
27 Aug 2014
distance2curve
Find the closest point on a (n-dimensional) curve to any given point or set of points
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
Comment only
31 Jan 2014
distance2curve
Find the closest point on a (n-dimensional) curve to any given point or set of points
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
Comment only