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$이어야 한다.

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

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

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

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

$$f(t)=(a^2-b^2)\cos t\sin t-ap\sin t+bq\cos t=0,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)을 구하면 수선의 연장선은 접촉원의 중심을 지난다.

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

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

$${\overrightarrow{r}}^{\prime }\approx \overrightarrow{c}+\frac{r}{q}\overrightarrow{q}$$

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

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

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

$$x^{″}=a\cos (t),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);
}

$m$개의 control point $\{{\mathbf{Q}}_i\}$가 주어지면  $p$차의 B-Spline curve는 basis 함수를 $N_{i,p}(t)$을 써서

$$\mathbf{B}(t)=\sum _{i=0}^{m-1}{N}_{i,p}(t){\mathbf{Q}}_i$$

로 표현할 수 있다. 이를 이용해서 일정한 순서로 샘플링된 평면상의 $N$개의 입력점 $\{{\mathbf{P}}_i\}$을 찾아보자. B-spline 곡선을 이용하면 이 문제는 control 점을 찾는 문제가 된다. 곡선이 입력점을 잘 표현하기 위해서는 곡선과 입력점과의 차이를 최소로 하는 control 점을 찾아야 한다:

$${\mathbf{Q}}^{\ast }=\text{argmin}(L),L:=\sum _{k=0}^{N-1}|\mathbf{B}(t_k)-{\mathbf{P}}_k{|}^2$$

여기서 $\{t_k|k=0,1,...,N-1\}$는 입력점이 얻어지는 sample time으로 $0=t_0\le t_1\le ...\le t_{N-1}=1$로 rescale 할 수 있다. 

행렬을 이용해서 식을 좀 더 간결하게 표현할 수 있다. 

$$\begin{gathered} L=\sum _{k=0}^{N-1}|\mathbf{B}(t_k)-{\mathbf{P}}_k{|}^2=\sum _{k=0}^{N-1}{|\sum _{i=0}^{m-1}{N}_{i,p}(t_k){\mathbf{Q}}_i-{\mathbf{P}}_k|}^2=\sum _{k=0}^{N-1}{|\sum _{i=0}^{m-1}{A}_{ki}{\mathbf{Q}}_i-{\mathbf{P}}_k|}^2, \\ {A}_{ki}=N_{i,p}(t_k) \end{gathered}$$

로 쓸 수 있으므로, $\hat{Q}=({\mathbf{Q}}_0,{\mathbf{Q}}_1,...,{\mathbf{Q}}_{m-1}{)}^t$, $\hat{P}=({\mathbf{P}}_0,{\mathbf{P}}_1,...,{\mathbf{P}}_{N-1}{)}^t$인 벡터, 그리고 $\hat{A}=(A_{ki})$인 행렬로 보면

$$L={|\hat{A}\cdot \hat{Q}-\hat{P}|}^2=({\hat{Q}}^t\cdot {\hat{A}}^t-{\hat{P}}^t)\cdot (\hat{A}\cdot \hat{Q}-\hat{P}).$$

위 식의 값을 최소로 하는 최소자승해는 ${\hat{Q}}^t$에 대한 미분 값이 0인 벡터를 찾으면 된다; $$\frac{\partial L}{\partial {\hat{Q}}^t}={\hat{A}}^t\cdot \hat{A}\cdot \hat{Q}-{\hat{A}}^t\cdot \hat{P}=0.$$ 이 행렬 방정식의 해는

$${\hat{Q}}^{\ast }=({\hat{A}}^t\cdot \hat{A}{)}^{-1}\cdot ({\hat{A}}^t\cdot \hat{P})$$ 로 표현된다. ${\hat{A}}^t\cdot \hat{A}$가 (banded) real symmetric ($m\times m$) 행렬이므로 Cholesky decomposion을 사용하면 쉽게 해를 구할 수 있다. $\hat{A}$가 banded matrix 형태를 가지는 이유는 basis가 local support에서만 0이 아닌 값을 취하기 때문이다. 

open b-spline(cubic)

