'Binarization'에 해당되는 글 5건

  1. 2012.02.04 Integral Image을 이용한 Adaptive Threshold
  2. 2010.01.28 Otsu-알고리즘의 새로운 해석
  3. 2008.07.14 Adaptive Binarization (1)
  4. 2008.07.01 EM : Binarization
  5. 2008.05.30 Otsu Algorithm (6)
Adaptive threshold 방법을 적용하는 데 있어서 윈도우 계산의 로드를 줄이는 방법은 integral image을 이용하면 된다(물론 메모리의 소요가 부가적으로 발생하지만, 현재의 기기에서는 메모리는 별로 문제가 안된다).

아래의 코드는 integral 이미지를 이용해서 주어진 moving 윈도우 내의 픽셀합을 구하고 그것의 평균 (= local average)을 기준으로 영상을 이진화 시키는 함수다 (정확히는 "평균값-3" 이다. 여기서 3은 바코드 인식 공개라이브러리인 zbar에서 쓰는 기준을 잡았다. zbar 라이브러리에서는 moving average를 구하여 임계값으로 사용하는데, 윈도우가 이동하면서 나가는 픽셀과 들어오는 픽셀을 업데이트 하는 과정이  이 라이브러리에서는 정확히 구현이 되어있지는 않다. 그러나, 근사적으로는 맞게 구현되어 있으므로 코드는 대부분의 경우에 원하는데로 잘 동작을 한다. 적분영상을 이용하면, 윈도우가 이동하면서 픽셀정보를 업데이트 하는 복잡한 과정이 필요없이 적분영상의 단순 합/차만 수행하면 된다)

"윈도우 평균-3" 대신 보다 합리적인 기준을 잡으려면, 보통은 윈도우의 표준편차을 이용할 수 있다. 그러나 이 경우에는 합의 제곱에 대한 적분영상이 하나 더 필요하고, 얼마의 편차를 허용할 것인지를 정해주어야 한다. 이 기준에 맞게 구현된 코드는 http://kipl.tistory.com/30 에서 찾을 수 있다.

2차원 바코드가 아닌 일차원 바코들 영상을 이진화할 때는 이 만큼 복잡한(?) 알고리즘을 쓸 필요가 없다. 일차원 바코드는 한 스캔라인의 정보만으로도 보통 인식이 가능하므로 라인단위의 이진화면 충분히다. 이 경우에도 이동평균을 사용하면 매우 간단하게, 그리고 adaptive한 임계값을 구할수 있는데, 라인기준이므로 적분영상이 따로 필요하지 않다.

