이미지에서 뭔가 유용한 정보를 얻어내기 위해서는 이미지가 담고 있는 객체들을 분리하는 작업을 해야한다. 이 작업중 가장 단순한 것 중의 하나가 바로 이진화(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

댓글을 달아 주세요

  1. 감사^^ 2008.09.10 01:21  댓글주소  수정/삭제  댓글쓰기

    세미나 준비하는 중에 막막했던 부분이었는데

    너무 잘 봤습니다.

    근데 한가지 궁금한점이 있어서요...

    두 클래스의 분산의 합에 대한 수식이 두 클래스의 평균의 차이로

    변하는 과정을 알고 싶은데요...혹시 알수 있을까요?

  2. helloktk 2008.09.11 06:22  댓글주소  수정/삭제  댓글쓰기

    전체평균=m=q1*m1 + q2*m2(정의로부터); q1 + q2=1;을 이용하면,
    q1*(m1-m)^2 + q2*(m2-m)^2;
    =q1*(m1-q1*m1-q2*m2) + q2*(m2-q1*m1-q2*m2)^2;
    =q1*((1-q1)*m1-q2*m2)^2 + q2*((1-q2)*m2-q1*m1)^2 ;
    =q1*(q2*m1-q2*m2)^2 + q2*(q1*m2-q1*m1)^2;
    =q1*q2^2*(m1-m2)^2 + q2*q1^2(m1-m2)^2;
    =q1*q2*(q1+q2)*(m1-m2)^2;
    =q1*q2*(m1-m2)^2;
    을 얻을 수 있습니다.

  3. 공대생 2010.05.29 22:58  댓글주소  수정/삭제  댓글쓰기

    혹시 보게 되시면요...

    저도 bmp 이미지를 받아

    히스토그램은..

    void DibHistogram(CDib& dib, float histo[256])
    {
    register int i,j;

    int w = dib.GetWidth();
    int h = dib.GetHeight();

    BYTE** ptr = dib.GetPtr();
    //Histogram 계산
    int temp[256];
    memset(temp, 0, sizeof(int)*256);
    for(j=0; j<h; j++)
    for(i=0; i<w; i++)
    {
    temp[ptr[j][i]]++;
    }
    //Histogram 정규화
    float area = (float)w*h;
    for(i=0; i<256; i++)
    histo[i] = temp[i]/area;
    }


    이진화는 현재,,,

    void DibBinarization(CDib& dib)
    {
    register int i, j;

    int w = dib.GetWidth();
    int h = dib.GetHeight();

    BYTE** ptr = dib.GetPtr();

    for(j=0; j<h; j++)
    for(i=0; i<w; i++)
    {
    ptr[j][i]=(ptr[j][i]>50) ? 255:0;
    }

    }


    이런식의 이진화 중이에요.

    근데 위에거 보고 따라하기 식으로 해도 안되더라구요...

    아직 초보자라서요

    소스 좀 가르쳐 주시면 감사하겠습니다(제꺼에 적용할수 있겠끔)

    부탁드립니다.

    E-mail : rockentop@cyworld.com

  4. helloktk 2010.06.02 11:48 신고  댓글주소  수정/삭제  댓글쓰기

    binirization함수시작부분에 아래 부분을 추가하고, 50 대신에 thresh를 기준으로 하면 됩니다.

    double chist[256], cxhist[256];
    memset(chist, 0, sizeof(chist));
    memset(cxhist, 0, sizeof(cxhist));

    chist[0] = histo[0];
    cxhist[0] = 0;
    for(i=1; i<256; i++){
    chist[i] = chist[i-1] + histo[i] ;
    cxhist[i] = cxhist[i-1] + double(i) * histo[i] ;
    };

    double cost_max = 0;
    int thresh = 0;
    double m = cxhist[255]; //total mean ;
    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 cost = q1*q2*(m1-m2)*(m1-m2) ;
    if(cost_max < cost) {
    cost_max = cost;
    thresh = i;
    };
    }

  5. 학생 2010.11.08 13:32  댓글주소  수정/삭제  댓글쓰기

    안녕하세요. 님의 귀중한 정보 덕분에 어느정도의 해결을 했습니다.
    저도 위에 공대생님처럼 bmp파일을 읽어서 이진화를 하려고하는데
    저는 bmp파일의 색상값이 아니라 다른 값으로 이진화를 하려고했지만
    히스토그램을 그려본 결과 봉우리가 3개가 나타났습니다. 원하는 이진화값은
    1,2번 봉우리 사이의 최소값이 나와야하는데 2,3번 봉우리 사이의 최소값이
    앞의 최소값보다 더 작아서 원하는 이진화가 안되고 있습니다.
    대략적으로 보면 앞의 최소값은 약 80정도이고 뒤의 최소값은 187인데
    뒤의 최소값을 없애던지 아니면 두개의 임계값을 찾을 방법은 없는것인지요..

    • helloktk 2010.11.11 09:53 신고  댓글주소  수정/삭제

      정확히 뭘 이야기하는지 잘 이해가 안되지만, 아마 님의 히스토그램에서 봉우리가 3개가 나타나고, 원하는 값은 1-2봉사이의 최소값 지점인 것 같은데. 이것이 2-3봉사이의 최소값보다도 커서 원하는 값이 안 되는(??) 상황같은데..OTSU방법은 이 최소값이 중요한 것이 아니라 전체적인 봉우리들의 모양이 중요합니다. 그리고 봉우리가 여럿있을 떄 반드시 봉우리 사이의 계곡에서 최소인지점에서 나타난다고 할 수 는 없습니다. 따라서 오수 알고리즘을 쓰면, 이 알고리즘이 주는 값에서 정해야 하고, 결과가 원하는 것이 아니면, 이 알고리즘을 쓰면 안되고 다른 알고리즘으로 구현하는 것이 더 적합할 것입니다. 이 알고리즘이 반드시 가장 좋은 결과를 주지는 않습니다.