3

Goal:

I am writing software to visualize 3-D objects in Python, using libraries such as sympy, numpy, and matplotlib.pyplot. I would like to fit the best surface to a small number of points. This is why I want to find the smallest distance between a point and a quadric surface

============================================================================== $$Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz + J = 0$$

==============================================================================

Questions:

1) How can I find a solution to this equation, given particular coefficients $A, B, ..., J$ ? This is for the purpose of solving question 2, which is the key question:

2) How can I find the smallest Euclidean distance between that surface and a point? I have access to a computer; I'm using python2 and sympy at the moment.

3) Is there software I can use to visualize such a surface?

4) Are these surfaces known as "manifolds" in proper mathematical terms? From brief reading, it sounds like a manifold is a more general mathematical object than the object the word "surface" I'm using denotes

I've looked here, here, here, here, and here.

Nathan majicvr.com
  • 111
platonicity
  • 51
  • 1
    You search for a line which passes the point of interest and which has the same direction as the gradient of the point on the surface which it also passes. – mathreadler Jun 29 '18 at 19:31
  • Some quadrics are manifolds, some aren’t. – amd Jun 29 '18 at 19:46
  • @amd I'm assuming a hyperboloid of 2 sheets isn't a manifold because of the 2 sheets – platonicity Jun 29 '18 at 22:34
  • @mathreadler Thanks! Is there an efficient (ideally simple) way to compute that? I don't even know how to easily parametrize the surface, so I'd need to spend quite a bit of time with paper to make progress – platonicity Jun 29 '18 at 22:37
  • 1
    Cones are the real problem. There’s no neighborhood of the vertex that looks like $\mathbb R^2$. – amd Jun 29 '18 at 23:13

2 Answers2

2

This solved the problem I cared about. Hope it's also useful to others. For convenience, once again the relevant equation is:

========================================================================= $$ Ax^2 + By^2 + Cz^2 + Dxy + Exz + Fyz + Gx + Hy + Iz + J = 0 $$

=========================================================================

python2 code:

from __future__ import division
from sympy import *

if __name__=='__main__':
    x = Symbol('x'); y = Symbol('y'); z = Symbol('z'); lambd = Symbol('lambd') 
    # lambda is a keyword; I didn't want to rename
    A = 0
    B = 0
    C = 0
    D = 0
    E = 0
    F = 0
    G = 1
    H = 1
    I = 1
    J = 0
    #x_e means x_external
    x_e = 1; y_e = 1; z_e = 1.1
    var_list            = [x,y,z,lambd]
    eq_list             = []
    eq_list.append(2*x - 2*x_e + 2*lambd*A*x + lambd*D*y + lambd*E*z + lambd*G)
    eq_list.append(2*y - 2*y_e + 2*lambd*B*y + lambd*D*x + lambd*F*z + lambd*H)
    eq_list.append(2*z - 2*z_e + 2*lambd*C*z + lambd*E*x + lambd*F*y + lambd*I)
    eq_list.append(   A*(x**2)  + B*(y**2)  + C*(y**2)
                    + D*x*y     + E*x*z     + F*y*z
                    + G*x       + H*y       + I*z
                    + J)
    print solve(eq_list, var_list, dict=True)

===================================================================

How to reuse my code for your problem:

x_e, y_e, and z_e should be the values of the relevant "external" point you care about. Similar for the values of A, B, ..., J. But this system of equations should work for solving the min dist between a general quadric and a point.

reference: Finding shortest distance between a point and a surface

===================================================================

Visualization:

I used GNU Octave to visualize. It's the free version of MATLAB, which you can use if you've paid for it. See the code from this post:

gv = linspace(-30,30,50); % adjust for appropriate domain
[xx yy zz]=meshgrid(gv, gv, gv);
F = A*xx.*xx + B*yy.*yy + C*zz.*zz+ ... etc

figure
isosurface(xx, yy, zz, F, 0)

For me, visualization is surprisingly unhelpful for checking how well my code is working. I really needed an extra "distance between point and surface" calculator to debug properly, which is not something easy to visualize, even with GNU Octave. Still, very cool to see the quadric graphed. Hope it helps someone else more

