Approximate Distance Between Ellipse and Point

타원은 베지어 곡선을 이용해서 근사적으로 표현할 수 있음을 잘 알고 있다(https://kipl.tistory.com/313). 표준형 타원의 경우 1사분면에서 표현만 알면 대칭성에 의해서 나머지 사분면에서는 바로 구할 수 있으므로 거리를 구하는 점이 1사분면에 있는 경우만 알아보면 충분하다. 장축 반지름이 $a$, 단축 반지름이 $b$인 타원은 다음과 같이 베이지 곡선으로 표현된다.

Bezier Curve Approximation of an Ellipse

$$\mathbf{B}(t)=(1-t^3){\mathbf{P}}_0+3t(1-t{)}^2{\mathbf{P}}_1+3t^2(1-t){\mathbf{P}}_2+t^3{\mathbf{P}}_3=\left(\begin{array}{c}ax(t)\\ by(t)\end{array}\right)$$

$$x(t)=3k(1-t{)}^{2}t+3(1-t)t^2+t^3$$

$$y(t)=x(1-t),0\le t\le 1$$

$k$ 값은 $t=1/2$일 때 $(a/\sqrt{2},b/\sqrt{2})$을 통과하는 조건을 부여하여 정하면

$$k=\frac{4}{3}(\sqrt{2}-1)=0.5522847498...$$임을 알 수 있다. 

 

한 점 $\mathbf{P}=(p,q)$에서 베지어 곡선 위의 한 점 $\mathbf{B}(t)$까지의 거리가 최단거리가 되기 위해서는 다음 직교 조건

$$f(t)=(\mathbf{B}(t)-\mathbf{P})\cdot {\mathbf{B}}^{\prime }(t)=0$$

$$f(t)=a^{2}x(t)x^{\prime }(t)-b^{2}y(t)y^{\prime }(t)-ap{x}^{\prime }(t)-bq{y}^{\prime }(t)=0$$

을 만족하는 근을 찾으면 된다. 풀어야 할 방정식이 5차이므로 직접적으로 근을 찾는 것은 불가능하므로 Newton-Raphson 방법을 사용하도록하자. 앞선 정확한 거리 계산에서와는 달리 초기 $t$ 값 설정에 민감하게 의존하지 않는다.

//cubic bezier_x: P0(0, 1), P1(k, 1),P2(1, k), P3(1, 0);
double B(double t) { //t in [0:1]
    const double k = 4.*(sqrt(2.)-1)/3;
    return t*(3*k + t*(3 - 6*k + (-2 + 3*k)*t));
}
// derivative of B(t);
double DB(double t) {
    const double k = 4.*(sqrt(2.)-1)/3;
    return 3*k + t*(6 - 12*k + (-6 + 9*k)*t);
}
// derivative of DB(t);
double D2B(double t) {
    const double k = 4.*(sqrt(2.)-1)/3;
    return 6 - 12*k + (-12 + 18*k)*t;
}
// ellipse radii=(a, b);
double dist2EllipseBezier3(double p, double q, double a, double b,
                           double& xt, double& yt) {
    if (a == b) return dist2Circle(p, q, a, xt, yt);
    double x = fabs(p), y = fabs(q);
    const double eps = 0.001 / max(a, b);
    double t = 0.5;  // mid
    while (1) {
        // Newton-Raphson;
        double f = a*a*B(t)*DB(t)-b*b*B(1-t)*DB(1-t)-a*x*DB(t)+b*y*DB(1-t);
        if (f == 0) break;
        double df = a*a*(DB(t)*DB(t)+B(t)*D2B(t))+b*b*(DB(1-t)*DB(1-t)+B(1-t)*D2B(1-t))
                    -a*x*D2B(t)-b*y*D2B(1-t);
        double dt = f / df;
        t -= dt;
        t = max(0, min(1, t));
        if (abs(dt) < eps ) break;
    }
 
    xt = a * B(t); yt = b * B(1-t);
    xt = p >= 0? xt: -xt;
    yt = q >= 0? yt: -yt;
    return hypot(p - xt, q - yt);
}

이심률이 큰 경우와 작은 경우