如何在C++中实现过指定点的椭圆法线计算函数？

问题背景 #

要计算某点到椭圆圆周的最短距离，核心是找到过该点的椭圆法线（及对应角度）——这条法线就是点到椭圆的最短路径所在直线。

数学推导 #

设椭圆的半长轴为a = ew/2，半短轴为b = eh/2（ew为椭圆宽度，eh为椭圆高度），外部点相对于椭圆中心的坐标为(x,y)。

椭圆上任意一点可表示为(a cos t, b sin t)，其中t为参数角。根据法线条件：该点的法线需经过外部点(x,y)，结合椭圆法线斜率公式，可推导出关于t的方程：

a² y cos t - b² x sin t = (a² - b²) sin t cos t

直接解析求解此方程需处理四次多项式，编程中更适合用牛顿迭代法进行数值求解，迭代收敛速度快且实现简单。

C++代码实现 #

#include <cmath>
#include <stdexcept>

/// Calculate the normal angle from ellipse circumference to the given point (relative to ellipse center)
/// @param ew The width of the ellipse
/// @param eh The height of the ellipse
/// @param x Relative x position (relative to ellipse center)
/// @param y Relative y position (relative to ellipse center)
/// @return The normal direction in radians (range: [0, 2π))
auto calculateNormal(float ew, float eh, float x, float y) -> double {
    const double a = ew / 2.0;
    const double b = eh / 2.0;
    const double a_sq = a * a;
    const double b_sq = b * b;
    const double diff_sq = a_sq - b_sq;

    // 特殊情况：点在椭圆中心
    if (std::fabs(x) < 1e-8 && std::fabs(y) < 1e-8) {
        return 0.0; // 任意方向均可，返回0弧度
    }

    // 初始值：取点(x,y)对应的极角作为迭代起点
    double t = std::atan2(y, x);
    const int max_iterations = 50;
    const double tolerance = 1e-8;

    for (int i = 0; i < max_iterations; ++i) {
        const double cos_t = std::cos(t);
        const double sin_t = std::sin(t);
        const double cos2t = cos_t * cos_t - sin_t * sin_t;

        // 计算f(t)
        const double f = a_sq * y * cos_t - b_sq * x * sin_t - diff_sq * sin_t * cos_t;
        // 计算f'(t)（导数）
        const double f_prime = -a_sq * y * sin_t - b_sq * x * cos_t - diff_sq * cos2t;

        // 导数接近0时，调整初始值重新迭代
        if (std::fabs(f_prime) < 1e-10) {
            t += M_PI / 4.0;
            continue;
        }

        // 牛顿迭代更新t
        const double delta_t = f / f_prime;
        t -= delta_t;

        // 确保t在[0, 2π)范围内
        t = std::fmod(t, 2 * M_PI);
        if (t < 0) t += 2 * M_PI;

        // 达到精度要求，停止迭代
        if (std::fabs(delta_t) < tolerance) {
            break;
        }
    }

    // 计算法线角度：取椭圆上点到外部点的连线方向
    const double ellipse_x = a * std::cos(t);
    const double ellipse_y = b * std::sin(t);
    return std::atan2(y - ellipse_y, x - ellipse_x);
}

代码说明 #

1. 参数处理：先将椭圆宽度、高度转换为半长轴a和半短轴b，计算平方值避免重复运算。
2. 特殊情况处理：当点在椭圆中心时，任意方向均为法线，直接返回0弧度。
3. 牛顿迭代：
   • 以点(x,y)的极角作为初始迭代值，保证快速收敛。
   • 每次迭代计算函数值f(t)和导数f’(t)，更新t直到满足精度要求。
   • 处理导数为0的异常情况，调整初始值避免死循环。
4. 角度计算：最终法线角度为椭圆上的点指向外部点(x,y)的方向角，用atan2直接计算得到弧度制结果。

内容的提问来源于stack exchange，提问作者Flovdis