Project 3D points to a surface

53 views (last 30 days)
psprinks
psprinks on 27 Apr 2017
Edited: DS on 17 Sep 2020
I have fit a surface through a set of x,y,z points. Image attached. I would like to be able to project points that lie within some distance threshold to that surface. Could anyone advise on a method for doing this?
  2 Comments
Sean de Wolski
Sean de Wolski on 28 Apr 2017
Is it the shortest euclidean distance to the surface of the distance to the surface at x, y for a specific point?
psprinks
psprinks on 30 Apr 2017
Yes , I think the minimum euclidean distance would be fine for this case.

Sign in to comment.

Accepted Answer

John D'Errico
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.
  2 Comments
psprinks
psprinks on 1 May 2017
John, Thanks for your time and effort in this response. This has become quite the rabbit hole I've fallen into. I thought this would be a way to improve some analyses I have already performed on one of my datasets for my M.S. thesis. I didn't realize it would be this involved. I have a few questions if you don't mind?
Since I didn't provide much information in my original post. My main objectives here are to be able to calculate minimum Euclidean distances from the earthquake locations to the surface of best fit, and additionally to be able to project the points in the 3d space onto their locations on the surface.
1)I used your gridfit tool to fit the surface in the image above. Should I use a tool like polyfitn to generate the coefficients for the polynomial that describes the surface since I have 2 independent variables?
2)I don't have much experience using the symbolic math toolbox in Matlab. Is that what you're referencing when you say you could have just as easily done the example with a 'call to roots' instead of the symbolic route? You wrote:
D2 = (P - xy(2))^2 + (x-xy(1))^2;
Is P (i.e. y(x) ) the equation for the polynomial describing the curve? How can this statement be evaluated if P isn't defined? This is likely something I'm missing from my lack of use with the symbolic toolbox.
3)Can you help me understand why the solution can't be in vectorized form?
I can provide code and data if needed to clear anything else up. Thanks again.
DS
DS on 17 Sep 2020
Edited: 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

Sign in to comment.

More Answers (1)

Joseph Cheng
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)),'.')
  2 Comments
psprinks
psprinks on 28 Apr 2017
Hi Joseph, Thanks for the reply. I just realized my image wasn't attached. It's there now. I used grift to fit the surface as I think it does a better job than just fit. Why do you load franke? I'm looking at the help for franke and it's not very informative. I'll give your method a try and see what I get.
What I'm really trying to do is find the orthographic/orthonormal projection for each point in the point cloud, but first I want to get rid of the points that are too far away.
Joseph Cheng
Joseph Cheng on 28 Apr 2017
Oh, i was just too lazy to generate a dummy curve and points to franke was just quickly available. now that i think about it just in the Z direction may not be proper thresholding if we're considering xyz jittering points. hmmm....

Sign in to comment.

Products

Community Treasure Hunt

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

Start Hunting!