How can I get a fast estimate for the distance between a point and a bicubic spline surface in Python?

Question

How can I get a fast estimate for the distance between a point and a bicubic spline surface in Python? Is there an existing solution that I could leverage in SciPy, NumPy, or some other package?

I've got surface defined by a bicubic interpolation as this:

import numpy as np
import scipy.interpolate

# Define regular grid surface
xmin,xmax,ymin,ymax = 25, 125, -50, 50
x = np.linspace(xmin,xmax, 201)
y = np.linspace(ymin,ymax, 201)
xx, yy = np.meshgrid(x, y)
z_ideal = ( xx**2 + yy**2 ) / 400
z_ideal += z_ideal + np.random.uniform(-0.5, 0.5, z_ideal.shape)
s_ideal = scipy.interpolate.interp2d(x, y, z_ideal, kind='cubic')   

and I've got some measured points of that surface:

# Fake some measured points on the surface
z_measured = z_ideal + np.random.uniform(-0.1, 0.1, z_ideal.shape)
s_measured = scipy.interpolate.interp2d(x, y, z_measured, kind='cubic')
p_x = np.random.uniform(xmin,xmax,10000)
p_y = np.random.uniform(ymin,ymax,10000)
p_z = s_measured( p_x, p_y )

I want to find the closest point on the surface s_ideal to each point in p. A general case could have multiple solutions for wildly varying splines, so I'm limiting the problem to surfaces that are known to have only one solution in the vicinity of the point's projection along z. This isn't a tiny number of measurement or surface definition points, so I'd like to optimize the speed even at the expense of accuracy to maybe 1E-5.

The method that comes to mind is to use a gradient descent approach and do something like for each measurement point p:

1. Use pt = [p_x, p_y, p_z] as the initial test point, where p_z = s_ideal(pt)
2. Calculate the slope (tangent) vector m = [ m_x, m_y ] at pt
3. Calculate the vector r from pt to p: r = p - pt
4. If the angle theta between r and m is within some threshold of 90 deg, then pt is the final point.
5. Otherwise, update pt as:

r_len = numpy.linalg.norm(r)
dx = r_len * m_x
dy = r_len * m_y
if theta > 90:
    pt = [ p_x + dx, p_y + dy ]
else:
    pt = [ p_x - dx, p_y - dy ]

I found this suggesting a method could produce fast results to a very high accuracy for the 1D case, but it's for a single dimension and might be too hard for me to convert to two.

Answer 1

The question seeks to minimize the Euclidian distance between a three-dimensional surface S(x,y,z) and another point x0,y0,z0. The surface is defined on a rectangular (x,y) mesh, where z(x,y) = f(x,y) + random_noise(x,y). The introduction of noise to the "ideal" surface adds considerable complexity to the problem, as it requires the surface to be interpolated using a 2-dimensional third-order spline.

It's not understood why the introduction of noise to the ideal surface is actually necessary. If the ideal surface were truly ideal, it should by understood well enough that a true polynomial fit in x and y could be determined, if not analytically, at least empirically. If the random noise were to simulate an actual measurement, then one only needs to record the measurement enough times until the noise is averaged out to zero. Similarly, the use of signal filtering can help eliminate noise and reveal the signal's true behavior.

To find the closest point on the surface to another point, the distance equation and its derivatives must be used. If the surface can truly only be described using a basis of splines, then one must reconstruct the spline representation and find its derivatives, which is non-trivial. Alternatively, the surface could be evaluated using a fine mesh, but here, one quickly runs into memory issues, which is why interpolation was used in the first place.

However, if we can agree that the surface can be defined using a simple expression in x and y, then the minimization becomes trivial:

For purposes of minimization, it is more convenient to look at the square of the distance d^2(x,y) (z is just a function of x and y) between two points, D(x,y), as it eliminates the square root. To find the critical points of D(x,y), we take its partial derivatives w.r.t x and y and find their roots by setting = 0: d/dx D(x,y) = f1(x,y) = 0 and d/dy D(x,y) = f2(x,y)=0. This is a non-linear system of equations, for which we can solve using scipy.optimize.root. We need only pass root a guess (the projection of the pt of interest onto the surface) and the Jacobian of the system of equations.

import numpy as np
import scipy.interpolate
import scipy.optimize

# Define regular grid surface
xmin,xmax,ymin,ymax = 25, 125, -50, 50
x = np.linspace(xmin,xmax, 201)
y = np.linspace(ymin,ymax, 201)
xx, yy = np.meshgrid(x, y)
z_ideal = ( xx**2 + yy**2 ) / 400

# Fake some measured points on the surface
z_measured = z_ideal + np.random.uniform(-0.1, 0.1, z_ideal.shape)
s_measured = scipy.interpolate.interp2d(x, y, z_measured, kind='cubic')
p_x = np.random.uniform(xmin,xmax,10000)
p_y = np.random.uniform(ymin,ymax,10000)

# z_ideal function
def z(x):
    return (x[0] ** 2 + x[1] ** 2) / 400

# returns the system of equations
def f(x,pt):
    x0,y0,z0 = pt
    f1 = 2*(x[0] - x0) + (z(x)-z0)*x[0]/100
    f2 = 2*(x[1] - y0) + (z(x)-z0)*x[1]/100
    return [f1,f2]

# returns Jacobian of the system of equations
def jac(x, pt):
    x0,y0,z0 = pt
    return [[2*x[0]+1/100*(1/400*(z(x)+2*x[0]**2))-z0, x[0]*x[1]/2e4],
    [2*x[1]+1/100*(1/400*(z(x)+2*x[1]**2))-z0, x[0]*x[1]/2e4]]

def minimize_distance(pt):
    guess = [pt[0],pt[1]]
    return scipy.optimize.root(f,guess,jac=jac, args=pt)

# select a random point from the measured data
x0,y0 = p_x[30], p_y[30]
z0 = float(s_measured(x0,y0))

minimize_distance([x0,y0,z0])

Output:

    fjac: array([[-0.99419141, -0.1076264 ],
       [ 0.1076264 , -0.99419141]])
     fun: array([ -1.05033229e-08,  -2.63163477e-07])
 message: 'The solution converged.'
    nfev: 19
    njev: 2
     qtf: array([  2.80642738e-07,   2.13792093e-06])
       r: array([-2.63044477, -0.48260582, -2.33011149])
  status: 1
 success: True
       x: array([ 110.6726472 ,   39.28642206])