Intersection point between line and surface in 3D python

Question

I want to be able to find the intersection between a line and a three-dimensional surface.

Mathematically, I have done this by taking the following steps:

Define the (x, y, z) coordinates of the line in a parametric manner. e.g. (x, y, z) = (1+t, 2+3t, 1-t)
Define the surface as a function. e.g. z = f(x, y)
Substitute the values of x, y, and z from the line into the surface function.
By solving, I would be able to get the intersection of the surface and the line

I want to know if there is a method for doing this in Python. I am also open to suggestions on more simple ways to solving for the intersection.

Welcome to Stack Overflow! Do you have any code which you have tried? Any code would be much appreciated and would allow the community to answer your question faster. — Jacob Lee
By "wasn't able to do it", do you mean you could not format your code into your question or that you were unable to come up with a solution? Even if it doesn't work fully, code examples help as an outline for others. Otherwise, they are further back than you are in their understanding of the issue. — Jacob Lee

cknoll cknoll · Accepted Answer · 2021-02-22T20:23:30

You can use the following code:

import numpy as np
import scipy as sc
import scipy.optimize
from matplotlib import pyplot as plt

def f(x, y):
    """ Function of the surface"""
    # example equation
    z = x**2 + y**2 -10
    return z

p0 = np.array([1, 2, 1]) # starting point for the line
direction = np.array( [1, 3, -1]) # direction vector

def line_func(t):
    """Function of the straight line.
    :param t:     curve-parameter of the line

    :returns      xyz-value as array"""
    return p0 + t*direction

def target_func(t):
    """Function that will be minimized by fmin
    :param t:      curve parameter of the straight line

    :returns:      (z_line(t) - z_surface(t))**2 – this is zero
                   at intersection points"""
    p_line = line_func(t)
    z_surface = f(*p_line[:2])
    return np.sum((p_line[2] - z_surface)**2)

t_opt = sc.optimize.fmin(target_func, x0=-10)
intersection_point = line_func(t_opt)

The main idea is to reformulate the algebraic equation point_of_line = point_of_surface (condition for intersection) into a minimization problem: |point_of_line - point_of_surface| → min. Due to the representation of the surface as z_surface = f(x, y) it is convenient to calculate the distance for a given t-value only on basis of the z-values. This is done in target_func(t). And then the optimal t-value is found by fmin.

The correctness and plausibility of the result can be checked with some plotting:


from mpl_toolkits.mplot3d import Axes3D
ax = plt.subplot(projection='3d')

X = np.linspace(-5, 5, 10)
Y = np.linspace(-5, 5, 10)
tt = np.linspace(-5, 5, 100)

XX, YY = np.meshgrid(X, Y)
ZZ = f(XX, YY)

ax.plot_wireframe(XX, YY, ZZ, zorder=0)

LL = np.array([line_func(t) for t in tt])
ax.plot(*LL.T, color="orange", zorder=10)
ax.plot([x], [y], [z], "o", color="red", ms=10, zorder=20)

Note that this combination of wire frame and line plots does not handle well, which part of the orange line should be below the blue wire lines of the surface.

Also note, that for this type of problem there might be any number of solutions from 0 up to +∞. This depends on the actual surface. fmin finds an local optimum, this might be a global optimum with target_func(t_opt)=0 or it might not. Changing the initial guess x0 might change which local optimum fmin finds.

Intersection point between line and surface in 3D python

1 Answers