728x90

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

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

$$ th= \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  (th[2]=255)
th[0]=88, th[1]=176, th[2]=255

 

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

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

댓글을 달아 주세요

728x90

영상의 히스토그램을 두 개의 가우시안 분포의 혼합으로 모델링하여 분리하려고 할 때 기준인 decision boundary 값 (threshold value)을 expectation maximization(EM) 알고리즘을 적용하여 구한다. 

 

 

E-step: compute responsibility of class 2; (for class 1, 1-γ_i)

 

 

M-step: compute the weighted means (μ1, μ2), variances (σ1, σ2) and mixing probability (π)

 

 

log-likelihood: 

decision boundary 값은 responsibility = 0.5인  bin 인덱스를 선택하면 된다.

 

아래 그림의 왼쪽은 히스토그램, 오른쪽은 최대우도 gaussian fitting 결과와 왼쪽 분포의 responsibility(1-γ_i)를 그린 것이다.

void estimGaussParams(double data[], int start, int end, double *mean, double *var) ;

더보기
void estimGaussParams(double data[], int start, int end, double *mean, double *var) {
    double s = 0, sx = 0, sxx = 0;
    for (int i = start; i < end; i++) {
        s += data[i];
        sx += data[i] * i;
        sxx += data[i] * i * i;
    }
    *mean = sx / s;
    *var = (sxx - sx * sx / s) / s;
};

void initGuess(double data[], int n, double mean[], double var[], double *mixprob);

더보기
void initGuess(double data[], int n, double mean[], double var[], double *mixprob) {
    int start = 0, end = n; 
    // trim null data;
    while (data[start] <= 0) start++;
    while (data[end - 1] <= 0) end--;
    // split given data into two equal size sets;
    int n2 = (end + start) / 2;
    // simple mean and variance;
    estimGaussParams(data, start, n2, &mean[0], &var[0]);
    estimGaussParams(data, n2, end, &mean[1], &var[1]);
    // initial guess for mixing probability;
    *mixprob = 0.5; 
};

#define PI (4.0 * atan(1.))

double gaussDist(double x, double mean, double var) ;  

더보기
double gaussDist(double x, double mean, double var) { 
    // N(mean, var);
    double arg = 0.5 * (x - mean) * (x - mean) / var;
    double factor = 1 / sqrt(2.* PI * var);
    return factor * exp(-arg); 
}

double responsibility2(double x, double mean[], double var[], double mixprob) ;   

더보기
double responsibility2(double x, double mean[], double var[], double mixprob) {   
    double a = (1 - mixprob) * gaussDist(x, mean[0], var[0]);
    double b = mixprob * gaussDist(x, mean[1], var[1]);  
    return b / (a + b); 
}

double weightedMeanVar(double data[], double gamma[], int n, double mean[], double var[]) ;

더보기
double weightedMeanVar(double data[], double gamma[], int n, double mean[], double var[]) { 
	// estimate new means;
    double s = 0, sx0 = 0, sx1 = 0, sg = 0;
    for (int i = 0; i < n; i++) {
        s   += data[i];
        sg  += data[i] * gamma[i]; 
        sx0 += data[i] * i * (1 - gamma[i]);
        sx1 += data[i] * i * gamma[i];
    }
    mean[0] = sx0 / (s - sg);
    mean[1] = sx1 / sg;
    // variances with new mean;
    double sv0 = 0, sv1 = 0;
    for (i = 0; i < n; i++) {
        sv0 += data[i] * (i - mean[0]) * (i - mean[0]) * (1 - gamma[i]);
        sv1 += data[i] * (i - mean[1]) * (i - mean[1]) * gamma[i];
    }
    var[0] = sv0 / (s - sg);
    var[1] = sv1 / sg;
    // return mixing probability = mixing ratio for class 2;
    return (sg / s);
};
#define EPSILON  1e-6
// Expectation Maximization algorithm applied to Two component Gaussian Mixture Model;
double emTwoCompGMM(double data[], int n) {
    double mean[2], var[2], mixprob;
    double *gamma = new double [n];     // responsibilities for class 2;
    initGuess(data, n, mean, var, &mixprob);
    // begin algorithm;
    while (1) {
        // E-step;
        for (int i = 0; i < n; i++) 
            gamma[i] = responsibility2(i, mean, var, mixprob);
        double old_mixprob = mixprob;
        // M-step;
        mixprob = weightedMeanVar(data, gamma, n, mean, var);
        TRACE("mixing probability= %f\n", mixprob);
        // check convergence(usually loglikelihood is tested);
        if (fabs(mixprob - old_mixprob) < EPSILON)
            break;
    }
    // estimate decision boundary;
    int k = n - 1;
    for ( ; k >= 0; k--) 
        if (gamma[k] < 0.5) break;
    delete [] gamma;
    return (2 * k + 1) / 2.; // = average of ;
};

