distancia de un punto dado a una elipse dada

12 minutos de lectura

Tengo una elipse, definida por el Punto central, el radio X y el radio Y, y tengo un Punto. Quiero encontrar el punto en la elipse que está más cerca del punto dado. En la siguiente ilustración, sería S1.

grafico1

Ahora ya tengo el código, pero hay un error lógico en algún lugar y parece que no puedo encontrarlo. Desglosé el problema en el siguiente ejemplo de código:

#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 siguiente imagen:

instantánea1

Puede ver que en realidad encuentra puntos “cercanos” en la elipse, pero no son los puntos “más cercanos”. Lo que quiero intencionalmente es esto: (disculpen mi pobre dibujo)

instantánea2

extenderías las líneas en la última imagen, cruzarían el centro de la elipse, pero este no es el caso de las líneas en la imagen anterior.
Espero que entiendas la imagen. ¿Alguien puede decirme qué estoy haciendo mal?

  • será más fácil si solo describe su método para encontrar el punto que el código real

    – rango1

    09/04/2014 a las 10:30

Considere un círculo límite alrededor del punto dado (c, d), que pasa por el punto más cercano de la elipse. A partir del diagrama, está claro que el punto más cercano es tal que una línea trazada desde él hasta el punto dado debe ser perpendicular a la tangente compartida de la elipse y el círculo. Cualquier otro punto estaría fuera del círculo y, por lo tanto, debe estar más lejos del punto dado.

ingrese la descripción de la imagen aquí

Así que el punto que estás buscando es no la intersección entre la línea y la elipse, sino el punto (x, y) en el diagrama.

Gradiente de tangente:

ingrese la descripción de la imagen aquí

Gradiente de línea:

ingrese la descripción de la imagen aquí

Condición para líneas perpendiculares – producto de gradientes = -1:

ingrese la descripción de la imagen aquí

ingrese la descripción de la imagen aquí

ingrese la descripción de la imagen aquí

Cuando se reorganiza y se sustituye en la ecuación de tu elipse…

ingrese la descripción de la imagen aquí

…esto dará dos desagradables ecuaciones cuarticas (polinomio de cuarto grado) en términos de x o y. AFAIK no hay general métodos analíticos (algebraicos exactos) para resolverlos. Puede probar un método iterativo: busque el algoritmo iterativo de búsqueda de raíces de Newton-Raphson.

Echa un vistazo a este muy buen artículo sobre el tema:
http://www.spaceroots.org/documents/distance/distance-to-ellipse.pdf

Perdón por la respuesta incompleta: culpo totalmente a las leyes de las matemáticas y la naturaleza …

EDITAR: Vaya, parece que tengo a y b al revés en el diagrama xD

  • joder, tienes razón. Eso resuelve el problema de mi implementación defectuosa, no al encontrar el error, sino al exponer mi enfoque como defectuoso en primer lugar. Gracias

    – usuario2950911

    9 de abril de 2014 a las 12:04

  • @ usuario2950911 Me alegro de ser de ayuda

    usuario3235832

    9 de abril de 2014 a las 12:05

  • @ mat69 Me gusta tu forma de pensar, pero no. La matriz de transformación en la que está pensando es una matriz de escala, que es una transformación lineal. Para un círculo, el punto más cercano sería de hecho la intersección que él asumió, pero cuando se aplica la transformación inversa, todos los puntos en esa línea se escalarían en la misma cantidad, lo que llevaría a la respuesta incorrecta que obtuvo. otra vez.

    usuario3235832

    31 de julio de 2015 a las 15:07

  • Los polinomios cuárticos tienen soluciones exactas (son el grado más alto para tenerlas, ver Teoría de Galois). Sin embargo, no son triviales, por lo que podría valer la pena investigar si el polinomio es reducible.

    – MSalters

    16 de febrero de 2016 a las 14:19

  • @MSalters tiene razón, pero hacerlo no vale la pena el esfuerzo computacional en comparación con la convergencia iterativa

    usuario3235832

    18/02/2016 a las 20:20

avatar de usuario
usuario364952

Existe un método numérico relativamente simple con mejor convergencia que el método de Newton. Tengo una publicación de blog sobre por qué funciona. http://wet-robots.ghost.io/simple-method-for-distance-to-ellipse/

Esta implementación funciona sin funciones trigonométricas:

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

    tx = 0.707
    ty = 0.707

    a = semi_major
    b = semi_minor

    for x in range(0, 3):
        x = a * tx
        y = b * ty

        ex = (a*a - b*b) * tx**3 / a
        ey = (b*b - a*a) * ty**3 / b

        rx = x - ex
        ry = y - ey

        qx = px - ex
        qy = py - ey

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

        tx = min(1, max(0, (qx * r / q + ex) / a))
        ty = min(1, max(0, (qy * r / q + ey) / b))
        t = math.hypot(ty, tx)
        tx /= t 
        ty /= t 

    return (math.copysign(a * tx, p[0]), math.copysign(b * ty, p[1]))

Convergencia

Crédito a Adrián Esteban Para el Optimización sin activación.

  • Este enfoque es excelente. Se puede realizar una optimización que evita por completo cualquier función trigonométrica (todo el crédito se debe al problema publicado en el repositorio de github para esta publicación de blog: github.com/0xfaded/ellipse_demo/issues/1 ). Hice una implementación que usa esta optimización en C# para Unity3D: gist.github.com/JohannesMP/777bdc8e84df6ddfeaa4f0ddb1c7adb3

    – Juan

    13 de julio de 2018 a las 6:16


  • El código no parece estar completo después de su edición. P.ej tx no está definido antes de su primer uso.

    – JHBonarius

    13 de julio de 2018 a las 7:15

  • Vaya, no copié la optimización correctamente. Gracias por revisar @JHBonarius

    – Juan

    14 de julio de 2018 a las 1:26

  • ¿Debería cambiarse esta línea de “for x in range(0, 3):” a “for i in range(0, 3):”?

    – quarkz

    7 ene a las 9:00


