최소 자승법 문제는 추정치와 실제의 측정치와의 차이를 최소로 만드는 parameter vector를 구하는 것이다.
$$\underset{\mathbf{x}} {\text{argmin}} ~|\mathbf{A}. \mathbf{x} - \mathbf{b}|^2.$$

여기서 design matrix $\mathbf{A}$는 추정에 사용이 된 basis 함수를 각각의 독립변수에서 계산한 값이고, $\mathbf{x}$는 구하고자 하는 parameter vector이며, $\mathbf{b}$는 측정값 vector이다. 예를 들면 주어진 측정값을 $n-1$차의 다항식을 이용하여서 피팅하려고 하는 경우에는 parameter는 다항식의 계수가 되며 $n$-차원의 vector space을 형성하게 되며, $\mathbf{A}$는

$$   A_{ij} = (x_i)^j ,  \quad j=0,..., n-1$$

가 일 것이다. 일반적으로 $A_{ij}$는 $x_i$에서 계산된 basis-함수의 값이 된다. 위의 식을 $\mathbf{x}$에 대해서 미분을 하면 극값을 취하면 최소자승 문제는 아래의 행렬식을 푸는 문제가 된다

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

$\mathbf{A}^T. \mathbf{A}$ 은 $n\times n$ matrix다. 이 행렬이 역행렬을 가지게 되면

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

를 하여서 원하는 parameter vector를 얻을 수 있다. 그러나 피팅 문제에 따라 행렬 $\mathbf{A}^T. \mathbf{A}$가 매우 singular 해져 역행렬을 구할 수 없게 되는 경우에 종종 생긴다. 예를 들면, 저주파의 신호를 고주파 기저 함수를 이용하여서 최소자승법을 사용하고자 하는 경우 등에 이러한 문제에 부딪히게 된다. 이런 경우에는 직접적으로 $\mathbf{A}^T. \mathbf{A}$의 역행렬을 구하는 방법을 이용할 수 없고

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

의 식을 $\mathbf{A}$의 SVD(Singular Value Decomposition)를 이용하여서 풀 수가 있다. $\mathbf{A}$를 SVD 하면 $\mathbf{A}_{m\times n}=\mathbf{U}_{m\times n} . \mathbf{w}_{n\times n}. \mathbf{V}_{n\times n}^T $의 형태로 분해할 수 있다. 여기서 $\mathbf{w}=\text{diag}(\underbrace{w_0, w_1,...}_{\text{nonzero}},0,..,0)$로 쓰여지는 대각행렬이다. matrix $\mathbf{U}$와 $\mathbf{V}$의 column vector를 사용하면
$$ \mathbf{A}  =\sum_{w_k \ne 0} w_k \mathbf{u}_k \otimes \mathbf{v}_k^T$$

의 형태로 쓰인다. $\mathbf{u}_k$는 $\mathbf{U}$의 $k$-번째 열벡터이고, $\mathbf{v}_k$는 $\mathbf{V}$의 $k$-번째 열벡터로 각각 orthonormal basis를 형성한다. parameter 벡터를 $\{ \mathbf{v}_k \}$ basis로 전개를 하면 영이 아닌 singularvalue에 해당하는 성분만 가지게 된다. 구체적으로 위의 $\mathbf{A}$ 분해와 $\mathbf{u}_j^T.\mathbf{u}_k=\delta_{jk}$, 그리고 $\sum_k \mathbf{v}_k \otimes \mathbf{v}_k^T= \mathbf{I}_{n\times n}$임을 이용하면,

\begin{gather}  \mathbf{v}_k^T . \mathbf{x} = \mathbf{u}_k^T . \mathbf{b} / w_k, \quad w_k \ne 0, \\                    \mathbf{v}_k^T . \mathbf{x} = 0, \quad w_k = 0, \\  \rightarrow ~~\mathbf{x} = \sum _{w_k \ne 0 } ( \mathbf{u}_k^T . \mathbf{b} / w_k)  \mathbf{ v} _k , \end{gather}

이어서 위의 해를 구할 수 있다. 이 해는 $|\mathbf{A} . \mathbf{x} -  \mathbf{b}|^2$를 최소화한다.

cubic polynomial fitting