Posted by helloktk

댓글을 달아 주세요

728x90

이미지를 이진화시키기 위해서 여러 알고리즘이 사용된다. 그중 이미지 전체에 대해 하나의 임계값으로 이진화시키는 전역 이진화 알고리즘은 간단하고 빠르기 때문에 많이 이용이 된다. 그러나 이미지를 형성할 때 조명 조건이 균일하지 않은 경우에는 전역 이진화는 원하는 결과를 얻기가 힘들다. 이런 경우에는 각각의 픽셀 주위의 그레이 값을 참조하여 임계치를 결정하는 국소적 이진화 방법을 사용한다. 국소적 이진화에서 임계값을 추출하는 간단한 방법은 윈도 내의 평균값을 이용하면 된다. 좀 더 개선된 알고리즘은 평균값($m(x, y)$)을 참조하되, 편차($\sigma(x, y)$)를 한번 더 고려해 주는 것이다. 이렇게 하여 잡은 국소적 임계값은 다음과 같이 표현된다: 

$$T_{(x, y)} = m_{(x, y)} [1+ \text{factor}(\sigma_{(x, y)}-128)]$$

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

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

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

더보기
void make_int_img12(BYTE *gray, int width, int height, *int intimage, int *intsqimg) {
    // first row accumulation;
    intimage[0] = gray[0];
    for (int x = 1; x < width; ++x) {
        int a = gray[x] ;
        intimage[x] = intimage[x - 1] + a;
        intsqimg[x] = intsqimg[x - 1] + a * a;
    }
    for (int y = 1, pos = y * width; y < height; ++y) {
        int linesum = 0, linesqsum = 0 ;
        for (int x = 0; x < width; ++x, ++pos) {
            int a = gray[pos];
            linesum   += a;
            linesqsum += a * a;
            intimage[pos] = intimage[pos - width] + linesum ;
            intsqimg[pos] = intsqimg[pos - width] + linesqsum;
        }
    }
};
#define integral_image(x, y) (intimage[(y) * width + (x)])
#define integral_sqimg(x, y) (intsqimg[(y) * width + (x)])
//
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, pos = 0; j < height; j++) {
        for (int i = 0; i < width; i++, pos++) {
            // 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) / area;
            double std  = sqrt((sqdiff - double(diff) * diff / area) / (area - 1));
            double threshold = mean * (1.0 + k * ((std / 128.0) - 1.));
            if (gray[pos] < threshold) bimage[pos] = 0;
            else                       bimage[pos] = 255;
        }
    }   
};

사용자 삽입 이미지

 

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

Fant's Resampling  (0) 2008.12.17
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
Otsu Algorithm  (6) 2008.05.30
Posted by helloktk

댓글을 달아 주세요

  1. 박영우 2011.08.26 13:01  댓글주소  수정/삭제  댓글쓰기

    글 잘 읽었습니다.
    좋은정보 주셔서 감사합니다.
    그런데 적응이진화를 적용하니 좌측상단에 정사각형 모양으로 점이 하니 찍히네요
    이것을 제가 할수는 없는건가요??

  2. 진만시 2020.04.13 13:04 신고  댓글주소  수정/삭제  댓글쓰기

    시간을 단축시킬만한 방법이 있을까요 ?

728x90

이미지의 히스토그램을 이용하여 전경과 배경을 분리하는 이진화는 가우시안 mixture model과 EM 알고리즘을 적용하기에 좋은 예다. 히스토그램에는 전경에 해당하는 픽셀 분포와 배경에 해당하는 픽셀 분포가 혼합되어 있다. 이를 두 가우시안의 혼합으로 모델링하고 EM 알고리즘을 사용해서 mixing parameter(πa), 각 클래스의 평균(μa) 과 표준편차(σa)를 추정한다. N개의 Gaussian mixture일 때,

Mixing parameter가 πa (a=1, 2,..., nclass)일 때 특정 픽셀 (=xi)이 클래스 a 소속일 posterior 확률은

 
로 쓸 수 있다. posterior 정보를 이용하면 mixing parameter, 평균 그리고 분산은 다음 식으로 주어진다. H[i]=Hi는 이미지의 히스토그램을 나타내고, bin 인덱스 i는 픽셀 값 xi를 나타낸다:
 

   

 

 

log-likelihood:

.


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

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

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

더보기

    double a = 1 / sqrt(2*M_PI * var);
    double b = 0.5*(x-mean)*(x-mean)/var;
    return a * exp(-b);
};

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

