2
$\begingroup$

I have a surface defined by $f(x,y)=ax^2+by^2+cxy+dx+ey+g$ and a 3D point $P=(x_1,y_1,z_1)$. How can I find the shortest Euclidean distance between P and the surface?

Perhaps I can use something like this: How to find distance between a point and 3D surface, solution to general quadric equation, and visualizing such a surface?

But how can I do this computation efficiently (on a computer) if I have a large number of such points P (eg., 1 million)?

$\endgroup$
4
  • $\begingroup$ Your surface is z = f(x,y) ? $\endgroup$ May 4, 2020 at 8:11
  • $\begingroup$ Do you need to find all point/paraboloid distances, or just the shortest ? $\endgroup$
    – user65203
    May 4, 2020 at 17:21
  • $\begingroup$ See this paper, among others. $\endgroup$
    – amd
    May 4, 2020 at 19:52
  • $\begingroup$ Just the shortest $\endgroup$ May 5, 2020 at 10:00

1 Answer 1

2
$\begingroup$

The surface is $$z(x, y) = a x^2 + b y^2 + c x y + d x + e y + f \tag{1}\label{AC1}$$ The axis of symmetry is parallel to the $z$ axis, and $z(x, y)$ has a stationary point where $(x, y)$ are at the axis of symmetry: $$\left\lbrace \begin{aligned} \displaystyle \frac{\partial z(x, y)}{\partial x} &= 0 \\ \displaystyle \frac{\partial z(x, y)}{\partial y} &= 0 \\ \end{aligned} \right. $$ This has exactly one solution, $$\left\lbrace \begin{aligned} \displaystyle x_0 &= \frac{2 b d - c e}{c^2 - 4 a b} \\ \displaystyle y_0 &= \frac{2 a e - c d}{c^2 - 4 a b} \\ \end{aligned}\right. \tag{2a}\label{AC2a}$$ Note that $$z_0 = z(x_0 , y_0) = f + \frac{a e^2 - c d e + b d^2}{c^2 - 4 a b} \tag{2b}\label{AC2b}$$ and that if $c^2 = 4 a b$, the surface is not a paraboloid, as it does not have an axis of symmetry.

To make the next steps clearer to follow, let's reparametrize the surface as $(u, v)$: $$\left\lbrace\begin{aligned} x(u, v) &= u \\ y(u, v) &= v \\ z(u, v) &= a u^2 + b v^2 + c u v + d u + e v + f \\ \end{aligned}\right. \tag{3}\label{AC3}$$ The distance squared between point $(x_i , y_i , z_i)$ and $(u, v) = (x(u,v), y(u,v), z(u,v))$ is $s(u, v)$, $$s(u, v) = (u - x_i)^2 + (v - y_i)^2 + (a u^2 + b v^2 + c u v + d u + e v + f - z_i)^2 \tag{4a}\label{AC4a}$$ i.e. $$\begin{array}{rl} s(x, y) & = (a^2) x^4 + (b^2) y^4 \\ ~ & + (2 a c) x^3 y + (2 b c) x y^3 \\ ~ & + (2 a d) x^3 + (2 b e) y^3 \\ ~ & + (c^2 + 2 a b) x^2 y^2 \\ ~ & + (2 a e + 2 c d) x^2 y \\ ~ & + (2 c e + 2 b d) x y^2 \\ ~ & + (1 + d^2 + 2 a (f - z_i)) x^2 \\ ~ & + (1 + e^2 + 2 b (f - z_i)) y^2 \\ ~ & + (2 d e + 2 c (f - z_i)) x y \\ ~ & + (2 d (f - z_i) - 2 x_i) x \\ ~ & + (2 e (f - z_i) - 2 y_i) y \\ ~ & + (x_i^2 + y_i^2 + (f - z_i)^2) \\ \end{array} \tag{4b}\label{AC4b}$$ Its partial derivatives define the gradient, $$\left\lbrace ~ \begin{array}{rl} \displaystyle \frac{\partial s(u,v)}{\partial u} = s_u(u,v) & = (4 a^2) u^3 ~ + (2 b c) v^3 \\ ~ & + (6 a c) u^2 v ~ + (2 c^2 + 4 a c) u v^2 \\ ~ & + (6 a d) u^2 ~ + (2 c e + 2 b d) v^2 \\ ~ & + (4 a e + 4 c d) u v \\ ~ & + (2 + 2 d^2 + 4 a f - 4 a z_i) u \\ ~ & + (2 d e + 2 c f - 2 c z_i) v \\ ~ & + (2 d f - 2 d z_i - 2 x_i) \\ \displaystyle \frac{\partial s(u,v)}{\partial v} = s_v(u,v) & = (2 a c) u^3 ~ + (4 b^2) v^3 \\ ~ & + (2 c^2 + 4 a b) u^2 v ~ + (6 b c) u v^2 \\ ~ & + (2 a e + 2 c d) u^2 ~ + (6 b e) v^2 \\ ~ & + (4 c e + 4 b d) u v \\ ~ & + (2 c f + 2 d e - 2 c z_i) u \\ ~ & + (2 + 2 e^2 + 4 b f - 4 b z_i) v \\ ~ & + (2 e f - 2 e z_i - 2 y_i) \\ \end{array}\right. \tag{5}\label{5}$$ Although the distance squared function has 15 coefficients (up to for $u^4$ and for $v^4$), it is a surprisingly smooth function. Starting at $u = x_i$ and $v = y_i$, gradient descent should find the true minimum.