int svd(double *A, int m, int n, double* w, double *V); // from cmath libary.
void fit_func(double x, double val[], int n) {          // polynomial fitting sample;
    val[0] = 1;
    for(int i = 1; i < n; ++i)
        val[i] = x * val[i - 1];
}
#define EPSILON 1.E-8
int svd_fit(const double x[], const double y[], const int m, const int n,
            void (*fit_func)(double , double [], int ),
            double params[],
            double *error)
{
    double *A = new double [m * n];
    double *w = new double [n];
    double *V = new double [n * n];
    // evaluate design matrix;
    for (int i = 0; i < m; ++i)
        fit_func(x[i], &A[i * n + 0], n) ;

    svd(A, m, n, w, V);
    // now A becomes U matrix;
    // truncate small singular values;
    double wmax = 0;
    for (int i = 0; i < n; ++i)
        if (w[i] > wmax) wmax = w[i];
    double thresh = wmax * EPSILON;
    for (int i = 0; i < n; ++i)
        if (w[i] < thresh) w[i] = 0;
    
    // back substitution;
    double *tmp = new double [n];
    for (int j = 0; j < n; ++j) {
        double s = 0;
        if (w[j]) {
            for (int i = 0; i < m; ++i)
                s += A[i * n + j] * y[i];
            s /= w[j];
        }
        tmp[j] = s;
    }
    for (int j = 0; j < n; ++j) {
        double s = 0;
        for (int jj = 0; jj < n; ++jj)
            s += V[j * n + jj] * tmp[jj];
        params[j] = s;
    };

    //estimate error;
    *error = 0;
    for (int i = 0; i < m; ++i) {
        fit_func(x[i], &A[i * n + 0], n); //use A as a tmp buffer;
        double sum = 0;
        for (int j = 0; j < n; ++j) sum += params[j] * A[i * n + j] ;
        double err = (y[i] - sum);
        *error += err * err ;
    }
    delete[] A; delete[] w; delete[] V;
    delete[] tmp;
    return 1;
}

 

728x90

'Image Recognition > Fundamental' 카테고리의 다른 글

Image rotation by FFT  (0) 2022.02.18
FFT를 이용한 영상의 미분  (0) 2022.02.12
Color Histogram Equalization  (4) 2022.02.07
Least Squares Fitting of Ellipses  (0) 2022.01.27
Circle Fitting: Pratt  (0) 2022.01.20
Posted by helloktk
,

$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}^* = \text{argmin}(L), \quad\quad 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 할 수 있다. 

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

$$L = \sum_{k = 0}^{N-1} | \mathbf{B}(t_k) - \mathbf{P}_k|^2 = \sum_{k = 0}^{N-1} \left| \sum_{i=0}^{m-1} N_{i, p}(t_k) \mathbf{Q}_i -  \mathbf{P}_k \right|^2 = \sum_{k = 0}^{N-1} \left| \sum_{i=0}^{m-1} A_{ki} \mathbf{Q}_i - \mathbf{P}_k \right|^2, \\ \quad A_{ki} = N_{i, p}(t_k) $$

로 쓸 수 있으므로, $\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= \left| \hat{A} \cdot \hat {Q} - \hat {P} \right|^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}^* = ( \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)

int BSplineFit_LS(std::vector<CPoint>& data,
                  int degree,             // cubic(3); 
                  int nc,                 // num of control points;
                  std::vector<double>& X,
                  std::vector<double>& Y) // estimated 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;

    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)

    // S = A^t * A; real-symmetric matrix;
    std::vector<double> Sx(nc * nc);
    std::vector<double> Sy(nc * nc);
    for (int i = 0; i < nc; i++) {
        for (int j = 0; j < nc; j++) {
            double s = 0;
            for (int k = data.size(); k-->0;)
                s += A[k * nc + i] * A[k * nc + j];
            Sx[i * nc + j] = s;
        }
    }
    //copy;
    for (int i = nc*nc; i-->0;) Sy[i] = Sx[i];
    // X = A^t * P.x;  Y = A^t * P.y
    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] = sx; Y[i] = sy;
    };
    // solve real symmetric linear system; S * x = X, S * y = Y;
    // solvps(S, X) destories the inputs;
    // ccmath-2.2.1 version;
    int res1 = solvps(&Sx[0], &X[0], nc);
    int res2 = solvps(&Sy[0], &Y[0], nc);
    return res1 == 0 && res2 == 0;
};

**네이버 블로그 이전;

 
728x90

'Computational Geometry' 카테고리의 다른 글

Approximate Distance Between Ellipse and Point  (0) 2024.03.08
Distance from a Point to an Ellipse  (0) 2024.03.06
Closest Pair of Points  (0) 2021.04.27
DDA Algorithm  (0) 2021.04.25
B-Spline  (0) 2021.04.25
Posted by helloktk
,

점집합을 일반적인 2차 곡선으로 피팅하는 경우에 방정식은

