distanza dal punto dato a un'ellisse dato

Question

Ho un'ellisse, definita da Punto centrale, raggioX e raggioY, e ho un punto. Voglio trovare il punto sull'ellisse più vicino al punto specificato. Nell'illustrazione seguente, sarebbe S1.

Ora ho già il codice, ma c'è un errore logico da qualche parte in esso, e mi sembra di essere in grado di trovarlo. Ho rotto il problema verso il basso per il seguente esempio di codice:

#include <vector> 
#include <opencv2/core/core.hpp> 
#include <opencv2/highgui/highgui.hpp> 
#include <math.h> 

using namespace std; 

void dostuff(); 

int main() 
{ 
    dostuff(); 
    return 0; 
} 

typedef std::vector<cv::Point> vectorOfCvPoints; 

void dostuff() 
{ 

    const double ellipseCenterX = 250; 
    const double ellipseCenterY = 250; 
    const double ellipseRadiusX = 150; 
    const double ellipseRadiusY = 100; 

    vectorOfCvPoints datapoints; 

    for (int i = 0; i < 360; i+=5) 
    { 
     double angle = i/180.0 * CV_PI; 
     double x = ellipseRadiusX * cos(angle); 
     double y = ellipseRadiusY * sin(angle); 
     x *= 1.4; 
     y *= 1.4; 
     x += ellipseCenterX; 
     y += ellipseCenterY; 
     datapoints.push_back(cv::Point(x,y)); 
    } 

    cv::Mat drawing = cv::Mat::zeros(500, 500, CV_8UC1); 

    for (int i = 0; i < datapoints.size(); i++) 
    { 
     const cv::Point & curPoint = datapoints[i]; 
     const double curPointX = curPoint.x; 
     const double curPointY = curPoint.y * -1; //transform from image coordinates to geometric coordinates 

     double angleToEllipseCenter = atan2(curPointY - ellipseCenterY * -1, curPointX - ellipseCenterX); //ellipseCenterY * -1 for transformation to geometric coords (from image coords) 

     double nearestEllipseX = ellipseCenterX + ellipseRadiusX * cos(angleToEllipseCenter); 
     double nearestEllipseY = ellipseCenterY * -1 + ellipseRadiusY * sin(angleToEllipseCenter); //ellipseCenterY * -1 for transformation to geometric coords (from image coords) 


     cv::Point center(ellipseCenterX, ellipseCenterY); 
     cv::Size axes(ellipseRadiusX, ellipseRadiusY); 
     cv::ellipse(drawing, center, axes, 0, 0, 360, cv::Scalar(255)); 
     cv::line(drawing, curPoint, cv::Point(nearestEllipseX,nearestEllipseY*-1), cv::Scalar(180)); 

    } 
    cv::namedWindow("ellipse", CV_WINDOW_AUTOSIZE); 
    cv::imshow("ellipse", drawing); 
    cv::waitKey(0); 
} 

Produce la seguente immagine:

Si può vedere che trova in realtà punti "vicino" sul ellisse, ma non sono i punti "più vicini". Quello che intenzionalmente voglio è questa: (scusate il mio povero disegno)

sarebbe si misura le righe l'ultima immagine, che avrebbe attraversato il centro dell'ellisse, ma questo non è il caso per le linee nell'immagine precedente.
Spero che tu abbia la foto. Qualcuno può dirmi cosa sto facendo di sbagliato?

Answer 1

consideri un cerchio di delimitazione attorno al punto dato (c, d), che passa attraverso il punto più vicino sull'ellisse. Dal diagramma è chiaro che il punto più vicino è tale che una linea tracciata da esso al punto dato deve essere perpendicolare alla tangente condivisa dell'ellisse e del cerchio. Qualsiasi altro punto sarebbe fuori dal cerchio e quindi deve essere più lontano dal punto specificato.

Quindi il punto che si sta cercando è non l'intersezione tra la linea e l'ellisse, ma il punto (x, y) nel diagramma.

gradiente di tangente:

gradiente della linea:

Condizione di linee perpedicular - prodotto di gradienti = -1:

Quando riarrangiati e sostituito nell'equazione del vostro dell'ellisse ...