If we call the point on the surface $(u_i, v_i) = (u + u_{\Delta}, v + v_{\Delta})$ closest to the point $(x_i , y_i, z_i)$, we have immediately limited the region of interest to $u_{\Delta}^2 + v_{\Delta}^2 \le s(u, v)$, because all other points on the surface are more distant from surface point $(u, v)$ than the point of interest $(x_i , y_i , z_i)$.

Finding the minimum distances for a million points $i$ should not take that long on current computers; a few seconds, I'd guess.

One way to speed up the calculation is to choose better starting points $(u, v)$ for the iterations (instead of simply $u = x_i$, $v = y_i$). A reasonable method would be to set up a regular rectangular grid, and find the $u, v$ coordinates for the closest point on the surface to the center of that cell. (You would still calculate the squared distance to both points on the surface, and start the gradient descent or other iterative method at the one with the smaller squared distance.) Note that you can easily set up a sentinel coordinate pair to indicate "unknown", and only find the coordinates when needed; this is important, if the points are clumped within the grid, with many grid cells without any points at all. Still, this method is only useful if each grid cell has many points inside it, on average.

When doing the gradient descent, remember that all but $x$ and $y$ are constants, and you can precalculate the coefficients to speed up the evaluation at each candidate point. (To evaluate the squared distance, and its partial derivatives, you'll need around 44 multiplications and around 32 additions or subtractions.)


It is possible to translate and rotate the coordinate system (around the $z$ axis, in the $x y$ plane) so that the surface is just $$z(x, y) = C_{20} x^2 + C_{02} y^2$$ and the squared distance between point $z(x, y)$ on the surface and point $(x_i , y_i , z_i)$ is $$s(x, y) = (x - x_i)^2 + (y - y_i)^2 + (C_{20} x^2 + C_{02} y^2 - z_i)^2$$ and its partial derivatives are $$\begin{aligned} s_x (x,y) &= (4 C_{20}^2) x^3 + (4 C_{20} C_{02}) x y^2 + (2 - 4 C_{20} z_i) x - 2 x_i \\ s_y (x,y) &= (4 C_{02}^2) y^3 + (4 C_{20} C_{02}) x^2 y + (2 - 4 C_{02} z_i) y - 2 y_i \\ \end{aligned}$$ For each candidate point $(x, y)$, you'll need to do seven multiplications and six additions or subtractions to get the squared distance; or a total of seventeen multiplications and fourteen additions or subtractions for the squared distance and the partial derivatives. This can be a significant speedup, if the gradient descent requires many steps to find the minimum.

If we use a translation by $(x_0, y_0, z_0)$ and a rotation of $\theta$ counterclockwise on the $x y$ plane, the original point $(x_i , y_i , z_i)$ is transformed to $$\left\lbrace ~ \begin{aligned} x_k &= x_0 + x_i \cos\theta - y_i \sin\theta \\ y_k &= y_0 + y_i \cos\theta + y_i \sin\theta \\ z_k &= z_0 + z_i \\ \end{aligned} \right. \iff \left\lbrace ~ \begin{aligned} x_i &= x_k \cos\theta + y_k \sin\theta - x_0 \\ y_i &= y_k \cos\theta - x_k \sin\theta - y_0 \\ z_i &= z_k - z_0 \\ \end{aligned} \right. $$

The surface in the transformed coordinates (subscript $k$) is $$z = C_{20} x^2 + C_{02} y_2$$ which in original coordinates corresponds to $$\begin{aligned} z + z_0 & = C_{20} \left( x_0 + x_i \cos\theta - y_i \sin\theta \right)^2 \\ & + C_{02} \left( y_0 + y_i \cos\theta + x_i \sin\theta \right)^2 \\ \end{aligned}$$ ie. $$\begin{aligned} z(x, y) & = ( C_{02} (\sin\theta)^2 + C_{20} (\cos\theta)^2 ) x^2 \\ ~ & + ( C_{20} (\sin\theta)^2 + C_{02} (\cos\theta)^2 ) y^2 \\ ~ & + ( 2 (C_{02} - C_{20}) \cos\theta\sin\theta ) x y \\ ~ & + ( 2 C_{20} x_0 \cos\theta - 2 C_{02} y_0 \sin\theta ) x \\ ~ & + ( 2 C_{02} y_0 \cos\theta - 2 C_{20} x_0 \sin\theta ) y \\ ~ & + ( C_{20} x_0^2 + C_{02} y_0^2 ) \\ \end{aligned}$$ Comparing the coefficients to the original form, we see that $$\begin{aligned} a &= C_{02} (\sin\theta)^2 + C_{20} (\cos\theta)^2 \\ b &= C_{20} (\sin\theta)^2 + C_{02} (\cos\theta)^2 \\ c &= ( 2 C_{02} - C_{20}) \cos\theta\sin\theta \\ d &= 2 C_{20} x_0 \cos\theta - 2 C_{02} y_0 \sin\theta \\ e &= 2 C_{02} y_0 \cos\theta - 2 C_{20} x_0 \sin\theta \\ f &= C_{20} x_0^2 + C_{02} y_0^2 \\ \end{aligned}$$

$\endgroup$
1
  • $\begingroup$ That's a brilliant answer! Thank you for your hard work! $\endgroup$ May 5, 2020 at 10:01

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .