'Least Square Method'에 해당되는 글 8건

  1. 2012.02.28 2차원 Savitzky-Golay Filters 응용
  2. 2012.02.15 Least Square Estimation of Perspective Transformation (4)
  3. 2010.03.24 Savitzky-Golay smoothing filter (2)
  4. 2010.01.20 Affine Transformation
  5. 2010.01.11 Isometric Transformation
  6. 2009.12.14 Similarity Transformation (1)
  7. 2008.07.21 RANSAC : Circle Fit
  8. 2008.07.08 Robust Line Fitting
Savitzky-Golay 필터는 일차원의 데이터에 대해서 일종의 이동평균을 취하는 경우와 동일하게 동작하는 필터이지만, 추정하는 지점의 주변의 모든 점에 동일한 가중치를 주는 방식(이동평균)을 택하지 않고, 그들을 보간하는 다항식을 최소자승법으로 찾아서 해당 지점의 값을 추정하는 방식을 택한다(frequency domain에서 분석을 하면 Savitzky-Golay 필터의 특성, 예를 들면, 피크의 위치나 등이 잘 유지되는 점등, 좀 더 다양하게 볼 수 있다). 이 필터를 쓰기 위해서는 몇 차의 다항식과 얼마의 윈도우 크기를 사용해야 하는지 설정을 해야 한다. (다항식의 찻수가 정해지면 최소의 윈도우 크기가 정해진다).
동일한 방식으로 이차원에 대해서도 Savitzky-Golay를 적용할 수 있다. 이 경우에는 다항식이 (x,y)의 2변수 함수 (2차원 평면에서의 정의되는 곡면)로 주어질 것이다. 이차원의 경우도 국소적인 필터로 사용하여서 영상의 smoothing 필터로 사용할 수 있지만, 필터의 윈도우를 영상전체로 잡아서, 즉 영상을 구성하는 픽셀값을 전 영역에서 보간하는 곡면을 찾을 수도 있다. 이렇게 찾은 곡면은 만들어진 영상의 배경조명이 균일하지 않는 경우에 이 추정된 곡면을 이용하면, 조명에 의한 효과를 예측할 수 있고, 배경조명이 보정된 영상을 만들어서 영상의 인식에 도움을 받을 수 있다. (문서인식에서 문서를 스캔할 때 생기는 균일하지 않은 배경이나, 2차원 코드 인식에서 배경의 추정등 다양한 부분에서 사용할 수 있다. 물론 간단한 경우에는 배경의 변화를 균일하게 기울어진 평면으로 근사를 하여서 추정할 수 있다)  

간단한 경우가 3차 다항식으로 영상을 보간하는 경우:

I(x, y) = a00
         + a
10*x + a01*y
         + a
20*x2 + a11*x*y + a02*y2
         + a
30*x3 + a21*x2*y + a12*x*y2 + a03*y3 
                                                      (x, y) ∈ image 


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

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

A.x = b

A는 N x 10 의 행렬로 각 행은 픽셀의 좌표로 구해진다:

A = [1, x0,  y0,  x02,  x0*y0,  y02,  x03,  x02*y0,  x0*y02,  y03]
     [1, x1,  y1,  x12,  x1*y1,  y12,  x13,  x12*y1,  x1*y12,  y13] 
     [1, x2,  y2,  x22,  x2*y2,  y22,  x23,  x22*y2,  x2*y22,  y23] 
     [1, .....................................................................] 
     [1, .....................................................................] 
                       ......................................
     [1, .....................................................................] 

여기서, 영상을 읽을 때, i-번째의 픽셀 위치가 (xi, yi) 로 주어진 경우다.

b 는 N-(열)벡터로 각각의 픽셀 위치에서 픽셀 값을 나타내는 벡터이다:

b = [ I(x0, y0) ]
     [ I(x1, y1) ] 
     [ I(x2, y2) ] 
     [ ............] 
     [ ............]
      ..............
     [.............]  


최소자승법을 적용하면, 추정된 다항식의 계수벡터 x는 

x = (A
T .A)-1.AT.b

임을 알 수 있다.

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

로컬버전 필터로 사용할 때는 1차원에서와 마찬가지로 필터계수를 lookup table로 만들어서 사용할 수 있으나, 전영역을 대상으로 할 때는 행렬의 크기가 매우 커져서 연산도 많아진다. 

""

""

""


저작자 표시 비영리 변경 금지
신고
Posted by helloktk

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

또는, 

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


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

또는 간단히 

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


에 대해서 미분을 하여서,

