101 views (last 30 days)

Show older comments

John D'Errico
on 30 Apr 2017

Edited: John D'Errico
on 30 Apr 2017

Minimum Euclidean distance to a general surface is a somewhat nasty problem to solve. With ONE independent variable, IF the function is a polynomial one, there is a solution, of sorts. Not trivial. But a solution. If the surface is non-polynomial, it gets nasty. And in multiple dimensions, things can get messy too. Lets see what happens with a one variable problem.

I'll pick a random cubic polynomial. Not even something where I'll pick the coefficients myself.

format long g

coef = randn(1,4)

coef =

0.318765239858981 -1.30768829630527 -0.433592022305684 0.34262446653865

I'll use these as coefficients of a cubic polynomial, to be evaluated by polyval.

ezplot(@(x) polyval(coef,x),[-3,5])

Now, suppose we choose some arbitrary point in the plane? What is the

xy = randn(1,2)*2

xy =

2.48288696432624 2.97939521557093

hold on

grid on

plot(xy(1),xy(2),'ro')

What is the point of closest approach? Be careful. This might change your mind. :)

axis equal

How would I solve for the projection in terms of minimum Euclidean distance? I would solve for the point that minimizes the function:

(y(x) - y0)^2 + (x-x0)^2

as a function of x. Do that by differentiating, and then searching for the roots of the polynomial. So the square of the Euclidean distance i the (x,y) plane is just:

D2 = (P - xy(2))^2 + (x-xy(1))^2;

xroots = double(solve(diff(D2,x)))

xroots =

0.167082630638436 + 0i

2.91225594186202 + 0i

4.72410209755873 + 0i

-0.48309085270852 - 1.29552403846027i

-0.48309085270852 + 1.29552403846027i

If we ignore the complex roots,

double(subs(D2,xroots))

ans =

12.8937799521223 + 0i

50.8355152126218 + 0i

5.09171123781745 + 0i

6.70892471381755 + 7.41419202614555i

6.70892471381755 - 7.41419202614555i

the real root with the minimum distance (squared) in the (x,y) plane is the third root. (Note that there will ALWAYS be at least one real solution to this problem. I'll leave the simple proof of that claim to the interested student.)

ezplot(@(x) polyval(coef,x),[-3,5])

hold on

grid on

plot([xy(1),xroots(3)],[xy(2),polyval(coef,xroots(3))],'r-o')

axis equal

As you can see, with the axes equal in spacing, the line is orthogonal to the curve, as it must be for a projected point of minimum distance.

So not trivial to solve for the point of minimum Euclidean distance, but not terribly difficult either, at least in one dimension. With very little extra thought, I could have written it as a call to roots, and not gone the symbolic route at all.

With two independent variables though, now we will end up with two multinomial equations, in two unknown. A solution will again always exist for the minimum distance, but we will not have recourse to a tool like roots at all. And there will again be in general multiple solutions, so we would need to find the one with the smallest distance. Solving for a zero of the gradient only yields a stationary point.

So possible, but not perfectly trivial. And you will not be able to compute a solution in a vectorized form.

DS
on 17 Sep 2020

Reply to 2): in case someone is interested in getting the complete set of results above from John's comment:

format long g

coef = [0.318765239858981 -1.30768829630527 -0.433592022305684 0.34262446653865];

ezplot(@(x) polyval(coef,x),[-3,5])

xy = [2.48288696432624 2.97939521557093];

hold on

grid on

plot(xy(1),xy(2),'ro')

axis equal

x = sym('x');

a = sym(coef(1));

b = sym(coef(2));

c = sym(coef(3));

d = sym(coef(4));

y = a*x^3 + b*x^2 + c*x + d;

D2 = (y - xy(2))^2 + (x-xy(1))^2;

xroots = double(solve(diff(D2,x)));

xroots

double(subs(D2,xroots))

ans =

12.8937799521223 + 0i

6.70892471381757 + 7.41419202614556i

6.70892471381757 - 7.41419202614556i

50.8355152126213 + 0i

5.09171123781739 + 0i

hold on

grid on

plot([xy(1),xroots(5)],[xy(2),polyval(coef,xroots(5))],'r-o')

axis equal

Joseph Cheng
on 28 Apr 2017

you could try something like this. my curve and fit method is probably different however what you can do is use find() to determine which ones are within some distance/error from the fitted curve for that point.

close all;

load franke

f = fit( [x, y], z, 'poly23' )

plot(f, [x,y], z)

withintol = find(f(x,y)-z<0.01);

hold on

ax=gca;

plot3(x(withintol),y(withintol),ax.ZLim(1)*ones(size(withintol)),'.')

Joseph Cheng
on 28 Apr 2017

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

Start Hunting!