Otsu 알고리즘은 이미지를 이진화시키는데 기준이 되는 값을 통계적인 방법을 이용해서 결정한다. 같은 클래스(전경/배경)에 속한 픽셀의 그레이 값은 유사한 값을 가져야 하므로 클래스 내에서 픽셀 값의 분산은 되도록이면 작게 나오도록 threshold 값이 결정되어야 한다. 또 잘 분리가 되었다는 것은 클래스 간의 거리가 되도록이면 멀리 떨어져 있다는 의미이므로 클래스 사이의 분산 값은 커야 함을 의미한다. 이 두 가지 요구조건은 동일한 결과를 줌을 수학적으로 보일 수 있다.

이미지의 이진화는 전경과 배경을 분리하는 작업이므로 클래스의 개수가 2개, 즉, threshold 값이 1개만 필요하다. 그러나 일반적으로 주어진 이미지의 픽셀 값을 임의의 개수의 클래스로 분리할 수도 있다. 아래의 코드는 주어진 이미지의 histogram을 Otsu의 아이디어를 이용해서 nclass개의 클래스로 분리하는 알고리즘을 재귀적으로 구현한 것이다. 영상에서 얻은 히스토그램을 사용하여 도수를 계산할 수 있는 0차 cumulative histogram($\tt ch$)과 평균을 계산할 수 있는 1차 culmuative histogram($\tt cxh$)을 입력으로 사용한다. 

$$ thresholds= \text {argmax} \left( \sigma^2_B = \sum_{j=0}^{nclass-1} \omega_j m_j^2 \right)$$

 

* Otsu 알고리즘을 이용한 이미지 이진화 코드: kipl.tistory.com/17

* 좋은 결과를 얻으려면 히스토그램에 적당한 필터를 적용해서 smooth하게 만드는 과정이 필요하다.

// 0 <= start < n;
double histo_partition(int nclass, double cxh[], int ch[], int n, int start, int th[]) {
    if (nclass < 1) return 0;
    if (nclass == 1) {
        int ws; double ms;
        if (start == 0) {
            ws = ch[n - 1];
            ms = cxh[n - 1] / ws;
        } else {
            ws = ch[n - 1] - ch[start - 1];             // start부터 끝까지 돗수;
            ms = (cxh[n - 1] - cxh[start - 1]) / ws;    // start부터 끝까지 평균값;
        }
        th[0] = n - 1;
        return ws * ms * ms;                            // weighted mean;
    }

    double gain_max = -1;
    int *tt = new int [nclass - 1];
    for (int j = start; j < n; j++) {
        int wj; double mj;
        if (start == 0) {
            wj = ch[j]; 
            mj = cxh[j];                    //mj = cxh[j] / wj;
        }
        else {
            wj = ch[j] - ch[start - 1];     //start부터 j까지 돗수;
            mj = (cxh[j] - cxh[start - 1]); //mj = (cxh[j] - cxh[start - 1]) / wj;
        }
        if (wj == 0) continue;
        mj /= wj;                           //start부터 j까지 평균;
        double gain = wj * mj * mj + histo_partition(nclass - 1, cxh, ch, n, j + 1, tt);
        if (gain > gain_max) {
            th[0] = j;
            for (int k = nclass - 1; k > 0; k--) th[k] = tt[k - 1];
            gain_max = gain;
        }
    }
    delete [] tt;
    return gain_max;
};

trimodal 분포의 분리;

class0: 0~th[0]

class1: (th[0]+1)~th[1],

class2: (th[1]+1)~th[2]=255

th[0]=103, th[1]=172&nbsp; (th[2]=255)
th[0]=88, th[1]=176, th[2]=255

더보기
// recursive histo-partition 테스트;
// 0--t[0];--t[1];...;--t[nclass-2];t[nclass-1]=255=n-1;
// nclass일 때 threshold 값은 0---(nclss-2)까지;
double GetThreshValues(int hist[], int n, int nclass, int th[]) {
    if (nclass < 1) nclass = 1;
    // preparing for 0-th and 1-th cumulative histograms;
    int *ch = new int [n];          // cdf;
    double *cxh = new double [n];   //1-th cdf;
    ch[0] = hist[0];
    cxh[0] = 0; // = 0 * hist[0]
    for (int i = 1; i < n; i++) {
        ch[i] = ch[i - 1] + hist[i] ;
        cxh[i] = cxh[i - 1] + i * hist[i];
    }
    // nclass=1인 경우도 histo_partition()내에서 처리할 수 있게 만들었다.
    double var_b = histo_partition(nclass, cxh, ch, n, 0, th);
    delete [] ch;
    delete [] cxh;
    return var_b;
}
 
