How to find min distance between a point and a surface

Question

I have equation of a surface (for example: x^2+y^2=400) and a point (for example: 1,2,3). I want to calculate the approx min. distance between that surface and a point.

I tried using following function for it:

def distance_to_surface(equation, point):  
    # Define a function that returns the gradient of the surface equation
    x, y, z = symbols('x y z')
    expr = eval(equation)
    dx = diff(expr, x)
    dy = diff(expr, y)
    dz = diff(expr, z)
    grad_expr = np.array([dx, dy, dz])
    grad_func = lambdify((x, y, z), grad_expr, modules=['numpy' , "sympy"])
    
    # Find a point on the surface by solving the equation for one of the variables
    # Here, we solve for z to find a point on the surface with the same x and y coordinates as the given point
    x0, y0, z0 = point
    z_func = lambdify((x, y), z - expr, 'numpy')
    z_surface = z_func(x0, y0)
    surface_point = np.array([x0, y0, z_surface])
    
    # Calculate the vector between the given point and the surface point
    v = point - surface_point
    
    # Calculate the normal vector to the surface at the surface point
    epsilon = 1e-6
    normal = np.array([
        grad_func(x0 + epsilon, y0) - grad_func(x0 - epsilon, y0),
        grad_func(x0, y0 + epsilon) - grad_func(x0, y0 - epsilon),
        -1
    ])
    normal = normal / np.linalg.norm(normal)
    
    # Calculate the distance between the given point and the surface
    distance = np.abs(np.dot(v, normal))
    
    return distance

Answer 1

Your example can be solved purely symbolically by SymPy with a Lagrange multiplier. We want to find the distance from the point (1, 2, 3) to the closest point on the cylinder x**2 + y**2 = 400.

In [86]: px, py, pz = 1, 2, 3  # The point

In [87]: from sympy import *

In [88]: x, y, z = symbols('x, y, z')

In [89]: l = symbols('lambda')  # Lagrange multiplier

In [90]: d2 = symbols('d^2')  # symbol for squared distance

In [91]: dist2 = (x - px)**2 + (y - py)**2 + (z - pz)**2 # squared distance

In [92]: eq = x**2 + y**2 - 400 # surface equation

In [93]: eqs = [d2 - dist2, L.diff(x), L.diff(y), L.diff(z), L.diff(l)]

In [94]: for eq in eqs: print(eq)
d^2 - (x - 1)**2 - (y - 2)**2 - (z - 3)**2
2*lambda*x + 2*x - 2
2*lambda*y + 2*y - 4
2*z - 6
x**2 + y**2 - 400

The solution:

In [95]: s1, s2 = solve(eqs, [x, y, z, l, d2], dict=True)

In [96]: print(s1)
{d^2: 405 - 40*sqrt(5), lambda: -1 + sqrt(5)/20, x: 4*sqrt(5), y: 8*sqrt(5), z: 3}

In [97]: print(s2)
{d^2: 40*sqrt(5) + 405, lambda: -1 - sqrt(5)/20, x: -4*sqrt(5), y: -8*sqrt(5), z: 3}

This gives the solutions for d^2 as 405 +- 40*sqrt(5). One of those is the squared distance to the nearest point and the other to the furthest point (furthest within the plane z=3).

This example is easy to solve symbolically because it is all polynomials. For more complicated surface equations you might want to use a numeric solver in which case just use the SciPy one rather than writing your own:

In [105]: from scipy.optimize import minimize, NonlinearConstraint

In [106]: f = lambdify([(x, y, z)], dist2)

In [107]: g = lambdify([(x, y, z)], eq)

In [108]: nlc = NonlinearConstraint(g, 0, 0)

In [114]: minimize(f, [1, 1, 1], constraints=nlc, method='SLSQP')
Out[114]: 
 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 315.55728090094
       x: [ 8.944e+00  1.789e+01  3.000e+00]
     nit: 17
     jac: [ 1.589e+01  3.178e+01  5.722e-05]
    nfev: 79
    njev: 17

The function value ~315 is approximately equal to the exact value computed before:

In [115]: N(405 - 40*sqrt(5))
Out[115]: 315.557280900008

Again that's the square of the distance so the actual distance is

In [116]: N(sqrt(405 - 40*sqrt(5)))
Out[116]: 17.7639320225002

That number makes sense because you have a point near the centre of a cylinder of radius 20.