If your surface is, as it seems, a hyperbolic paraboloid, you can parametrize a point s
on it as:
s = p0 + u * (p1 - p0) + v * (p3 - p0) + u * v * (p2 - p3 - p1 + p0)
Doing things this way, the line p0p3
has equation u = 0
, p1p2
is u = 1
, p0p1
is v = 0
and p2p3
is v = 1
. I haven't been able to figure out a way of coming up with an analytical expression for the closest point to a point p
, but scipy.optimize
can do the job numerically for you:
import numpy as np
from scipy.optimize import minimize
p0 = np.array([1.15, 0.62, -1.01])
p1 = np.array([1.74, 0.86, -0.88])
p2 = np.array([1.79, 0.40, -1.46])
p3 = np.array([0.91, 0.79, -1.84])
p = np.array([1.17, 0.94, -1.52])
def fun(x, p0, p1, p2, p3, p):
u, v = x
s = u*(p1-p0) + v*(p3-p0) + u*v*(p2-p3-p1+p0) + p0
return np.linalg.norm(p - s)
>>> minimize(fun, (0.5, 0.5), (p0, p1, p2, p3, p))
status: 0
success: True
njev: 9
nfev: 36
fun: 0.24148810420527048
x: array([ 0.16446403, 0.68196253])
message: 'Optimization terminated successfully.'
hess: array([[ 0.38032445, 0.15919791],
[ 0.15919791, 0.44908365]])
jac: array([ -1.27032399e-06, 6.74091280e-06])
The return of minimize
is an object, you can access the values through its attributes:
>>> res = minimize(fun, (0.5, 0.5), (p0, p1, p2, p3, p))
>>> res.x # u and v coordinates of the nearest point
array([ 0.16446403, 0.68196253])
>>> res.fun # minimum distance
0.24148810420527048
Some pointers on how to go about finding a solution without scipy... The vector joining the point s
parametrized as above with a generic point p
is p-s
. TO find out the closest point you can go two different ways that give the same result:
- Compute the length of that vector,
(p-s)**2
, take its derivatives w.r.t. u
and v
and equate them to zero.
- Compute two vectors tangent to the hypar at
s
, which can be found as ds/du
and ds/dv
, and impose that their inner product with p-s
be zero.
If you work these out, you'll end up with two equations that would require a lot of manipulation to arrive at something like a third or fourth degree equation for either u
or v
, so there is no exact analytical solution, although you could solve that numerically with numpy only. An easier option is to work out those equations until you get these two equations, where a = p1-p0
, b = p3-p0
, c = p2-p3-p1+p0
, s_ = s-p0
, p_ = p-p0
:
u = (p_ - v*b)*(a + v*c) / (a + v*c)**2
v = (p_ - u*a)*(b + u*c) / (b + u*c)**2
You can't come up with an analytical solution for this easily, but you can hope that if you use those two relations to iterate a trial solution, it will converge. For your test case it does work:
def solve_uv(x0, p0, p1, p2, p3, p, tol=1e-6, niter=100):
a = p1 - p0
b = p3 - p0
c = p2 - p3 - p1 + p0
p_ = p - p0
u, v = x0
error = None
while niter and (error is None or error > tol):
niter -= 1
u_ = np.dot(p_ - v*b, a + v*c) / np.dot(a + v*c, a + v*c)
v_ = np.dot(p_ - u*a, b + u*c) / np.dot(b + u*c, b + u*c)
error = np.linalg.norm([u - u_, v - v_])
u, v = u_, v_
return np.array([u, v]), np.linalg.norm(u*a + v*b +u*v*c + p0 - p)
>>> solve_uv([0.5, 0.5], p0, p1, p2, p3, p)
(array([ 0.16446338, 0.68195998]), 0.2414881041967159)
I don't think this is guaranteed to converge, but for this particular case it seems to work pretty fine, and only needs 15 iterations to get to the requested tolerance.