How do I implement the math function to calculate the ellipse normal that crosses a point into C++?

Question

To calculate the shortest distance from a point to the circumference of an ellipse, I have to calculate the normal/angle that passes through this point. See the following illustration:

The math to calculate this normal is explained here. In short:

Where E is the ellipse and L is the line between the point p and the circumference that creates the normal.

In C++, I have the following code fragment, where I need to insert this formula - solving it to variable t:

#include <cmath>

/// Calculate the normal of a line from ellipse circumference to point.
///
/// @param ew The width of the ellipse
/// @param eh The height of the ellipse
/// @param x The relative x position (relative to ellipse center)
/// @param y The relative y position (relative to ellipse center)
/// @return The normal in radians.
///
auto calculateNormal(float ew, float eh, float x, float y) -> double {
   // ???
   return std::atan(...);
}; 

How can I implement the formula in C++ to get the normal as shown?

Answer 1

What you currently have is not yet solved system of equations which as I originally taught from experience with similar but slightly simpler problem the algebraic solution leads to horrible stuff that would be too slow and inaccurate to compute (math solver throw about 2 pages of equations for roots which I throw away right away) so I would attack this with search approach instead.

If you are worrying about speed due to using goniometrics you can avoid that completely by simply computing normal as numeric derivation of the circumference points and just compare the slopes between normal (nx,ny) and found_point(x,y) - input_point(px,py).

The algo is like this:

1. search x in range <-a,+a>

2. compute y=f(x) and y1=f(x+epsilon)

   I simply rewrite the implicit ellipse equation into something like this:

   float ellipse_y(float rx,float ry,float x)   // y = f(x)
        {
        if (fabs(x)>rx) return 0.0
        return sqrt((rx*rx)-(x*x))*ry/rx;
        }
   

   of coarse the result should be +/- depending on the quadrant so if py<0 use negative values...

   The epsilon should be some small value but not too small, I used 0.001*rx where rx,ry are sizes of the ellipse half axises.

3. compute normal (nx,ny)

   so simply take two consequent points (x,y) and (x+epsilon,y1) substract them and rotate by 90deg by exchanging their coordinates and negate one of them. Once put together I got this:

   void ellipse_n(float rx,float ry,float &nx,float &ny,float x,float &y)   // x',y',y = f(x)
        {
        if (x<=-rx){ y=0.0; nx=-1.0; ny=0.0; return; }
        ny=x+(0.001*rx);                // epsilon
        if (ny>=+rx){ y=0.0; nx=+1.0; ny=0.0; return; }
        y=ellipse_y(rx,ry,x);           // first point
        nx=y-ellipse_y(rx,ry,ny);       // second point
        ny=ny-x;
    /*
        // normalize
        x=divide(1.0,sqrt((nx*nx)+(ny*ny)));
        nx*=x;
        ny*=x;
    */
        }
   

   The normalization is optional (I commented it out for speed as its not needed for the search itself).

4. compute error e for the search

   Simply the slopes (x-px,y-py) and (nx,ny) should be equal so:

   e=fabs(((y-py)*nx)-((x-px)*ny));
   

   The x search should minimize the e towards zero.

Do not forget to handle py<0 by negating y. Putting all togetherusing my approx search leads to:

//---------------------------------------------------------------------------
float ellipse_y(float rx,float ry,float x)  // y = f(x)
    {
    if (fabs(x)>rx) return 0.0;
    return sqrt((rx*rx)-(x*x))*ry/rx;
    }
//---------------------------------------------------------------------------
void ellipse_pn(float rx,float ry,float &nx,float &ny,float x,float &y) // x',y',y = f(x) if (py>=0)
    {
    if (x<=-rx){ y=0.0; nx=-1.0; ny=0.0; return; }
    ny=x+(0.001*rx);                // epsilon
    if (ny>=+rx){ y=0.0; nx=+1.0; ny=0.0; return; }
    y=ellipse_y(rx,ry,x);           // first point
    nx=y-ellipse_y(rx,ry,ny);       // second point
    ny=ny-x;
    }
//---------------------------------------------------------------------------
void ellipse_nn(float rx,float ry,float &nx,float &ny,float x,float &y) // x',y',y = f(x) if (py<=0)
    {
    if (x<=-rx){ y=0.0; nx=-1.0; ny=0.0; return; }
    ny=x+(0.001*rx);                // epsilon
    if (ny>=+rx){ y=0.0; nx=+1.0; ny=0.0; return; }
    y=-ellipse_y(rx,ry,x);          // first point
    nx=y+ellipse_y(rx,ry,ny);       // second point
    ny=ny-x;
    }
//---------------------------------------------------------------------------
void this_is_main_code()
    {
    float rx=0.95,ry=0.35;  // ellipse
    float px=-0.25,py=0.15; // input point
    float x,y,nx,ny;
    approx ax; double e;
    if (py>=0.0)
        {
        for (ax.init(-rx,+rx,0.25*rx,3,&e);!ax.done;ax.step())
            {
            x=ax.a;
            ellipse_pn(rx,ry,nx,ny,x,y);
            e=fabs(((y-py)*nx)-((x-px)*ny));
            }
        x=ax.aa; y=+ellipse_y(rx,ry,x);
        }
    else{
        for (ax.init(-rx,+rx,0.25*rx,3,&e);!ax.done;ax.step())
            {
            x=ax.a;
            ellipse_nn(rx,ry,nx,ny,x,y);
            e=fabs(((y-py)*nx)-((x-px)*ny));
            }
        x=ax.aa; y=-ellipse_y(rx,ry,x);
        }
    // here (x,y) is found solution and (nx,ny) normal
    }
//---------------------------------------------------------------------------

You still can adjust the search parameters (step and recursions) to have desired precision and speed. Also you can optimize by inlinening and removing heap thrashing and moving the range checks before the iterations...

Also you should adjust the normal computation for range <rx-epsilon,rx> as it truncates it to x=+rx (you could use negative epsilon for that).

[Ed1t1]

If I see it right the point with matching normal is also closest points so we can search directly that (which is even faster):

// search closest point
if (py>=0.0)
    {
    for (ax.init(-rx,+rx,0.25*rx,3,&e);!ax.done;ax.step())
        {
        x=ax.a;
        y=ellipse_y(rx,ry,x);
        x-=px; y-=py; e=(x*x)+(y*y);
        }
    x=ax.aa; y=+ellipse_y(rx,ry,x);
    }
else{
    for (ax.init(-rx,+rx,0.25*rx,3,&e);!ax.done;ax.step())
        {
        x=ax.a;
        y=-ellipse_y(rx,ry,x);
        x-=px; y-=py; e=(x*x)+(y*y);
        }
    x=ax.aa; y=-ellipse_y(rx,ry,x);
    }

Also I feel there still might be some better solution using graphic approach like rescale to circle solve for circle and then scale back to ellipse +/- some corrections however too lazy to try...