c++ - propiedades - lado recto de la elipse

distancia del punto dado a la elipse dada

Question

Tengo una elipse, definida por Center Point, radiusX y radiusY, y tengo un punto. Quiero encontrar el punto en la elipse más cercano al punto dado. En la ilustración de abajo, eso sería S1.

Ahora ya tengo el código, pero hay un error lógico en algún lugar y parece que no puedo encontrarlo. Rompí el problema con 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:

Puedes ver que realmente encuentra puntos "cercanos" en la elipse, pero no son los puntos "más cercanos". Lo que intencionalmente quiero es esto: (disculpe mi pobre dibujo)

Si ampliara 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?

Answer 1

Solo necesita calcular la intersección de la línea [P1,P0] con su elipse que es S1 .

Si la equeación de línea es:

y la equesion del elipse es:

que los valores de S1 serán:

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

Answer 2

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 se ha 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;
    }

Answer 3

Considere un círculo delimitador alrededor del punto dado (c, d), que pasa por el punto más cercano en la elipse. 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 alejado del punto dado.

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

Gradiente de tangente:

Gradiente de línea:

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

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

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

Eche 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: uy, parece que tengo a y b al revés en el diagrama xD

Answer 4

El siguiente código python implementa las ecuaciones descritas en " Distancia de un punto a una elipse " y usa el método de Newton para encontrar las raíces y, desde allí, el punto más cercano de la elipse hasta el punto.

Desafortunadamente, como se puede ver en el ejemplo, parece ser solo exacto fuera de la elipse. Dentro de la elipse ocurren 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 una elipse y la distancia desde el límite de la elipse como un mapa de calor. Como se puede ver, en el límite la distancia es cero (azul profundo).

Answer 5

Hay un método numérico relativamente simple con una mejor convergencia que el Método de Newton. Tengo una publicación en el blog sobre por qué funciona 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]))