Geometry & Recognition

물체의 형상은 폴리곤이나 폴리곤의 집합으로 근사적으로 표현할 수 있다. 예를 들면 snake나 active shape model (ASM) 등에서 손 모양이나 얼굴의 윤곽, 또는 의료 영상 등에서 장기의 모양 등을 표현할 때 사용이 된다. 이러한 응용에서 주어진 형상을 기준으로 주어진 형상에 정렬을 시켜야 필요가 생긴다. 일반적으로 카메라를 써서 얻은 각 영상에서 추출한 정보들 사이에는 서로 사영 변환의 관계로 연결된다. 그러나 많은 경우에는 in-plane 변형만 고려해도 충분할 때가 많다. 이 경우에 가장 일반적인 형상의 변형은 affine 변환으로 표현된다. 회전(rotation), 평행 이동(translation), 크기 변환(scale transformation) 그리고 층 밀림(shear)을 허용하는 변환이다. 물론, 간단한 경우로는 shear를 제외할 수도 있고 (similarity transformation), 더 간단하게는 크기 변환을 제외할 수도 있다 (isometric transformation).

$N$개의 꼭짓점을 갖는 두 개의 형상 $S=\{(x_1, y_1), (x_2, y_2),..., (x_N, y_N) \}$, $S'=\{(x'_1, y'_1), (x'_2, y'_2),..., (x'_N, y'_N) \}$이 affine 변환에 의해서 연결이 되는 경우에 각 꼭짓점 사이의 관계는

\begin{align} x'_i &= a x_i  + b y_i + t_x \\ y'_i &= c x_i + d y_i + t_y, \quad (i=1,2,..., N);\end{align}

의 6개의 매개변수$(a, b, c, d, t_x, t_y)$에 의해서 기술이 된다(평행 이동: $x/y$축 방향 2개, 회전: 1개, shear: 1개, 스케일: $x/y$축 방향 2개). Affine 변환에 의해서 평행인 두 직선은 변환 후에도 평행인 관계를 유지한다.

꼭짓점 위치는 실제로 다양한 영상처리 과정에 의해서 얻어지므로 필연적으로 노이즈를 포함하게 되어서 일종의 랜덤 변수로 생각해야 한다. 주어진 랜덤 변수에서 최적으로 매개변수를 추출하기 위해 최소자승법을 이용한다. Affine 변환된 좌표와 실제 측정된 좌표 사이의 거리 차이를 최소화하는 매개변수를 찾도록 하자:

$$L=\sum_i \big| x'_i - a x_i - b y_i - t_x \Big|^2 + \big| y'_i - c x_i -d y_i - t_y\big|^2 $$

