从给定点到给定椭圆的距离

Question

我有一个椭圆,由中心点,radiusX和radiusY定义,我有一个Point.我想在椭圆上找到最接近给定点的点.在下图中,那将是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);
}

它产生以下图像:

您可以看到它实际上在椭圆上找到"近"点,但它不是"最近"的点.我故意想要的是:(请原谅我糟糕的画作)

你会扩展最后一个图像中的线条,它们会穿过椭圆的中心,但是前一个图像中的线条不是这种情况.
我希望你能得到这张照片.谁能告诉我我做错了什么？

Answer 1

考虑给定点(c,d)周围的边界圆,它通过椭圆上的最近点.从图中可以清楚的是,最近点是这样的,从中得出给定的点的线必须垂直于椭圆和圆的切线共享.任何其他点都在圆圈之外,因此必须远离给定点.

因此,您要寻找的不是线与椭圆之间的交点,而是图中的点(x,y).

切线渐变:

线条渐变:

perpedicular line的条件 - 梯度的乘积= -1:

重新排列并替换为椭圆方程时......

...这将给出两个令人讨厌的四次(四次多项式)方程,无论是x还是y.AFAIK没有通用的分析(精确代数)方法来解决它们.您可以尝试迭代方法 - 查找Newton-Raphson迭代根寻找算法.

看看在这个主题非常好纸: http://www.spaceroots.org/documents/distance/distance-to-ellipse.pdf

对不完整的答案感到抱歉 - 我完全责怪数学和自然法则......

编辑:哎呀,我似乎在图xD中有一个错误的方法

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

感谢阿德里安·斯蒂芬斯的触发-免费的优化。