Isometric Transformation

기준 좌표계에 대해서 원점을 이동하고 좌표축을 회전시킨 새로운 좌표계에서 점의 좌표는 바뀐다. 원래의 좌표와 바뀐 좌표값 사이의 관계를 주는 변환이 Isometric transformation (isometry)이다. 평면에서 이 변환은 평행이동을 나타내는 파라미터 2개, 그리고 1개의 회전각 파라미터에 의해서 결정이 된다. 회전각이 $\theta$고, 평행이동이 $(t_x,t_y)$인 isometry에 의해서 두 점 $(x,y)$가 $(u,v)$로 연결이 되는 경우에, 아래의 식으로 표현이 된다:

$$u=\cos (\theta )x-\sin (\theta )y+t_x$$

$$v=\sin (\theta )x+\cos (\theta )y+t_y;$$

따라서 isometry로 연결이 되는 두 점의 조합 $\{(x_1,y_1)\to (u_1,v_1),(x_2,y_2)\to (u_2,v_2)\}$ 만 있으면 이들 파라미터를 정확히 결정할 수 있다. 그러나 변환에 필요한 점 정보를 얻는 과정은 필연적으로 노이즈의 영향을 받게 되므로 주어진 모든 점을 정확히 연결하는 변환을 일반적으로 구할 수 없다. 이 경우에는 isometry 파라미터는 일반적으로 최소자승법에 의해서 결정될 수 있다. 

최소자승법을 사용하기 위해서는 회전각 $\theta$보다는 $a=\cos \theta$, $b=\sin \theta$로 정의된 새로운 파라미터로 식을 표현하는 것이 더 편리하다. 그러나 이 경우에 파라미터 $a,b$는 서로 독립적이 아니고 $a^2+b^2=1$의 제한 조건을 만족시켜야 한다.  

평행이동 파라미터는 질량중심의 isometry 관계로 해결이 되므로, 이 전체 계산을 각각의 질량중심을 원점으로 하는 좌표로 옮겨서 적용하면 더 이상 평행이동을 고려할 필요 없이 회전만 계산하면 된다.

최소자승법의 원리에 따라 입력점의 isometry 결과와 대응점 사이의 거리의 제곱 합 $L$을 주어진 제약조건 내에서 최소화시키는 파라미터 $a,b,\lambda$를 찾으면 된다:

$$L=\sum _i[(a{x}_i-b{y}_i-u_i{)}^2+(b{x}_i+a{y}_i-v_i{)}^2]+\lambda (a^2+b^2-1);$$

여기서 $\lambda$는 제한 조건 $a^2+b^2=1$를 넣기 위한 Lagrange multiplier이다. 극값을 찾기 위해서 $L$를 각각 $a,b,\lambda$에 대해서 미분해서 다음 조건을 얻는다:

$$\sum _i(a{x}_i-b{y}_i-u_i)x_i+(b{x}_i+a{y}_i-v_i)y_i+\lambda a=0$$

$$\sum _i(a{x}_i-b{y}_i-u_i)(-y_i)+(b{x}_i+a{y}_i-v_i)x_i+\lambda b=0$$

$$a^2+b^2=1$$

이 식들을  $a,b,\lambda$에 대해서 풀면 다음의 관계식을 얻는다:

$$a=\sum (x_i{u}_i+y_i{v}_i)/\sum (x_i^2+y_i^2+\lambda )$$

$$b=\sum (x_i{v}_i-y_i{u}_i)/\sum (x_i^2+y_i^2+\lambda )$$
또한, Lagrange 멀티플라이어 $\lambda$는

$$A=\sum (x_i{u}_i+y_i{v}_i)$$

$$B=\sum (x_i{v}_i-y_i{u}_i);$$

로 놓으면, $a^2+b^2=1$ 에서

$$\sum (x_i^2+y_i^2+\lambda )=\sqrt{A^2+B^2};$$

임을 쓰면 된다. 따라서 회전각은

$$\cos \theta =a=\frac{A}{\sqrt{A^2+B^2}}$$

$$\sin \theta =b=\frac{B}{\sqrt{A^2+B^2}};$$

로 주어진다.

질량중심을 빼기 전 좌표 $(x,y)$의 질량중심과 $(u,v)$의 질량중심은 서로 isometry에 의해서 연결이 되므로, 이 관계에서 평행이동 파라미터 $(t_x,t_y)$는 결정이 된다:
$$(x_c,y_c)\to (u_c,v_c)$$

$$u_c=a{x}_c-b{y}_c+t_x$$

$$v_c=b{x}_c+a{y}_c+t_y;$$

참고:
** affine transformation = similarity transformation + shear;
** similarity transformation = isometry transformation + overall scaling;

/* struct CfPt { double x, y;};
*      u = T[0] * x + T[1] * y +T[4] ;
*      v = T[2] * x + T[3] * y + T[5] ; 
*/
BOOL IsometryTransform(std::vector<CfPt> &A, std::vector<CfPt> &U, double T[6]) {
    // A.size()==U.size();
    double cx = 0, cy = 0;
    double ux = 0, uy = 0;
    for (int i = A.size(); i-->0;) {
        cx += A[i].x ;  cy += A[i].y ;
        ux += U[i].x ;  uy += U[i].y ;
    };
    //center of mass ;
    cx /= A.size(); cy /= A.size();
    ux /= A.size(); uy /= A.size();

    //centering 된 좌표계에서 계산;
    double dot = 0 , cross = 0;
    for (int i = A.size(); i-->0;) {
        double x = A[i].x - cx, y = A[i].y - cy;
        double u = U[i].x - ux, v = U[i].y - uy;
        dot += (x * u + y * v);
        cross += ( x * v - y * u) ;
    };
    double norm = sqrt(dot * dot + cross * cross) ;
    double a = dot / norm ;
    double b = cross / norm ;

    T[0] = a ; T[1] = -b ; T[2] = b; T[3] = a; 
    T[4] = ux - (a * cx - b * cy) ;
    T[5] = uy - (b * cx + b * cy) ;
    return 1;
} ;