avatar de usuario
juanml1135

Aquí está el código traducido a C# implementado a partir de este documento para resolver la elipse:
http://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf

Tenga en cuenta que este código no está probado; si encuentra algún error, hágamelo saber.

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

El siguiente código de Python implementa las ecuaciones descritas en “Distancia de un punto a una elipse” y utiliza el método de Newton para encontrar las raíces y, a partir de ahí, el punto más cercano de la elipse al punto.

Desafortunadamente, como se puede ver en el ejemplo, parece que solo es preciso fuera de la elipse. Dentro de la elipse suceden cosas raras.

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

Aquí hay un ejemplo:

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

que genera elipse girada una elipse y la distancia desde el límite de la elipse como un mapa de calor. Como puede verse, en el límite la distancia es cero (azul oscuro).

avatar de usuario
nwellnhof

Dada una elipse mi en forma paramétrica y un punto PAGS

ecuación

el cuadrado de la distancia entre PAGS y mi

ecuación

El mínimo debe satisfacer

ecuación

Usando las identidades trigonométricas

ecuación

y sustituyendo

ingrese la descripción de la imagen aquí

produce la siguiente ecuación de cuarto grado:

ingrese la descripción de la imagen aquí

Aquí hay un ejemplo de función C que resuelve la cuarta directamente y calcula pecado

void nearest(double a, double b, double x, double y, double *ecos_ret, double *esin_ret) {
    double ax = fabs(a*x);
    double by = fabs(b*y);
    double r  = b*b - a*a;
    double c, d;
    int switched = 0;

    if (ax <= by) {
        if (by == 0) {
            if (r >= 0) { *ecos_ret = 1; *esin_ret = 0; }
            else        { *ecos_ret = 0; *esin_ret = 1; }
            return;
        }
        c = (ax - r) / by;
        d = (ax + r) / by;
    } else {
        c = (by + r) / ax;
        d = (by - r) / ax;
        switched = 1;
    }

    double cc = c*c;
    double D0 = 12*(c*d + 1);      // *-4
    double D1 = 54*(d*d - cc);     // *4
    double D  = D1*D1 + D0*D0*D0;  // *16

    double St;
    if (D < 0) {
        double t = sqrt(-D0);             // *2
        double phi = acos(D1 / (t*t*t));
        St = 2*t*cos((1.0/3)*phi);        // *2
    } else {
        double Q = cbrt(D1 + sqrt(D));    // *2
        St = Q - D0 / Q;                  // *2
    }

    double p    = 3*cc;                          // *-2
    double SS   = (1.0/3)*(p + St);              // *4
    double S    = sqrt(SS);                      // *2
    double q    = 2*cc*c + 4*d;                  // *2
    double l    = sqrt(p - SS + q / S) - S - c;  // *2
    double ll   = l*l;                           // *4
    double ll4  = ll + 4;                        // *4
    double esin = (4*l)    / ll4;
    double ecos = (4 - ll) / ll4;

    if (switched) {
        double t = esin;
        esin = ecos;
        ecos = t;
    }

    *ecos_ret = copysign(ecos, a*x);
    *esin_ret = copysign(esin, b*y);
}

¡Pruébelo en línea!

  • ¿Existen límites inherentes a los valores que puede tomar? check_dist(4.0, 2.0, 0, 6); resultados en solo nan por ejemplo.

    – Juan

    13 de julio de 2018 a las 4:07

  • @Johannes Probablemente olvidé agregar un caso especial para x = 0.

    – nwellnhof

    13 de julio de 2018 a las 9:06

  • Devoluciones NaN por (1.0, 0.75, -0.128053, -0.743808)

    – brad

    19 oct 2018 a las 15:43

  • @Brad Correcto, en algunos casos, mi código elige la solución incorrecta del cuartico.

    – nwellnhof

    19 oct 2018 a las 17:09

Solo necesitas calcular la intersección de la línea. [P1,P0] a tu elipse que es S1.

Si la ecuación lineal es:

ingrese la descripción de la imagen aquí

y la ecuación de elipse es:

Ecuación de elipse

que los valores de S1 estarán:

ingrese la descripción de la imagen aquí

Ahora solo necesitas calcular la distancia entre S1 a P1 la fórmula (para A,B puntos) es:
distancia

  • ¿Existen límites inherentes a los valores que puede tomar? check_dist(4.0, 2.0, 0, 6); resultados en solo nan por ejemplo.

    – Juan

    13 de julio de 2018 a las 4:07

  • @Johannes Probablemente olvidé agregar un caso especial para x = 0.

    – nwellnhof

    13 de julio de 2018 a las 9:06

  • Devoluciones NaN por (1.0, 0.75, -0.128053, -0.743808)

    – brad

    19 oct 2018 a las 15:43

  • @Brad Correcto, en algunos casos, mi código elige la solución incorrecta del cuartico.

    – nwellnhof

    19 oct 2018 a las 17:09

avatar de usuario
Roberto Hudjakov

He resuelto el problema de la distancia a través de puntos focales.

Para cada punto de la elipse

r1 + r2 = 2*a0

dónde

r1 – Distancia euclidiana desde el punto dado al punto focal 1

r2 – Distancia euclidiana desde el punto dado hasta el punto focal 2

a0 – longitud del semieje mayor

También puedo calcular r1 y r2 para cualquier punto dado, lo que me da otra elipse en la que se encuentra este punto que es concéntrica a la elipse dada. Entonces la distancia es

d = Abs((r1 + r2) / 2 – a0)

¿Ha sido útil esta solución?