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
,

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 + 2 a_1 + 4 a_2 -d_4|^2$$

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

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

$$\text{필터 출력} =a_0 = (-3d_0 + 12 d_1 + 17 d_2 + 12 d_3 - 3 d_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.a = d}$ 형태로 쓸 수 있다. $\mathbf{A}$가 정방행렬이 아니므로 역행렬을 바로 구할 수 없지만, $|\mathbf{ A \cdot a - d} |^2$을 최소로 하는 최소제곱해는
$$\bf (A^T\cdot A)\cdot a = A^T \cdot d$$를 만족시켜야 하므로 

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

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

 


2차 다항식일 때, 같은 방식으로 다양한 윈도 크기에 따른 계수를 구할 수 있다.
*크기($n$)에 따른 필터값 결정계수 (중앙에 대해 좌우대칭이다);
\begin{align} n=5;\quad &W[~] = \{-3, 12, \mathbf{17}, 12, -3\};\\ n=7;\quad &W[~] = \{-2, 3,  6, \mathbf{7}, 6, 3, -2\};\\ n=9;\quad &W[~] = \{-21, 14, 39, 54, \mathbf{59}, 54, 39, 14, -21\}; \\n=11;\quad&W[~]=\{-36,9,44,69,84,\mathbf{89},84,69,44,9,-36\};\\ n=13;\quad &W[~]=\{13,-11,0,9,16,21,24,\mathbf{25},24,21,16,9,0,-11,13\}\end{align} 

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

int SavitzkyGolayFilter(std::vector<double>& data, double W[], int wsz, 
                        std::vector<double>& out) {
    int hwsz = wsz >> 1;
    wsz = (hwsz<<1) + 1;
    std::vector<double> data_pad(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++) { 
        data_pad[i]                  = data[hwsz-i];
        data_pad[i+data.size()+hwsz] = data[data.size()-2-i];
    }
    for (int i = data.size(); i-->0;) data_pad[i+hwsz] = data[i];
    //
    out.resize(data.size());
    double wsum = 0;
    for (int i = 0; i < wsz; i++) wsum += W[i]; 
    for (int i = data.size(); i-->0;) {
        double *pdata_pad = &data_pad[i];
        double fsum = 0;
        for (int k = 0; k < wsz; k++) fsum += pdata_pad[k] * W[k];
        out[i] = fsum / wsum;
    }
    return 1;
};

wsz=25, green=Savitzky-Golay, red=Moving average

 

728x90

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

Adaboost  (0) 2010.12.28
Blur Detection  (0) 2010.05.25
Retinex 알고리즘  (11) 2010.02.03
Gaussian Mixture Model & KMeans  (4) 2010.01.30
Image Morphing  (0) 2010.01.24
Posted by helloktk
,