$$ a x^2 + by^2 + cxy +d x + ey +f = 0$$

의 계수를 주어진 데이터를 이용하여서 구해야 한다. 실제 문제에서는 타원, 포물선 쌍곡 선등의 타입에 따라 몇 가지 제약 조건을 넣어 피팅을 한다. 원은 타원의 특별한 경우로 일반적으로 $a = b$, $c = 0$의 제약 조건이 필요하다. 그러나 보다 엄밀하게 제약을 하게 되면 $a = b = 1$의 추가 조건을 줄 수 있다. 이 경우는 점들이 모두 일직선에 있는 경우를 ($a = b = 0$) 취급할 수 없게 된다. 이 예외적인 경우를 제외하고는 최소자승법을 사용하면 계수를 매우 쉽게 구할 수 있기 때문에 많이 이용된다.

 

문제: 주어진 데이터를 fitting 하는 이차곡선(원)

$$x^2  + y^2 + A x + B  y + C = 0$$

의 계수 $A, B, C$를 최소자승법을 사용해서 구하라. 

 

주어진 점집합이 원 위의 점이면 우변이 0이 되어야 하나, 실제 데이터를 얻는 과정에서 여러 노이즈에 노출되므로 일반적으로 0이 되지 않는다. 최소자승법은 주어진 점들이 원에서 벗어나는 정도의 제곱 합이 최소가 되도록 하는 계수 $A, B, C$를 결정한다.  원과 점의 편차의 제곱합
$$ L=\sum_ i   \left |x_i^2 + y_i^2 + A x_i + B y_i + C \right|^2 , $$

의 극값을 찾기 위해서 $A, B,$ 그리고 $C$에 대해 미분을 하면

$$\frac{\partial L}{\partial A} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) x_i = 0, $$

$$\frac{\partial L}{\partial B} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) y_i = 0, $$

$$\frac{\partial L}{\partial C} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) = 0. $$

이 연립방정식을 풀면  $A, B, C$를 구할 수 있다. 우선 세 번째 식에서 

$$ CN = -S_{x^2} - S_{y^2} - AS_x - BS_y ,$$

을 얻고, 이를 첫번째와 두 번째 식에 각각 대입하면

$$A ( NS_{x^2} - S_x^2) + B ( N S_{xy} - S_x S_y ) =-N S_{x^3} - N S_{xy^2} + S_{x^2} S_x + S_{y^2} S_x, $$

$$A ( NS_{xy} - S_x S_y ) + B ( N S_{y^2} - S_y^2) = -N S_{y^3} - N S_{x^2 y}  +S_{x^2} S_y +S_{y^2} S_y, $$

을 얻을 수 있다. 다시 정리하면 두 개의 연립방정식

$$-a_1 A - a_2 B = 2c_1,$$

$$-b_1 A - b_2 B = 2c_2,$$

을 얻는다. $a_1, a_2 =  b_1, b_2, c_1, c_2$는 코드에서 정의되어 있다. 그리고

따라서, 추정된 원의 중심 $(c_x, c_y)$는 

$$ c_x = - \frac{A}{2} = \frac{c_1 b_2  - c_2 a_2}{ a_1 b_2 - a_2 b_1},$$

$$ c_y = - \frac{B}{2} = \frac{-c_1 b_1 + c_2 a_1}{a_1 b_2 -a_2 b_1},$$

로 주어지고, 반지름은 

$$r^2 =  c_x^2 +c_y^2 - C = c_x^2 + c_y^2 + \frac{1}{N}( S_{x^2}+S_{y^2}- 2c_x S_x - 2 c_y S_y)$$

로 주어진다.

Note: 질량중심 좌표계로 옮긴 후 moment를 계산하면 $S_x = S_y =0$이어서 식이 더 간결해지고, 수치적인 안정성도 좋아질 수 있다.

Ref: I. Kasa, A curve fitting procedure and its error analysis. IEEE Trans. Inst. Meas., 25:8-14, 1976