Affine변환을 규정하는 매개변수를 구하기 위해서는 L을 각 매개변수에 대해서 미분해서 극값을 가질 조건을 구하면 된다:

        ∂L/∂a = -2 * ∑ (x'_i - a * x_i - b * y_i - t_x) * x_i ;
        ∂L/∂b = -2 * ∑ (x'_i - a * x_i - b * y_i - t_x) * y_i ;
        ∂L/∂c = -2 * ∑ (y'_i - c * x_i - d * y_i - t_y) * x_i ;
        ∂L/∂d = -2 * ∑ (y'_i - c * x_i - d * y_i - t_y) * y_i ; 
        ∂L/∂t_x = -2 * ∑ (x'_i - a * x_i - b * y_i - t_x) ;
        ∂L/∂t_y = -2 * ∑ (y'_i - c * x_i - d * y_i - t_y); 

각 식을 0으로 놓아서 얻어지는 연립방정식을 행렬식으로 다시 정리하면,

$$\left[\begin{array}{ccc} S_{xx} & S_{xy} & S_x \\ S_{xy} & S_{yy} & S_y \\ S_x & S_y & N \end{array}\right]\left[ \begin{array}{ll} a & c \\ b & d\\ t_x & t_y \end{array} \right] = \left[\begin{array}{cc} S_{xx'} & S_{x y'} \\ S_{y x'} & S_{yy'} \\ S_{x'} & S_{y'}\end{array} \right]$$

여기서,
\begin{align} & S_{xx}= ∑ x^2, ~S_{yy} = ∑ y^2, ~S_{xy} = ∑ xy, \\ &S_x = ∑ x, ~S_y = ∑ y, ~S_{x'} = ∑ x', ~S_{y'} = ∑ y' \\ & S_{xx'} = ∑ xx', ~S_{xy'} = ∑ xy', ~S_{yx'} =∑ yx' \end{align} 이다.

// dst = (A,T)src;
//  [u]  = [ A0 A1 ][x] + A4
//  [v]  = [ A2 A3 ][y] + A5
//
BOOL GetAffineParameter(const std::vector<CPoint> &srcPts, 
                        const std::vector<CPoint> &dstPts, 
                        double AT[6]) 
{
    double Sx, Sy, Sxx, Sxy, Syy;
    double Su, Sv, Sxu, Sxv, Syu, Syv ;
    double A[9], invA[9];
    Sx = Sy = Sxx = Sxy = Syy = 0;
    Su = Sv = Sxu = Sxv = Syu = Syv = 0;
    for (int i = srcPts.size(); i-->0;) {
        double x = srcPts[i].x, y = srcPts[i].y ;
        double u = dstPts[i].x, v = dstPts[i].y ;
        Sx += x;        Sy += y ;
        Sxx += (x * x); Sxy += (x * y); Syy += (y * y);
        Su += u;        Sv += v ;
        Sxu += (x * u); Sxv += (x * v); Syu += (y * u); Syv += (y * v);
    }
    A[0] = Sxx; A[1] = Sxy; A[2] = Sx;
    A[3] = Sxy; A[4] = Syy; A[5] = Sy;
    A[6] = Sx ; A[7] = Sy ; A[8] = srcPts.size() ;
    double det = (A[0]*(A[4]*A[8]-A[5]*A[7])-\
                  A[1]*(A[3]*A[8]-A[5]*A[6])+\
                  A[2]*(A[3]*A[7]-A[4]*A[6]));
    if (det != 0.) {
        det = 1. / det; 
        invA[0] = (A[4]*A[8] - A[5]*A[7]) * det;
        invA[1] = (A[2]*A[7] - A[1]*A[8]) * det;
        invA[2] = (A[1]*A[5] - A[2]*A[4]) * det;
        invA[3] = (A[5]*A[6] - A[3]*A[8]) * det;
        invA[4] = (A[0]*A[8] - A[2]*A[6]) * det;
        invA[5] = (A[2]*A[3] - A[0]*A[5]) * det;
        invA[6] = (A[3]*A[7] - A[4]*A[6]) * det;
        invA[7] = (A[1]*A[6] - A[0]*A[7]) * det;
        invA[8] = (A[0]*A[4] - A[1]*A[3]) * det;
    }
    else return FALSE;

    AT[0] = invA[0] * Sxu + invA[1] * Syu + invA[2] * Su;
    AT[1] = invA[3] * Sxu + invA[4] * Syu + invA[5] * Su;
    AT[4] = invA[6] * Sxu + invA[7] * Syu + invA[8] * Su;
    AT[2] = invA[0] * Sxv + invA[1] * Syv + invA[2] * Sv;
    AT[3] = invA[3] * Sxv + invA[4] * Syv + invA[5] * Sv;
    AT[5] = invA[6] * Sxv + invA[7] * Syv + invA[8] * Sv;
    return TRUE ;
};

아래의 그림은 지문에서 얻은 특징점을 가지고 변환을 한 것이다. 밑에 그림이 기준 template (붉은 점)이고 윗 그림은 이 기준  template와 입력된 지문의 특징점(노란 점+ 녹색점) 사이에 서로 메칭이 되는 특징점(노란색)을 찾고, 그것을 기준으로 두 지문 영상 간의 affine 파라미터를 찾아서 기준 template을 변환시킨 것이다. 이렇게 하면 새로 찾은 특징점 중에서 기준 template에 없는 특징점(녹색점)을 발견할 수 있고, 이 특징점을 기준 template에 추가하여서 좀 더 넓은 범위를 커버할 수 있는 template을 만들 수 있다. 물론 추가된 녹색점이 신뢰할 수 있는 것인가에 대한 판단을 하기 위해서는 추가적인 정보가 더 요구된다.

기준 좌표계에 대해서 원점을 이동하고 좌표축을 회전시킨 새로운 좌표계에서 점의 좌표는 바뀐다. 원래의 좌표와 바뀐 좌표값 사이의 관계를 주는 변환이 Isometric transformation (isometry)이다. 평면에서 이 변환은 평행이동을 나타내는 파라미터 2개, 그리고 1개의 회전각 파라미터에 의해서 결정이 된다. 회전각이 $θ$고, 평행이동이 $(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) \rightarrow(u_1, v_1), (x_2, y_2)\rightarrow(u_2, v_2)\}$ 만 있으면 이들 파라미터를 정확히 결정할 수 있다. 그러나 변환에 필요한 점 정보를 얻는 과정은 필연적으로 노이즈의 영향을 받게 되므로 주어진 모든 점을 정확히 연결하는 변환을 일반적으로 구할 수 없다. 이 경우에는 isometry 파라미터는 일반적으로 최소자승법에 의해서 결정될 수 있다. 

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

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

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

$$L = \sum_i \big [ (a  x_i - b  y_i - u_i)^2 + (b  x_i + a  y_i - v_i)^2 \big] + λ  (a^2 + b^2 - 1) ;$$

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

$$\sum_i  ( a  x_i - b  y_i - u_i) x_i + ( b  x_i + a  y_i - v_i) y_i + λ a = 0 ;\\ \sum_i  ( a  x_i - b  y_i - u_i) (-y_i) + ( b  x_i + a  y_i - v_i) x_i + λ b = 0;\\ a^2 + b^2 = 1; $$

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

$$a = ∑(x_i u_ i + y_ i v_ i) / ∑ (x_ i^2 + y_i^2 + λ) ; \\ b = ∑ (x_i v_ i - y_i u_i) / ∑ (x_i^2 + y_i^2 + λ); $$
또한, Lagrange 멀티플라이어 $λ$는

$$A  = ∑ (x_i u_i + y_i v_i); \\  ~B =  ∑ (x_i v_i - y_i u_i);$$

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

$$∑ ( x_i^2 + y_i^2 + λ ) = \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) \rightarrow (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;
} ;

2차원 이미지의 기하학적인 변형 중에서 평행이동, 회전 및 전체적인 크기의 변화를 주는 변환이 similarity transformation이다. 이 변환은 두 직선이 이루는 각을 보존하고 길이 비를 유지한다. 따라서 similarity 변환 후 물체의 모양은 변환 전과 같은 형태를 가진다. 이 변환보다도 더 일반적인 2차원의 기하학적인 변환은 affine transformation이다. Affine 변환은 한쪽 방향으로의 밀림(sheer)도 허용한다. 평행한 두 직선은 affine 변환 후에도 여전히 평행하다.

Similarity transformation은 전체적인 크기를 바꾸는 scale parameter($s$) 1개와 회전각($θ$) 1개, 그리고 $x, y$축으로의 평행이동을 나타내는 parameter ($t_x$, $t_y$) 2 개를 합해서 총 4개가 있어야 한다. 이 parameter에 의해서 원본 이미지의 픽셀 $(x, y)$가 변환된 이미지의 픽셀 $(u, v)$에 대응한다고 하면, 이들 간의 관계는 다음식으로 주어진다. $$u =  s\cos (θ) x - s \sin (θ) y + t_x;$$ $$v =  s \sin (θ) y + s \cos (θ) y + t_y;$$ 따라서 원본 영상의 2점에 대응하는 정보만 주어지면 파라미터 $(s, θ, t_x, t_y)$를 유일하게 결정할 수 있다.     $$(x_1, y_1) \rightarrow  (u_1, v_1),\\ (x_2 , y_2)  \rightarrow (u_2, v_2) $$그러나 많은 경우에는 기준점을 잡는데 에러 등을 고려하여서 일반적으로 원본 영상의 $N(\ge 2)$ 개의 점에 대응하는 정보를 주게 되는데, 이 경우에 변환 관계식은 overdetermined 되어서 해를 구할 수 없는 경우도 있다. 이 경우에는 최소자승법을 써서 변환점과 변환식에 의해서 의해서 주어지는 값의 차이를 최소화시키는 파라미터를 구해서 쓰면 된다.$$L =  \sum_{i} | u_i - (s\cos(θ) x_i - s \sin(θ) y_i + t_x)|^2 + |v_i - (s \sin(θ) x_i + s \cos(θ) y_i + t_y)|^2, \\ (s, \theta, t_x, t_y) =\text {argmin}(L);$$

이 식을 최소화시키는 파라미터는 $(a= s \cos(θ), b=s \sin(θ)$로 놓으면) $a, b, t_x, t_y$에 대해서 극값을 가질 조건에서 얻을 수 있다. $$\frac {\partial L}{\partial a}=0: \quad \sum_{i} (u_i - (ax_i - by_i + t_x))(-x_i) + (v_i - (bx_i + ay_i + t_y))(-y_i) = 0,\\ \frac {\partial L}{\partial b}=0:\quad \sum _{i} (u_i - (ax_i - by_i + t_x))(y_i) + (v_i - (bx_i + a y_i + t_y))(-x_i) = 0, \\ \frac {\partial L}{\partial t_x}=0: \quad \sum_{i} (u_i - (ax_i - by_i + t_x)) = 0, \\ \frac {\partial L}{\partial t_y}=0: \quad  \sum_{i} (vi - (bx_i + ay_i + t_y)) = 0.$$

따라서, $S_u = \sum_i  u_i$, $S_v = \sum_i v_i$, $S_{ux} = \sum _i  u_i x_ i$, $S_{uy} = \sum _i  u_iy_i$, $S_{vx} = \sum_i  v_i x_i$, $S_{vy} = \sum _i v_i y_i$, $S_x = \sum  x_i$, $S_y=\sum _i y_i$, $S_{xx} = \sum_i  x_i^2$, $S_{xy} = \sum_i x_iy_i$, $S_{yy}=\sum_i y_i^2$라고 하면,$$-S_{ux}  + a   S_{xx} + t_x  S_x - S_{vy} + a S_{yy} + t_y S_y = 0; \\ S_{uy} + b  S_{yy} - t_x  S_y -S_{vx} + b S_{xx}  + t_y S_x = 0;\\ S_u - a S_x + bS_y - t_x  N = 0; \\ S_v - b S_x - aS_y - t_y  N = 0;$$의 4개의 식을 얻으므로 $(a, b, t_x, t_y)$에 대한 1차 연립방정식을 풀면 된다.

$$\left [\begin {array}{cccc} S_x&-S_y&N&0\\S_y &S_x&0&N\\ S_{xx}+S_{yy}&0&S_x&S_y\\0 &S_{xx}+S_{yy}&-S_y&S_x\end {array} \right]\left [\begin {array}{c} a\\b\\t_x\\t_y \end {array}\right]=\left [\begin {array}{c} S_u\\S_v\\S_{ux} +S_{vy}\\S_{vx}-S_{uy}\end {array}\right] \\   \text{or}~~ \mathbf{A}\cdot \mathbf{x} = \mathbf{b}$$ $4\times 4$ 행렬 $\mathbf{A}$의 역행렬은 다음과 같이 쉽게 구해진다. 

$$ \mathbf{A}^{-1} = \frac{1}{S_x^2+S_y^2 -N(S_{xx}+S_{yy})}\begin{bmatrix}  S_x & S_y & -N & 0 \\ -S_y & S_x & 0 & -N \\ -(S_{xx}+S_{yy}) & 0 & S_x & -S_y \\ 0 & -(S_{xx}+S_{yy}) & S_y & S_x  \end{bmatrix} $$

아래의 코드는 이것을 구현한 것이다. 물론, $N=2$개인 경우에는 파라미터는 유일하게 정해지고 이보다도 더 간단한 식으로 주어진다.

// dst = (S|T)(src)
BOOL SimilarTransParams(std::vector<CPoint>& src, std::vector<CPoint>& dst, double ST[4]) {
    double Sx = 0, Sy = 0, Sxx = 0, Syy = 0;
    double Su = 0, Sv = 0, Sxu = 0, Sxv = 0, Syu = 0, Syv = 0;
    for (int i = srcPts.size(); i-->0;) {
        double x = src[i].x, y = src[i].y;
        double u = dst[i].x, v = dst[i].y;
        Sx  += x;        Sy  += y;
        Sxx += (x * x);  Syy += (y * y);
        Su  += u;        Sv  += v;
        Sxu += (x * u);  Syv += (y * v);
    }
    double Z = Sxx + Syy;
    double denorm = Sx * Sx + Sy * Sy - src.size() * Z;
    // det = (denorm)^2;
    if (denorm == 0) return FALSE;
    invA[16] = { Sx, Sy, -src.size(),           0,
                -Sy, Sx,           0, -src.size(),
                -Z,   0,          Sx,         -Sy,
                 0,  -Z,          Sy,          Sx};
    for (int i = 0; i < 16; i++) invA[i] /= denorm;
    //
    double b[4] = {Su, Sv, Sxu + Syv, Sxv - Syu};
    for (int i = 0; i < 4; i++) {
    	double s = 0;
        for (int j = 0; j < 4; j++) 
            s += invA[i * 4 + j] * b[j];
        ST[i] = s;
    }
    return TRUE ;
};

InvertMatrix4x4()는 4x4행렬의 역행렬을 구한다(OpenCV에서)

BOOL InvertMatrix4x4_d(double* srcMatr, double* dstMatr) {
    double di = srcMatr[0];
    double d = 1.0 / di;

    dstMatr[0] = d;
    dstMatr[4] = srcMatr[4] * -d;
    dstMatr[8] = srcMatr[8] * -d;
    dstMatr[12] = srcMatr[12] * -d;
    dstMatr[1] = srcMatr[1] * d;
    dstMatr[2] = srcMatr[2] * d;
    dstMatr[3] = srcMatr[3] * d;
    dstMatr[5] = srcMatr[5] + dstMatr[4] * dstMatr[1] * di;
    dstMatr[6] = srcMatr[6] + dstMatr[4] * dstMatr[2] * di;
    dstMatr[7] = srcMatr[7] + dstMatr[4] * dstMatr[3] * di;
    dstMatr[9] = srcMatr[9] + dstMatr[8] * dstMatr[1] * di;
    dstMatr[10] = srcMatr[10] + dstMatr[8] * dstMatr[2] * di;
    dstMatr[11] = srcMatr[11] + dstMatr[8] * dstMatr[3] * di;
    dstMatr[13] = srcMatr[13] + dstMatr[12] * dstMatr[1] * di;
    dstMatr[14] = srcMatr[14] + dstMatr[12] * dstMatr[2] * di;
    dstMatr[15] = srcMatr[15] + dstMatr[12] * dstMatr[3] * di;
    di = dstMatr[5];
    dstMatr[5] = d = 1.0 / di;
    dstMatr[1] *= -d;
    dstMatr[9] *= -d;
    dstMatr[13] *= -d;
    dstMatr[4] *= d;
    dstMatr[6] *= d;
    dstMatr[7] *= d;
    dstMatr[0] += dstMatr[1] * dstMatr[4] * di;
    dstMatr[2] += dstMatr[1] * dstMatr[6] * di;
    dstMatr[3] += dstMatr[1] * dstMatr[7] * di;
    dstMatr[8] += dstMatr[9] * dstMatr[4] * di;
    dstMatr[10] += dstMatr[9] * dstMatr[6] * di;
    dstMatr[11] += dstMatr[9] * dstMatr[7] * di;
    dstMatr[12] += dstMatr[13] * dstMatr[4] * di;
    dstMatr[14] += dstMatr[13] * dstMatr[6] * di;
    dstMatr[15] += dstMatr[13] * dstMatr[7] * di;
    di = dstMatr[10];
    dstMatr[10] = d = 1.0 / di;
    dstMatr[2] *= -d;
    dstMatr[6] *= -d;
    dstMatr[14] *= -d;
    dstMatr[8] *= d;
    dstMatr[9] *= d;
    dstMatr[11] *= d;
    dstMatr[0] += dstMatr[2] * dstMatr[8] * di;
    dstMatr[1] += dstMatr[2] * dstMatr[9] * di;
    dstMatr[3] += dstMatr[2] * dstMatr[11] * di;
    dstMatr[4] += dstMatr[6] * dstMatr[8] * di;
    dstMatr[5] += dstMatr[6] * dstMatr[9] * di;
    dstMatr[7] += dstMatr[6] * dstMatr[11] * di;
    dstMatr[12] += dstMatr[14] * dstMatr[8] * di;
    dstMatr[13] += dstMatr[14] * dstMatr[9] * di;
    dstMatr[15] += dstMatr[14] * dstMatr[11] * di;
    di = dstMatr[15];
    dstMatr[15] = d = 1.0 / di;
    dstMatr[3] *= -d;
    dstMatr[7] *= -d;
    dstMatr[11] *= -d;
    dstMatr[12] *= d;
    dstMatr[13] *= d;
    dstMatr[14] *= d;
    dstMatr[0] += dstMatr[3] * dstMatr[12] * di;
    dstMatr[1] += dstMatr[3] * dstMatr[13] * di;
    dstMatr[2] += dstMatr[3] * dstMatr[14] * di;
    dstMatr[4] += dstMatr[7] * dstMatr[12] * di;
    dstMatr[5] += dstMatr[7] * dstMatr[13] * di;
    dstMatr[6] += dstMatr[7] * dstMatr[14] * di;
    dstMatr[8] += dstMatr[11] * dstMatr[12] * di;
    dstMatr[9] += dstMatr[11] * dstMatr[13] * di;
    dstMatr[10] += dstMatr[11] * dstMatr[14] * di;
    return TRUE;
}

2개의 대응점만 주어진 경우 $(x_1, y_1), (x_2, y_2) \rightarrow (u_1, v_1), (u_2, v_2)$;

bool SimilarTransParams(double x1, double y1, double x2, double y2, 
                        double u1, double v1, double u2, double v2,
                        double ST[4]) {
    double x21 = x2 - x1, y21 = y2 - y1;
    double u21 = u2 - u1, v21 = v2 - v1;
    double det = x21 * x21 + y21 * y21;
    if (det == 0.) return false;
    double a = (x21 * u21 + y21 * v21) / det ;
    double b = (x21 * v21 - y21 * u21) / det ;
    double tx = u1 - a * x1 + b * y1;
    double ty = v1 - b * x1 - a * y1;
    ST[0] = a; ST[1] = b; ST[2] = tx; ST[3] = ty;
    return true;
};

얼굴인식용 training data set을 만들기 위해서 얼굴을 정렬시키는 데 사용한 예:
- 양 눈의 위치 변환: (70,93), (114, 84) --> (30,45), (100,45)로 변환( linear interpolation사용)
- 실제로 사용되는 변환은 정해진 dst영역으로 매핑하는 src영역을 찾아야 하므로, 역변환이 필요하다.
- 필요한 역변환은 src와 dst의 역할만 바꾸면 쉽게 구할 수 있다.