를 만족시키는 극값 x*을 구하면 된다. 는 8x8의 대칭행렬이어서 대각화가 가능하므로 역행렬을 구할 수 있다 (주어진 점들 중 한 직선 위에 놓이지 않는 점이 4개 이상이 있어야 한다). 따라서, 최소제곱해는 다음과 같이 주어진다:

.

저작자 표시 비영리 변경 금지
신고

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

2차원 Savitzky-Golay Filters 응용  (0) 2012.02.28
webcam용 QR code detector  (0) 2012.02.19
Least Square Estimation of Perspective Transformation  (4) 2012.02.15
Perspective Transformation  (2) 2012.02.14
Integral Image을 이용한 Adaptive Threshold  (0) 2012.02.04
Peak Finder  (1) 2012.02.02
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, d0), (-1, d1), (0, d2), (1, d3), (2, d4)}

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

p[x] = a0 + a1*x + a2*x2

이다. 다항식의 계수들은 다항식의 값과 실제 데이터의 값과의 차이를 최소화시키는 것들로 잡아야 한다.. 즉, 최소자승의 원리를 적용하여서 구하면 된다. 계산된 다항식의 값과 실제 데이터 값 사이의 차의 제곱을 구하면:

L = |a0-2*a1+4*a2-d0|2+|a0 -a1+ a2-d1|2+|a0-d2|2+|a0+a1+a2-d3|2
                          
+|a0+2*a1+4*a2-d4|2

이 식의 a0, a1, a2에 대한 극값을 가질 조건은

5*a0 + 10*a2 = d0 + d1 + d2 + d3 + d4

10*a1 = -2d0– d1 + d3+ 2d4

10*a0 + 34*a2 = 4d0 + d1+ d3 + 4d4        

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

필터출력 = a0 = (-3*d0 + 12*d1 + 17*d2 + 12*d- 3*d4) / 35;

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

""

좌변의 5행3열 행렬을 A, a = [a0, a1, a2]T, d = [d0, d1, d2, d3, d4]T로 놓으면, 이 행렬방정식은 A.a = d 형태로 쓸 수 있고, 양변에 AT를 곱해서 왼쪽을 대칭행렬로 바뀌면

(AT.A).a = AT.d

따라서 해는,

a = (AT.A)-1.(AT.d)

이 식은 임의의 k-차 다항식을 이용한 경우에도 사용할 수 있다. 이 경우에 행렬 (ATA) 는 (k+1)x(k+1)의 대칭행렬이 된다. 행렬 A는 다항식의 찻수와 피팅에 사용이 될 데이터의 구간의 크기가 주어지면 정해지므로,  윗 식에서 (AT.A)-1.AT의 첫행 (a0을 d로 표현하는 식의 계수들)을 구해서 프로그램상에서는 lookup table를 만들어서 사용할 수 있다. 아래표는 mathematica 를 이용해서 윈도우 크기가 7 (7개 점)인 경우 2차 다항식을 사용할 때 계수를 구하는 과정이다.

 

""

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

n=5;  W[] = {-3, 12, 17, 12, -3};
n=7;  W[] = {-2, 3,  6, 7, 6, 3. -2};
n=9;  W[] = {-21, 14. 39, 54, 59, 54, 39, 14, -21};
n=11; W[] = {-36, 9, 44, 69, 84, 89, 84, 69, 44, 9, -33};

필터출력 =   ∑ W[i]d[i] / ∑ W[i]

저작자 표시 비영리 변경 금지
신고

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

Adaboost  (0) 2010.12.28
Blur Detection  (0) 2010.05.25
Savitzky-Golay smoothing filter  (2) 2010.03.24
Watershed Algorithm Code  (0) 2010.03.19
Retinex 알고리즘  (10) 2010.02.03
Gaussian Mixture Model & KMeans  (4) 2010.01.30
Posted by helloktk

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

N개의 꼭지점을 갖는 두개의 형상 S={(x1,y1),...(xN,yN)}, S'={(x'1,y'1),...,(x'N,y'N)}이 affine 변환에 의해서 연결이 되는 경우에 각각의 꼭지점들 간에는

          x'i = a * xi + b * yi + tx ;
          y'i = c * xi + d * yi + ty ;         i=1,...,N; 

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