더보기
    double marginal = 0;
    for (int c = 0; c < nclass; c++) {
        marginal += mclass[c].prob * gauss1d(x, mclass[c].mean, mclass[c].var) ;
    };
    // Bayes 공식 = prior * PDF;
    return mclass[cid].prob * gauss1d(x, mclass[cid].mean, mclass[cid].var) / marginal;
}
// posterior (class_prob[i][c]) table 만들기;
void update_class_prob(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
더보기
        for (int i = 0; i < nbins; i++) {
            for (int c = 0; c < nclass; c++) {
                class_prob[i][c] = classprob(double(i), nclass, mclass, c);
            }
        }
};
// 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) {
더보기
        double ntot = 0;
        for (int i = 0; i < nbins; i++) ntot += hist[i];
        for (int c = 0; c < nclass; c++) {
            double s = 0;
            for (int i = 0; i < nbins; i++) s += hist[i] * class_prob[i][c];
            mclass[c].prob = s / ntot;
        }
};
// mu[c]; 클래스의 평균;
void update_mean(int nbins, double * hist, int nclass, mixclass* mclass,  double ** class_prob) {
더보기
        double ntot = 0;
        for(int i=0; i<nbins; i++) ntot += hist[i];
        for (int c = 0; c < nclass; c++) {
            double sx = 0.0;
            for (int i = 0; i < nbins; i++) sx += hist[i] * i * class_prob[i][c];
            mclass[c].mean = sx / (ntot * mclass[c].prob);
        }
};
// var[c]; 클래스의 분산;
void update_var(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
더보기
    double ntot = 0;
    for (int i = 0; i < nbins; i++) ntot += hist[i];
    for (int c = 0; c < nclass; c++) {
        double m= mclass[c].mean ;
        double sxx = 0;
        for (int i = 0; i < nbins; i++) sxx += hist[i] * SQR(i - m) * class_prob[i][c];
        mclass[c].var = sxx / (ntot * mclass[c].prob);
    }
};
// M-step; 
void update_parameters(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
더보기
    // mixture 파라미터를 갱신;
    update_prob(nbins, hist, nclass, mclass, class_prob);
    // 각 클래스의 평균을 갱신;
    update_mean(nbins, hist, nclass, mclass, class_prob);
    // 각 클래스의 분산을 갱신;
    update_var(nbins, hist, nclass, mclass, class_prob);
};
// initialization;
void init_em(int nbins, double * hist, int nclass, mixclass* mclass) {
더보기
        srand(unsigned(time(0)));
        double mean1 = 0, var1 = 0, ntot = 0;
        int i, k;
        for (k = 0; k < nbins; k++) ntot += hist[k];
        for (i = 0; i < nbins; i++) mean1 += hist[i] * i;
        mean1 /= ntot;
        for (i = 0; i < nbins; i++) var1 += hist[i] * SQR(i - mean1);
        var1 /= ntot;
        for (int c = 0; c < nclass; c++) {
            mclass[c].prob = 1.0 / nclass;          //same mixing parameter;
            mclass[c].mean = rand() % nbins; // random mean;
            mclass[c].var = var1;                     // same standard deviation;
        }
};
// calculate log-likelihood;
double mixLLK(int nclass, mixclass* mclass) { 
더보기
    double llk = 0;
    for (int i = 0; i < nbins; i++) {
        double s = 0 ;
        for (int c = 0; c < nclass; c++) 
            s += mclass[c].prob * gauss1d(double(i), mclass[c].mean, mclass[c].var);
        llk+= log(s);
    }
    return llk;
};
// check termination condition;
bool check_tol(double llk, double llk_p, double  eps) {
더보기
    return (fabs(llk - llk_p) / fabs(llk)) > 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 정도임.

  • Otsu Threshold source code: kipl.tistory.com/17

 

사용자 삽입 이미지

 

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

KMeans Algorithm  (0) 2008.07.19
Robust Line Fitting  (0) 2008.07.08
EM: Binarization  (0) 2008.07.01
EM Algorithm: Line Fitting 예  (0) 2008.06.29
Gaussian Mixture Model  (2) 2008.06.07
Rasterizing Voronoi Diagram  (0) 2008.05.26
Posted by helloktk

댓글을 달아 주세요

728x90

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

사용자 삽입 이미지

'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
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방법은 이 최소값이 중요한 것이 아니라 전체적인 봉우리들의 모양이 중요합니다. 그리고 봉우리가 여럿있을 떄 반드시 봉우리 사이의 계곡에서 최소인지점에서 나타난다고 할 수 는 없습니다. 따라서 오수 알고리즘을 쓰면, 이 알고리즘이 주는 값에서 정해야 하고, 결과가 원하는 것이 아니면, 이 알고리즘을 쓰면 안되고 다른 알고리즘으로 구현하는 것이 더 적합할 것입니다. 이 알고리즘이 반드시 가장 좋은 결과를 주지는 않습니다.