... questo vi darà due equazioni brutto quarto grado (4 ° grado polinomiale) in termini di x o y. AFAIK non esistono i metodi analitici (esatti algebrici) generali per risolverli. Potresti provare un metodo iterativo: cerca l'algoritmo iterativo di individuazione delle radici di Newton-Raphson.

Date un'occhiata a questo ottimo documento sul tema: http://www.spaceroots.org/documents/distance/distance-to-ellipse.pdf

Ci scusiamo per la risposta incompleta - Sono assolutamente colpa le leggi della matematica e della natura ...

EDIT: oops, mi sembra avere aeb turno modo errato nel diagramma xD

Answer 2

Hai solo bisogno di calcolare l'intersezione della linea [P1,P0] per l'eliporto che è S1.

Se la linea equeation è:

e l'equesion ellisse è:

rispetto ai valori di S1 saranno:

Ora solo bisogno di calcolare la distanza tra S1-P1, la formula (per A,B punti) è:

Answer 3

Ecco il codice tradotto a C# attuato da questo documento per risolvere per l'ellisse: http://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf

Nota che questo codice non è stato testato, se trovi qualche errore fammelo sapere.

//Pseudocode for robustly computing the closest ellipse point and distance to a query point. It 
    //is required that e0 >= e1 > 0, y0 >= 0, and y1 >= 0. 
    //e0,e1 = ellipse dimension 0 and 1, where 0 is greater and both are positive. 
    //y0,y1 = initial point on ellipse axis (center of ellipse is 0,0) 
    //x0,x1 = intersection point 

    double GetRoot (double r0 , double z0 , double z1 , double g) 
    { 
     double n0 = r0*z0; 
     double s0 = z1 - 1; 
     double s1 = (g < 0 ? 0 : Math.Sqrt(n0*n0+z1*z1) - 1) ; 
     double s = 0; 
     for (int i = 0; i < maxIter; ++i){ 
      s = (s0 + s1)/2 ; 
      if (s == s0 || s == s1) {break; } 
      double ratio0 = n0 /(s + r0); 
      double ratio1 = z1 /(s + 1); 
      g = ratio0*ratio0 + ratio1*ratio1 - 1 ; 
      if (g > 0) {s0 = s;} else if (g < 0) {s1 = s ;} else {break ;} 
     } 
     return s; 
    } 
    double DistancePointEllipse(double e0 , double e1 , double y0 , double y1 , out double x0 , out double x1) 
    { 
     double distance; 
     if (y1 > 0){ 
      if (y0 > 0){ 
       double z0 = y0/e0; 
       double z1 = y1/e1; 
       double g = z0*z0+z1*z1 - 1; 
       if (g != 0){ 
        double r0 = (e0/e1)*(e0/e1); 
        double sbar = GetRoot(r0 , z0 , z1 , g); 
        x0 = r0 * y0 /(sbar + r0); 
        x1 = y1 /(sbar + 1); 
        distance = Math.Sqrt((x0-y0)*(x0-y0) + (x1-y1)*(x1-y1)); 
        }else{ 
         x0 = y0; 
         x1 = y1; 
         distance = 0; 
        } 
       } 
       else // y0 == 0 
        x0 = 0 ; x1 = e1 ; distance = Math.Abs(y1 - e1); 
     }else{ // y1 == 0 
      double numer0 = e0*y0 , denom0 = e0*e0 - e1*e1; 
      if (numer0 < denom0){ 
        double xde0 = numer0/denom0; 
        x0 = e0*xde0 ; x1 = e1*Math.Sqrt(1 - xde0*xde0); 
        distance = Math.Sqrt((x0-y0)*(x0-y0) + x1*x1); 
       }else{ 
        x0 = e0; 
        x1 = 0; 
        distance = Math.Abs(y0 - e0); 
      } 
     } 
     return distance; 
    } 

Answer 4

Il seguente codice Python implementa le equazioni descritte in "Distance from a Point to an Ellipse" e usa il metodo di Newton per trovare le radici e da quel punto più vicino sull'ellisse al punto.

