Geometry & Recognition

두 영상 사이의 perspective 변환은 8개의 매개변수 $(a,b,c,d,e,f,g,h)$에 의해서 다음 식처럼 기술이 된다. (see, http://kipl.tistory.com/86)

또는, 

따라서, 매개변수를 찾기 위해서는 두 영상에서 서로 대응하는 점이 4개 이상 주어져야 한다. N개의 대응점들이 주어진 경우

각각의 대응점을 위의 식에 대입해서 정리하면 아래의 행렬식을 얻을 수 있다.(좌변 행렬의 마지막 열은 전부 - 부호가 들어가야 한다) 
 

또는, 간단히 

$$\mathbf{A}\cdot \mathbf{x}\mathbf{=}\mathbf{b}$$

로 쓸 수 있다. 그러나 대응점을 찾을 때 들어오는 noise로 인해서 실제 데이터를 이용하는 경우에는 정확히 등호로 주어지지 않는다. 따라서, 실제 문제에서는 좌변과 우변의 차이의 제곱을 최소로 만드는 $\mathbf{x}$를 찾아야 할 것이다.

$${\mathbf{x}}^{\ast }=\underset{\mathbf{x}}{\text{argmin }}|\mathbf{A}\cdot \mathbf{x}-\mathbf{b}{|}^2.$$

최소자승해를 찾기 위해 ${\mathbf{x}}^{\mathbf{T}}$에 대해 미분을 하면

$$\mathbf{(}{\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{A}\mathbf{)}\cdot \mathbf{x}\mathbf{=}{\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{b}\mathbf{,}$$

를 얻고, 이 식을 풀어서 ${\mathbf{x}}^{\ast }$을 구하면 된다. ${\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{A}$는 $8\times 8$의 대칭 행렬로 역행렬을 구할 수 있다 (주어진 점들 중 한 직선 위에 놓이지 않는 점이 4개 이상이 있어야 한다). 따라서 최소자승해는 다음과 같이 쓸 수 있다:

$${\mathbf{x}}^{\mathbf{\ast }}\mathbf{=}\mathbf{(}{\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{A}{\mathbf{)}}^{\mathbf{-}\mathbf{1}}\cdot \mathbf{(}{\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{b}\mathbf{)}\mathbf{.}$$

The Savitzky–Golay method essentially performs a local polynomial regression (of degree $k$) on a series of values (of at least $k+1$ points which are treated as being equally spaced in the series) to determine the smoothed value for each point.

Savitzky–Golay 필터는 일정한 간격으로 주어진 데이터들이 있을 때(이들 데이터는 원래의 정보와 노이즈를 같이 포함한다), 각각의 점에서 주변의 점들을 가장 잘 피팅하는 $k$-차의 다항식을 최소자승법으로 찾아서 그 지점에서의 출력값을 결정하는 필터이다. 이 필터는 주어진 데이터에서의 극대나 극소, 또는 봉우리의 폭을 상대적으로 잘 보존한다.(주변 점들에 동등한 가중치를 주는 Moving Average Filter와 비교해 볼 수 있다).

간단한 예로, 2차의 다항식과 5개의 데이터 점

$$\{(-2,d_0),(-1,d_1),(0,d_2),(1,d_3),(2,d_4)\}$$

을 이용해서 중앙에서의 값을 결정하는 방법을 살펴보자. 사용하려고 하는 다항식은

$$p(x)=a_0+a_{1}x+a_2{x}^2$$

이다. 다항식의 계수는 다항식의 값과 실제 데이터의 값과의 차이를 최소화시키도록 선택해야 한다. 즉, 최소자승의 원리를 적용하여서 구하면 된다. 계산된 다항식의 값과 실제 데이터 값 사이의 차의 제곱을 구하면:

$$L=|a_0-2a_1+4a_2-d_0{|}^2+|a_0-a_1+a_2-d{|}^2+|a_0-d_0{|}^2$$ $$+|a_0+a_1+a_2-d_3{|}^2+|a_0+2a_1+4a_2-d_4{|}^2$$

이 식의 $a_0,a_1,a_2$에 대한 극값을 가질 조건은($\partial L/\partial a_i=0$) $$5a_0+10a_2=d_0+d_1+d_2+d_3+d_4$$$$10a_1=-2d_0–d_1+d_3+2d_4$$$$10a_0+34a_2=4d_0+d_1+d_3+4d_4$$

이 식을 만족시키는 $a_0$를 구하면, 필터의 출력(원도 중앙에서 값)이 결정된다.

$$\text{필터 출력:}{a}_0=(-3d_0+12d_1+17d_2+12d_3-3d_4)/35$$

위에서 계수 $a_0$, $a_1$, $a_2$를 결정하는 방정식은 행렬로 정리하면 아래의 식과 같이 표현할 수 있다. 

좌변의 5행 3열 행렬을 $\mathbf{A}$, $\mathbf{a}=[a_0,a_1,a_2{]}^T$, $\mathbf{d}=[d_0,d_1,d_2,d_3,d_4{]}^T$로 놓으면, 이 행렬방정식은 $\mathbf{A}\mathbf{.}\mathbf{a}\mathbf{=}\mathbf{d}$ 형태로 쓸 수 있다. $\mathbf{A}$(design matrix)가 정방행렬이 아니므로 역행렬을 바로 구할 수 없지만, $|\mathbf{A}\cdot \mathbf{a}\mathbf{-}\mathbf{d}{|}^2$을 최소로 하는 최소제곱해는
$$\mathbf{(}{\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{A}\mathbf{)}\cdot \mathbf{a}\mathbf{=}{\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{d}$$를 만족시켜야 하므로 

$$\mathbf{a}\mathbf{=}\mathbf{(}{\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{A}{\mathbf{)}}^{\mathbf{-}\mathbf{1}}\cdot \mathbf{(}{\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{d}\mathbf{)}$$로 주어짐을 알 수 있다.

이 식은 임의의 $k$-차 다항식을 이용한 경우에도 사용할 수 있다. 이 경우 행렬(scattering matrix) ${\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{A}$는 $(k+1)\times (k+1)$의 대칭행렬이 된다. 행렬 $\mathbf{A}$는 다항식의 찻수와 피팅에 사용이 될 데이터의 구간의 크기가 주어지면 정해지므로, 윗 식에서 $({\mathbf{A}}^T\cdot \mathbf{A}{)}^{-1}\cdot {\mathbf{A}}^T$의 첫 행 ($a_0$을 $d$로 표현하는 식의 계수들)을 구하면 코드 내에서 결과를 lookup table로 만들어서 사용할 수 있다. 아래 표는 mathematica 를 이용해서 윈도 크기가 7 (7개 점)인 경우 2차 다항식을 사용할 때 계수를 구하는 과정이다.

2차 다항식일 때, 같은 방식으로 다양한 윈도 크기에 따른 계수를 구할 수 있다.
*크기($n$)에 따른 필터값 결정계수 (중앙에 대해 좌우대칭이다);

$$\begin{array}{rl}n=5;& W[]=\{-3,12,\mathbf{17},12,-3\};\\ n=7;& W[]=\{-2,3,6,\mathbf{7},6,3,-2\};\\ n=9;& W[]=\{-21,14,39,54,\mathbf{59},54,39,14,-21\};\\ n=11;& W[]=\{-36,9,44,69,84,\mathbf{89},84,69,44,9,-36\};\\ n=13;& W[]=\{-11,0,9,16,21,24,25,24,21,16,9,0,-11\}\\ n=15;& W[]=\{-78,-13,42,87,122,147,162,167,162,147,122,87,42,-13,-78\}\\ n=25;& W[]=\{-253,-138,-33,62,147,222,287,342,387,422,447,462,467,\\ & 462,447,422,387,342,287,222,147,62,-33,-138,-253\}\end{array}$$

$$\text{필터 출력}=\frac{\sum _{i}W[i]d[i]}{\sum _{i}W[i]}$$

std::vector<double> SavitzkyGolayFilter(std::vector<double>& data, double W[], int wsz) {
    const int hwsz = wsz >> 1;
    wsz = (hwsz<<1) + 1;
    std::vector<double> padded(data.size() + 2 * hwsz);
    // reflective boundary conditions;   
    // [hw]...[1][0][1][2]...[hw]....[n-1+hw][n-2+hw][n-3+hw]....[n-1];
    for (int i = 0; i < hwsz; i++) { 
        padded[i]                  = data[hwsz-i];
        padded[i+data.size()+hwsz] = data[data.size()-2-i];
    }
    for (int i = data.size(); i-->0;) padded[i+hwsz] = data[i];
    //
    std::vector<double> smoothed(data.size());
    double wsum = 0;
    for (int i = 0; i < wsz; i++) wsum += W[i]; 
    for (int i = data.size(); i-->0;) {
        double *ppad = &padded[i];
        double fsum = 0;
        for (int k = 0; k < wsz; k++) fsum += ppad[k] * W[k];
        smoothed[i] = fsum / wsum;
    }
    return smoothed;
};

$$f(x)=1-(|x|+0.01{)}^{0.01}+N(0,0.002),-2\le x\le 2$$

일반적인 찻수(order)에 대한 필터 계수 계산;

std::vector<double> sgFilterCoeffs(int width, int order/* = 2*/) {
    // ASSERT(width >= order+1);
    int hwidth = width >> 1;
    width = (hwidth << 1) + 1; // width = odd number;
    // design matrix;
    std::vector<double> D((order+1) * width);
    for (int i = 0; i < width; i++) {
        double x = double(i) - hwidth;  // -hwidth <= x <= hwidth;
        double *pD = &D[i*(order+1)];
        pD[0] = 1;
        for (int k = 1; k <= order; k++) 
            pD[k] = pD[k-1] * x;
    };
    // scattering matrix;
    std::vector<double> S((order+1)*(order+1)); // S=~D.D;
    for (int i = 0; i <= order; i++)
        for (int j = i; j <= order; j++) {
            double s = 0;
            for (int k = 0; k < width; k++)
                s += D[k*(order+1) + i] * D[k*(order+1) + j];
            S[i*(order+1) + j] = s;
        }
    for (int i = 0; i <= order; i++) 
        for (int j = 0; j < i; j++) 
            S[i*(order+1) + j]= S[j*(order+1) + i];
			
    // ccmath library;
    psinv(&S[0], order+1);
    // filter coeffs = 0-th row fo S.~D;
    std::vector<double> filter;
    for (int i = 0; i < width; i++) {
        double s = 0;
        for (int k = 0; k <= order; k++)
            s += S[k] * D[i*(order+1) + k];
        filter.push_back(s);
    }
    return filter;
}

물체의 형상은 폴리곤이나 폴리곤의 집합으로 근사적으로 표현할 수 있다. 예를 들면 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^{\prime }=\{(x_1^{\prime },y_1^{\prime }),(x_2^{\prime },y_2^{\prime }),...,(x_N^{\prime },y_N^{\prime })\}$이 affine 변환에 의해서 연결이 되는 경우에 각 꼭짓점 사이의 관계는

$$\begin{array}{rl}{x}_i^{\prime }& =a{x}_i+b{y}_i+t_x\\ y_i^{\prime }& =c{x}_i+d{y}_i+t_y,(i=1,2,...,N);\end{array}$$

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

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

$$L=\sum _i|x_i^{\prime }-a{x}_i-b{y}_i-t_x{|}^2+|y_i^{\prime }-c{x}_i-d{y}_i-t_y{|}^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}_{x{x}^{\prime }}& S_{x{y}^{\prime }}\\ S_{y{x}^{\prime }}& S_{y{y}^{\prime }}\\ S_{x^{\prime }}& S_{y^{\prime }}\end{array}\right]$$

여기서,
$$\begin{array}{rl}& S_{xx}=\sum x^2,S_{yy}=\sum y^2,S_{xy}=\sum xy,\\ & S_x=\sum x,S_y=\sum y,S_{x^{\prime }}=\sum x^{\prime },S_{y^{\prime }}=\sum y^{\prime }\\ & S_{x{x}^{\prime }}=\sum x{x}^{\prime },S_{x{y}^{\prime }}=\sum x{y}^{\prime },S_{y{x}^{\prime }}=\sum y{x}^{\prime }\end{array}$$ 이다.

// 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을 만들 수 있다. 물론 추가된 녹색점이 신뢰할 수 있는 것인가에 대한 판단을 하기 위해서는 추가적인 정보가 더 요구된다.