Intersection of two sfit planes?

I fitted two planes to two datasets, and the planes are outputted as 1x1 sfits (see figure).

How do I get the intersection between these two planes? I want to be able to extract the xyz info along this line.

I'm not used to working with symbolic equations in MATLAB, so I tried this, without much success:

%
%Extract surface formula
cliff=formula(fit_cliff);
platform=formula(fit_platform);
%Extract formula Coefficients
coef_cliff=coeffvalues(fit_cliff);
coef_platform=coeffvalues(fit_platform);
c1=coef_cliff(1);
c2=coef_cliff(2);
c3=coef_cliff(3);
c4=coef_cliff(4);
c5=coef_cliff(5);
p1=coef_platform(1);
p2=coef_platform(2);
p3=coef_platform(3);
p4=coef_platform(4);
p5=coef_platform(5);
%Write equation and solve for intersection
syms x y 
z1 = c1+(c2*x)+(c3*y)+(c4*(x^2))+(c5*(x*y));
z2 = p1+(p2*x)+(p3*y)+(p4*(x^2))+(p5*(x*y));
intersect = solve(z1==z2)

Answer 1

[Star Strider]

I would subtract one from the other and then use the contour function on the result to determine the zero value. See the contour documentation for Display Single Contour Line but use:

v = [0, 0];

instead to find the intersection of the result of the subtraction. See the documentation on ContourMatrix to return these values to your workspace so you can use them in the rest of your code.

Answer 2

[John D'Errico]

So what was the problem? :)

Let me guess. This problem will come down to solving what is effectively a 4th degree polynomial. The results will look nasty. Lots of roots. Lots of things that say root of. So the solution you got looked like a mess. (Because it is, and it must be so.) The fact is, there is no simple formula for the solution. Certainly not anything that will be easy to look at.

However, there is a simple solution. You know the range in x and y that you care about. So pick a sequence of values for x (or y.) Just loop over this set of values. So for any value of x, now use vpasolve to solve for the solution to your problem, with x substituted into that expression. Or you can use fzero.

For any x, you get the value of y, to a reasonably high degree of precision. You will actually get a set of probably 4 values for y for any x. Pick the REAL root that falls in the region of interest. This gives you a closely spaced set of solutions to the problem. Easy to do. It gives you what you want. And the nice thing is, all of these points on the curve are essentially exact, to within the ability of the solver to give you high accuracy.

You can use a spline to approximate the curve if you want something smooth. No matter what, you WON'T get a simple function, as I said before. So just use a spline in the end.

[Julia Howe]

I ended up going with the plug and chug method - the point data is what I wanted in the end anyway. Thanks for the help!

Answer 3

[Teja Muppirala]

This works for me:

%%1. Making some test coefficients

c1=1; c2=2; c3=3; c4=4; c5=5;
p1=-1.1; p2=2.2; p3=-3.3; p4=4.3; p5=-5.5;

syms x y  

z1 = c1+(c2*x)+(c3*y)+(c4*(x^2))+(c5*(x*y))
z2 = p1+(p2*x)+(p3*y)+(p4*(x^2))+(p5*(x*y))

%%2. Do The solving

y_of_x = solve(z1==z2,y)
z_of_x = simplify(subs(z1,y,y_of_x))

%%3. Plot the result

figure
h1 = fsurf(z1,[-2 2]);
set(h1,'facecolor','r','facealpha',0.5);
hold  all
h2 = fsurf(z2,[-2 2]);
set(h2,'facecolor','b','facealpha',0.5);
hp = fplot3(x,y_of_x,z_of_x,[-2 2]);
set(hp,'LineWidth',2,'Color','g');
ylim([-2 2]);

This gives in the command window:

z1 =

2*x + 3*y + 5*x*y + 4*x^2 + 1

z2 =

(11*x)/5 - (33*y)/10 - (11*x*y)/2 + (43*x^2)/10 - 11/10

y_of_x =

(3*x^2 + 2*x - 21)/(105*x + 63)

z_of_x =

(x*(87*x + 44))/21

Replace "fsurf" with "ezsurf" and "fplot3" with "ezplot3" if your MATLAB is old and doesn't have these functions.

[Julia Howe]

Hi,

So conceptually this makes sense and I think I was on the right track before. It's when I get the solution and go to plot that everything goes wrong. This is the ezplot of z1:

The axes are completely wrong (see axes limits from original figure), and when I go to edit them, the plane isn't even located in the original coordinate system. How can this be? I pulled the equation and coefficients directly from the sfit planes. Until I figure this out, the intersection line has no meaning. It's the elevation data specifically that I need, so if the coordinates are off then the line is useless.