728x90

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

Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Binary Image에서 Convex Hull  (0) 2021.01.06
Kuwahara Filter  (2) 2020.12.28
Moving Average을 이용한 Thresholding  (0) 2020.11.26
Posted by helloktk
,

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

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

$$B(i,j) = \left\{\begin{array}{ll} m_0, & \text{if}~A(i,j) \le T\\ m_1, &\text{otherwise}\end{array}\right. $$

로 나타난다. 여기서 $m_0$는 배경 픽셀의 평균값이고, $m_1$은 전경 픽셀의 평균값이다. 이 값은 임계값 $T$에 따라 달라진다. 임계값이 높으면 $m_0$는 커지고, 반대로 $m_1$은 작아진다

 

임계값이 $T$일 때 배경 픽셀 비를 $p$, 전경 픽셀 비를 $q(=1- p)$라 하면 segmented된 영상 B는 각 영역에서의 픽셀 값을 평균으로 대체했으므로 원본 영상의 평균과 같다. 또한, 원본 영상의 분산은 임계값에 무관하게 일정한 값을 유지한다. 이를 정리하면,

$$E(A)=E(B)=m=\text{pixel mean}=p m_0 + q m_1$$

$$V(A)=\text{variance} =T\text{-independent} = \text{const}$$

$$V(B)=pm_0^2 + q m_1^2 - m^2 = pq (m_0 - m_1)^2$$

$$E(A,B)= p m_0^2 + q m_1^2 $$

$$E(A,B) - E(A) E(B) = V(B)$$ 이므로, 

\begin{align}\text{Correlation}(A,B) &=\frac{ {E(A,B)-E(A)E(B)} }{\sqrt{V(A)V(B)} } \\ &=\frac{\sqrt{pq(m_0 - m_1)^2 } }{\sqrt{V(A)} }\\ &\propto \sqrt{pq(m_0 -m_1)^2 }\\ &=\sqrt{\text{interclass variance}}\end{align}

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

 

참고: Otsu Algorithm 구현 예.

728x90

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

Is Pow of 2  (0) 2012.02.13
Fixed-point RGB2Gray  (0) 2012.01.25
Object Orientation  (1) 2010.01.17
Bicubic Interpolation  (1) 2010.01.14
Bezier Curve을 이용한 Histogram Smoothing  (0) 2010.01.10
Posted by helloktk
,

이미지에서 어떤 유용한 정보를 추출하기 위해서는 이미지가 담고 있는 객체들을 분리하는 작업을 해야 한다. 가장 단순한 것 방법 중의 하나가 이진화(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-번째 cumulative histogram과 1-번째 cumulative histogram을 사용하면 각 클래스의 가중치와 평균을 쉽게 계산할 수 있다. 
    cumulative_histogram_0th[k] = 0...k 까지 값이 나타날 확률.
    cumulative_histogram_1th[k]/
cumulative_histogram_0th[k] = 0...k까지 값의 평균.


Otsu 알고리즘은 2개의 클래스뿐만 아니라 여러 클래스로 히스토그램을 분리하도록 확장할 수 있다. 재귀 알고리즘을 이용하면 쉽게 구현할 수 있다. (see: kipl.tistory.com/258)

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

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

    // make 0-th and 1-st cumulative histogram;
    chist[0] = hist[0];
    cxhist[0] = 0;
    for (int 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-st cumulative histogram ;
    };
    
    double gain_max = 0.;
    int thresh = 0;   
    double m = cxhist[255];                     //total mean = q1 * m1 + q2 * m2;
    int mul_count = 1;                          //number of degenerate maxima;
    for (int 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 (int i = 0; i < ntot; i++) dst[i] = (src[i] >= thresh) ? 0xFF : 0x00 ;

    return thresh;
}

 

사용자 삽입 이미지

 

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

사용자 삽입 이미지

728x90

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

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