동일한 물체나 풍경이라도 촬영에 사용된 센서의 종류나, 전반적인 조명조건때문에 얻어진 사진이 더 어둡거나 더 밝게 표현될 수 있다. 이 경우에 사진을 잘 찍혀진 사진과 비교하여 보정하려면 어떤 기준를 가지고 작업을 해야 할까?  동일한 내용의 영상이라면 단순히 히스토그램을 늘리거나 줄이는 과정(histogram stretching)을 사용할 수 있다. 그러나 같은 물체나 풍경을 담지 않은 기준영상의 히스토그램과 같은 내용의 히스토그램을 같도록 사진을 조작하고 싶으면 단순히 늘리거나 줄여서는 안된다. 

영상의 주어진 그레이 값이 기준영상의 어떤 그레이 값으로 바뀌어야 하는가를 결정해야 하는데, 이때 기준은 주어진 영상의 i번째 그레이 값까지 픽셀이 나타날 확률합과 같은 확률합을 주는 기준영상의 그레이 값이 k일 때 

i(영상) -> k(기준영상)

으로 보내면 될것이다. 이 과정은 히스토그램 스트레칭에도 그대로 적용이 된다. 누적히스토그램(cumulative histogram)을 사용하면 주어진 주어진 그레이 값까지의 확률합을 쉽게 구한다. i-번째 cumulative histogram은 그레이값 0부터 i까지의 히스토그램의 합을 의미한다. 따라서 히스토그램 매칭(histogram matching, histogram specification)은 같은 cumulative histogram의 값을 갖는 그레이값 사이의 변환을 의미한다 (두  영상의 크기가 같지 않은 경우에는, cumulative histogram를 쓸 수 없고, 이 경우에는 전체 픽셀수로 정규화시킨 cumulative distribution function(cdf)을 써야 한다)



히스토그램 평탄화(histogram equalization)은 각 그레이값이 나타날 확률이 균일한 기준영상을 사용한 히스토그램 매칭에 해당한다. 따라서 기준영상이 필요하지 않는다.


//map을 lookup table로 이용해서 영상을 변환할 수 있다.

void matchHistogram(int srcHist[256], int refHist[256], int map[256]) {

    double cumSrc[256], cumRef[256];

    cumSrc[0] = srcHist[0];

    cumRef[0] = refHist[0];

    for (int i = 1; i < 256; i++) {

        cumSrc[i] = cumSrc[i - 1] + srcHist[i];

        cumRef[i] = cumRef[i - 1] + refHist[i];

    }

    // normalize ==> make CDF; 두 비교영상이 같은 크기가 아닐 수 있으므로 필요하다;

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

        cumSrc[i] /= cumSrc[255];

        cumRef[i] /= cumRef[255];

    }

    // 두 비교영상의 cumulant histogram이 가장 가까운(같은) 픽셀끼리 변환한다;

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

        int k = 255;  

        while (cumRef[k] > cumSrc[i]) k--;

        map[i] = k < 0 ? 0 : k;

        // 좀 더 smooth한 interpolation이 필요함;

    }

};

실제 사용은 grey 이미지인 경우에;

