ROUGH DRAFT authorea.com/48063
Main Data History
Export
Show Index Toggle 2 comments
  •  Quick Edit
  • Note: Closest point on ellipsoid

    Equations

    The equation for an ellipsoid with axises \(e_1, e_2\) and \(e_3\): \[E(x) := \frac{x_1^2}{e_1^2} + \frac{x_2^2}{e_2^2} + \frac{x_3^2}{e_3^2} - 1 = 0\] Let \(y = \begin{bmatrix} y_1 & y_2 & y_3 \end{bmatrix}^T\) be a point outside the ellipsoid. For any point \(x = \begin{bmatrix} x_1 & x_2 & x_3 \end{bmatrix}^{T}\) on the ellipsoid we have1: \[y - x = \frac{1}{2} t \nabla E(x),\] or equivalently in component form, \[y_k = \left( 1 + \frac{t}{e_k^2} \right) x_k, \qquad k=1,2,3.\] The values \(x_k\) can be expressed as a function of \(t\): \[x_k(t):=\frac{e_k^2}{t+e_k^2} y_k, \qquad k=1,2,3,\] and their derivatives are given by, \[x_k'(t) = -\frac{e_k^2}{(t+e_k^2)^2} y_k = -\frac{x_k(t)}{t+e_k^2}, \qquad k=1,2,3.\] The equation of the ellipsoid is reduced to a scalar equation in \(t\): \[E(t)=\frac{x_1^2(t)}{e_1^2} + \frac{x_2^2(t)}{e_2^2} + \frac{x_3^2(t)}{e_3^2} - 1.\] The derivative \(E'(t)\) is simply: \[E'(t)=\frac{2 x_1(t) x_1'(t)}{e_1^2} + \frac{2 x_2(t) x_2'(t)}{e_2^2} + \frac{2 x_3(t) x_3'(t)}{e_3^2}.\]


    1. http://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf

    Analytical solution

    When \(e_1=e_2=e_3=R\) we can solve \(E(t)=0\) analytically. \[y_1^2+y_2^2+y_3^2-(t+R^2)^2=0\] \[t = - R^2 + \sqrt{y_1^2+y_2^2+y_3^2}\] Note: We consider the positive \(t\). Negative \(t\) corresponds to the point \(\tilde{x}\) on the other side of the ellipsoid.

    Numerical algorithm

    1. Find a starting guess \(t_0\).

    2. Iterate a few Newton Raphson steps: \[t_{n+1} = t_n - \frac{E(t_n)}{E'(t_n)}, \qquad n=0,1,\dots.\]

    3. Until convergence \(|\Delta t| < TOL\) and \(|E(t)| < RTOL\).

    Code

    int EllipsoidSolver(Vec3 E, Vec3 y,
                        float& error, float& resid, Vec3& x, Vec3& n,
                        float tol=1e-4, float rtol=1e-4,
                        int maxiter=10, bool useX=false)
    {
        const float Rmin = 1e-4;
        const Vec3 One = Vec3(1.0f);
    
        const Vec3 y2 = y*y;
        const Vec3 E2 = E*E;
        const float yn = norm(y);
        const float R = min(E);
    
        float t;
        int iter = 0;
    
        // Check if y is too close to center of the ellipsoid.
        if (yn < Rmin)
            return -maxiter;
    
        // Start guess.
        if (useX)
            t = norm(y-x)/norm(x/E2);
        else
            t = -R + norm(y);
    
        while (iter < maxiter)
        {
            x = (E2/(Vec3(t)+E2))*y;
            Vec3 xp = -x/(Vec3(t)+E2);
            float F = dot(One, x*x/E2) - 1.0f;
            float Fp = dot(One, 2.0f*x*xp/E2);
            float dt = F/Fp;
            t = t - dt;
            iter = iter + 1;
            error = std::fabs(dt);
            resid = std::fabs(F);
            if (error < tol && resid < rtol)
            {
                n = x/E2;
                float rn = 1.0f/norm(n);
                n = rn*n;
                return iter;
            }
        }
        
        return -maxiter;
    }