Geometry & Recognition

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

$$\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} a x(t) \\ b y(t) \end {array}\right) $$

$$ x(t) = 3k (1-t)^2 t + 3 (1-t) t^2 + t^3 \\ y(t) = x(1-t) ,  \quad 0 \le t \le 1$$

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

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

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

$$f(t)= ({\bf B}(t) - {\bf P})\cdot {\bf B}'(t) = 0 \\ f(t)=a^2x(t)x'(t)-b^2y(t)y'(t)-ap x'(t)-bqy'(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);
}

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