Geometry & Recognition

타원  외부의 한 지점에서 타원까지의 최단거리를 구하자. 타원 피팅 문제에서 피팅 에러를 예측하는데 유용하게 사용할 수 있다. 일반적인 타원은 타원의 중심과 회전각을 구해서 원점으로 평행이동 후 다시 반대로 회전시켜서 표준형으로 만든 후 구하면 되기 때문에 표준형에 대해서만 고려해도 충분하다.(외부점도 동일하게 평행이동 시키고 회전시켜야 한다). 또한 표준형 타원의 경우도  대칭성에 의해서 외부점이 1 사분면에 있는 경우만 살펴보면 충분하다. 

외부점 $(p, q)$에 타원 상의 한 점 $(x, y)$까지의 거리가 최단 거리가 되기 위해서는  두 점을 잇는 선분과 타원에서 접선이 서로 직교해야 한다는 사실에서 $(x,y)$ 위치를 정할 수 있다.

우선 타원 상의 한 점 $(x,y)$에서  접선의 기울기는 

$$ \frac{dy}{dx} = - \frac{b^2 x}{a^2 y}$$

이고, $(p,q)$까지 선분의 기울기가 $$ m =\frac{ y - q}{x-p}$$이다. 서로 직교하기 위해서는 이 둘의 곱이 $-1$이어야 한다.

$$ \left(-\frac{b^2 x}{a^2 y} \right)  \frac{y-q}{x-p} = -1 \quad \to\quad (a^2 -b^2 )xy + b^2 qx - a^2 y p = 0$$

접점 $(x, y)$는 타원 위에 있다는 점을 고려하면 $x$나 $y$에 대한 4차 방정식을 풀면 답을 analytic하게  쓸 수 있다. 하지만 직접적으로 이용하기에는 너무 복잡한 형태이므로 대부분의 경우는 Newton-Raphson 방법과 같은 반복적인 수치해석 기법을 이용해서 답을 얻는다. 이를 위해 타원을 

$$ x= a \cos t, \quad y = b \sin t$$

로 매개화하면 최단거리를 주는 점을 찾는 문제는 다음 $f(t)$의 근을 찾는 문제로 환원이 된다.

$$ f(t) =  (a^2 -b^2) \cos t \sin t - ap \sin t + bq \cos t = 0, \quad 0 \le t  \le\frac{\pi}{2}$$

$z= \tan(t/2)$로 치환하면 이 식의 근을 구하는 문제는 $z$에 대한 4차 방정식의 근을 구하는 문제와 동일함을 보일 수 있다.

/*(a, b)=radii of ellipse; 타원에 포함되는 점에 대해서는 초기조건 설정에 주의;*/
double dist2Ellipse(double p, double q, double a, double b, double& xt, double& yt) {
    if (a == b) return dist2Circle(p, q, a, xt, yt); //원의 경우 별도 처리;
    if (x == 0 && y == 0) //원점의 경우 별도처리;
        if (a > b) { xt = 0,  yt = b; return b;} 
        else { xt = a, yt = 0; return a;}
    double x = fabs(p), y = fabs(q);
    const double eps = 0.01 / max(a, b);
    const int maxiters = 10;
    double t;
    if ((x*x/a/a + y*y/b/b) < 1) {// 타원 내부의 경우 초기조건의 조정이 필요함;
        double t1 = atan2(y/b, x/a), t2 = atan(1.);
        if (a > b) t = max(t1, t2);
        else t = min(t1, t2);
    } else t = atan2(y/b, x/a);
    for (int iter = 0; iter < maxiters; iter++) {
        double c = cos(t), s = sin(t);
        double f = (a*a - b*b)*c*s - a*x*s + b*y*c;
        if (f == 0) break;
        double Df = (a*a - b*b)*(c*c - s*s) - a*x*c - b*y*s;
        double dt = -f/Df;
        t += dt;
        t = max(0, min(2*atan(1.), t)); //[0:pi/2];
        if (abs(dt) < eps ) break;
    }
    xt = a * cos(t); yt = b * sin(t);
    xt = p >= 0? xt: -xt;
    yt = q >= 0? yt: -yt;
    return hypot(p-xt, q-yt);
}