Sfortunatamente, come si può vedere dall'esempio, sembra essere preciso solo al di fuori dell'ellisse. All'interno dell'ellisse succedono cose strane.

from math import sin, cos, atan2, pi, fabs 


def ellipe_tan_dot(rx, ry, px, py, theta): 
    '''Dot product of the equation of the line formed by the point 
    with another point on the ellipse's boundary and the tangent of the ellipse 
    at that point on the boundary. 
    ''' 
    return ((rx ** 2 - ry ** 2) * cos(theta) * sin(theta) - 
      px * rx * sin(theta) + py * ry * cos(theta)) 


def ellipe_tan_dot_derivative(rx, ry, px, py, theta): 
    '''The derivative of ellipe_tan_dot. 
    ''' 
    return ((rx ** 2 - ry ** 2) * (cos(theta) ** 2 - sin(theta) ** 2) - 
      px * rx * cos(theta) - py * ry * sin(theta)) 


def estimate_distance(x, y, rx, ry, x0=0, y0=0, angle=0, error=1e-5): 
    '''Given a point (x, y), and an ellipse with major - minor axis (rx, ry), 
    its center at (x0, y0), and with a counter clockwise rotation of 
    `angle` degrees, will return the distance between the ellipse and the 
    closest point on the ellipses boundary. 
    ''' 
    x -= x0 
    y -= y0 
    if angle: 
     # rotate the points onto an ellipse whose rx, and ry lay on the x, y 
     # axis 
     angle = -pi/180. * angle 
     x, y = x * cos(angle) - y * sin(angle), x * sin(angle) + y * cos(angle) 

    theta = atan2(rx * y, ry * x) 
    while fabs(ellipe_tan_dot(rx, ry, x, y, theta)) > error: 
     theta -= ellipe_tan_dot(
      rx, ry, x, y, theta)/\ 
      ellipe_tan_dot_derivative(rx, ry, x, y, theta) 

    px, py = rx * cos(theta), ry * sin(theta) 
    return ((x - px) ** 2 + (y - py) ** 2) ** .5 

Ecco un esempio:

rx, ry = 12, 35 # major, minor ellipse axis 
x0 = y0 = 50 # center point of the ellipse 
angle = 45 # ellipse's rotation counter clockwise 
sx, sy = s = 100, 100 # size of the canvas background 

dist = np.zeros(s) 
for x in range(sx): 
    for y in range(sy): 
     dist[x, y] = estimate_distance(x, y, rx, ry, x0, y0, angle) 

plt.imshow(dist.T, extent=(0, sx, 0, sy), origin="lower") 
plt.colorbar() 
ax = plt.gca() 
ellipse = Ellipse(xy=(x0, y0), width=2 * rx, height=2 * ry, angle=angle, 
        edgecolor='r', fc='None', linestyle='dashed') 
ax.add_patch(ellipse) 
plt.show() 

che genera un'ellisse e la distanza dal confine dell'ellisse come una mappa di calore. Come si può vedere, al confine la distanza è zero (blu profondo).

Answer 5

Esiste un metodo numerico relativamente semplice con una convergenza migliore rispetto al metodo Newton.Ho un post sul blog sul perché funziona http://wet-robots.ghost.io/simple-method-for-distance-to-ellipse/

def solve(semi_major, semi_minor, p): 
    px = abs(p[0]) 
    py = abs(p[1]) 

    t = math.pi/4 

    a = semi_major 
    b = semi_minor 

    for x in range(0, 3): 
     x = a * math.cos(t) 
     y = b * math.sin(t) 

     ex = (a*a - b*b) * math.cos(t)**3/a 
     ey = (b*b - a*a) * math.sin(t)**3/b 

     rx = x - ex 
     ry = y - ey 

     qx = px - ex 
     qy = py - ey 

     r = math.hypot(ry, rx) 
     q = math.hypot(qy, qx) 

     delta_c = r * math.asin((rx*qy - ry*qx)/(r*q)) 
     delta_t = delta_c/math.sqrt(a*a + b*b - x*x - y*y) 

     t += delta_t 
     t = min(math.pi/2, max(0, t)) 

    return (math.copysign(x, p[0]), math.copysign(y, p[1]))