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

I have two surfaces, and I would like to plot the line in three-dimensional space where the two surfaces intersect.

Answer 1

[MathWorks Support Team]

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.