# How do I plot the line of intersection between two surfaces?

487 views (last 30 days)
MathWorks Support Team on 16 Nov 2009
Commented: arthur poquet on 12 Jun 2020
I have two surfaces, and I would like to plot the line in three-dimensional space where the two surfaces intersect.

MathWorks Support Team on 16 Nov 2009
There are three different scenarios to consider depending on how the two surfaces are defined. In general there are two different ways to define a surface: explicitly or implicitly. An explicitly defined surface is one in which the height of the surface (z) can be written as a function of the location (x and y). For example:
z = f(x, y)
z = 2*x + y^2
An implicitly defined surface is one in which z cannot be written as a function of x and y. The surface can instead be defined as the points which satisfy an equation of three variables (x, y, and z). The classic example of an implicitly defined surface is a sphere:
f(x, y, z) = 0
x^2 + y^2 + z^2 - 1 = 0
Case 1: The intersection of two explicitly defined surfaces.
In the case of two explicitly defined surfaces, we must find the difference between the two surface heights at each point and then trace the contour where that difference is zero. This will give us the x- and y-locations of points on the line of intersection. To find the z-locations of the points, we can interpolate one of the original surfaces at those points (it should not matter which surface is interpolated since the surface heights are equal at those points). Here is an example:
% Define the input grid
[x, y] = meshgrid(linspace(-1, 1));
% Calculate the two surfaces
z1 = y.^2 + 2*x;
z2 = 2*y.^3 - x.^2;
% Visualize the two surfaces
surface(x, y, z1, 'FaceColor', [0.5 1.0 0.5], 'EdgeColor', 'none');
surface(x, y, z2, 'FaceColor', [1.0 0.5 0.0], 'EdgeColor', 'none');
view(3); camlight; axis vis3d
% Take the difference between the two surface heights and find the contour
% where that surface is zero.
zdiff = z1 - z2;
C = contours(x, y, zdiff, [0 0]);
% Extract the x- and y-locations from the contour matrix C.
xL = C(1, 2:end);
yL = C(2, 2:end);
% Interpolate on the first surface to find z-locations for the intersection
% line.
zL = interp2(x, y, z1, xL, yL);
% Visualize the line.
line(xL, yL, zL, 'Color', 'k', 'LineWidth', 3);
Case 2: The intersection of an explicitly defined surface with an implicitly defined surface.
In this case, we must express the two surfaces as f1(x,y,z) = 0 and f2(x,y,z) = 0. We compute f1 and f2 over some region of space and compute the difference between these two fields (f3 = f1 - f2). However, this will include information about all of the level sets of f1 and f2. For the surfaces of interest, we only need to focus on the level set at 0. Therefore we must compute the heights of the explicitly defined surface over the input area, and interpolate the difference field (f3) on this surface. Then we can find the contour on this surface where the difference is zero and proceed as in Case 1 to find the x-, y-, and z-locations of the line of intersection. Here is an example:
% Define the input grid (in 3D)
[x3, y3, z3] = meshgrid(linspace(-1,1));
% Compute the implicitly defined function x^2 + y^2 + z^3 - 0.5^2 = 0
f1 = x3.^2 + y3.^2 + z3.^2 - 0.5^2;
% This surface is z = 2*y - 6*x^3, which can also be expressed as
% 2*y - 6*x^3 - z = 0.
f2 = 2*y3 - 6*x3.^3 - z3;
% Also compute z = 2*y - 6*x^3 in the 'traditional' way.
[x2, y2] = meshgrid(linspace(-1,1));
z2 = 2*y2 - 6*x2.^3;
% Visualize the two surfaces.
patch(isosurface(x3, y3, z3, f1, 0), 'FaceColor', [0.5 1.0 0.5], 'EdgeColor', 'none');
patch(isosurface(x3, y3, z3, f2, 0), 'FaceColor', [1.0 0.5 0.0], 'EdgeColor', 'none');
view(3); camlight; axis vis3d;
% Find the difference field.
f3 = f1 - f2;
% Interpolate the difference field on the explicitly defined surface.
f3s = interp3(x3, y3, z3, f3, x2, y2, z2);
% Find the contour where the difference (on the surface) is zero.
C = contours(x2, y2, f3s, [0 0]);
% Extract the x- and y-locations from the contour matrix C.
xL = C(1, 2:end);
yL = C(2, 2:end);
% Interpolate on the first surface to find z-locations for the intersection
% line.
zL = interp2(x2, y2, z2, xL, yL);
% Visualize the line.
line(xL,yL,zL,'Color','k','LineWidth',3);
Case 3: The intersection of two implicitly defined surfaces.
There is no direct way to compute the line of intersection between two implicitly defined surfaces. You can try solving the equation f1(x,y,z) = f2(x,y,z) for y and z in terms of x either by hand or using the Symbolic Math Toolbox.

Antone Chacartegui on 27 Jan 2020
I used this method (for case 1 above) and it worked to graph the lines of intersection. Is there a way that I can extract the curve function so that I can use it to calculate curvature at a point on the surface that the curve passes through?
J AI on 30 May 2020
is there a way to do it for multiple intersecting planes? more specifically, let's say multiple planes intersect, is there a way to find the vertices due to intersection? I have attached a figure as an example. The right side figure is a top-view of the left one. The envelope/top surface has 8 vertices and I am looking to find the coordinates of the planes. The equations for the three planes in the diagram are:
[X,Y] = meshgrid(0:0.01:1,0:0.01:1);
Z1 = -4.3587*X-2.9209*Y+88.4568;
Z2 = 91.4443*X-34.8074*Y+88.1613;
Z3 = 4.7371*X+70.0061*Y+84.3237; arthur poquet on 12 Jun 2020
I am also interrested with multiple intersections