위 알고리즘은 Newton-Raphson iteration이 수렴하지 않을 수도 있다. 좀 더 좋은 방법은 외부점에서 타원에 내린 수선의 발에서 접촉원(osculating circle)을 구하면 수선의 연장선은 접촉원의 중심을 지난다.

$$ \vec{c} =\left( \frac{ a^2 - b^2 }{a} \cos ^3 t , \frac{b^2 - a^2}{b} \sin^3 t\right)  $$

만약 타원의 한 지점에서 접촉원 중심과 외부점을 연결하는 직선이 접선에 수직이 아닌 경우에는 그 지점은 최단점이 아니다. 이때는 이 직선과 타원의 교점까지의 거리가 좀 더 짧아진다. 다시 이 교점에서 접촉원을 구한 후 외부점과 중심을 연결하는 선분이 교차하는 새로운 교점을 찾는 과정을 반복하면 점점 교점은 최단점에 접근함을 알 수 있다. 이를 이용하면 좀 더 효율적인 점과 타원사이의 최단거리를 구하는 알고리즘을 만들 수 있다. 타원 위의  한 점 $\vec{e}$에서 접촉원의 중심을 구하면 중심  $\vec{c}$와  외부점$\vec{p}$을 연결하는 직선과 타원의 교점은 대략적으로 $\vec c$와 $\vec p$를 $|\vec{r}=\vec{e}-\vec{c}|: |\vec{q}=\vec{p}-\vec{c}| - |\vec{e}-\vec{c}|=r:q-r$로 내부하는 지점이다. 따라서 새로운 교점의 위치는

$$ \vec{r}' \approx \vec{c} + \frac{r}{q} \vec{q} $$

이다. 물론 이 지점이 정확히 타원 위에 있지는 않다. 타원 위에 있는 점은 각 축의 반지름으로 각성분을 나누면 단위원 위에 있어야 한다는 사실을 이용하면 $\vec{r}'$에 가장 가까이 있는 타원 상의 한 지점을 얻을 수 있다.  즉, 

$$  \cos t = \frac{x'/a}{\sqrt{(x'/a)^2+(y'/b)^2}},\quad    \sin t = \frac{y'/b}{\sqrt{(x'/a)^2 + (y'/b)^2 }}$$

인 $t$을 구하면 새로운 타원 위의 지점은

$$ x''= a \cos (t), \quad y'' = b \sin t$$

로 된다.

double dist2Ellipse2(double p, double q, double a, double b, double& xt, double& yt) {
    if (a==b) return dist2Circle(p, q, a, xt, yt);//원의 경우에는 별도로;
    const double px = abs(p), double py = abs(q);
    const double eps = 0.01/max(a,b);
    const double cc = a*a - b*b;
    double cost, sint; 
    //근이 여러개 존재하므로 초기값 설정에 주의해야 한다;
    if ((px*px/a/a + py*py/b/b) < 1) //타원에 포함되는 경우;
    	cost = sint = 1/sqrt(2.); //  pi/4에서 출발;
    else { //타원바깥에 있는 경우;
        cost = cos(atan2(py/b, px/a));
    	sint = sin(atan2(py/b, px/a));
    }
    double t = 1, told = t;
    while (1) {
        // 곡률중심;evoute
        double cx = cc * pow(cost, 3.) / a;
        double cy = - cc * pow(sint, 3.) / b;
        // 곡률중심에 대한 타원 위의 점(교점)
        double x = a*cost - cx, y = b*sint - cy;
        // 곡률중심에 대한 외부점 변위
        double qx = px - cx, qy = py - cy;
        //거리;
        double r = hypot(x, y), q = hypot(qx, qy);
        cost = min(1, max(0, (qx * r / q + cx) / a));
        sint = min(1, max(0, (qy * r / q + cy) / b));
        double t = hypot(cost, sint); //1이 되면 좋지만 안되는 경우;
        cost /= t; sint /= t;
        if (abs(t-told) < eps) break;
        told = t;
    }
    xt = abs(a *cost), yt = abs(b *sint);
    xt = p >= 0 ? xt: -xt;
    yt = q >= 0 ? yt: -yt;
    return hypot(p - xt, q - yt);
}