꼭지점의 위치는 실제로 다양한 영상처리 방법에 의해서 주어지므로 필연적으로 노이즈를 포함하게 되어서 일종의 램덤변수로 생각해야 한다. 따라서 주어진 랜덤변수에서 최적으로 매개변수를 추출하기 위해서 최소자승방법을 이용한다. 즉, affine변환된 좌표값과 실제값의 차이를 최소화하는 매개변수를 찾는 것이다. 

         L = ∑ ( |x'i - a * xi - b * yi - tx |2 + |y'i - c * xi - d * yi - ty|2 ),

affine변환을 규정하는 매개변수들을 구하기 위해서는 L을 각각의 매개변수들에 대해서 미분을 하여 그 극값을 구하면 된다.

        ∂L/∂a = -2 * ∑ (x'i - a * xi - b * yi - tx) * xi ;
        ∂L/∂b = -2 * ∑ (x'i - a * xi - b * yi - tx) * yi ;
        ∂L/∂c = -2 * ∑ (y'i - c * xi - d * yi - ty) * xi ;
        ∂L/∂d = -2 * ∑ (y'i - c * xi - d * yi - ty) * yi ; 
        ∂L/∂tx = -2 * ∑ (x'i - a * xi - b * yi - tx) ;
        ∂L/∂ty = -2 * ∑ (y'i - c * xi - d * yi - ty); 

각각의 식들을 0으로 놓아서 얻어지는 연립방정식을 행렬의 형태로 다시 정리하여서 쓰면 그 해는 쉽게 구할 수 형태로 된다.

             [ Sxx   Sxy   Sx ][ a  c ]     [  Sxx'   Sxy' ]
             [ Sxy   Syy   Sy ][ b  d ]  =  [  Syx'   Syy' ]
             [ Sx     Sy    N   ][tx, ty ]     [  Sx'     Sy'  ]
 

여기서,
Sxx = ∑ x*x, Syy = ∑ y*y, Sxy = ∑ x*y,
Sx = ∑ x, Sy = ∑ y,
Sxx' = ∑ x*x', Sxy' = ∑
x*y', Syx' = ∑ y*x',
Sx' = ∑ x', Sy' = ∑ y'

이다.

// dst = (A,T)src;

//  [u]  = [ A0 A1 ][x] + A4
//  [v]  = [ A2 A3 ][y] + A5
//
BOOL GetAffineParameter(POINT *srcPts, POINT *dstPts, int n, double AT[6]) {
    double Sx, Sy, Sxx, Sxy, Syy;
    double Su, Sv, Sxu, Sxv, Syu, Syv ;
    double A[9], invA[9] ;
    double det ;
    Sx = Sy = Sxx = Sxy = Syy = 0;
    Su = Sv = Sxu = Sxv = Syu = Syv = 0;
    for (int i = 0; i < n; i++) {
        double x = srcPts[i].x ;
        double y = srcPts[i].y ;
        double u = dstPts[i].x ;
        double 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] = n  ;
    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을 만들 수 있다. 물론 녹색점이 신뢰할 수 있는 것인가에 대한 것은 다른 확률적인 추정을 할 수 있는 정보들이 더 필요하다.

저작자 표시 비영리 변경 금지
신고

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

Image Morphing  (0) 2010.01.24
Fant's Algorithm  (0) 2010.01.22
Affine Transformation  (0) 2010.01.20
Color Counting  (0) 2010.01.18
Isometric Transformation  (0) 2010.01.11
Active Shape Model (3)  (0) 2009.12.30
Posted by helloktk
기준 좌표계에 대해서 원점이 이동하고 회전한 경우에 이 새로운 좌표계에서 점들의 좌표는 바뀐다. 원래의 좌표와 바뀐 좌표값 사이의 관계를 주는 변환이 Isometric transformation (isometry)이다. 평면상에서 이변환은 평행이동 2개 회전각도 1개에 의해서 결정이 된다. 다수의 점과 그 대응점들이 주어진 경우에 두 둘간의 isometry 파라미터는 일반적으로 최소자승법에 의해서 결정될 수 있다.
회전각도가 θ이고, 평행이동이 (tx, ty)인 isometry에 의해서 두 점 (x,y)가 (u,v)로 연결이 되는 경우에, 아래의 식으로 표현이 된다:

u = cosθ * x - sinθ * y + tx ;
v = sinθ * y + cosθ * x + ty ;

최소자승법을 사용하기 위해서는 회전각도 θ보다는 a = cosθ, b = sinθ로 정의된 새로운 파라미터로 식을 표현하는 것이 더 편리하다. 그러나 이 경우에 파라미터 a, b는 서로 독립적이 아니고  a^2 + b^2 = 1의 제한조건을 만족시켜야 한다.  
평행이동 파라미터는 질량중심의 isometry 관계로 해결이 되므로, 이  전체계산을 각각의 질량중심을 기준으로 하는 좌표로 옮겨서 하면 더 이상 평행이동을 고려할 필요가 없고 회전만 생각하면 된다. 

최소자승법의 원리에 따라서 아래의 라그랑지안을 최소화시키는 파라미터 a, b, λ를 찾으면 된다

L = ∑ ( a * xi - b * yi - ui)^2 + ( b * xi + a * yi - vi)^2 + λ*(a^2 + b^2 - 1) ;

여기서 λ는 제한조건 a^2 + b^2 = 1를 넣기 위한 Lagrange multiplier이다. 라그랑지안 L를 각각 a, b, λ에 대해서 미분하여서 다음의 조건을 얻는다:

∑ ( a * xi - b * yi - ui)*xi + ( b * xi + a * yi - vi)*yi + λ*a = 0
∑ ( a * xi - b * yi - ui)*(-yi) + ( b * xi + a * yi - vi)*xi + λ*b = 0
a^2 + b^2 = 1;

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

a = (∑ xi*ui + yi*vi) / (∑ xi^2 + yi^2 + λ) ;
b = (∑ xi*vi - yi*ui) / (∑ xi^2 + yi^2 + λ) ;

또한, 라그랑지 멀티플라이어 λ는

A  = (∑ xi*ui + yi*vi);   B =  (∑ xi*vi - yi*ui);

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

∑ xi^2 + yi^2 + λ = sqrt(a^2 + b^2) 

임을 쓰면 된다. 따라서 우리가 구하는 회전각도는

a = A / sqrt(A*A + B*B);
b = B / sqrt(A*A + B*B);

로 주어진다.

질량중심을 빼기전의 좌표 (x,y)의 질량중심과 (u,v)의 질량중심은 서로 isometry에 의해서 연결이 되므로, 이것에서  평행이동 파라미터(tx, ty) 는 결정이 된다;
(xc, yc) --->(uc, vc) ;
uc = a * xc - b * yc + tx ;
vc = b * yc + a * yc + ty ;

참고:

** 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(CfPt A[], CfPt U[], int n, double T[6]) {
    double cx = 0, cy = 0;
    double ux = 0, uy = 0;
    for(int i = 0; i < n ; i++) {
        cx += A[i].x ;  cy += A[i].y ;
        ux += U[i].x ;  uy += U[i].y ;
    };
    //center of mass ;
    cx /= n; cy /= n;
    ux /= n; uy /= n;

    //centering 된 좌표계에서 계산;
    double dot = 0 , cross = 0;
    for(i = 0; i < n; i++) {
        double x = A[i].x - cx ;
        double y = A[i].y - cy ;
        double u = U[i].x - ux ;
        double 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;
} ;
저작자 표시 비영리 변경 금지
신고

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

Affine Transformation  (0) 2010.01.20
Color Counting  (0) 2010.01.18
Isometric Transformation  (0) 2010.01.11
Active Shape Model (3)  (0) 2009.12.30
Eigenface (2)  (0) 2009.12.28
Active Shape Model :: Face  (1) 2009.12.27
Posted by helloktk

2차원 이미지의 기하학적인 변형중에서 일정하게 이동을 시키고, 일정각도 회전 및 전체적인 크기의 변화를 주는 변환이 similarity transformation이다. 이 변화은 두 직선이 이루는 각도를 보존하고 길이비를 변경시키지 않는다. 이 변환보다도 더 일반적인 2차원의 기하학적인 변환은 affine transformation이다. 이것은 한쪽 방향으로의 밀림(sheer)도 허용한다. 그러다 두 직선의 각도는 변환후에도 변하지 않는다.
similarity transformation은 전체적인 크기를 바꾸는 스케일 파리미터(s) 1개와 회전각도(θ) 1개 그리고 x, y축으로의 이동을 나타내는 파라미터 2 (tx, ty)개를 합해서 총 4개의 파라미터가 필요하다. 이 파라미터에 의해서 원본 이미지의 (x,y)점이 변환된 이미지의 (u,v)에 대응한다고 하면, 이들간의 관계는 아래의 식으로 주어진다.

    u =  s*cos(θ)*x - s*sin(θ)*y + tx ;
    v =  s*sin(θ)*y + s*cos(θ)*y + ty;

 이 따라서 원본영상의 2 점에 대응하는 정보만 주어지면 파라미터(s, θ, tx, ty)를 유일하게 결정할 수 있다. 

    (x1, y1) -----> (u1, v1)
    (x2, y2) -----> (u2, v2)

그러나 많은 경우에는 기준점을 잡는데 에러등을 고려하여서 일반적으로 원본 영상의 N(>=2)개의 점에 대응하는 정보를 주게되는데, 이 경우에 변환관계식은 overdetermined되어서 해를 구할 수 없는 경우도 있다. 이 경우에는 최소자승법을 써서 변환점과 변환식에 의해서 의해서 주어지는 값의 차이를 최소화시키는 파라미터를 구해서 쓰면 된다.

L =  ∑ |ui - (s*cos(θ)*xi - s*sin(θ)*yi + tx)|^2 + |vi - (s*sin(θ)*xi + s*cos(θ)*yi + ty)|^2,           (i=1,...,N)

                        (s, θ, tx, ty) = arg min (L)

이 식을 최소화사키는 파라미터는  (a= s*cos(θ), b=s*sin(θ)로 놓으면)  a, b, tx, ty에 대해서 미분하여서 0인 조건을 만족시키면 된다.

a에 대한 미분 :   ∑ (ui - (a*xi - b*yi + tx))*(-xi) + (vi - (b*xi + a*yi + ty))*(-yi) = 0,   
b에 대한 미분 :   ∑ (ui - (a*xi - b*yi + tx))*(yi) + (vi - (b*xi + a*yi + ty))*(-xi) = 0,   
tx에 대한 미분 :  ∑ (ui - (a*xi - b*yi + tx)) = 0.
ty에 대한 미분 :  ∑ (vi - (b*xi + a*yi + ty)) = 0

따라서, Su = ∑ ui, Sv = ∑ vi, Sux = ∑ ui*xi, Suy = ∑ ui*yi, Svx = ∑ vi*xi, Svy = ∑ vi*yi , Sx = ∑ xi, Sy=∑ yi, Sxx = ∑ xi*xi, Sxy = ∑ xi*yi, Syy=∑ yi*yi 라고 하면,

-Sux  + a * Sxx - b * Sxy + tx * Sx - Svy + b*Sxy + a*Syy + ty *Sy = 0; 
  Suy - a * Sxy + b * Syy - tx * Sy -Svx + b*Sxx + a*Sxy + ty*Sx = 0;
  Su - a* Sx + b*Sy - tx * N = 0;
  Sv - b* Sx - a*Sy - ty * N = 0;

의 4개의 식을 얻으므로 4원(a,b,tx,ty) 1차 연립방정식을 풀면 된다. 이식의 답은 쉽게 구할 수 있고, 아래의 코드는 이것을 구현한 것이다. 물론, 점이 2개인 경우에는 파라미터는 유일하게 정해지고, 이보다도 더 간단한 식으로 주어진다.

//dstPt = (ST)(srcPt)
BOOL SimilarTransParams(POINT *srcPts, POINT *dstPts, int n, double ST[4]) {
    double Sx, Sy, Sxx, Syy;
    double Su, Sv, Sxu, Sxv, Syu, Syv ;
    Sx = Sy = Sxx = Syy = 0;
    Su = Sv = Sxu = Sxv = Syu = Syv = 0;
    for (int i = 0; i < n; i++) {
        double x = srcPts[i].x ;
        double y = srcPts[i].y ;
        double u = dstPts[i].x ;
        double v = dstPts[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 C1 = Sxu + Syv ;
    double C2 = Sxv - Syu ;
    double A[16] , invA[16] ;
    A[0]  = Sx; A[1]  = -Sy;  A[2]  =   n; A[3] = 0.;
    A[4]  = Sy; A[5]  =  Sx;  A[6]  =  0.; A[7] = n ;
    A[8]  = Z ; A[9]  =  0.;  A[10] =  Sx; A[11] = Sy;
    A[12] = 0.; A[13] =  Z;   A[14] = -Sy; A[15] = Sx;
    InvertMatrix4x4_d(A, invA) ;
    double R[4] ;
    R[0] = Su ; R[1] = Sv; R[2] = C1; R[3] = C2 ;
    // ax = scale * cos(angle) ;
    double ax = invA[0]*R[0]  + invA[1]*R[1]  + invA[2]*R[2]  + invA[3]*R[3];
    // ay = scale * sin(angle) ;
    double ay = invA[4]*R[0]  + invA[5]*R[1]  + invA[6]*R[2]  + invA[7]*R[3];
    // x-translation ;
    double tx = invA[8]*R[0]  + invA[9]*R[1]  + invA[10]*R[2] + invA[11]*R[3];
    // y-translation ;
    double ty = invA[12]*R[0] + invA[13]*R[1] + invA[14]*R[2] + invA[15]*R[3];
    ST[0] = ax ;
    ST[1] = ay ;
    ST[2] = tx ;
    ST[3] = ty ;

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

더보기

2개의 대응점만 주어진 경우 (x1,y1), (x2, y2)-->(u1,v1), (u2,v2);
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;
    double y21=y2-y1;
    double u21=u2-u1;
    double 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의 역할만 바꾸면 쉽게 구할 수 있다.

신고

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

Active Shape Model :: Face  (1) 2009.12.27
Active Shape Model (ASM)  (2) 2009.12.25
Similarity Transformation  (1) 2009.12.14
Eigenface  (0) 2009.12.12
Retinex 알고리즘 관련 자료.  (1) 2009.04.29
Spline Based Snake  (0) 2008.08.15
Posted by helloktk

RANSAC 알고리즘을 이용하여서 원을 추정한다. 원을 만들기 위해서는 최소한 3개의 일직선상에 있지 않은 점을 잡아야 하면, 이때 원은 세점의 외접원으로 주어진다. 이 원을 기본으로 해서 다시 원의 중심과 반지름을 추정한다. 추정된 원들 중에서 inlier 들의 벗어남의 편차가 가장 작은 것을 결과로 이용한다.
// 참고: http://en.wikipedia.org/wiki/RANSAC 
//

// [0,max_num) 사이에서 3개의 숫자를 무작위로 선택함;
void GetRandomTriplet(int max_num, int triplet[3]) {

더보기

// 세 점의 외접원을 구함;
int CircumCircle(double x1, double x2, double x3,
                 double y1, double y2, double y3,
                 double *cx, double *cy, double *rad) {

더보기

// 3x3 선형방정식을 푼다: A * x = b;
bool SolveLinearEQ3x3(double A[9], double bb[3], double x[3]) {

더보기

// x^2 + y^2 + A*x + B*y + C = 0;
// Least square fit for A, B, C;
int CircleFit_LS(int N, double xp[], double yp[],
                   double *cx, double *cy, double *rad) {  

더보기

// 데이터 중에서 주어진 원에서 거리 임계값 이내의 데이터만 골라냄;
int findInlier(double xp[], double yp[], int N,
                double cx, double cy, double rad,
                double dist_th,
                double consensusx[], double consensusy[], double *var) {

더보기

int RansacCircleFit(int N, double xp[], double yp[],
                    int sample_th,      //# of inliers; >= 50% of data(N), 66%;
                    double dist_th,     // 25% of radius;   |dist-rad|/rad< dist_th
                    int max_iter,
                    double *centerx,
                    double *centery,
                    double *radius)
{    
    double pr = double(sample_th) / double(N);
    double trials99 = log(1. - 0.99)/log(1. - pow(pr, 3));
    double tx[3], ty[3];
    int tri[3];
    int found = 0;  
    int iter = min(max_iter, trials99);
    //inlier buffer;
    std::vector<double> consensusx(N), consensusy(N) ;
    double min_dev = 1.e+10, var, sdev;
    if (sample_th < 3) sample_th = 3;
    while (iter) {
        GetRandomTriplet(N, tri);
        for (int i = 0; i < 3; i++) {
            tx[i] = xp[tri[i]];
            ty[i] = yp[tri[i]];
        }
        double cx, cy, rad;
        if (!CircumCircle(tx[0], ty[0], tx[1], ty[1], tx[2], ty[2], &cx, &cy, &rad)) {
            // if tree points are degenerate or on the sample line, discard them!
            continue ;
        }
        int ninlier = findInlier(xp, yp, N, cx, cy, rad, dist_th, &consensusx[0], &consensusy[0], &var);
        if (ninlier >= sample_th) {          
            // estimate model params using maybe-inliers;
            if (!CircleFit_LS(ninlier, &consensusx[0], &consensusy[0], &cx, &cy, &rad))
                continue; //least-square fitting fails;
            // collect real-inliers;
            ninlier = findInlier(xp, yp, N, cx, cy, rad, dist_th / 2, &consensusx[0], &consensusy[0], &var);
            if (ninlier < sample_th)
                continue;
            sdev = sqrt(var / ninlier);
            // estimate model params using inliers again;
            if (!CircleFit_LS(ninlier, &consensusx[0], &consensusy[0], &cx, &cy, &rad))
                continue;            
            TRACE("cx = %f, cy = %f, rad=%f\n", cx, cy, rad);
            if (sdev < min_dev) {
                *centerx = cx;
                *centery = cy;
                *radius  = rad;
                min_dev = sdev;
                found = 1;
                // we need more elaborate termination conditions;
#if _DEBUG
#endif
            }
        }
        --iter;
    }
    return found;
};

  • sample_th = 2 * N / 3;
  • dist_deviate_th = 0.25;
  • 파란색: 최소자승법을 이용한 피팅 (outlier의 영향을 많이 받음);
  • 붉은색: RANSAC 피팅 결과 (2017.01.04 수정)


see also http://blog.naver.com/helloktk/80029898273


신고

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

Chamfer Match  (0) 2008.08.01
Retinex Algorithm  (1) 2008.07.26
RANSAC : Circle Fit  (0) 2008.07.21
KMeans Algorithm  (0) 2008.07.19
Robust Line Fitting  (0) 2008.07.08
Bayesian Spam Filtering  (0) 2008.07.03
Posted by helloktk

이미지에서 관찰된 점들이 {(xi, yi) | i = 1, 2,..., N}이 있다. 이 점들을 직선 y = a + b*x 로 피팅을 하고 싶을 때, 보통 최소 자승법을 이용하는 데,  원리는 직선의 방정식이 예측한 y값과 실제 관찰한 y값의 차이의 제곱(=square deviation)을 최소화시키는 직선의 파라미터 a, b를 찾는 것과 동일하다:


확률적으로는, 데이터를 얻는 과정인 측정에서 측정값 yi 는 랜덤에러를 가지게 되고, 이 에러는 참값 y(x) 근방에서 정규분포의 모양을 가지게 된다고 가정을 한다. 만약에 모든 측정의 에러가 동일한 표준편차 σ를 갖게 된다면, N개의 관측데이터가 나타날 확률(likelihood)는 (개별측정은 IID조건을 만족한다고 전제)
     

의 형태가 된다. 따라서 최소자승방법은 이 likelihood를 최대화시키는 파라미터를 찾는 방법이 된다. 이 최소 자승법은 피팅 파라미터를 주어진 데이터를 이용해서 표현할 수 있는 장점이 있지만, outliers에 대해서 매우 취약하다 (아래의 결과 그림을 참조). 이것은 적은 수의 outlier도 χ2 에 큰 기여를 할 수 있기 때문이다. 따라서 피팅을 좀더 robust하게 만들기 위해서는 이들 outlier들이 likelihood에 기여하는 정도를 줄여야 한다. 그러기 위해서는 likelihood의 지수인자가 큰 에러에서 덜 민감하게 반응하는 꼴로 바뀌어야 한다. 이런 형태 중에서 가장 간단한 것 중의 하나가 square deviation 대신에 absolute deviation을 이용하는 것이다:


그러나 이 식을 이용하는 경우에는 최소자승법과는 달리, 파라미터 a, b를 주어진 데이터 {(xi, yi)}로 표현할 수 없고, 반복적인 방법을 이용해서 찾는 수 밖에 없다. 


주어진 수열 {c
i}와 a 와의 차이의 절대값의 합
i |ci - a| 는 a가 이 수열의 median 값인 경우에 최소된다는 사실을 이용하면 (증명: 극값을 구하기 위해서 a에 대해서 미분하면, 0 = (i 1) - (i 1) : 합은 a가 ci 보다 큰 경우와 작은 경우로 분리=> 따라서 0이 되기 위해서는 작은 경우와 큰 경우의 수가 같아야 한다 => 즉, a = median{ci}), 고정된 b값에 대해서 위 absolute deviation을 최소화 시키는 a값은

 


임을 알 수 있다. 그리고 absolute deviation 식을 b에 대해서 미분을 하여서


을 얻는데 위의 a 값을 대입해서 위의 식을 만족시키는 b값을 bracketing and bisection 방법을 이용해서 얻을 수 있다(불연속함수여서, 일반적으로 근을 구하는 방법을 사용하는 것은 위험하다). 아래의 코드는 이를 구현한 것이다.


참고 : Numerical Recipes. 

// qsort-comparator; 

static int comp(const void *a, const void * b) { 

    double d = *((double *)a) - *((double *)b);

    return d > 0 ? 1 : d < 0 ? -1 : 0; 

// 기울기(bb)가 주어진 경우에 y-절편(median = aa)값을 추정하고, 이 때 AD값을 얻는다. 

double RhoFunc(double *x, double* y, int n, 

               double bb, double *aa, double *abdev) { 

    double *h = new double [n];

    for (int i = 0; i < n; i++)  h[i] = y[i] - bb * x[i];

    qsort(h, n, sizeof(double), comp); 

    // median; 

    *aa = (n & 1) ? h[n / 2] : (h[n / 2] + h[n / 2 - 1]) / 2;

    

    double sum = 0;  

    *abdev = 0; 

    for (i = 0; i < n; i++) { 

        double d = y[i] - (bb * x[i] + (*aa)); 

        *abdev += fabs(d); 

        // sign-함수의 원점에서 모호함을 없애기 위해서 증폭을 시킴; 

        if (y[i] != 0.) d /= fabs(y[i]);

        if (fabs(d) > DBL_EPSILON) // sign 함수의 모호함에서 벗어나면,

            sum += (d >= 0 ? x[i] : -x[i]); 

    } 

    delete [] h;

    return sum; // return sum{xi * sign(yi - (b * xi + a))} 

}; 

// 최소자승법을 이용한 라인추정: 

// return (sigma[dy] / sigma[x]);

double FitLine_LS(double *x, double *y, int n, double *a, double *b) { 

    double sx = 0, sy = 0, sxx = 0, sxy = 0; 

    for (int i = 0; i < n; i++) { 

        sx += x[i]; 

        sy += y[i]; 

        sxx += x[i] * x[i] ; 

        sxy += x[i] * y[i] ; 

    }; 

    double det = n * sxx - sx * sx; 

    if (det == 0.) return -1;                   // vertical line;

    *a = (sxx * sy - sx * sxy) / det; 

    *b = (n * sxy - sx * sy) / det;  

    double chi2 = 0;  

    for (i = 0; i < n; i++) { 

        double t = y[i] - (*a + *b * x[i]); 

        chi2 += t * t; 

    }  

    det /= n; // det -> var(x) * n; 

    // chi2 = var(dy) * n;

    // (dy vs x의 편차비)

    return  sqrt(chi2 / det);  ; 

}

// y = a + b * x ; 

// Least Absolute Deviation: 

void FitLine_MAD (double *x, double *y, int n, 

                  double *a, double *b, double *mad) 

    // least square esimates for (aa, bb);      

    double aa, bb; 

    double sigb = FitLine_LS(x, y, n, &aa, &bb); 

    double b1 = bb; 

    double abdev; 

    double f1 = RhoFunc(x, y, n, b1, &aa, &abdev); 

    /* bracket 3-sigma away in the downhill direction; */ 

    double b2 = fabs(3 * sigb); 

    b2 = bb + (f1 < 0 ? -b2 : b2); 

    double f2 = RhoFunc(x, y, n, b2, &aa, &abdev); 

    

    // if conditional added to take care of the case of a

    // line input into this function. It is needed to avoid an

    // infinite loop when (b1 == b2) (within floating point error)

    if (fabs(b2 - b1) > (sigb + 0.005)) { 

        // bracketing;

        while ((f1 * f2) > 0) { 

            bb = 2 * b2 - b1; 

            b1 = b2;

            b2 = bb;

            f1 = f2;  

            f2 = RhoFunc(x, y, n, b2, &aa, &abdev) ; 

        } 

    } 

    // refine until the error is a negligible number of std; 

    sigb *= 0.01; 

    while (fabs(b2 - b1)> sigb) { 

        // bisection; 

        bb = (b1 + b2) / 2.; 

        if ((bb == b1) || (bb == b2)) break ; 

        double f = RhoFunc(x, y, n, bb, &aa, &abdev) ; 

        if ((f * f1) >= 0) { 

            f1 = f; 

            b1 = bb; 

        } 

        else { 

            f2 = f; 

            b2 = bb; 

        } 

    } 

    *a = aa; 

    *b = bb; 

    *mad = abdev / n; 

}

// 붉은 선--> 최소자승법에 의한 피팅.: outlier에 매우 취약한 구조.
// 파란 선--> least absolute deviation을 이용한 피팅: outlier에 매우 robust하다.

"사용자

신고

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

RANSAC : Circle Fit  (0) 2008.07.21
KMeans Algorithm  (0) 2008.07.19
Robust Line Fitting  (0) 2008.07.08
Bayesian Spam Filtering  (0) 2008.07.03
EM : Binarization  (0) 2008.07.01
EM Algorithm : Line Fitting 예  (0) 2008.06.29
Posted by helloktk