platonicity
  • 51
  • 1
    If I’m reading your code correctly, all of the coefficients of the second-degree terms are zero, so what you’re really doing is computing the closest point on a plane, and one that passes through the origin, no less. That’s a much simpler problem that can be solved via a straightforward projection. For a general quadric, the system of equations that you have might have two solutions: both the nearest and farthest points. – amd Jun 29 '18 at 23:28
  • @amd yeah, I jumped the gun. It's now trying to compute a real quadric, and it's taking awhile – Nathan majicvr.com Jun 29 '18 at 23:31
  • 1
    You can eliminate $\lambda$ from the calculations by using the equations that you get from $(x-x_e,y-y_e,z-z_e)\times\nabla f(x,y,z)=0$ instead. This condition says that the line is parallel to the normal, which for quadrics works just as well. – amd Jun 29 '18 at 23:36
  • @amd Perhaps more a coding question than a math question, but why does sympy often spit out values with capital I in them? I'm getting {x: -99.0000000000000, z: 1.11122448979592, y: 99.0*I, lambd: -1.01010101010101}. Also, appreciate all the help :) – platonicity Jun 29 '18 at 23:42
  • Haven’t used sympy myself, but I’d guess that means the solution to the system of equations involves complex numbers. I assume that there’s some way to tell it that you only want real solutions. – amd Jun 29 '18 at 23:44
  • @amd Forget asking you to debug my code. I think that's a waste of our time. Real question: is there an intuitive reason your perfectly good constraint outperforms the Lagrange Multipliers method I was using? I'll definitely try yours, but I did a calculation, and your method gives me a system of equations including $0 = (x-x_e)(2By + Dx + Fz +H) - (y-y_e)(2Ax + Dy + Ez + G) = 0$, which has both an $x^2$ term and a $y^2$ term, as well as many cross terms – platonicity Jun 30 '18 at 00:44
  • 1
    Off the top of my head, you’ve done some of sympy’s work for it by combining equations and reducing the number of variables. If you eliminate $\lambda$ from the Lagrange multiplier-generated stuff, you end up with similar-looking second-degree equations. You also end up with similar equations if you eliminate $\lambda$ from $(x,y,z)+\lambda\nabla f(x,y,z)=(x_e,y_e,z_e)$. – amd Jun 30 '18 at 02:05
1

Here is a sketch of a partial solution. Hopefully I can expand a bit later.

Let $$f(x,y,z) = Ax^2+By^2+Cz^2+Dxy+Exz+Fyz+Gx+Hy+Iz+J=0$$

The gradient is $$\nabla(f) = \begin{bmatrix}2xA+Dy+Ez+G\\2yB+Dx+Fz+H\\2zC+Ex+Fy+I\end{bmatrix}$$

now you want to find $c$ which minimizes $$\underset{c}{\min}\{\|c({\bf p}-[x,y,z]^T) -\nabla f(x,y,z)\|\}$$

Where $\bf p$ is the point outside which you want to find shortest line to. We are not interested in the $c$, but we are interested in the error vector of the solution. If the error is zero vector (and if $f(x,y,z)=0$) then we have found solution. So now we have a stopping criteria for an algorithm to solve the problem.

mathreadler
  • 25,824
  • Is there a fast way to compute this? I've ended up implementing coordinate descent that basically 1) initializes at a point 2) calculates distance between $p$ and that point, 3) checks points on the surface near $p$ by adding or subtracting $\Delta p$ 4) moves to the closer point and repeats steps 2-4 – Nathan majicvr.com Jul 04 '18 at 18:58
  • Also, please correct me if I'm wrong, but it sounds like there are no analytic solutions that take substantially less time to calculate. I hope I'm using the word "analytic" correctly; by that I mean solutions where I don't have to try any points on the surface and can just directly extract the answer – Nathan majicvr.com Jul 04 '18 at 19:01
  • 1
    This one is probably analytically solvable, but a general case I would not hope for it @frank. We want to find $[x,y,z]^T + k \nabla f(x,y,z) = {\bf p}$. Everything is actually linear, (except for the $f(x,y,z)=0$ constraint). – mathreadler Jul 04 '18 at 20:52