int histogramSpec(BYTE *src, int ws, int hw, BYTE *ref, int wr, int hr, BYTE *out) {


Posted by helloktk

삼각형의 외접원을 기하학적인 방법이 아닌 대수적인 방법을 이용해서 구해보자. 중심이 (xc, yc)이고 반지름이 R인 원의 방정식이 


으로 주어진다. 이 식에서 xc, yc, R 세개의 변수를 결정해야는데, 각각의 변수에 대해서 2차 방정식의 형태로 주어져 있으므로 좀 더 쉬운 1차 식의 형태가 되도록 변형을 해보자. 방정식을 전개하면



이제 원의 방정식은 xc , yc, d에 대한 1차 식으로 정리되었다. 삼각형의 외접원의 방정식은 세 꼭지점을 통과하므로 세 꼭지짐 A(x1, y1), B(x2, y2), C(x3, y3)을 위의 식에 넣으면 3개의 방정식을 얻는다. 이 방정식을 행렬로 표현하면


따라서, Cramer의 규칙을 적용하면 다음과 같이 외접원의 중심 (또한 d을 구해서 외접원의 반지름을 구할 수 있다)을 얻을 수 있다.


determinant를 정리하면,


을 얻을 수 있다. 이 식들의 분모는 세 꼭지점에 동일직선 위에 있지 않으면 0이 아니다. 한 꼭지점과 다른 꼭지점간의 좌표의 차이 4개 (x2-x1, y2-y1, x3-x1, y3-y1)와 좌표의 합 4개(x1+x2, y1+y2, x1+x3, y1+y3)를 계산하여서 중심점의 좌표를 얻을 수 있다. 

//

구현된 코드는 http://kipl.tistory.com/113; 또는

int circumcenter2(CPoint P, CPoint Q, CPoint R, double *xc, double *yc) {

double A = Q.x - P.x;

double B = Q.y - P.y;

double C = R.x - P.x;

double D = R.y - P.y;

double E = A * (P.x + Q.x) + B * (P.y + Q.y);

double F = C * (P.x + R.x) + D * (P.y + R.y);

double G = 2. * ( A * D - B * C);

if (G != 0.) { // 세 점이 일직선에 놓이지 않는 경우; 이 경우만 외접원이 정의된다;

*xc = (D * E - B * F) / G;

*yc = (A * F - C * E) / G;

return 1;

}

else return 0;

}

분모 G의 계산을 G = 2 * ( A * (R.y - Q.y) - B * (R.x - Q.x)) 로 쓰는 경우가 많다. 위의 행렬식을 보면 같은 결과이다: G = 2.*(A * (D - B) - B * (C - A)) = 2 * (A * D - B * C).  단지,식이 보다 대칭적으로 보일 뿐이다.

Posted by helloktk

평면 위의 세 점 A, B, C가 만드는 삼각형을 외접하는 외접원의 중심을 구해보자. 점 C을 원점으로 하면, A, B는 각각


의 두 벡터로 표현이 된다 (원점을 C점으로 평행이동하였다고 생각하면 된다). 이 두 벡터에 각각 수직이고 삼각형의 평면에 놓인 (시계방향으로 회전하여서 만든) 두 벡터 m, n를 다음과 같이 정의하자.


그러면, 변 CA의 중점을 지나면서 수직인 직선의 방정식은


또, CB의 중점을 지나면서 수직인 직선의 방정식은 

이 두 직선의 교점이 외접원의 중심이 된다. 매개변수 s, t를 구하기 위해서 두 식을 빼면


 매개변수 s는 b와 n이 수직인 조건을 이용하면


따라서, 외접원의 중심은  



여기서, 3개 벡터의 외적의 성질을 이용해서 (평면에서 외적은 숫자이다)


가 항등식으로 성립하고, 또한 a, m이 서로 수직이면서 길이가 같고, m 도 b에 수직이면 길이가 같으므로


따라서, 

여기서, 윗식의 분모를 보면


이므로, 이 값이 0이 아니기 위해서는 세 꼭지점이 일직선상에 놓여있지 않으면 된다.


그런데 이것은 점 C를 원점으로 하여서 계산을 한 것이므로, 원래의 좌표계에 대한 식으로 바꾸려면 (Cx, Cy)를 더해 주어야 한다.


또는, 성분으로 표현하면



이 설명은 공간상에 놓인 삼각형에 대해서도, 삼각형이 놓인 평면을 xy평면으로 생각하면 적용된다. 단지 조심해야 할 것은 m, n벡터를 구하는 것이다.


int circumcenter(CPoint A, CPoint B, CPoint C, double *xc, double *yc){

double ax = A.x - C.x ;

double ay = A.y - C.y ;

double bx = B.x - C.x ;

double by = B.y - C.y ;

double asq = ax * ax + ay * ay;

double bsq = bx * bx + by * by;

double ccw = ax * by - ay * bx;

if (ccw != 0.) { // 세 점임 일직선 위에 있지 않는 경우; 이 경우만 외접원이 정의됨;

*xc = C.x + (by * asq - ay * bsq) / (2. * ccw) ;

*yc = C.y + (-bx * asq + ax * bsq) / (2. * ccw) ;

return 1;

}

else return 0;

}



Posted by helloktk

세 개의 점 A, B, C가 만드는 삼각형의 외접원의 반지름을 구해보자. 


우선 세점이 일직선상에 놓이지 않기 위해서는 


이어야 한다. 

그림의 삼각형에서 각(BOC)가 각(BAC)의 두배이므로,


이 각도를 이용하면 삼각형의 면적은 세 변의 길이(a, b, c)와 외접원의 반지름(R)로 표현된다:


또한, 삼각형의 면적이 세 점이 직선상에 있는가를 테스트 하는 값으로 아래 식처럼 쓰여지므로,


외접원의 반지름은 다음과 같이 구할 수 있다:


#define SQR(x) ((x)*(x))

double circumradius(CPoint A, CPoint B, CPoint C) {

double ax = C.x - B.x; 

double ay = C.y - B.y;

double bx = A.x - C.x; 

double by = A.y - C.y;

double crossab = ax * by - ay * bx;

if (crossab != 0.) { 

double a = sqrt(SQR(ax) + SQR(ay));

double b = sqrt(SQR(bx) + SQR(by)); 

double cx = B.x - A.x; 

double cy = B.y - A.y;       

double c = sqrt(SQR(cx) + SQR(cy));

return (0.5 * a * b * c/fabs(crossab));

else

return ( -1.0 ); 

}


Posted by helloktk

conic section은 직교좌표계에서 (x, y)에 대한 2차식으로 쓰여진다;



이 conic section의 파라미터 {a,b,c,d,e,f}를 이용해서 타원의 중심과 장축 및 단축의 길이, 그리고 회전각에 대한 공식을 구해보자. 2차 항을 행렬을 써서 표현하면



따라서, 적절한 회전변환을 하면 두 좌표의 곱으로 주어지는 xy-항을 없앨 수 있다.  회전변환은 행렬의 determinant를 보존하고, 그 값은 행렬의 두 eigenvalue의 곱으로 쓰여진다.



따라서, 회전후의 방정식이 타원의 방정식(원점이 이동된)을 기술하기 위해서는 det > 0 이어야 한다 (회전시킨 후 식에서 x^2항과 y^2의 항의 계수는 각각 eigenvalue로 주어지므로 같은 부호를 가져야 한다.)



conic section F가 타원을 기술한다면, 다음과 같이 평행이동 및 회전(rigid transformation)변환을 해서 타원의 중심과 회전각도를 얻을 수 있다



이 식을 대입해서  xy-항이 없어지는 조건과 1차 항이 없어지는 조건을 사용하면


그리고, x와 y의 제곱의 계수는


이 eigenvalue는 위에서 구한 회전각도를 이용해서 표현이 가능하다:



나머지 상수항은



로 주어짐을 알 수있다. 따라서 타원의 각축의 길이는



로 주어진다.

// conic_params(a, b, c, d, e, f);

//ellipse_params(major_axis, minor_axis, center_x, center_y, tilt_angle);

bool conic_to_ellipse(double conic_param[6], double ellipse_param[5]) {

    const double a = conic_param[0];

    const double b = conic_param[1];

    const double c = conic_param[2];

    const double d = conic_param[3];

    const double e = conic_param[4];

    const double f = conic_param[5];

    //get ellipse orientation

    const double theta = 0.5 * atan2(b, a - c);


    //get scaled major/minor axes

    const double ct = cos(theta);

    const double st = sin(theta);

    const double ap = a * ct * ct + b * ct * st + c * st * st;

    const double cp = a * st * st - b * ct * st + c * ct * ct;


    //get translations

    const double cx = (2 * c * d - b * e) / (b * b - 4 * a * c);

    const double cy = (2 * a * e - b * d) / (b * b - 4 * a * c);


    //get scale factor

    const double val = a * cx * cx + b * cx * cy + c * cy * cy;

    const double scale_inv = val - f;


    if (scale_inv / ap <= 0 || scale_inv / cp <= 0) {

        //TRACE("Error! ellipse parameters are imagina\n");

        return 0;

    }


    ellipse_param[0] = sqrt(scale_inv / ap);

    ellipse_param[1] = sqrt(scale_inv / cp);

    ellipse_param[2] = cx;

    ellipse_param[3] = cy;

    ellipse_param[4] = theta;

    return 1;

};

Posted by helloktk
32비트 머신에서 float 형의 변수는 4바이트의 메모리 공간을 차지한다. 즉, int 와 같은 메모리가 할당이 된다. 그리고, 4바이트의 최상위 비트가 1로 세팅이 되면, 이 float형의 변수는 음수를 의미한다. 따라서, float 형의 변수의 절대값을 구하고 싶으면, 메모리에 접근해서 4바이트를 얻어 온 후에, 최상위비트만 0으로 만들고 나머지는 그대로 두어야 한다. 따라서 비트마스크를 이용해서 이 과정을 수행하고 싶으면, and 연산의 비트마스크를 

01111111 11111111 11111111 11111111 = 2^31 - 1 = 0x7F FF FF FF 

처럼 잡으면 된다.

float  fast_abs(float fx) {
    int  ix = *(int *)&fx ;   // 4바이트 얻어 오기;
    ix &= 0x7FFFFFFF;    // MSB 지우기
    return *(float *)&ix ;   // float 형으로 되돌려 줌;
}

 

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

삼각형의 외접원 : 외접원의 반지름  (0) 2012.10.13
Ellipse Parameters  (0) 2012.10.13
float 타입 변수의 절대값은?  (0) 2012.02.17
x 보다 크거나 같은 가장 작은 2^n ?  (0) 2012.02.13
Is Pow of 2  (0) 2012.02.13
Fixed-point RGB2Gray  (0) 2012.01.25
Posted by helloktk
양의 정수  x가 주어질 때, 이보다 크거나 같은 가장 작은 2^n의 숫자를 찾아보자. 왜? FFT 때문이다. 물론,  n = (int)ceil(log(x) / log(2.)) 로 계산하거나, 
                    int a = 1;
                    while (a < x) a <<= 1; 
로 찾을 수 있다;

2^(n-1) < x <= 2^n 사이에 있으면, 원하는 답은 2^n이다. x=2^n 인 경우를 제외하고 모두 n번째 비트가 1로 세팅이 된다. 따라서, 통일시키기 위해서 1을 빼면, x-1 은 n 번째 비트가 항상 1로 주어진다.  이제, 남은 일은 n-1 번째에서 0번째까지 모든 비트를 1로 채우면 된다. 그 결과에 1을 더하면 원하는 2^n을 얻는다.  n-번째 비트가 1이고 이를 이용해서 하위비트를 채우기 위해서는 차례로 >> 연산을 한 후에 or 연산을 하면 된다.

a = x - 1
a                   = 1xxxxxxxxxxxxxxxxxxxxxxxx
a >> 1             = 01xxxxxxxxxxxxxxxxxxxxxxx
a >> 2             = 001xxxxxxxxxxxxxxxxxxxxxx 
.......

31번의 >> 연산을 한 후에 or 연산을 하면 n이하의 모든 비트가 1로 채워진다(31은 최대로의 경우다. 사전에 답을 모르므로 31번 >>까지 해야 한다). 이것은 너무 낭비다.  잘 살펴보면,

a | (a>>1) = 11xxxxxxxxxxxxxxxxxxxxxxxxx

형태이어서, 상위 두 자리 비트가 이미 세팅이 되어 있으므로, 이 결과를 다시  >> 2 하여서 이전 결과와 or연산을 하면 상위 4자리의 비트가가 채워지고, 또다시 이 결과를 >>4 하고 난 후에 직전 결과와 or 연산을 하면 상위 8자리의 비트가 1로 채워진다. 이런식으로 하면 >>16까지 하면 32자리를 모두 커버할 수 있다.

따라서, 전체적인 알고리즘은 (x > 0)
 

   x-- ; 
   x |= (x >> 1);  //상위 2 자리 채움
   x |= (x >> 2);  //상위 4자리
 채움 
   x |= (x >> 4);  //상위 8자리
 채움 
   x |= (x >> 8);  //상위 16자리
 채움 
   x |= (x >> 16);//상위 32자리
 채움 

   x++;
   return x;


단, 32비트 머신에서만이다. 64비트 프로그래밍에서는 >>32 도 추가로 해주어야 한다.

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

Ellipse Parameters  (0) 2012.10.13
float 타입 변수의 절대값은?  (0) 2012.02.17
x 보다 크거나 같은 가장 작은 2^n ?  (0) 2012.02.13
Is Pow of 2  (0) 2012.02.13
Fixed-point RGB2Gray  (0) 2012.01.25
A Fast 2D Otsu Thresholding Algorithm Based on Improved Histogram  (0) 2010.11.17
Posted by helloktk
FFT를 수행하려면 이미지의 폭이나 높이가 2의 지수승인 경우가 가장 단순하게 주어진다. 따라서  주어진 정수 x가 2의 지수승인가 판별할 수 있어야 한다. x가 양의 정수이고, 2^n으로 표현이 된다면, 2진수로 나타낼 때, (n+1) 번째 비트만 1이고 (0부터 센다), 나머지 비트는 모두 0이다. 그리고 x-1은  n 번째에서 0번째까지 모든 비트가 1이 된
다.  
    x                    x(2진수)            x-1
    1                      1                        0
    2                    10                         1
    4                   100                       11
    8                  1000                     111
     ................................................................
이 표를 보면, x와 x-1사이에는 겹치는 비트가 없다. 따라서 두 수를  and 연산을 하면 0 이 되는 경우에는 2의 지수승이고, 그 이외의 경우에는 0이 아님을 알 수 있다. x = 2의 지수승  판별은    
 

                           return  x & (x - 1) == 0         /* x != 0 인 정수*/
 

인가를 보면 된다. 

그런데 이 판별식은 x = 0 인 경우에는 성립이 안된다. x-1 = -1 이므로 32비트 자리 전부가 1로 채워지므로 x & (x-1) = 0 이어서 2의 지수승으로 판별한다. 따라서 0을 제외하는 방법을 찾아야 한다. (물론 함수 인자에서 양수로 제한을 하면 되지만 폼이 안난다). 음수는 최상위비트가 1로 채워지므로 이 사실을 이용하자. 최상위 비트를 1로 만들려면,

~0U                    =111111111111111..11111111(32개) 
~0U>>1              = 011111111111111..11111111
~(~0U>>1)          =100000000000000..00000000 
~(~0U>>1)|x       = x의 최상위 비트를 항상 1로 채워준다(음수 일때는 자동으로 만족)
                             나머지 비트는 그대로 둔다.

따라서 이 값과 x-1을 and 연산을 하면 0 이하인 수가 들어오면 연산이 결과를 항상 0이 아니게 된다.


                          return  !((~(~0U>>1)|x) & (x - 1))          /*x = 정수 */    
 
0 하나를 예외처리하기 위해서 너무 많은 과정을 거치는 것이 아닌가? 
Posted by helloktk

matlab에서 rgb 컬러이미지를 gray 이미지로 바꿀 때, 호출하는 함수가 rgb2gray이다. 이 함수는 주어진 RGB값에 다음과 같은 weight를 주어서 명암값을 얻는다:

gray = 0.2989 * r + 0.5870 * g + 0.1140 * b ;

이 계산을 fixed-point 버전으로 바꾸어 보자. 정밀도를 유지하기 위해서는 소수점 이하 4자리까지 유지해야 하는데, 이것은 weight에 각각 10000을 곱한 값을 사용한 후에 다시 10000으로 나누면 된다 (10000을 곱하더라도 r, g, b가 0-255사이의 값이므로 최대로 255 * 10000 보다 작은 값이 나와서 32-비트 정수 범위 내에 있게 된다)

gray = (2989 * r + 5870 * g + 1140 * b) / 10000;

여기서 좀 더 개선을 할 수 있다. 10000으로 나누는 과정을 shift 연산으로 바꾸면 된다. 정밀도를 보존하기 위해서 10000보다 큰 2의 power의 수를 찾으면 2^14 = 16384 이 주어진다. 따라서 10000 대신에 2^14을 곱한 weight를 쓰고, 계산 결과를 오른쪽으로 14만큼 shift 연산을 하면 된다.

0.2989 * 2^14 = 4897.1776;
0.5870 * 2^14 = 9617.408;
0.1140 * 2^14 = 1867.776;

gray = (4897 * r + 9617 * g + 1868 * b) >> 14;

weight의 총 합 = 4897 + 9617 + 1868 < 2^14 이어서 gray 값은 255를 넘지 않고, 중간 계산은 32-비트 정수 범위 내에서만 이루어 진다.

정말 빠를까?
Posted by helloktk

A Fast 2D Otsu Thresholding Algorithm Based on Improved Histogram
Ningbo Zhu, Gang Wang, Gaobo Yang, Weiming Dai
Pattern Recognition Letters 16 (1995) 653-666

Based on the gray-level histogram of an image, Otsu’s thresholding method takes the variance between clusters as the criterion to select the optimal threshold. This method has been cited as an effective thresholding technique and widely used in real thresholding tasks. However, there are always uncertain
disturbing factors in practical applications. These disturbing factors probably make one-dimensional (1D) histogram fail to possess obvious vally between two peaks. As a result, Otsu method can not obtain optimal segmentation results. To solve this problem, Liu et. al.[5] extended 1D histogram into
two-dimensional(2D) histogram, which utilizes not only gray value distribution of an image, but also the dependency of pixel in its neighborhood. This method greatly improved the segmentation results, specially in low SNR condition.

그레이 이미지의 2차원 히스토그램은 원래 이미지의 픽셀값 분포를 x-축으로, 주변 픽셀값의 평균을 이용해서 만든 이미지에서의 픽셀값의 분포를 y-축으로 하여서 만든 일종의 2차원 영상이다(256x256). 
원본이미지에서 전경 픽셀(검정)은 주변 픽셀을 평균한 값도 전경에 비슷한 값을 가져야 한다. 그러나 그 픽셀이 배경(흰색)에 나타난 검정색의 노이즈라면, 주변 픽셀의 평균값은 흰색에 가까운 값을 가지게 된다. 반대로, 배경픽셀이면, 주변 픽셀의 평균값도 배경에 비슷한 값을 가져야 하고, 만약에 그것이 전경에 놓인 흰색잡음이었다면, 주변픽셀의 평균값은 검정으로 나타날 것이다(전경과 배경의 경계 근처에 놓인 픽셀들도 마찬가지다). 따라서 이 2차원 히스토그램을 보면 잡음이 많은 경우에 어떤 픽셀값들을 노이즈로 처리할 수 있는가를 보다 명확하게 알 수 있다.
위 2차원 히스토그램의 x-축은 원본, y-축은 주변픽셀을 평균이미지의 픽셀값 분포(0~255=L)를 표현한다. A 영역에 해당하는 픽셀은 전경(검정)픽셀로 생각할 수 있고, C-영역에 해당하는 픽셀은 배경(흰색)으로 생각할 수 있다. B영역에 해당하는 픽셀은 배경에 놓인 검정색 노이즈이고, D-영역은 전경에 놓인 흰색 노이즈 또는 전경과 배경의 경계지점이라고 생각할 수 있다(전경이 배경보다도 작은 영역을 차지하는 경우에)

이미지크기 : 512x512;

픽셀평균을 구하는 윈도우 크기를 5x5로 할 때, 2차원 히스토그램(돗수를 log-스케일로 바꾸어서 표현한 것임. 붉은색일수록 돗수가 많음을 의미한다):
Posted by helloktk
"Optimizing Connected Component Labeling Algorithms"
Kesheng Wu, Ekow Otoo and Arie Shoshani
Lawrence Berkeley National Laboratory, Berkeley, CA, USA

ABSTRACT
This paper presents two new strategies that can be used to greatly improve the speed of connected component labeling algorithms. To assign a label to a new object, most labeling algorithms use a scanning step that examines some of its neighbors. The first strategy exploits the dependencies among the neighbors to reduce the number of neighbors examined. When considering 8-connected components in a 2D image, this can reduce the number of neighbors examined from four to one in many cases. The second strategy uses an array to store the equivalence information among the labels. This replaces
the pointer based rooted trees used to store the same equivalence information. It reduces the memory required and also produces consecutive final labels. Using an array based instead of the pointer based rooted trees speeds up the connected component labeling algorithms by a factor of 5 100 in our tests on random binary images.

Keywords: Connected component labeling, Union-Find, optimization


Posted by helloktk

제대로 세그먼트된 그레이 영상은 원래의 영상이 나타내고자 하는 전경이 제대로 표현이 된 것이다. 이것은 원래 영상과 세그먼트된 영상의 상관관계가 높아야 함을 의미한다. 따라서 세그먼트를 위한 임계값의 기준으로 이 상관계수를 최대로 하는 임계값을 찾는 것도 좋은 기준중의 하나가 될 수 있다.

 

여기서 사용할 상관계수는 원래의 영(A)과 전경과 배경을 그들의 평균 그레이값으로 대체한 세그먼트된 영상(B)간의 상관계수를 사용한다. 세그먼트된 영상 B는 임계값이 T인 경우에

B(i, j) =  m0   if A(i, j) <= T

            m1    otherwise ;

로 정의된다. 여기서 m0는 배경부분의 평균 그레이값이고, m1은 전경부분의 평균 그레이값을 의미한다. 이 값은 임계값 T에 따라서 달라진다. 임계값이 높으면 m0는 커지고, 반대로 m1은 작아진다.

임계값이 T일 때, 배경의 픽셀수 비를 p 라고 하고, 전경의 픽셀수 비를 q(=1- p) 라고 하면

E(A) = E(B) = m = 그레이_평균 = p*m0 + q*m1;

V(A) = 그레이_분산 = T-값에 무관하게 일정;

V(B) = p*m02 + q*m12 – m2 = p*q*(m0 – m1)2;

E(A, B) = p*m02 + q*m12 ;

이므로,

Correlation(A, B)

= (E(A,B) – E(A)*E(B))/sqrt(V(A)*V(B))

= sqrt(p*q*(m0 – m1)2) / sqrt(V(A));

= sqrt(p*q*(m0 – m1)2) up to constant factor;

= interclass variance;

= Otsu’s criteria

, 원래의 그레이영상 A와 전경과 배경을 각각의 평균값으로 대체한 영상간의 상관계수는 전경과 배경 두 클래스간의 분산이 최대일 때, 가장 크게 나타난다. 이 기준은 Otsu 알고리즘에서 사용한 기준과 같다.

참고: Otsu Algorithm 구현 예.

Posted by helloktk

영상에서 전경물체가 어떤 방향을 정렬이 되어있는가를 찾는 문제는 여러 응용에서 나타난다. 예를 들면, 영상에서 사람의 머리가 어떤 자세를 취하고 있는가를 묻는 것에 대한 답이나, 손바닥인식에서 손이 가리키는 방향에 대한 정보를 제공한다. 

물체가 있을 때 이 물체의 정렬방향(orientation)의 의미는 물체를 구성하는 각각의 픽셀에서 정렬방향을 정의하는 직선까지의 거리의 합이 최소인 직선의 방향/기울기를 의미한다. 물론 이 직선은 물체의 질량중심 (xc, yc)을 지나야 한다.

아래 그림에서 물체의 중심을 지나고 각도가 Θ만큼 기울어진 직선에 대해서 생각해보면, 임의의 픽셀 (i, j)에서 직선까지의 거리는

으로 주어진다. (직선에 수직한 단위벡터 (-sinΘ, cosΘ)에 대한 (i-xc, j-yc) 정사영임을 생각하면 쉽게 이해할 수 있다) 

따라서, 우리가 구하고자 하는 orientation은 최소자승법의 의미에서  거리의 제곱을 전체 픽셀에 대해서 더한 값을 최소화 시키는 각도 Θ를 구하는 걸로 귀결된다.

 

따라서, S(Θ)를 Θ에 대해서 미분을 한 후에 정리하면,

로 주어짐을 알 수 있다. μpq는 영상의 (p, q)차원 central momemt이다. 

물론, 이들 central moment을 이용해서 만든 공분산행렬(covariance matrix)의 

eigenvector / eigenvalue 두 개를 구하면 큰 eigenvalue을 갖는 eigenvector의 방향이 물체의 정렬방향이 된다.


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

Optimizing Connected Component Labeling Algorithms  (0) 2010.10.23
Otsu-알고리즘의 새로운 해석  (0) 2010.01.28
Object Orientation  (1) 2010.01.17
Bicubic Interpolation  (1) 2010.01.14
Histogram Smoothing via Bezier Curve  (1) 2010.01.10
Running Median Filter  (0) 2010.01.07
Posted by helloktk

이미지처리에서 픽셀좌표는 간격이 1인 2차원의 그리드 교차점들로 볼 수 있다. 이미지를 확대하거나 축소할 때, 픽셀과 픽셀사이의 값이 필요한 경우가 발생한다. 이러한 경우에 보통은 주변의 가장 가까은 픽셀의 값을 가져다 쓰거나, 아니면 주변의 픽셀값을 선형보간하여서 쓴다. 2차원 그리드에서 선형보간은 bilinear interpolation이다. 이러한 보간법들은 속도면에서 빠르지만, 픽셀값의 분포가 부드럽게 이어지지 않는 단점이 있다. 픽셀값의 분포가 부드럽게 이어지기 위해서는 1차 미분이 적어도 연속인 보간법을 사용해야하는데, 이러한 조건을 만족시키는 가장 낮은 찾수의 다항식 보간법이 bicubic interpolation이다.

네점 (0,0), (0,1), (1,0),(1,1)을 꼭지점으로하는 정사각형내의 임의의 지점(x, y) 0<=x<=1, 0<= y <=1에서의 픽셀값을 주변의 16개 지점 (i,j) -1<= i <=2, -1 <= j <= 2에서의 픽셀값을 이용해서 구하도록하자. 픽셀값 f(x,y)은 x와 y의 3차 함수로 주어지진다

이제, 문제는 3차 다항식의 계수 aij = a(i,j) 16 개를 구하는 것이다.  따라서 16개의 조건이 필요한데 이것을 주변의 16개의 픽셀값을 이용해서 만들 수 있다.

1. f(x, y) : 4 꼭지점에서 원래의 값:
    f(0,0) = a00;
    f(1,0) = a00 + a10 + a20 + a30;
    f(0,1) = a00 + a01 + a02 + a03;
    f(1,1) = a00 + a01 + a02 + a03 + a10 + a11 + a12 + a13 + a20 + a21 + a22 + a23 + a30 + a31 + a32 + a33;
2. 4 꼭지점에서 x/y-미분값:
    fx(0,0) = a10;
    fx(1,0) = a10 + 2*a20 + 3*a30;
    fx(0,1) = a10 + a11 + a12 + a13;
    fx(1,1) = a10 + 2*a20 + 3*a30 + a11 + 2*a21 + 3*a31 + a12 + 2*a22 + 3*a32 + a13 + 2*a23 + 3*a33;
    fy(0,0) = a01;
    fy(1,0) = a01 + a11 + a21 + a31;
    fy(0,1) = a01 + 2*a02 + 3*a03;
    fy(1,1) = a01 + a11 + a21 + a31 + 2*a02 + 2*a12 + 2*a22 + 2*a32 + 3*a03 + 3*a13 + 3*a23 + 3*a33;
3. 4 꼭지점에서 교차미분값:
    fxy(0, 0) = a11;
    fxy(1, 0) = a11 + 2*a21 + 3*a31;
    fxy(0, 1) = a11 + 2*a12 + 3*a13;
    fxy(1, 1) = a11 + 2*a21 + 3*a31 + +2*a12 + 4*a22 + 6*a32 + 3*a13 + 6*a23 + 9*a33;

(x,y) = (0,0), (1,0), (0,1), (1,1) 지점에서 f-값이나 f-미분값은 계수 a(i,j)로 표현이 되고, 또 이것은 꼭지점에서 픽셀값 및 인접픽셀값과의 차이로 표현(미분의 경우)되므로 16개의 다항식의 계수를 16차원 벡터 v로
    v=[a00, a10, a20, a30, a01, a11, a21, a31, a02, a12, a22, a32, a03, a13, a23, a33];
각각의 꼭지점에서의 픽셀값, 미분값을 f 벡터로 나타내면
    f =[f00, f10, f01, f11, fx00, fx10, fx01, fx11, fy00, fy10, fy01, fy11, fxy00, fxy10, fxy01, fxy11];
방정식으로 해는 
    v = A.f 
로 쓰여진다. 이때 16x16 행렬 A는 
   A= {1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0},
        {-3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0},
        {2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0},
        {0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0},
        {0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1},
        {0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1},
        {-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0},
        {9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2},
        {-6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2},
        {2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0},
        {0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0},
        {-6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1},
        {4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1}
이 결과는 mathematica나 maple등의 프로그램을 이용하면 쉽게 답을 구할 수 있고, 출력을 C-형태로 얻을 수 있다. 픽셀에서의 미분값은 인접픽셀값의 평균변화율로 정의할 수 있는데, 예를 들면
fx00 = fx(0,0) = ( p(1,0) - p(-1,0))/2,.......
fxy00 = fxy(0,0) = (p(1,0) - p(1,-1) - p(-1, 1) + p(-1,-1)) / 4....
와 같다.

/* 픽셀의 위치가 grid의 교차점이 되도록 표현하였다*/


/* [ix, ix+1][iy+iy+1] 영역에서 보간을 하는 경우에, 주변픽셀의 정보는
p[a][b]에 들어가고, p[a][b]= Image (ix - 1 + a, iy - 1 + b)*/
double BicubicInterpolate(double p[4][4], double x, double y) { /* 0 <= x <= 1, 0 <= y <= 1 */
    double a00 = p[1][1];
    double a01 = (-p[1][0] + p[1][2])/2.;
    double a02 = p[1][0] - (5*p[1][1])/2. + 2*p[1][2] - p[1][3]/2.;
    double a03 = (-p[1][0] + 3*p[1][1] - 3*p[1][2] + p[1][3])/2.;
    double a10 = (-p[0][1] + p[2][1])/2.;
    double a11 = (p[0][0] - p[0][2] - p[2][0] + p[2][2])/4.;
    double a12 = (-2*p[0][0] + 5*p[0][1] - 4*p[0][2] + p[0][3] + 
                         2*p[2][0] - 5*p[2][1] + 4*p[2][2] - p[2][3])/4.;
    double a13 = (p[0][0] - 3*p[0][1] + 3*p[0][2] - p[0][3] -
                        p[2][0] + 3*p[2][1] - 3*p[2][2] + p[2][3])/4.;
    double a20 = p[0][1] - (5*p[1][1])/2. + 2*p[2][1] - p[3][1]/2.;
    double a21 = (-2*p[0][0] + 2*p[0][2] + 5*p[1][0] - 5*p[1][2] -
                         4*p[2][0] + 4*p[2][2] + p[3][0] - p[3][1])/4.;
    double a22 = (4*p[0][0] - 10*p[0][1] + 8*p[0][2] - 2*p[0][3] - 10*p[1][0] + 
                      25*p[1][1] - 20*p[1][2] + 5*p[1][3] + 8*p[2][0] -
                      20*p[2][1] + 16*p[2][2] - 4*p[2][3] - 2*p[3][0] + 7*p[3][1] - 6*p[3][2] + p[3][3])/4.;
    double a23 = (-2*p[0][0] + 6*p[0][1] - 6*p[0][2] + 2*p[0][3] + 5*p[1][0] - 15*p[1][1] +
                        15*p[1][2] - 5*p[1][3] - 4*p[2][0] + 12*p[2][1] - 12*p[2][2] +
                          4*p[2][3] + p[3][0] - 4*p[3][1] + 4*p[3][2] - p[3][3])/4.;
    double a30 = (-p[0][1] + 3*p[1][1] - 3*p[2][1] + p[3][1])/2.;
    double a31 = (p[0][0] - p[0][2] - 3*p[1][0] + 3*p[1][2] + 
                     3*p[2][0] - 3*p[2][2] - p[3][0] + p[3][1])/ 4.;
    double a32 = (-2*p[0][0] + 5*p[0][1] - 4*p[0][2] + p[0][3] + 6*p[1][0] - 15*p[1][1] + 
                        12*p[1][2] - 3*p[1][3] - 6*p[2][0] + 15*p[2][1] -
                        12*p[2][2] + 3*p[2][3] + 2*p[3][0] -  7*p[3][1] + 6*p[3][2] - p[3][3])/4.;
    double a33 = (p[0][0] - 3*p[0][1] + 3*p[0][2] - p[0][3] - 3*p[1][0] + 9*p[1][1] - 9*p[1][2] +
                     3*p[1][3] + 3*p[2][0] - 9*p[2][1] + 9*p[2][2] - 3*p[2][3] -
                        p[3][0] + 4*p[3][1] - 4*p[3][2] + p[3][3])/4 ;
        
    double x2 = x * x;
    double x3 = x2 * x;
    double y2 = y * y;
    double y3 = y2 * y;
    //
    return a00 + a01 * y + a02 * y2 + a03 * y3 +
        a10 * x + a11 * x * y + a12 * x * y2 + a13 * x * y3 +
        a20 * x2 + a21 * x2 * y + a22 * x2 * y2 + a23 * x2 * y3 +
        a30 * x3 + a31 * x3 * y + a32 * x3 * y2 + a33 * x3 * y3;
}
*원본 16x16 이미지:

bicubic interpolation 결과 (256x256): (경계지점에서 P를 구성하는데 필요한 값은 경계값을 반복해서 적용하였다. 히스토그램 스트레치한 후의 결과이다. 필셀 간격을 확대한 것이 아니라 픽셀 수를 확대한 것이므로 확대비율이 큰 경우에 오른쪽과 하단은 2 픽셀을 두 중복해서 사용한 결과가 나온다. 폭/높이가 큰 이미지에서는 문제가 안되는 작은 이미지에서는 문제가 된다. 이것을 방지하려면, 픽셀의 값을 그리드의 중앙에서의 값으로 생각하고 보간을 하면 된다) 


*simple stretch 결과 (256x256) :




코드 구현 일부:

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

Otsu-알고리즘의 새로운 해석  (0) 2010.01.28
Object Orientation  (1) 2010.01.17
Bicubic Interpolation  (1) 2010.01.14
Histogram Smoothing via Bezier Curve  (1) 2010.01.10
Running Median Filter  (0) 2010.01.07
Fast Median Filter 3x3  (0) 2009.11.27
Posted by helloktk
이미지의 히스토그램에서 피크를 찾는 작업을 하기 위해서 보통 히스토그램을 평균필터나 가우시안 필터로 컨볼류선을 한 후에 사용한다. 여기서는 히스토그램의 각각의 레벨에서의 값(=(레벨, 히스토그램값))을 베지어 컨트롤 포인트로 하여서 만든 베지어 커브로 히스토그램을 근사하려고 한다. 따라서 베지어곡선은 255 찻수의 곡선이 된다. 높은 차수의 베지어커브를 계산하기 위해서 Berstein 함수를 그대로 사용하는 경우에는 수치에러가 발생하여 값이 불안정해진다. 이런 경우에는 De Casteljau's algorithm을 이용하여서 반복적으로 계산을 하면 된다.
//
int SmoothenHistogram (int hist[], int numLevels/*=255*/) {
    // deCasteljau's algorithm
    std::vector<double> p(numLevels+1);
    std::vector<double> hist2(numLevels+1);
    // back-up;
    for(int k = 0; k <= numLevels; k++)
        hist2[k] = hist[k];
    for (int j = 0; j <= numLevels; j++) {
        double t = j*1.0/(numLevels+1);
        for(int i=0; i <= numLevels; i++) p[i] = hist2[i];
        // numLevels-order;
        for(int k = 1; k <= numLevels; k++) {
            for(int i = 0; i < (numLevels + 1 - k); i++)
                p[i] = (1 - t) * p[i] + t * p[i+1] ;
        }
        // final interploation --> p[0];
        hist[j] = (int) p[0];
    }
    return 1;
};

빨간색 프로파일이 원래 데이터이고, 파란색  경계가  Bezier smoothing 을 적용한 결과이다.

베지어곡선 위키피디아: http://en.wikipedia.org/wiki/B%C3%A9zier_curve

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

Object Orientation  (1) 2010.01.17
Bicubic Interpolation  (1) 2010.01.14
Histogram Smoothing via Bezier Curve  (1) 2010.01.10
Running Median Filter  (0) 2010.01.07
Fast Median Filter 3x3  (0) 2009.11.27
Fant's Resampling  (0) 2008.12.17
Posted by helloktk

영상처리에서 한 픽셀의 수정을 위해서 주변픽셀의 정보를 요구하는 윈도우에 기반한 필터들은 일반적 연산비용이 큰 편이다. 한변의 길이가 W인 윈도우를 사용할 때  W^2의 번의 픽셀을 참조해야 하므로 윈도우가 클수록 그 비용이 제곱으로 커진다. 
선형 필터 중에는 x방향과 y 방향으로 연산이 분리가 가능해서 이 연산비용이 필터의 크기에만 비례하도록 만들 수 있다. 메디안 필터는 이런 분리가 안되는 비선형 필터중 하나다. 물론 근사적으로 x방향으로 메디안 필터링을 하고, 그 결과를 y방향으로 다시 메디안 필터링을 하여서 근사적인 메디안 필터를 구성하여 사용하기도 한다.  
그러나 윈도우가 움직이면서 윈도우내의 모든 픽셀이 다 바뀌는 것이 아니라 움직이는 방향에 수직인 가장자리 픽셀만 바뀐다는 사실을 이용하면 각각의 픽셀위의 윈도우에서 메디안을 찾는 작업을 할 필요가 없다.
예를 들면 스캔라인 방향(x-방향)으로 윈도우를 움직이면서 메디안 필터를 작용할 때, 윈도우가 오른쪽으로 1픽셀 움직이면 윈도우의 왼쪽 가장자리의 수직픽셀들은 새 윈도우에서 사라지고, 기존윈도우의 오른쪽 수직가장자리 앞의 픽셀들이 새로이 들어온다. 따라서 처음 한번만 윈도우에 대해서 메디안을 찾으면, 이 나가는 픽셀과 새로 들어오는 픽셀들만 고려하면 (2*W개)되므로 연산시간이 윈도우에 선형적으로만 늘어나게 되어서 연산비용을 많이 줄일 수 있다. 이 기법은 윈도우가 사각형인 필터들(평균, 최대, 최소...)에 적용이 가능하다.

아래의 코드는 이것을 구현한 것이다. 
-영상의 가장자리에서 윈도우폭의 절반에 해당하는 픽셀들은 따로 취급 해주어야 한다;
-y값의 변화는 영상을 일차원배열로 취급하였기 때문에 항상 y*width의 형태로만 나타난다. 따라서 곱셈을 없애기 위해서 모든 y의 변화는 width를 곱한 형태로 구현되게 했다.
-좌표값은 윈도우의 top-left의 위치를 기준으로 만들었고, 적용 픽셀은 윈도우의 중심픽셀이다.

// 경계에서의 픽셀처리가 없다;
// 윈도우의 topleft를 기준으로 작동하도록 만들어졌다(픽셀값은 윈도우 중심에서);
// y --> y * width 단위로 변경하여서 곱셈연산을 줄였다;
BOOL RunningMedianFilter(BYTE* image, int width, int height, int wSize/*=7*/, BYTE* out) {
    int hist[256] ;
    int hwSize = wSize/2;               // half-windows-size;
    wSize  = (hwSize<<1)+1;             // wSize=odd-number;
    int nW2     = wSize * wSize / 2;        // half number of pixels in the windows;

    int xoff = hwSize ;                 //outpixel의 x-출발점;
    int xEnd = width - wSize;           //window topleft의 한계;
    int yoff = hwSize * width ;         //outpixel의 y-출발점;
    int yEnd = (height - wSize) * width;

    //window의 topleft 위치를 옮기면서,  윈도우 중심위치의 픽셀의 median값을 찾는다;
    for (int y = 0 ; y < yEnd; y += width) {
        memset(hist, 0, sizeof(hist));
        //각각이 행에서 첫번째 윈도우에서의 메디안값은 정상적으로 구한다;
        int iyEnd = y + wSize * width;          //window-bottom;
        for (int iy = y ; iy < iyEnd; iy += width) {
            for (int ix = 0; ix < wSize; ix++)
                hist[image[iy + ix]]++ ;
        }
        //첫 윈도우에서 메디안 찾기;
        int median = 0;
        int ltmed = 0;                  //메디안 값보다 작은 픽셀의 수; <= nW2 이어야 한다;
        while(ltmed < nW2) {
            ltmed += hist[median];
            median++ ;
        };//==>ltmed > nW일 때도 생김-->이것을 없애기 위해서;
        while(ltmed > nW2) {
            median--;
            ltmed -= hist[median];
        };
        out[y + yoff + xoff] = median;
        //
        for (int x = 1; x < xEnd; x++) {
            for (iy = y; iy < iyEnd ; iy += width) {
                //왼편에서 기존의 윈도우에서 나가는 픽셀 제외;
                int v = image[iy + x - 1];
                hist[v]-- ;
                if(v < median) ltmed--;
                //오른편에서 들어오는 픽셀 추가;
                v = image[iy + x + wSize - 1] ;
                hist[v]++ ;
                if(v < median) ltmed++ ;
            }
            //update median;
            //ltmed <= nW2;
            while(ltmed < nW2) {
                ltmed += hist[median];
                median++ ;
            }; //ltmed>nW일 수 있으므로;
            while(ltmed > nW2) {
                median--;
                ltmed -= hist[median];
            }
            //윈도우 중심에서의 픽셀값;
            out[y + yoff + x + xoff] = median ;
        }
    }
    //결과이미지에서 경계영역이 처리안되어 있으므로 이부분에 대한 처리가 필요함;
    return TRUE ;
};
wSize = 11 결과(가장자리 5픽셀은 원본을 그대로 복사한 것이다)

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

Bicubic Interpolation  (1) 2010.01.14
Histogram Smoothing via Bezier Curve  (1) 2010.01.10
Running Median Filter  (0) 2010.01.07
Fast Median Filter 3x3  (0) 2009.11.27
Fant's Resampling  (0) 2008.12.17
Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Posted by helloktk
REF:  http://graphics.cs.williams.edu/papers/MedianShaderX6/
/*
3x3 Median
Morgan McGuire and Kyle Whitson
http://graphics.cs.williams.edu
*/

// Input texture
uniform sampler2D T;

#ifndef G3D
// vec2(1/width, 1/height) of the texture
    uniform vec2 Tinvsize;
#else
#    define Tinvsize g3d_sampler2DInvSize(T))
#endif 


// Change these 2 defines to change precision,
#define vec vec3
#define toVec(x) x.rgb

//#define vec vec4
//#define toVec(x) x.rgba

#define s2(a, b)		temp = a; a = min(a, b); b = max(temp, b);
#define mn3(a, b, c)		s2(a, b); s2(a, c);
#define mx3(a, b, c)		s2(b, c); s2(a, c);

#define mnmx3(a, b, c)		mx3(a, b, c); s2(a, b);                         // 3 exchanges
#define mnmx4(a, b, c, d)	s2(a, b); s2(c, d); s2(a, c); s2(b, d);         // 4 exchanges
#define mnmx5(a, b, c, d, e)	s2(a, b); s2(c, d); mn3(a, c, e); mx3(b, d, e); // 6 exchanges
#define mnmx6(a, b, c, d, e, f) s2(a, d); s2(b, e); s2(c, f); \
                                mn3(a, b, c); mx3(d, e, f);                // 7 exchanges

void main() {

  vec v[9];

  // Add the pixels which make up our window to the pixel array.
  for(int dX = -1; dX <= 1; ++dX) {
    for(int dY = -1; dY <= 1; ++dY) {		
      vec2 offset = vec2(float(dX), float(dY));
		    
      // If a pixel in the window is located at (x+dX, y+dY), 
      // put it at index (dX + R)(2R + 1) + (dY + R) of the
      // pixel array. This will fill the pixel array, 
      // with the top left pixel of the window at pixel[0] and the
      // bottom right pixel of the window at pixel[N-1].
      v[(dX + 1) * 3 + (dY + 1)] = toVec(texture2D(T, gl_TexCoord[0].xy + offset * Tinvsize));
    }
  }

  vec temp;

  // Starting with a subset of size 6, remove the min and max each time
  mnmx6(v[0], v[1], v[2], v[3], v[4], v[5]);
  mnmx5(v[1], v[2], v[3], v[4], v[6]);
  mnmx4(v[2], v[3], v[4], v[7]);
  mnmx3(v[3], v[4], v[8]);
  toVec(gl_FragColor) = v[4];

}

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

Histogram Smoothing via Bezier Curve  (1) 2010.01.10
Running Median Filter  (0) 2010.01.07
Fast Median Filter 3x3  (0) 2009.11.27
Fant's Resampling  (0) 2008.12.17
Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (1) 2008.07.14
Posted by helloktk


// 배열첨자(dj)와 픽셀의 실제위치(srcindex, destindex)를 따로 분리하여서
// 열방향을 따라서 작업을 하더라도 메모리 복사가 필요없이 처리가 가능하도록 하였음.
BOOL
resampleRGB(BYTE *source, BYTE* dest, int length,    // 
               int stride,                                                   //RGB=3(행방향), 한행의 총 바이트수(열방향) 
               double outpos[length+1])                            //forward map:source->dest; 
{
    double *inpos = new double[length+1];
    //si=source index;
    //dj=dest index;
    // Find inverse map; dest--> source;
    for (int si = 0, dj = 0; dj < length; dj++) {
        while (si < (length-1) && outpos[si+1] < dj) si++;       //max(si)=length-1;
        //
        if(si < length-1) {
        inpos[dj] = si + (dj - outpos[si]) / (outpos[si + 1] - outpos[si]);
        }else{
            //max(inpos[dj])=lenght-1;
            inpos[dj]=si+1 ;
        }
    }
    inpos[length] = length;

    double inseg  = 1.0;                //첫입력픽셀은 완전히 사용가능;
    double outseg = inpos[1];       //첫출력을 위해서 필요한 입력의 양은
                                                //inpos[1]-inpos[0]이나 inpos[0]=0으로 함;
    double sizefac = outseg ;
    int srcindex = 0;
    int destindex = 0;
    //첫번째 픽셀;
    int b = source[0] ;
    int g = source[1] ;
    int r = source[2] ;
    srcindex += stride ;
    //두번째 픽셀;
    int nextb = source[srcindex+0] ;
    int nextg = source[srcindex+1] ;
    int nextr = source[srcindex+2] ;
    srcindex += stride ;
    //
    double bsum = 0, gsum = 0, rsum = 0;
    for(dj = 1; dj < length; ) {
        //linear interpolation;
        double bintensity = inseg* b + (1-inseg)*nextb ;
        double gintensity = inseg* g + (1-inseg)*nextg ;
        double rintensity = inseg* r + (1-inseg)*nextr ;
        //
        if(inseg < outseg) {//output cycle;
            //accumulation of weighted contrib;
            bsum += inseg * bintensity ;
            gsum += inseg * gintensity ;
            rsum += inseg * rintensity ;
            //기존에 현재의 픽셀을 추가하므로, 출력을 위한 입력픽셀의 소요량이 줄어든다.
            outseg -= inseg ;
           
            //현재의 현재의 픽셀값을 갱신한다.
            b = nextb ;
            g = nextg ;
            r = nextr ;
            //새로 들어올 입력픽셀을 준비한다.
            inseg = 1.0 ;
            //srcindex < endindex-2;
            nextb = source[srcindex+0];
            nextg = source[srcindex+1];
            nextr = source[srcindex+2];
            srcindex += stride ;
        }else{ //input cycle;
            //accumulation of weighted contrib;
            bsum +=outseg* bintensity ;
            gsum +=outseg* gintensity ;
            rsum +=outseg* rintensity ;
            //hack;
            if(sizefac==0) sizefac = 1;
            //set dest pixel if inpos[dj]>=0; dj=1,2,3,....;
            //source가 dest 내부로 들어가는 경우임;
            if(inpos[dj-1]>=0){
                //x=0,1,2...
                // 출력픽셀을 만드는데 sizefac만큼의 입력픽셀이 들어갔으므로 나누어 주어야 한다.
                dest[destindex+0] = (BYTE)min(bsum/sizefac, 0xFF);
                dest[destindex+1] = (BYTE)min(gsum/sizefac, 0xFF);
                dest[destindex+2] = (BYTE)min(rsum/sizefac, 0xFF);
            }
            destindex += stride ;

            // reset accumulator for next output;
            bsum = gsum = rsum = 0;
            // source has been consumed outseg fraction;
            // 현재의 입력픽셀이 다음 출력픽셀에 기여할 수 있는 남아 있는 양;
            inseg -= outseg ;
            // set new outseg; = 다음 출력픽셀을 완성하는데 필요한 입력픽셀의 양;
            outseg = inpos[dj+1] - inpos[dj] ;
            // 출력픽셀을 완성하는데 필요한 입력픽셀의 양(outseg는 다음 출력시까지 계속 변하므로
            // 처음 세팅할 때 기억해 두어야 한다;
            sizefac = outseg ;
            dj++ ;
        }
    }

    delete[] inpos ;
    return TRUE;
}

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

Running Median Filter  (0) 2010.01.07
Fast Median Filter 3x3  (0) 2009.11.27
Fant's Resampling  (0) 2008.12.17
Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (1) 2008.07.14
Histogram Equalization  (0) 2008.06.22
Posted by helloktk

Histogram equalization(HE)은 주어진 이미지의 픽셀분포가 모든 픽셀값에서 같은 확률로 나타나도록 픽셀값의 변경시켜서 이미지를 보다 잘 인식되도록하는 처리과정이다. 그리고, 이러한 균일한 분포는 또한 엔트로피의 관점에서 살펴볼 때 최대의 엔트로피값을 주는 분포이기도 하다. 그러나, HE는 원본이미지의 밝기를 유지하지 않는다. 따라서, 보다 현실적으로 원본이미지의 밝기를 유지하면서, 엔트로피를 최대화시키는 히스토그램 분포를 찾아서, 그것으로 변환을 시도하는 것이 방법을 고려해 볼 수 있다. 그러한 목표 히스토그램을 f(s) (연속적인 확률밀도함수로 생각함)라고 하면,

히스토그램일 조건 ;  

밝기 보존: 

이러한 제한조건에서 아래의 최대 엔트로피 조건을 만족시켜야 한다.

위의 조건를 만족시키는 f(s)는 lagrange-multiplier를 이용하고, 변분방법을 쓰면 쉽게 구해진다.

그리고 λ는 

의 근이다. 오르쪽은 λ에 대해서 단조증가함수이므로 근은 유일하게 하나가 존재한다.
원본이미지에서 픽셀의 평균값을 계산하고, 위의 방정식에 넣어서 λ를 구하면, 목표로하는 히스토그램을 얻을 수 있으면, 이것의 cumulative 히스토그램과, 원본이미지의 cumulative히스토그램을 구해서, 그 차이를 최소로하는 픽셀값의 대응관계를 찾으면 된다(histogram specification).

참고: Bright Preserving Histogram Equalization with Maximum Entropy: A Variational Perspective.
        C. Wang and Z.Ye, IEEE Trans. Consumer Electronics. V51. No4. (2005);

// BPHEME의 cumulative histogram(continuous version)::integral of f(s) over s;

double cdf(double s, double mu, double lambda) {

// histogram specification;1==>2
void hist_spec(double chist1[],  //source cumulative histogram;
                      double chist2[], // target cumulative histogram;
                      int n,                 //=256;
                      int lut[])             // resultant mapping(1->2); {

//

void BPHEME(BYTE* src, int width, int height, BYTE *dst) {
    double hist[256], hist2[256];           //src(dst) histogram;

    double chist[256], chist2[256];        //src(dst) cumulative histogram;
    int lut[256];                                //histogram specification mapping;

    int n = width * height;
    make_hist(src, width, height, hist);

    normalize_hist(hist, 256);

    //cumulative histogram;

    make_cumulative_hist0(hist, 256, chist);
    TRACE("chist[255]=%f\n", chist[255]);
    // gray-mean;

    double mu = hist_mean(hist, 256);
    TRACE("mean=%f(%d)\n", mu, int(mu * 255+.5));
    // determin lambda;
    double lambda, th=1.e-15;
    FindRoot(mu, th, lambda);
    TRACE("lambda=%f\n", lambda);
    // entropy of src;

    double entropy1 = hist_entropy(hist, 256)

     // dst-cumulative

    for (i = 0; i < 256; i++) chist2[i] = cdf(double(i)/255., mu, lambda);
    TRACE("chist2[255]=%f\n", chist2[255]);

    // histogram-specification;
    hist_spec(&chist[0], &chist2[0], 256, &lut[0]);

    // make dist image;
    for (i=0; i<n; i++) dst[i] = lut[src[i]];

    // make histogram of dst;

    make_hist(dst, width, height, hist2);

    // normalize hist2;

    normalize_hist(hist2, 256);

    // mean of dst;

    double mu2 = hist_mean(hist2, 256);
    TRACE("new mean=%f(%d)\n", mu2, int(mu2 + .5));

    // entropy of dst;

    double entropy2 = hist_entropy(hist2, 256);
    TRACE("total entropy(before, after)= (%f, %f)\n", entropy1, entropy2);
};


// Root finding procedure ::
// define F(s, mu);
double F(double s, double mu) {

//derivative of F(s, mu) w.r.t. s;
double DF(double s, double mu) {   

// Find root of F(s,mu) based on Newton-Raphson method.
int FindRoot(double mu, double threshold, double& s) {

사용자 삽입 이미지

histogram equaliztion 결과;

사용자 삽입 이미지


BPHEME  결과:

사용자 삽입 이미지

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

Fast Median Filter 3x3  (0) 2009.11.27
Fant's Resampling  (0) 2008.12.17
Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (1) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Posted by helloktk

이미지를 이진화시키기 위해서 여러가지 알고리즘이 사용된다. 이미지 전체에 대해서 하나의 임계값을 이용하여서 이진화 시키는 알고리즘은 간단하고 빠르기 때문에 많이 이용이 된다. 그러나 이미지를 형성할 때 조명조건이 균일하지 않은 경우에는 이 방법으로는 원하는 결과를 얻을 수 없다. 이런 경우에는 각각의 픽셀 주위의 그레이값을 참조하여서 임계치를 결정하여야 한다. 간단한 방법은 윈도우내의 그레이값의 평균을 참조하는 것이다. 이보다 좀더 개선된 알고리즘은 평균값을 참조하되, 그레이값의 편차를 한번 더 고려해 주는 것이다. 이렇게 하여서 잡은 국소적인 임계치는 

128 은 그레이값이 갖을 수 있는 최대편차를 의미한다. 이 경우에는 단순 평균값으로 취한다는 의미이다. 그 외의 경우는 표준편차와 128의 차이(항상 음수다)에 비례하는 값으로 윈도우 평균값을 오프셋한 값을 임계치로 잡는다. factor는 일반적으로 정해지지 않고, 실험적으로 [0.2, 0.5]사이의 값이 취해진다. (문서처럼 배경이 흰색인 경우는 factor>0이지만, 만약에 반대의 경우, 즉 검정 배경에 흰색글씨를 처리하는 경우는 음수의 값을 취하는 것이 맞다)
 
국소적인 이진화 알고리즘은  매 픽셀마다 윈도우를 잡아서 계산해야하는 과정으로 인해서 계산비용이 많이 든다. 만약에 충분한 메모리를 갖춘 시스템의 경우에는 적분이미지(integral image)를 이용하면 윈도우 계산에 소요되는 비용을 줄일 수는 있다.

윈도우 크기와 factor는 결정하는 기준은 무었일까? 이것은 해결하고자 하는 문제의 특성에, 예를 들면 스캔된 문서를 이진화시키는 경우에는 윈도우에 충분한 글자가 들어 있어야 한다..등등, 많이 의존한다.

#define integral_image(x, y) (intimage[(y)*width+(x)])
#define integral_sqimg(x, y) (intsqimg[(y)*width+(x)])

void make_int_img12(BYTE *gray, int width, int height, *int intimage, int *intsqimg) {

//
void adap_binariztion(BYTE *gray, int width, int height,
          int w       /*window size = 15*/,
          double k  /*factor           = 0.2*/,
          BYTE *bimage)
{
    int whalf=w>>1; //half of adaptive window;
    int diff, sqdiff;
    // make integral image && square integral image;
    // if image is sufficiently large, use int64 or floating point number;
    std::vector<int> intimage(width * height) ;
    std::vector<int> intsqimg(width * height) ;

    //make integral image and its square integral image;
    make_int_img12(gray, width, height, &intimage[0], &intsqimg[0]);  
    //algorithm main;
    for(int j=0; j<height; j++){
        for(int i=0; i<width; i++){
            //clip windows
            int xmin = max(0,i-whalf);
            int ymin = max(0,j-whalf);
            int xmax = min(width-1,i+whalf);
            int ymax = min(height-1,j+whalf);
            int area = (xmax-xmin+1)*(ymax-ymin+1);
            // calculate window mean and std deviation;
            if(!xmin && !ymin){     // origin
                diff   = integral_image(xmax,ymax);
                sqdiff = integral_sqimg(xmax,ymax);
            }
            else if(!xmin && ymin){ // first column
                diff   = integral_image(xmax,ymax) - integral_image(xmax,ymin-1);
                sqdiff = integral_sqimg(xmax,ymax) - integral_sqimg(xmax,ymin-1);
            }
            else if(xmin && !ymin){ // first row
                diff   = integral_image(xmax,ymax) - integral_image(xmin-1,ymax);
                sqdiff = integral_sqimg(xmax,ymax) - integral_sqimg(xmin-1,ymax);
            }
            else{ // rest of the image
                int diagsum    = integral_image(xmax,ymax) + integral_image(xmin-1,ymin-1);
                int idiagsum   = integral_image(xmax,ymin-1) + integral_image(xmin-1,ymax);
                diff           = diagsum - idiagsum;
                int sqdiagsum  = integral_sqimg(xmax,ymax) + integral_sqimg(xmin-1,ymin-1);
                int sqidiagsum = integral_sqimg(xmax,ymin-1) + integral_sqimg(xmin-1,ymax);
                sqdiff         = sqdiagsum - sqidiagsum;
            }
            // threshold = window_mean *( 1 + factor * (std_dev/128.-1));
            // 128 = max_allowed_std_deviation in the gray image;
            double mean = double(diff)/double(area);
            double std  = sqrt((sqdiff - double(diff)*diff/double(area))/double(area-1));
            double threshold = mean*(1.0 + k*((std/128.0) - 1.));
            if(gray[j * width + i] < threshold)
                bimage[j * width + i] = 0;
            else
                bimage[j * width + i] = 255;
        }
    }  
}

사용자 삽입 이미지

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

Fant's Resampling  (0) 2008.12.17
Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (1) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Otsu Algorithm  (6) 2008.05.30
Posted by helloktk
영상에는 보통 우리가 관심이 있는 물체와 배경을 동시에 담고 있다. 그런데 어떤 경우에는 전체영상의 픽셀값의 다이나믹 레인지가 매우 좁아서, 달리 이야기하면 영상의 히스토그램을 살펴볼 때, 픽셀값의 분포가 매우 좁은 영역에 한정이 되어있는 경우에 물체와 배경이 잘 구별이 안되는 영상을 접하게 된다. 이러한 영상의 경우에 몇가지 영상처리방법을 동원하여서 좋은 콘트라스트를 갖는 영상을 만들 수 있다.

Histogram Equalization은 이러한 기법중의 하나이다. 이것은 좁은 영역에 분포하는 픽셀의 그레이 값을 가능한 모든 그레이 값의 범위에 골고루 분포하도록 재배치하여서 영상의 constrast 를 증대시키는 방법이다. Histogram Equalization에서 사용하는 그레이 값의 mapping은 출력영상의 픽셀 값은 모든 그레이레벨의 값을 다 가져야 한다는 것과, 각각의 그레이 레벨에 해당하는 픽셀의 수가 거의 일정하여야 한다(픽셀 값이 이산적이기 때문에 완벽하게 상수가 될 수 없다)는 조건에서 결정이 된다.



입력영상에서 그레이 값이 k 인 픽셀은 출력영상에서 어떤 그레이 값을 가져야 할까? 조건에 따르면 입력영상의 k 그레이레벨이 출력영상의 j 그레이 레벨로 매핑이 된다면, 출력영상에서는 부터 까지의 그레이레벨에 해당하는 픽셀의 수는 모든 그레이 레벨에서 픽셀분포가 균일하다는 조건 때문에 

로 주어진다. 이 값과 입력영상에서 부터 까지의 그레이 레벨을 갖는 픽셀의 수가 같으면 된다.


부터 k 까지의 그레이 값을 갖는 입력영상의 픽셀의 수는 


또, 총 픽셀 수는 cumulative histogram[255]로 주어진다. 따라서,  원하는 매핑 j = T[k]



이 매핑에 의하면 입력영상의 픽셀이 가지는 가장 큰 그레이 값은 출력영상에는 255에 해당하게 된다. 출력영상은 정의에 의해서 cumulative histogram이 직선적으로 증가하는 형태를 가지게 된다. 

큰 영상의 경우에 이 매핑을 Lookup 테이블로 저장하여서 연산의 횟수를 줄일 수 있다.


note)
1. 0차 cumulative histogram을 이용하면 0 부터 k 까지의 histogram의 합은 한 번의 호출로 얻어진다.
2. 0차 cumulative histogram의 마지막 bin의 값은 전체 픽셀수의 합이 된다.


// histogram의 합은 cumulative histogram을 구성하면 된다;

void HistogramEqualize(BYTE *src, int width, int height, BYTE* dst)

{

    int histogram[256], cumhist[256];

    memset(histogram, 0, sizeof(histogram)) ;

    for (int y = 0, k = 0; y < height; y++){

        for (int x = 0; x < width; x++, k++){

            histogram[src[k]]++ ;

       }

    }

    // make 0-st order cumulative histogram;

    cumhist[0] = histogram[0];

    for (int i = 1; i < 256; i++){

        cumhist[i] = cumhist[i-1] + histogram[i] ;

    };


    // make output;

    // cumhist[255] = number of pixels;

    double v = 255.0/cumhist[255];

    for (y = 0, k = 0; y < height; y++){

        for (int x = 0; x < width; x++, k++){
            int  a = int(v * cumhist[src[k]]);

            dst[k] = a < 0 ? 0 : a > 255 ? 255 : a;  //clamp to [0,255] 

        }

    }

};


사용자 삽입 이미지


평탄화 전의 히스토그램(붉은색) 밎 누적히스토그램(연두색)


사용자 삽입 이미지



평탄화 후의 히스토그램(붉은색) 밎 누적히스토그램(연두색)

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

Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (1) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Otsu Algorithm  (6) 2008.05.30
Hough Transform  (2) 2008.05.22
Posted by helloktk

사용자 삽입 이미지

사용자 삽입 이미지

/* 1차원 FFT:: dir=1 for forward, -1 for backward;
** x = real-component, and y = imaginary componet.
** nn = length of data: 2의 제곱승이어야 한다.
*/
int fft1d(int dir, int nn, double *x, double *y) {
    int  i,i1,j,k,i2,l,l1,l2;
    double c1,c2,tx,ty,t1,t2,u1,u2,z;
    /* Calculate the number of points */
    int m = 0;
    while(nn > 1) {
        nn >>= 1;  m++;
    }
    nn = 1 << m;

    /* Do the bit reversal */
    i2 = nn >> 1;
    j = 0;
    for (i=0;i<nn-1;i++) {
        if (i < j) {
            tx = x[i];    ty = y[i];
            x[i] = x[j]; y[i] = y[j];
            x[j] = tx;    y[j] = ty;
        }
        k = i2;
        while (k <= j) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
   
    /* Compute the FFT */
    c1 = -1.0;
    c2 = 0.0;
    l2 = 1;
    for (l=0;l<m;l++) {
        l1 = l2;
        l2 <<= 1;
        u1 = 1.0;
        u2 = 0.0;
        for (j=0;j<l1;j++) {
            for (i=j;i<nn;i+=l2) {
                i1 = i + l1;
                t1 = u1 * x[i1] - u2 * y[i1];
                t2 = u1 * y[i1] + u2 * x[i1];
                x[i1] = x[i] - t1;
                y[i1] = y[i] - t2;
                x[i] += t1;
                y[i] += t2;
            }
            z =  u1 * c1 - u2 * c2;
            u2 = u1 * c2 + u2 * c1;
            u1 = z;
        }
        c2 = sqrt((1.0 - c1) / 2.0);
        if (dir == 1)
            c2 = -c2;
        c1 = sqrt((1.0 + c1) / 2.0);
    }
   
    /* Scaling for forward transform */
    if (dir == 1) {
        for (i=0;i<nn;i++) {
            x[i] /= (double)nn;
            y[i] /= (double)nn;
        }
    }
   
    return(1);
}
/* 2차원 FFT. data= (2*nn) x mm 크기의 버퍼. 각각의 행은 nn개의 실수성분다음에 nn개의 허수성분이 온다.
** nn, mm은 은 각각 행과 열은 나타내며, 2의 제공승의 숫자이어야 한다.
** isign=1 for forward, -1 for backward.
** copy = buffer of length 2*mm ;
*/
void fft2d(double* data, int nn, int mm, int isign, double* copy)
{
    int i, j, index1, index2;
   
    /* Transform by ROWS for forward transform */
    if (isign == 1)
    {
        index1 = 0;
        for (i = 0; i < mm; i++)
        {
            fft1d(isign, nn, &data[index1], &data[index1+nn]);
            index1 += (nn << 1);
        }
    }
   
    /* Transform by COLUMNS */
    for (j = 0; j < nn; j++)
    {
        /* Copy pixels into temp array */
        index1 = j ;
        index2 = 0;
        for (i = 0; i < mm; i++)
        {
            copy[index2] = data[index1];
            copy[index2+mm] = data[index1 + nn];
            index2++;
            index1 += (nn << 1);
        }
       
        /* Perform transform */
        fft1d(isign, mm, &copy[0], &copy[mm]);
       
        /* Copy pixels back into data array */
        index1 = j;
        index2 = 0;
        for (i = 0; i < mm; i++)
        {
            data[index1] = copy[index2];
            data[index1 + mm] = copy[index2+mm];
            index2++;
            index1 += (nn << 1);
        }
    }
   
    /* Transform by ROWS for inverse transform */
    if (isign == -1)
    {
        index1 = 0;
        for (i = 0; i < mm; i++)
        {
            fft1d(isign, nn, &data[index1], &data[index1+nn]);
            index1 += (nn << 1);
        }
    }
}
void butterworth( double * data, int Xdim, int Ydim,
                 int Homomorph, int LowPass,
                 double Power, double CutOff, double Boost)
{
    int x, y, x1, y1, halfx, halfy;
    double Filter;
    double CutOff2 = CutOff * CutOff;
   
    /* Prepare to filter */
    halfx = Xdim / 2;
    halfy = Ydim / 2;
    double *rdata = &data[0];
    double *idata = &data[Xdim];
    /* Loop over Y axis */
    for (y = 0; y < Ydim; y++) {
        if (y < halfy)  y1 = y;
        else             y1 = y - Ydim;
       
        /* Loop over X axis */
        for (x = 0; x < Xdim; x++) {
            if (x < halfx)  x1 = x;
            else             x1 = x - Xdim;
           
            /* Calculate value of Butterworth filter */
            if (LowPass)
                Filter = (1 / (1 + pow((x1 * x1 + y1 * y1) / CutOff2, Power)));
            else if ((x1 != 0) || (y1 != 0))
                Filter = (1 / (1 + pow( CutOff2 / (x1 * x1 + y1 * y1), Power)));
            else
                Filter = 0.0;
            if (Homomorph)
                Filter = Boost + (1 - Boost) * Filter;
           
            /* Do pointwise multiplication */
            rdata[x] *= Filter;
            idata[x] *= Filter;
        };
        rdata += (Xdim << 1) ; //go to next-line;
        idata += (Xdim << 1) ;
    }
}
Butterworth Low Pass filter 적용의 예(order=2, cutoff=20)

사용자 삽입 이미지

PowerSpectrum 변화
사용자 삽입 이미지

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

Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (1) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Otsu Algorithm  (6) 2008.05.30
Hough Transform  (2) 2008.05.22
Posted by helloktk
TAG FFT, fft2d

이미지에서 뭔가 유용한 정보를 얻어내기 위해서는 이미지가 담고 있는 객체들을 분리하는 작업을 해야한다. 이 작업중 가장 단순한 것 중의 하나가 바로 이진화(binarization)이다. 이진화는 이미지의 그레이값에 따라서 픽셀을 0 과 1(또는 255)로 바꾸는 작업이다. 이러한 이진화 작업을 거치면 이미지가 담고 있는 객체를 배경에서 분리할 수 있다. 그런데 이때 어떤 그레이값으로 픽셀을 분리해야하는가라는 문제가 생긴다. 임계값(threshold value)의 설정에 대한 여러가지 알고리즘이 알려져 있는데, 그 중에서 통계적인 방법을 이용한 Otsu의 알고리즘이 가장 자연스러운 임계값 설정의 기준을 준다.

Otsu 알고리즘은 classification 기법을 이용하고 있다. 임계값을 설정하는데 있어서 비용함수를 설정하고 그 비용함수의 최소값을 주는 값으로 임계값을 취하는 방식이다. 그럼 어떻게 비용함수를 설정할 것인가?  이미지상에서 나타나는 그레이값을 2개의 클래스로 분리할 때, 좋은 분리는 각각의 클래스에 속한 픽셀의 그레이값의 분포가 유사해야 한다. 즉, 같은 클래스내에  들어있는 픽셀의 그레이값의 분산이 작아야 한다는 의미다. 따라서 비용함수는 가중치를 갖는 클래스내의 분산의 합이 될것이고, 임계값은 이 비용함수을 최소화하는 그레이값이다.

           비용함수 = (가중치1 * 분산1) + (가중치2 * 분산2)  <= 2개 클래스로 분리시
                        =   q1 * V1 + q2 * V2 ;      

                     q1 =  전체 이미지에서 클래스1에 해당하는 픽셀이 나타날 확률
                     q2 =  전체 이미지에서 클래스2에 해당하는 픽셀이 나타날 확률
                     V1 = 클래스1 내의 픽셀들의 그레이값의 분산.
                     V2 = 클래스2 내의 픽셀들의 그레이값의 분산.

             임계값  -->  MIN ( 비용함수 )

이미지상이 그레이값의 분포는 히스토그램으로 표현이 되므로, 임계값은 히스토그램으로 분리하는 레벨값이고, 클래스1은 그 값보다도 작은 부분, 클래스2는 큰 부분을 의미한다. 비용함수의 의미를 좀더 살펴보기 위해서 비용함수에 간단한 계산을 하면

비용함수 = 전체의 분산 - (가중치1*(전체평균-평균1)^2 + 가중치2*(전체평균-평균2)^2);
             = V - ( q1 * (m1 - m)^2  + q2 * (m2 - m)^2) ;
                         
                          V = 전체분산;
                          m = 전체 평균,
                          평균1 (m1) = 클래스1의 평균,
                          평균2 (m2) = 클래스2의 평균,

그런데 q1*(m-m1)^2 + q2*(m-m2)^2는 클래스들의 분산이다. 전체분산은 어떤식으로 클래스를 분리하더라도 항상 일정한 값을 가지므로, 비용함수를 최소화하는 것은 클래스들의 분산을 최대화하는 것과 같다. 그래서 새로운 비용함수(엄밀한 의미로 이득함수이다)를 이것으로 잡으면
 
             이득함수 = q1 * (m1 - m)^2 + q2 * (m2 - m)^2;
             임계값 --> MAX ( 이득함수)
 
새로운 이득함수는 약간의 계산을 하면 그 의미가 더 명확한 표현인
             
             이득함수 = q1 * q2 * (m1 - m2)^2 ;

로 쓸 수 있다. 즉, 클래스 분리하는 값은 두 클래스의 평균값의 차이(가중치를 갖는)를 최대화시키는 값으로 주어진다.

이 알고리즘의 구현은 히스토그램의 각각의 레벨에 대해세 좌 우를 각각 클래스 1과 2로 설정하고, 이들 클래스에 대해서 이득함수의 최대값을 업데이트하는 방식을 취하면 된다. 0차 및 1차 cumulative 히스토그램을 이용하면 확률과 평균값은 바로 구할 수 있다.

Otsu 알고리즘은 2개의 클래스 뿐만 아니라 여러 개의 클래스로 히스토그램을 분리하도록 확장할 수 있다. 그리고 구현은 재귀알고리즘으로 사용하면 쉽다. 또한, 0번째 cumulative histogram과 1번째 cumulative histogram을 사용하면, 각 클래스의 가중치와 평균값을 쉽게 계산할 수 있다. 
        cumulative_histogram_0th[k] = 0...k까지 값이 나타날 확률.
        cumulative_histogram_1th[k]/
cumulative_histogram_0th[k] = 0...k까지 값의 평균.         


/*  Otsu임계화 예: 설명을 위해 최적화는 하지 않은 것이다. 반드시 cumulative 히스토그램을 이용할
** 필요는 없다: 이득함수의 최대값을 주는 지점이 여러개일 때는 평균값을 취하도록 한다(2016.04.26)
*/
int OtsuThreshold(BYTE *src, int width, int height, BYTE *dst) {
    double hist[256], chist[256], cxhist[256];
    memset(hist, 0, sizeof(hist));
    memset(chist, 0, sizeof(chist));
    memset(cxhist, 0, sizeof(cxhist));
    int ntot = width * height;

    // make histogram ;
    for (int i=0; i<ntot; i++) hist[src[i]] += 1;
   
    // normalize;
    for (i = 0; i < 256; i++) hist[i] /= ntot;

    // make 0- and 1-st cumulative histogram;
    chist[0] = hist[0];
    cxhist[0] = 0;
    for (i = 1; i < 256; i++) {
        chist[i] = chist[i-1] + hist[i] ;                    //0-th cumulative histogram ;
        cxhist[i] = cxhist[i - 1] + double(i) * hist[i] ; //1-th cumulative histogram ;
    };
   
    double gain_max = 0;
    int thresh = 0;   
    double m = cxhist[255];                     //total mean ;
    int mul_count = 1;                            //number of degenerate maxima;
    for (i = 0; i < 256; i++) {
        if (chist[i] == 0.) continue ;
        double q1 = chist[i] ;                  //weight1;
        double q2 = 1 - q1;                     //weight2;
        if (q2 == 0.) break;
        double m1 = cxhist[i] / q1;             //mean1 ;
        double m2 = (m - cxhist[i]) / q2;       //mean2 ;
        double gain = q1 * q2 * (m1 - m2) * (m1 - m2) ;
        if (gain_max < gain) {
            gain_max = gain;
            thresh   = i;
            mul_count = 1;                         //reset mul_count=1;
        } else if (gain_max == gain)             //degenerate case;
            mul_count ++;
    }
    if (mul_count > 1) thresh = thresh + (mul_count - 1) / 2;    //2016.04.26;

    // threshold image;
    for (i = 0; i < ntot; i++) dst[i] = (src[i] >= thresh) ? 0xFF : 0x00 ;

    return thresh;
}

사용자 삽입 이미지

히스토그램 (계산된 임계치는 100이다)

사용자 삽입 이미지

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

Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (1) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Otsu Algorithm  (6) 2008.05.30
Hough Transform  (2) 2008.05.22
Posted by helloktk

주어진 영상에서 선분을 찾고자 하면 어떻게 하여야 할까. 우선 에지영상을 만들어 선분을 강조하고, 선분이 하나만 있는 경우에는 에지데이터를 이용하여 line-fitting을 이용하여 해당선분을 기술하는 파라미터를 추출할 수 있다. 그러나 영상에 선분이 여러개 있는 경우에는 이 방법을 이용하기 어렵다(불가능한 것은 아니지만 노이즈등에 의한 영향등을 고려해야하기 때문에 복잡한 과정을 거쳐야 한다). 평면상에서 선분은 원점에서의 거리와 그것의 수직(수평)인 방향만 주어지면 하나의 직선이 결정이 된다. 원점에서의 직선까지의 거리를 r, 법선이 x-축과 이루는 각도를 θ 라고 하면

                          r = cosθ * x + sinθ * y;

의 형태로 주어진다. 즉, (r,θ) 한쌍이 주어지면 (x,y)평면에서 직선이 하나가 정의된다 (하나의 (r, θ) 값이 주어지는 경우, 만들어지는 직선을 a * x + b * y = c  꼴로 쓰면, 계수들은 a = cosθ, b = sinθ, c = r 로 주어진다).

이 직선의 방정식은 rθ-평면의 stripe을  xy-평면으로 보내는 변환으로도 생각할 수있다. 이 변환은 rθ-평면의 한점을 xy-평면의 직선으로 보낸다. 그럼 역으로 xy-평면의 한점은 rθ-평면의 곡선으로 변환이 된다. xy-평면상의 각각의 점들은 rθ-평면의 곡선을 하나씩 만드는데 이 곡선들이 공통으로 만나는 교점이 생기게 된다. rθ-평면에서의 이 교점의 좌표가 바로 xy-평면에서의 직선을 기술한다.

Hough Transform은 이 변환의 특성을 이용한 것이다. 이미지에서 에지에 해당하는 점들을 위의 변환에 의해서 rθ-평면의 곡선으로 보내는 것이다 (프로그램적으로는, 에지점 좌표 (x,y)가 주어지면 θ=0, π까지 일정 간격으로 증가시키면서 곡선 r=x*cos(θ) + y*sin(θ)의 값을 구해서 메모리상의 [r,θ]지점의 돗수를 1씩 증가시킨다). 만약 에지에 해당되는 점들이 직선의 관계를 가지게 되면, 각 점들에 해당하는 rθ-평면에서 곡선들은 한 교점을 형성하게 될 것이다. 긴 선분은 rθ-평면의 한점에서 더 많은 곡선들의 교점을 형성하게 된다. 따라서, rθ-평면에서 이러한 교점의 히스토그램을 구하여, 교점수가 일정 이상 누적이 된 경우를 취하면, xy- 이미지 평면에서 일정한 길이 이상을 갖는 선분을 골라 낼 수 있다.

그러면 왜 rθ-평면으로 변환을 생각하는 것일까? 직선이  y= a*x+b의 형태로 주어 질수 있으므로 기울기와 y축과의 교점을 기술하는 ab-공간을 이용할 수 도 있으나, 이 경우에 y-축에 평행인 직선의 경우에 기울기 a 값이 무한히 커지는 경우가 발생하여 저장공간 할당에 문제가 생긴다. 실제적인 문제에서 되도록이면 compact한 공간으로 변환을 고려해야 하는 것이다. rθ-공간으로의 변환은 유한한 이미지의 경우 항상 유한한 rθ-공간의 영역으로 변환된다.

아래의 그림은 x-y=-5인 직선상의 세 점 (5,10), (10,15), (15,20) 에 해당하는 rθ-평면상의 세 곡선을 보인것이다: r = 5*cosθ + 10*sinθ, r = 10*cosθ + 15 * sinθ, r = 15*cosθ + 20*sinθ. 이 세 곡선은 (r,θ) = (5/sqrt(2), 3*Pi/4) 지점에서 만남을 알 수 있다. 이 만나는 지점을 구하면, 거꾸로 세 점이  x-y=-5 인 직선 위의 점들이었다는 것을 추정할 수 있다.

사용자 삽입 이미지

BOOL HoughTransformLine(BYTE* image, int width, int height, /*background=0*/
                                       double rho/*=2*/,
                                       double theta/*=Pi/180.*/,
                                       int threshold/*=20*/,
                                       std::vector<Line>& vecLine) {
    /* use cosine and sine look-up tables*/
    double irho = 1 / rho;
    int numangle = (int) (Pi / theta);
    int numrho = (int) (((width + height) * 2 + 1) / rho);
    int *accum = (int*)malloc( sizeof(accum[0]) * (numangle+2) * (numrho+2) );
    memset( accum, 0, sizeof(accum[0]) * (numangle+2) * (numrho+2) );
    double ang;
    int n, rwidth=(numrho+2);
    for(int  y=0; y<height; y++ ){
        for(int x=0; x<width; x++ ){
            if( image[y * width+ x] != 0 ){//only for foreground pixels;
                for(ang = 0, n = 0; n < numangle; ang+= theta, n++ ){
                    double rho=(x * cos(ang) + y * sin(ang)) * irho ;
                    int r=rho>=0? int(rho+.5) : (-int(-rho+.5));   //round to nearest integer;
                    //accum을 이차원배열로 생각하고, 1픽셀의 border를 두는 형태로 한다.
                    //이는 아래의 local-maxima를 찾을때 경계를 고려할 필요가 없어서 유리하다.
                    r += (numrho - 1) / 2;
                    accum[(n+1) * rwidth + r+1]++;
                }
            }
        }
    }
    //find local maxima;
    for(int  r = 0; r < numrho; r++ ){
        for(int  n = 0; n < numangle; n++ ){
            int base = (n+1) * rwidth + r+1;
            //test if it is local-maximum(4-way);
            if( accum[base] > threshold &&
                accum[base] > accum[base - 1] && accum[base] > accum[base + 1] &&
                accum[base] > accum[base - rwidth] && accum[base] > accum[base + rwidth] )
            {              
                Line line;
                line.rho = (r - (numrho - 1) *.5) * rho;
                line.angle = n * theta;
                vecLine.push_back( line );
            }
        }
    }
    free(accum);
    return TRUE;
}

아래 그림에서 첫번째는 소스이미지이고, 이 이미지에서 rθ-평면으로 HoughTransform한 것이 두번째 것이다(rho=2, theta=Pi/360). 원본이미지에는 8개의 선분이 있고, 4개씩 같은 방향을 갖고 있다. HoughTransform된 이미지에서 이러한 특징을 볼수 있다(가로 축이 θ이고 세로축이 r이다. r축은 가운데가 r=0 이다). 누적이 많이 된 부분이 두 군데의 동일한 θ에서 나타나고 있다. 그러나 결과에서 8개의 피크를 분리하기는 쉽지 않다. 위의 코드는 각각의 점에서 4방향으로 체크하여 극대값을 찾고 있으나, 항상 잘 동작하는 것은 아니다.

local -maxima를 좀 더 잘 잡고 싶으면, 위처럼 주변의 4점만 체크할 것이 아니라, 윈도우를 좀더 크게 잡고, 그 윈도우내에서 최대값을 찾아 보는 것도 한 방법이 된다.

사용자 삽입 이미지

사용자 삽입 이미지


/**
** http://blog.naver.com/helloktk/80051779331
*/

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

Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (1) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Otsu Algorithm  (6) 2008.05.30
Hough Transform  (2) 2008.05.22
Posted by helloktk