/* 구현 코드:*/
double circleFit_LS(std::vector<CPoint>& Q, double& cx, double& cy, double& radius) {
    if (Q.size() < 3) return -1;
    double sx  = 0.0,  sy = 0.0;
    double sx2 = 0.0, sy2 = 0.0, sxy  = 0.0;
    double sx3 = 0.0, sy3 = 0.0, sx2y = 0.0, sxy2 = 0.0;
    double mx = 0, my = 0;            /* center of mass;*/
    for (int k = Q.size(); k-->0;)
        mx += Q[k].x, my += Q[k].y;
    mx /= Q.size(); my /= Q.size();
    /* compute moments; */
    for (int k = Q.size(); k-->0;) { /* offset (mx, my)*/
        double x = Q[k].x - mx, xx = x * x;
        double y = Q[k].y - my, yy = y * y;
        sx   += x;       sy   += y;
        sx2  += xx;      sy2  += yy;      sxy  += x * y;
        sx3  += x * xx;  sy3  += y * yy;
        sx2y += xx * y;  sxy2 += yy * x;
    }
    /* compute a's,b's,c's; */
    const int N = Q.size();
    double a1 = 2.0 * (sx * sx - sx2 * N);
    double a2 = 2.0 * (sx * sy - sxy * N);
    double b1 = a2;
    double b2 = 2.0 * (sy * sy - sy2 * N);
    double c1 = (sx2 + sy2) * sx - (sx3 * N + sxy2) * N;
    double c2 = (sx2 + sy2) * sy - (sy3 * N + sx2y) * N;
    double det = a1 * b2 - a2 * b1;
    if (fabs(det) < 1.e-10) return -1;    /*collinear한 경우임;*/
    /* center; */
    cx = (c1 * b2 - c2 * b1) / det;
    cy = (a1 * c2 - a2 * c1) / det;
    /* radius squared */
    double radsq = cx * cx + cy * cy + (sx2 + sy2 - 2 * (sx * cx  + sy * cy)) / N;
    radius = sqrt(radsq);
    cx += mx; cy += my; /* recover offset; */
    return FitError(Q, cx, cy, radius);
}

 

 
 
 
 
 
 
 
 
728x90

'Image Recognition > Fundamental' 카테고리의 다른 글

PCA Line Fitting  (0) 2020.11.12
Histogram Equalization  (0) 2020.11.12
Integer Sqrt  (0) 2020.11.11
Parabolic Interpolation in Peak Finding  (3) 2020.11.10
Histogram Matching  (0) 2012.11.03
Posted by helloktk
,

Savitzky-Golay 필터는 일차원의 데이터에 대해 이동평균을 취하는 경우와 같은 방식으로 동작하는 필터이지만, 윈도의 모든 점에 동일한 가중치를 주는 이동평균과 다르게 윈도 픽셀 값을 보간하는 다항식을 최소자승법으로 찾아서 해당 지점의 값으로 할당하는 방식을 택한다(frequency domain에서 분석하면 Savitzky-Golay 필터의 특성, 예를 들면, 피크의 위치가 잘 유지되는 점과 같은 특성을 좀 더 다양하게 볼 수 있다). 이 필터를 쓰기 위해서는 다항식의 찾수와 윈도 크기를 정해야 한다. (다항식의 찾수가 정해지면 최소 윈도 크기는 정해진다).

동일한 방식으로 이차원에 대해서도 Savitzky-Golay를 적용할 수 있다. 이 경우 다항식은 $(x, y)$의 2 변수 함수로 2차원 평면에서 정의되는 곡면으로 나타낸다. 2차원 영상의 경우도 국소 필터를 사용할 수 있지만, 필터 윈도를 영상 전체로 잡아서 전 영역을 보간하는 곡면을 찾을 수도 있다. 배경 조명이 균일하지 않는 영상의 경우 이 곡면을 이용하면 조명에 의한 효과를 예측할 수 있고, 이를 보정한 영상을 이용하면 인식에 도움을 받을 수 있다. (문자 인식에서 문서를 스캔할 때 생기는 균일하지 않은 배경이나, 2차원 바코드 인식에서 배경의 추정 등 다양한 부분에서 사용할 수 있다. 좀 더 간단하게는 배경의 변화를 균일하게 기울어진 평면으로 근사를 하여 추정할 수 있다) 

3차 다항식으로 영상을 보간하는 경우: \begin{align} I(x, y)&= a_{00}\\ &+a_{10} x + a_{01} y \\ &+a_{20} x^2 + a_{11} xy + a_{02} y^2\\ &+a_{30} x^3+a_{21} x^2y+a_{12} xy^2+a_{03} y^3, \quad (x, y)\in \mbox {image} \end{align}

다항식은 $x= [a_{00}, a_{10},..., a_{03}]^T$ 의 10개의 필터 계수를 추정하면 얻어진다. 추가적으로 Savitzky-Golay을 이용하면 영상의 미분 값을 쉽게 구할 수 있다. 로컬 버전의 필터인 경우에 필터 적용 값은 윈도의 중심인 $(x, y) = (0, 0)$에서 다항식 값인 $a_{0}$이다. 이 지점에서 $x$-방향의 편미분 값은 $a_{10}$, $y$-방향의 편미분 값은 $a_{01}$로 주어진다.

