From: Casey Muratori
Distance from plane to axis-aligned ellipsoid  
2002-10-24 01:45

 I'm new to this list and I noticed the tail end of a thread about the
 distance from a plane to an ellipse.  It appeared as if the person's
 question hadn't been answered, or perhaps it hadn't be answered
 algebraically, it was unclear - people were talking about transforms and
 what you can and can't do in non-uniformly scaled spaces and so on, which
 seemed somewhat irrelevant to the basic question that he was asking.  So
 at the risk of redundency if one has been posted before (since I didn't
 see the whole thread), here is a computational solution.
 
 A plane is
 
     p(x, y, z) = p_x x + p_y y + p_z z + p_d = 0
 
 An ellipsoid is
 
     e(x, y, z) = e_x x^2 + e_y y^2 + p_z z^2 - e_r^2 = 0
 
 I've put the ellipsoid at the origin to make the equations slightly
 shorter to type, since obviously the problem moves to the origin
 trivially regardless of where it is in space (and it does not require
 rotation, since the problem statement was for an axis-aligned ellipsoid).
 
 Using LaGrange multipliers, you get
 
     grad e(x, y, z) = l grad p(x, y, z)
 
 which evaluates to
 
     | 2 e_x x |     | p_x |
     | 2 e_y y | = l | p_y |
     | 2 e_z z |     | p_z |
 
 solving for the point, you get
 
          l p_x
      x = -----
          2 e_x
 
          l p_y
      y = -----
          2 e_y
 
          l p_z
      z = -----
          2 e_z
 
 The only unknown is l, the LaGrange multiplier.  So you just plug in the
 point to the original equation of the object on which you want the point
 (ie., the plane or the ellipsoid), and you have your closest point on that
 object to the other one.
 
 In the case of the plane, that is
 
      p_x x + p_y y + p_z z + p_d = 0                   Plane equation
 
      l_p   p_x^2   p_y^2   p_z^z
      --- ( ----- + ----- + ----- ) + p_d = 0           Substitute x, y, z
       2    e_x     e_y     e_z
 
                     -2 p_d
      l_p = -------------------------                   Solve for lambda
              p_x^2   p_y^2   p_z^z
            ( ----- + ----- + ----- )
               e_x     e_y     e_z
 
 In the case of the ellipse, that is
 
      e_x x^2 + e_y y^2 + e_z z^2 - e_r^2 = 0           Ellipsoid equation
 
      l_e^2   p_x^2   p_y^2   p_z^z
      ----- ( ----- + ----- + ----- ) - e_r^2 = 0       Substitute x, y, z
        4      e_x     e_y     e_z
 
                                      1
      l_e = +/- 2 e_r sqrt(-------------------------)   Solve for lambda
                             p_x^2   p_y^2   p_z^z
                           ( ----- + ----- + ----- )
                              e_x     e_y     e_z
 
 Now you can just plug l_p or l_e in for l in the equation of the point,
 and out pops the closest point on that object.  Note that since an
 ellipsoid has two extremum, there is a +/- in the solution for l_e, which
 means you will get both the closest AND the furthest point from the plane,
 so you're kind of getting one point "for free" (for some definition of
 free).
 
 You can pull out the common term in the denominators of both equations and
 this ends up being pretty easy to compute in code:
 
     float c = (p_x*p_x / e_x) + (p_y*p_y / e_y) + (p_z*p_z / e_z);
     float c_inv = 1.0f / c;
 
     float l_p = (-2 * p_d) * c_inv;
     float l_e = SquareRoot(4 * e_r*e_r * c_inv);
 
     float x_term = p_x / (2 * e_x);
     float y_term = p_y / (2 * e_y);
     float z_term = p_z / (2 * e_z);
 
     float PointOnPlane[3];
     PointOnPlane[0] = l_p * x_term;
     PointOnPlane[1] = l_p * y_term;
     PointOnPlane[2] = l_p * z_term;
 
     float PointOnEllipse[2][3];
     PointOnEllipse[0][0] = l_e * x_term;
     PointOnEllipse[0][1] = l_e * y_term;
     PointOnEllipse[0][2] = l_e * z_term;
     PointOnEllipse[1][0] = -l_e * x_term;
     PointOnEllipse[1][1] = -l_e * y_term;
     PointOnEllipse[1][2] = -l_e * z_term;
 
 That is the "complete" solution to the problem, in the sense that in those
 equations you now have all the information necessary to do any proximity
 oriented operation you might want between a plane and an ellipsoid.
 Since I missed too much of the thread to know what was actually being
 attempted with this distance formulation, I'm not able to make
 recommendations about what should be going on, so I'll simply leave it at
 this and hope that whoever asked the question can take it from there and
 construct what he needs from the equations.
 
 But it should be noted that if the only thing you wanted this for was to
 compute a measure of separation between the two objects, then I would
 recommend solving only for l_p and ignoring l_e entirely.  That way
 all you do is plug-in l_p, get x, y, and z, then plug those into the
 ellipsoid equation, which "is its own distance function", since it's in
 implicit form.  Which is to say:
 
     // / -> 7
     // * -> 13
     // + -> 5
     // - -> 1
     //      26
 
     float c = (p_x*p_x / e_x) + (p_y*p_y / e_y) + (p_z*p_z / e_z);
 
     float l_p = -p_d / c;
     float x = (l_p * p_x) / e_x;
     float y = (l_p * p_y) / e_y;
     float z = (l_p * p_z) / e_z;
 
     float d = e_x*x*x + e_y*y*y + e_z*z*z - e_r*e_r
 
 And hey, if this is a rigid bounding volume, you could store 1/e_x,
 1/e_y, and 1/e_z along with the bounding volume, and you'd only have
 a single divide... you'd just have:
 
     // / -> 1
     // * -> 19
     // + -> 5
     // - -> 1
     //      26
 
     float c = p_x*p_x*inv_e_x + p_y*p_y*inv_e_y + p_z*p_z*inv_e_z;
 
     float l_p = -p_d / c;
     float x = l_p * p_x * inv_e_x;
     float y = l_p * p_y * inv_e_y;
     float z = l_p * p_z * inv_e_z;
 
     float d = e_x*x*x + e_y*y*y + e_z*z*z - e_r*e_r
 
 So, 26 ops ain't bad.  The divide is not great, but I'm not sure if it can
 be eliminated.  I would have to give it more thought.
 
 Hope that helped.  Apologies in advance for any typos, it is really
 difficult to transcribe even simple math into ASCII.
 - Casey