任意点到斜椭圆最小距离

最近由于项目需求，要计算任意点到斜椭圆的距离，方程相信大家都知道怎么建立，但是它是一个4次方程，matlab都很不好解，写到程序中也很难，搜了很多帖子，都是只给方程，不给解的，其实原理大家都知道，如何用程序解才是大家关心的，所以一直都没有找到合理的答案，关键时刻还是google给力，度娘一般的东西还撮合，可学术性的、英文的确是一塌糊涂，网络上收集到这样一篇论文,我也是用这个算法最终解决我的问题的，它只使用与标准方程，斜椭圆只是多了平移和旋转操作而已。

###第一象限中的点到标准椭圆的距离###

double DistancePointEllipseSpecial(const double e[2], const double y[2], double x[2])
{
    double distance;
    if (y[1] > 0)
    {
        if (y[0] > 0)
        {
            // bisect to compute the root of F(t) for t >= -e1*d1
            double esqr[2] = {e[0] * e[0], e[1] * e[1]};
            double ey[2] = {e[0] * y[0], e[1] * y[1]};
            double t0 = -esqr[1] + ey[1];
            double t1 = -esqr[1] + sqrtf(ey[0]*ey[0] + ey[1]*ey[1]);
            double t = t0;
            for (int i = 0; ; ++i)
            {
                t = 0.5*(t0 + t1);
                if (t == t0 || t == t1)
                {
                    break;
                }
                double r[2] = { ey[0]/(t + esqr[0]), ey[1]/(t + esqr[1]) };
                double f = r[0]*r[0] + r[1]*r[1] - 1;
                if (f > 0)
                {
                    t0 = t;
                }
                else if (f < 0)
                {
                    t1 = t;
                }
                else
                {
                    break;
                }
            }
            x[0] = esqr[0]*y[0]/(t + esqr[0]);
            x[1] = esqr[1]*y[1]/(t + esqr[1]);
            double d[2] = { x[0] - y[0], x[1] - y[1] };
            distance = sqrtf(d[0]*d[0] + d[1]*d[1]);
        }
        else // y0 == 0
        {
            x[0] = 0;
            x[1] = e[1];
            distance = fabs(y[1] - e[1]);
        }
    }
    else // y1 == 0
    {
        double denom0 = e[0]*e[0] - e[1]*e[1];
        double e0y0 = e[0]*y[0];
        if (e0y0 < denom0)
        {
            // y0 is inside the subinterval.
            double x0de0 = e0y0/denom0;
            double x0de0sqr = x0de0*x0de0;
            x[0] = e[0]*x0de0;
            x[1] = e[1]*sqrtf(fabs(1 - x0de0sqr));
            double d0 = x[0] - y[0];
            distance = sqrtf(d0*d0 + x[1]*x[1]);
        }
        else
        {
            // y0 is outside the subinterval. The closest ellipse point has
            // x1 == 0 and is on the domain-boundary interval (x0/e0)^2 = 1.
            x[0] = e[0];
            x[1] = 0;
            distance = fabs(y[0] - e[0]);
        }
    }
    return distance;
}

任意象限的点到标准椭圆的距离###

double DistancePointEllipse (const double e[2], const double y[2], double x[2])
{
    // Determine reflections for y to the first quadrant.
    bool reflect[2];
    int i, j;
    for (i = 0; i < 2; ++i)
    {
        reflect[i] = y[i] < 0;
    }
    // Determine the axis order for decreasing extents.
    int permute[2];
    if (e[0] < e[1])
    {
        permute[0] = 1; permute[1] = 0;
    }
    else
    {
        permute[0] = 0; permute[1] = 1;
    }
    int invpermute[2];
    for (i = 0; i < 2; ++i)
    {
        invpermute[permute[i]] = i;
    }
    double locE[2], locY[2];
    for (i = 0; i < 2; ++i)
    {
        j = permute[i];
        locE[i] = e[j];
        locY[i] = y[j];
        if (reflect[j])
        {
            locY[i] = -locY[i];
        }
    }
    double locX[2];
    double distance = DistancePointEllipseSpecial(locE, locY, locX);
    // Restore the axis order and reflections.
    for (i = 0; i < 2; ++i)
    {
        j = invpermute[i];
        if (reflect[i])
        {
            locX[j] = -locX[j];
        }
        x[i] = locX[j];
    }
    return distance;
}

任意点到斜椭圆的距离

double dist_point_ellipse( const POINT &point1, const ELLIPSE &ellipse )
{
    double distance;
    double yin[2], xout[2];
    double ein[2] = {ellipse.a, ellipse.b};
    double cost = cos(ellipse.theta * _PI / 180);
    double sint = sin(ellipse.theta * _PI / 180);

    POINT point = point1;
    point = point - ellipse.center;
    yin[0] = cost * point.x + sint * point.y;
    yin[1] = -sint * point.x + cost * point.y;

    distance = DistancePointEllipse(ein, yin, xout);

    POINT pt;
    pt.x = cost * xout[0] - sint * xout[1];
    pt.y = sint * xout[0] + cost * xout[1];
    pt = pt + ellipse.center;

    distance = sqrtf((pt.x - point1.x) * (pt.x - point1.x) + (pt.y - point1.y) * (pt.y - point1.y));

    return distance;
}