필터의 계수 $x$는 최소자승법을 적용하면 얻을 수 있다. 위의 다항식에 $N(= width\times height)$개의 픽셀로 구성된 영상의 각 픽셀에서 좌표와 픽셀 값을 대입하면, $N$개의 식을 얻는다. 이를 행렬로 표현하면, 

$$\bf A\cdot x = b$$

$\bf A$는 $N\times10$ 의 행렬로 각 행은 픽셀의 좌표로 구해진다: 

$${\bf A} =\left[ \begin{array}{cccccccccc} 1&x_0&y_0&x_0^2&x_0y_0&y_0^2&x_0^3&x_0^2y_0&x_0y_0^2&y_0^3\\ 1&x_1&y_1&x_1^2& x_1y_1& y_1^2& x_1^3& x_1^2 y_1 & x_1 y_1^2 & y_1^3\\ 1& x_2& y_2 &x_2^2 & x_2 y_2& y_2^2 & x_2^3 & x_2^2 y_2 & x_2 y_2^2 & y_2^3 \\ &&&&\vdots \end{array} \right]$$

여기서, $i$-번째의 픽셀 위치가 $(x_i, y_i)$로 주어진 경우다. $\bf b$는 $N$-(열) 벡터로 각 픽셀 위치에서 픽셀 값을 나타내는 벡터다: 

$${\bf b}=\left[\begin{array}{c} I(x_0, y_0)\\I(x_1,y_1)\\I(x_2, y_2)\\ \vdots \end{array}\right]$$

최소자승법을 적용하면, 추정된 다항식의 계수 벡터 $\bf x$는 $|\bf A\cdot x - b|^2$을 최소로 하는 벡터로,

$$\bf x = (A^T \cdot A)^{-1} \cdot A^T \cdot b$$

로 주어짐을 알 수 있다. $\bf A^T \cdot A$는 $10\times 10$의 대칭 행렬로 역행렬은 쉽게 구할 수 있다.

이렇게 추정된 2차원 곡면은 영상에서 추정된 배경의 픽셀 값 분포를 의미한다. 문자인식의 예를 들면, 보통 경우에 흰 배경에 검은색 활자를 인식한다. 스캔된 영상에 검은색 활자 때문에 추정된 곡명은 일반적으로 주어진 픽셀이 만드는 곡면보다도 낮게 된다. 픽셀 값이 추정된 곡면보다 더 낮은 픽셀들은 보통 검은색 문자들을 의미하므로, 이 차이의 평균값을 구하면, 대략적으로 어떤 픽셀이 배경에 속하는지 (곡면과 차이가 평균보다 작고, 또한 픽셀 값이  곡면의 아래에 놓인 경우), 아니면 문자 영역인지(곡면과 차이가 평균보다 크고, 픽셀 값이 곡면의 아래에 놓인 경우)를 구별할 있게 된다.   

이제 이 정보들을 이용해서 추정을 다시 하는데 이번에는 1차 추정에서 글씨 영역으로 분류된 픽셀을 제외하고 배경을 추정하면 좀 더 정확한 배경을 기술하는 곡면을 얻을 수 있다.
로컬 필터로 사용할 때는 1차원에서와 마찬가지로 필터 계수를 lookup table로 만들어서 사용할 수 있으나, 전 영역을 대상으로 할 때는 행렬의 크기가 매우 커져서 연산량도 많아진다. 

영상:

 

1차 추정 배경 영상:

 

2차 추정 배경 영상:

 

728x90

'Image Recognition' 카테고리의 다른 글

Statistical Region Merging  (2) 2012.03.25
Local Histogram Equalization  (0) 2012.03.10
webcam용 QR code detector  (0) 2012.02.19
Least Squares Estimation of Perspective Transformation  (4) 2012.02.15
Perspective Transformation  (2) 2012.02.14
Posted by helloktk
,

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

또는, 

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

 

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

 

 

 

 

또는, 간단히 

$$ \bf A \cdot x = b$$

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

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

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

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

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

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

728x90

'Image Recognition' 카테고리의 다른 글

2차원 Savitzky-Golay Filters 응용  (0) 2012.02.28
webcam용 QR code detector  (0) 2012.02.19
Perspective Transformation  (2) 2012.02.14
Integral Image을 이용한 Adaptive Threshold  (0) 2012.02.04
Peak Finder  (1) 2012.02.02
Posted by helloktk
,