void makeIntegralImage(BYTE *image, int width, int height, int* intImage) {

더보기

/*
** moving window의 중심에 해당픽셀을 놓을 필요는 없다; 
*/
void thresholdByIntegralImage(BYTE *image, int width, int height, int wsz,
                              BYTE *matrix) { 
    std::vector<int> intImage(width * height);
    makeIntegralImage(image, width, height, &intImage[0]);
    const int winArea = wsz * wsz ;
    /*const int wsz = 10;*/
    for (int y = 0, offset = 0; y < height; y++, offset += width) {
        int top = y - (wsz>>1) ;
        if (top < 0 ) top = 0;
        else if (top > height - wsz) top = height - wsz;
        int bottom = top + wsz - 1;
        // y-range = [top, bottom];
        for (int x = 0; x < width; x++) {
            int left = x - (wsz>>1);
            if (left < 0) left = 0;
            else if (left > width - wsz) left = width - wsz;
            int right = left + wsz - 1;
            // xrange = [left, right];
            //
            int sum1 = (left > 0  && top > 0) ? intImage[(top - 1) * width + left - 1] : 0;
            int sum2 = (left > 0) ? intImage[bottom * width + left - 1] : 0;
            int sum3 = (top > 0) ? intImage[(top - 1) * width + right] : 0;
            //
            int graySum = intImage[bottom * width + right] - sum3 - sum2 + sum1;
            // overflow ? 
            // Threshold T = (window_mean - 3); why 3?
            if ((image[offset + x] + 3) * winArea <= graySum) {
                matrix[offset + x] = 0xFF; //inverted!
            } else {
                matrix[offset + x] = 0x00;
            }
        }
    }
}



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

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

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
QR-code : decoder  (0) 2012.01.26
QR-code : detector  (0) 2012.01.12
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

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

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

이미지의 히스토그램을 이용하여서 전경과 배경을 분리하는 이진화 과정은 가우시안 mixture model을 이용하여서 EM 알고리즘을 적용하기에 좋은 예다. 전경에 해당하는 픽셀값의 분포와 배경에 해당하는 픽셀값의 분포는 히스토그램상에 섞여서 나타나는데. 이것을 두 가우시안의 혼합으로 생각할 때  em 알고리즘은 각각의 믹싱정도와 각 가우시안의 평균 및 표준편차를 추정한다. mixing parameter πa (a=1, 2,..., nclass)로 표시하는 경우에 특정한 픽셀값에 대한 posterior는


로 쓸수 있다. 여기서 f(x; θ)는 정규분포를 의미한다. 

posterior 정보를 이용하면 mixing parameter πa와 평균, 분산은 다음식으로 갱신이 된다. H[i]는 이미지의 히스토그램을 나타낸다.


   


log-likelihood:

여기서 a는 클래스 인덱스이고, i는 히스토그램 bin 인덱스이다.



// mixing 클래스를 기술하는 클래스;

struct mixclass {
    double prob ;               //mixing parameter;
    double mean ;               // mean
    double var ;                //variance;
} ;
// N(mean, var);

double gauss1d(double x, double mean, double var) {

더보기

// posterior; Pr(Zi = c | xi, Theta);
// 주어진 관측값 x이 클래스 cid에 속할 posterior;
double classprob(double x, int nclass, mixclass*  mclass, int cid) {

더보기

// posterior (class_prob[i][c]) table 만들기;
void update_class_prob(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {

더보기

// E-step;  pi[c] = mixture parameter for class c;
// posterior를 이용해서 특정클래스의 mixing 정도를 계산;==> next prior;
void update_prob(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {

더보기

// mu[c]; 클래스의 평균;
void update_mean(int nbins, double * hist, int nclass, mixclass* mclass,  double ** class_prob) {

더보기

// var[c]; 클래스의 분산;
void update_var(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {

더보기

// M-step; 
void update_parameters(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {

더보기

// initialization;
void init_em(int nbins, double * hist, int nclass, mixclass* mclass) {

더보기

// calculate log-likelihood;
double mixLLK(int nclass, mixclass* mclass) { 

더보기

// check termination condition;
bool check_tol(double llk, double llk_p, double  eps) {

더보기

// 입력은 이미지의 히스토그램;
double em(int nbins/*=256*/, double hist[/*256*/],
    int nclass/*=2*/, mixclass mclass[/*=2*/], double eps/*=1.e-10*/) {
    double llk = 0, prev_llk = 0;
    // allocate memory buffers for the posterior information;
    double ** class_prob = (double**)malloc(sizeof(double*) * nbins);
    class_prob[0] = (double*)malloc(sizeof(double) * nbins * nclass) ;
    for (int i = 1; i < nbins; i++) class_prob[i] = class_prob[i - 1] + nclass;

    // initialization of algorithm;
    init_em(nbins, hist, nclass, mclass);
    //
    do {
        prev_llk = llk;
        // E-step ;
        update_class_prob(nbins, hist, nclass, mclass, class_prob);
        // M-step;
        update_parameters(nbins, hist, nclass, mclass, class_prob);
        llk = mixLLK(nclass, mclass);
        // TRACE("mean1=%f, mean2=%f\n", mclass[0].mean, mclass[1].mean);
        TRACE("log-likelihood=%e\n", llk);
    } while (!check_tol(llk, prev_llk, eps));
    // clean ;
    free(class_prob[0]);
    free(class_prob) ;
    return llk;
};
  • 적색 : 히스토그램 
  • 청색, 녹색 : posterior(membership); 
  • Otsu 알고리즘을 쓰는 경우에 100에서 threshold 값이 결정되고 EM은 110 정도임.

사용자 삽입 이미지

신고

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

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
Shuffling  (0) 2008.06.21
Bayesian Decision Theory  (1) 2008.06.17
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차 및 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