Abstand von gegebenem Punkt zu gegebener Ellipse

Question

Ich habe eine Ellipse, definiert durch Mittelpunkt, RadiusX und RadiusY, und ich habe einen Punkt. Ich möchte den Punkt auf der Ellipse finden, der dem gegebenen Punkt am nächsten ist. In der Abbildung unten wäre das S1.

Jetzt habe ich bereits Code, aber es ist ein logischer Fehler irgendwo drin, und ich scheine nicht in der Lage sein, es zu finden. Ich brach das Problem auf das folgende Codebeispiel:

#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); 
} 

Es ergibt folgendes Bild:

Sie können sehen, dass es „in der Nähe“ Punkte auf der Ellipse tatsächlich findet, aber es sind nicht die "nächsten" Punkte. Was ich absichtlich will, ist dies: (entschuldigen Sie meine schlechte Zeichnung)

würden Sie weit die Linien im letzten Bild, sie das Zentrum der Ellipse überqueren würden, aber dies ist nicht der Fall für die Linien im vorherigen Bild.
Ich hoffe, Sie bekommen das Bild. Kann mir jemand sagen, was ich falsch mache?

Answer 1

einen Begrenzungskreis Betrachten um den gegebenen Punkt (c, d), der den nächsten Punkt durchläuft auf der Ellipse. Aus dem Diagramm wird klar, dass der nächste Punkt so ist, dass eine Linie, die von ihm zu dem gegebenen Punkt gezogen wird, senkrecht zu der geteilten Tangente der Ellipse und des Kreises sein muss. Alle anderen Punkte befinden sich außerhalb des Kreises und müssen daher weiter von dem gegebenen Punkt entfernt sein.

So ist der Punkt, den Sie für Ihre Suche ist nicht der Schnittpunkt zwischen der Leitung und der Ellipse, aber der Punkt (x, y) in dem Diagramm.

Gradient der Tangenten:

Gradient der Linie:

Zustand für perpedicular Linien - Produkt von Gradienten = -1:

Wenn neu geordnet und in die Gleichung der Ellipse ersetzt ...

... diese zwei böse quartic (4. Grades Polynom) Gleichungen geben in Bezug auf entweder x oder y. AFAIK gibt es keine allgemeine analytische (genaue algebraische) Methoden, um sie zu lösen. Sie könnten eine iterative Methode ausprobieren - sehen Sie sich den iterativen Newton-Raphson-Root-Finding-Algorithmus an.

in diesem sehr guten Papier zum Thema Werfen Sie einen Blick: http://www.spaceroots.org/documents/distance/distance-to-ellipse.pdf

Sorry für die unvollständige Antwort - ich die Gesetze der Mathematik und Natur völlig schuld ...

EDIT: oops, ich scheine a und b falsch herum im Diagramm xD

Answer 2

Sie müssen nur die Schnittmenge der Linie [P1,P0] zu Ihrer Elipse berechnen, die S1 ist.

Wenn die Leitung equeation ist:

und die elipse equesion ist:

als die Werte von S1 wird:

Jetzt müssen Sie nur den Abstand zwischen S1-P1 berechnen, die Formel (für A,B Punkte) ist:

Answer 3

Hier zu haben ist der Code in C# aus diesem Papier umgesetzt übersetzt für die Ellipse zu lösen: http://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf

Beachten Sie, dass dieser Code nicht getestet wurde - wenn Sie Fehler finden, lassen Sie es mich wissen.

//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

Das folgende Python-Code implementiert die beschriebenen Gleichungen auf „Distance from a Point to an Ellipse“ und verwendet Newtons Methode, die Wurzeln zu finden und dass der nächste Punkt auf der Ellipse auf den Punkt.

Leider, wie aus dem Beispiel ersichtlich, scheint es nur außerhalb der Ellipse genau zu sein. Innerhalb der Ellipse passieren seltsame Dinge.

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 

Hier ist ein Beispiel:

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() 

Welche eine Ellipse und den Abstand von der Grenze der Ellipse als Heatmap erzeugt. Wie zu sehen ist, ist der Abstand an der Grenze Null (tiefblau).

Answer 5

Es gibt eine relativ einfache numerische Methode mit besserer Konvergenz als die Newtons-Methode.Ich habe einen Blog-Post über, warum es funktioniert 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]))