std::vector<CfPt> BSplineFit_LS(const std::vector<CfPt>& data,
                                const int degree,    // cubic(3); 
                                const int nc)        // num of control points;
{
    // open b-spline;
    std::vector<double> knot((nc - 1) + degree + 2);
    for (int i = 0; i <= nc + degree; i++) knot[i] = i;
    
    std::vector<double> t(data.size());       // parameter;
    double scale = (knot[nc] - knot[degree]) / (data.size()-1);
    for (int i = data.size(); i-->0;) 
        t[i] = knot[degree] + scale * i;

    // design matrix;
    std::vector<double> A(ndata * nc); 
    for (int i = data.size(); i-->0;)
        for (int j = 0; j < nc; j++)
            A[i*nc + j] = Basis(j, degree, &knot[0], t[i]); //A(i,j)=N_j(t_i)

    // scattering matrix: S = A^t * A; real-symmetric matrix;
    std::vector<double> S(nc * nc); 
    for (int i = 0; i < nc; ++i)
        for (int j = i; j < nc; ++j) {
            double s = 0;
            for (int k = data.size(); k-->0;)
                s += A[k*nc + i] * A[k*nc + j];
            S[i*nc + j] = s;
        }
    // fill lower triangle;
    for (int i = 0; i < nc; ++i) 
        for (int j = 0; j < i; j++) 
            S[i*nc + j] = S[j*nc + i];
    
    std::vector<CfPt> X(nc);
    // X = A^t * P;
    for (int i = 0; i < nc; i++) {
        double sx = 0, sy = 0;
        for (int k = data.size(); k-->0;) {
            sx += A[k*nc + i] * data[k].x;
            sy += A[k*nc + i] * data[k].y;
        };
        X[i] =CfPt(sx, sy);
    };
    // psinv(A, n) = inverse of real symmetric matrix A;
    // ccmath-2.2.1 version;
    if (0 == psinv(&S[0], nc)) {
        std::vector<CfPt> Q(nc);
        for (int i = 0; i < nc; ++i) {
            souble sx = 0, sy = 0;
            for (int k = 0; k < nc; ++k) {
            	sx += S[i*nc + k] * X[k].x;
                sy += S[i*nc + k] * X[k].y;
            }
            Q[i] = CfPt(sx, sy);
        }
        return Q;  // estimated control points;
    }
    return std::vector<CfPt> ();
};

**네이버 블로그 이전;

Brute-force : $O(n^2)$,

Optimal: $O(n\log n)$

참고: people.csail.mit.edu/indyk/6.838-old/handouts/lec17.pdf

참고: arxiv.org/pdf/1911.01973.pdf

//  points should be sorted in order of increasing x-component.
//
#define DIST2(p, q) ((p).x-(q).x)*((p).x-(q).x) + ((p).y-(q).y)*((p).y-(q).y)
// 수정: index range -> [left, right]; ClosestPair(points, 0, points.size()-1, closest);
int ClosestPair(std::vector<CPoint>& points, int left, int right, int closest[2]) {
    if (right - left >= 3) {
        int lpair[2], rpair[2], min_dist2;
        int mid = left + (right - left) / 2;                    // half index;
        int ldist = ClosestPair(points, left, mid, lpair);     // left side;
        int rdist = ClosestPair(points, mid+1, right, rpair); // right side;
        if (ldist < rdist) {
            closest[0] = lpair[0]; closest[1] = lpair[1];
            min_dist2 = ldist;
        } else {
            closest[0] = rpair[0]; closest[1] = rpair[1];
            min_dist2 = rdist;
        }
        // find points which lies near center strip(2d-strip);
        // Note our distance is the squar of actual distance;
        int width = int(sqrt(double(min_dist2))+ .5);
        int ll = left;
        while (points[ll].x < (points[mid].x - width - 1)) ll++;
        int rr = idx2;
        while (points[rr].x > (points[mid + 1].x + width + 1)) rr--;
        for (int i = ll; i < rr; i++) {
            for (int j = i + 1; j <= rr; j++) {
                int dist2 = DIST2(points[i], points[j]);
                if (min_dist2 > dist2) {
                    min_dist2 = dist2;
                    closest[0] = i; closest[1] = j;
                }
            }
        }
        return min_dist2;
    } 
    else if (right == left + 2) {
        return ClosestPair3(points, left, closest);
    }
    else if (right == left + 1) {
        closest[0] = left; closest[1] = right;
        return DIST2(points[left], points[right]);
    }
    else return INT_MAX;
};