Расстояние от заданной точки до заданного эллипса

У меня есть эллипс, определяемый центральной точкой, радиусом X и радиусом Y, и у меня есть точка. Я хочу найти точку на эллипсе, ближайшем к данной точке. На приведенном ниже рисунке это будет S1.

Теперь у меня уже есть код, но в нем есть логическая ошибка, и я, похоже, не могу его найти. Я нарушил проблему до следующего примера кода:

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

Он создает следующее изображение:

Вы можете видеть, что на самом деле он находит "близкие" точки на эллипсе, но это не "ближайшие" точки. Я намеренно хочу, чтобы это: (извините мой плохой рисунок)

Вы протянете линии в последнем изображении, они пересекут центр эллипса, но это не относится к строкам предыдущего изображения.
Надеюсь, вы получите картину. Может ли кто-нибудь сказать мне, что я делаю неправильно?

Ответы

Ответ 1

Рассмотрим ограниченную окружность вокруг данной точки (c, d), которая проходит через ближайшую точку на эллипсе. Из диаграммы видно, что ближайшая точка такова, что линия, взятая от нее к данной точке, должна быть перпендикулярна общей касательной эллипса и круга. Любые другие точки будут вне круга и поэтому должны быть дальше от данной точки.

Итак, точка, которую вы ищете, - это не пересечение линии и эллипса, а точка (x, y) на диаграмме.

Градиент касательной:

Градиент строки:

Условие для перкудикулярных линий - произведение градиентов = -1:

При перестановке и подстановке в уравнение вашего эллипса...

... это даст два неприятных четвертичных (4-градусных многочлена) уравнения в терминах либо x, либо y. AFAIK нет общих аналитических (точных алгебраических) методов их решения. Вы можете попробовать итеративный метод - найти алгоритм итеративного корневого поиска Newton-Raphson.

Взгляните на эту очень хорошую статью по этому вопросу: http://www.spaceroots.org/documents/distance/distance-to-ellipse.pdf

Извините за неполный ответ - я полностью обвиняю законы математики и природы...

ИЗМЕНИТЬ: oops, я, кажется, имеет a и b неправильный путь на диаграмме xD

Ответ 2

Существует относительно простой численный метод с лучшей сходимостью, чем метод Ньютона. У меня есть запись в блоге о том, почему это работает 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])

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

Благодарим Адриана Стивенса за оптимизацию без триггеров.

Ответ 3

Вот код, переведенный на С#, реализованный из этой статьи, чтобы решить для эллипса: http://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf

Обратите внимание, что этот код не проверен - если вы обнаружите какие-либо ошибки, сообщите мне.

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

Ответ 4

Следующий код python реализует уравнения, описанные в "" Расстояние от точки до эллипса" и использует метод newton для поиска корней и от нее ближайшая точка на эллипсе до точки.

К сожалению, как видно из примера, он кажется точным только вне эллипса. Внутри эллипса происходят странные вещи.

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

Вот пример:

rx, ry = 12, 35  # major, minor ellipse axis
x0 = y0 = 50  # center point of the ellipse
angle = 45  # ellipse 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()

Что генерирует эллипс и расстояние от границы эллипс как тепловая карта. Как видно, на границе расстояние равно нулю (темно-синий).

Ответ 5

Вам просто нужно вычислить пересечение линии [P1,P0] с вашим elipse, который равен S1.

Если выравнивание линии:

и равенство elipse равно:

чем значения S1 будут:

Теперь вам просто нужно рассчитать расстояние между S1 до P1, формула (для A,B точек):

Ответ 6

Дан эллипс E в параметрической форме и точка P

квадрат расстояния между P и E (t) равен

Минимум должен удовлетворять

Использование тригонометрических тождеств

и заменяя

дает следующее квартическое уравнение:

Вот пример функции C, которая напрямую решает задачу и вычисляет sin (t) и cos (t) для ближайшей точки на эллипсе:

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

Попробуйте онлайн!

Ответ 7

Я решил проблему расстояния с помощью фокусных точек.

Для каждой точки на эллипсе

r1 + r2 = 2 * a0

где

r1 - евклидово расстояние от заданной точки до фокуса 1

r2 - евклидово расстояние от заданной точки до фокуса 2

а0 - длина большой полуоси

Я также могу вычислить r1 и r2 для любой заданной точки, которая дает мне еще один эллипс, на котором лежит эта точка, который является концентрическим по отношению к данному эллипсу. Так что расстояние

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