Geometry & Recognition

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

1. 정규화된  조건 (픽셀 값을 $[0,255] \rightarrow [0,1]$ 로 rescale 함)
   $$f(s) \ge 0, \quad \int_0^1 f(s)ds = 1.$$
2. 밝기 보존: 
   $$\int_0^1 sf(s) ds = \mu=\text{pixel mean of input image}.$$

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

$$\underset{f}{\text{argmax}}\int_0^1 \Big[ -f(s) \ln f(s)\Big] ds.$$

위 조건를 만족시키는 $f(s)$는 Lagrange-multiplier를 이용하고, 변분 방법을 쓰면 쉽게 구해진다: 다음의 범함수 $J [f]$

$$J=-\int_0^1  f(s) \ln f(s) ds + \lambda_1  \left( \int_0^1 f(s) ds-1\right) + \lambda \left(\int_0^1 s f(s) ds -\mu\right)$$

의 극값을 구하면 된다.

$$\frac{\delta J}{\delta f}=0$$

에서

$$ -\log f - 1 + \lambda_1 f +\lambda  s = 0$$

이를 정리하면,

$$f(s)=e^{\lambda s } e^{\lambda_1-1}$$

을 얻고, 제약 조건을 적용하면 ($ e^{\lambda_1 -1}=\lambda/(e^{\lambda}-1))$

$$f(s) = \left\{\begin{array}{ll} 1,& \mu = 0.5 \\ \frac{\lambda e^{\lambda s}}{e^\lambda -1},& \mu \ne 0.5 \end{array} \right.$$

그리고 $λ$는 

$$ \mu = \frac{\lambda e^\lambda - e^\lambda +1}{\lambda (e^\lambda -1)},\quad(\mu \ne 0.5)$$

의 근이다. 오른쪽은 $λ$에 대해서 단조함수이므로 근이 하나만 존재한다.  원본 이미지에서 픽셀 평균을 계산하여 위 방정식에 대입하면 $λ$를 구할 수 있고, 목표 히스토그램을 얻을 수 있다. 목표 히스토그램의 cdf와 ($\text {chist}(x)=\int_0^x sf(s) ds$), 원본 이미지의 cdf를 구해서, 그 차이를 최소로 하는 픽셀 값의 대응관계를 찾으면 된다 (histogram specification):

$$x \rightarrow y = F(x) = \int_0^x f(s) ds$$

참고: Bright Preserving Histogram Equalization with Maximum Entropy: A Variational Perspective.
        C. Wang and Z.Ye, IEEE Trans. Consumer Electronics. V51. No4. (2005);
// 알고리즘 적용 단계에서 픽셀 값의 범위를 [0,255] -> [0,1]로 변환해야 한다.
// BPHEME의 cumulative histogram(continuous version)::integral of f(s) over s;

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

더보기

    if (mu == 0.5) return s ;
    else return (exp(lambda * s) - 1) / (exp(lambda) - 1);
}

// 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); {

더보기

    for (int i = 0; i < n; i++) {
        int j = n - 1;
        while (chist2[j] > chist1[i]) j--;
        lut[i] = j < 0 ? 0 : j
    }
}

void BPHEME(BYTE* src, int width, int height, BYTE *dst) {
    double hist[256] = {0};                     // src(dst) histogram;
    double chist[256], chist2[256];             // src(dst) cumulative histogram;
    int lut[256];                               // histogram specification mapping;
    const int n = width * height;
    make_hist(src, width, height, hist);
    normalize_hist(hist, 256);
    //cumulative histogram;
    make_cumulative_hist0(hist, 256, chist);
    // gray-mean; note, pixel range should be changed [0,255] -> [0,1]
    double mu = hist_mean(hist, 256);
    // determine lambda;
    double lambda, th =1.e-15;
    FindRoot(mu, th, lambda);
    // entropy of src;
    double entropy1 = hist_entropy(hist, 256);
     // dst-cumulative
    for (int i = 0; i < 256; i++) chist2[i] = cdf(double(i) / 255., mu, lambda);
    // histogram-specification;
    hist_spec(&chist[0], &chist2[0], 256, &lut[0]);
    // make dst image;
    for (int i = 0; i < n; i++) dst[i] = lut[src[i]];
};

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

더보기

    if (s == 0)
        return (mu - 0.5);
    double ev = exp(s);
    return mu - (s * ev - ev + 1.) / (s * (ev - 1.)) ;
};

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

더보기

    if (s == 0)
        return -1.0 / 12.0;
    double ev = exp(s);
    double ss = s * s;
    return -(ev * (ev - ss - 2.) + 1.) / (ss * (ev - 1.)*(ev - 1.));
}

// Find root of F(s,mu) based on Newton-Raphson method.
int FindRoot(double mu, double threshold, double& s) {     
     int iter = 0, max_iter = 100;
     s = 0. ; // initial guess; problem specific!
     for (iter = 0; iter < max_iter; ++iter) {
          double sold = s;
          // ASSERT(DF)
          s = s - F(s, mu) / DF(s, mu) ;
          if (fabs(s - sold) < threshold)
              break ;
     } ;
     TRACE("terminating iters=%d\n", iter);
     return (iter < max_iter) ;
};

histogram equaliztion 결과;

BPHEME  결과:

영상에는 보통 우리가 관심이 있는 물체와 배경을 동시에 담고 있다. 그런데 어떤 영상의 경우 전체 픽셀 값의 다이내믹 레인지(dynamic range)가 매우 좁아서, 달리 이야기하면 영상의 히스토그램을 살펴볼 때 픽셀 값의 분포가 매우 좁은 영역에 한정이 되어있으면 물체와 배경의 구별이 잘 안되는 영상이 되게 된다. 이러한 영상의 경우 몇 가지 영상 처리 방법을 이용하면 좋은 contrast를 갖는 영상을 만들 수 있다.

Histogram Equalization(HE)은 이러한 기법 중의 하나로 좁은 영역에 분포하는 픽셀 값을 가능한 모든 범위에 골고루 분포하도록 재배치하여서 영상의 contrast를 증대시키는 방법이다. Histogram Equalization에서 사용하는 픽셀 값의 mapping은 

(1) 출력 영상의 픽셀 값 범위는 전 그레이 레벨이어야 하고 

(2) 각각의 그레이 레벨에 해당하는 픽셀의 수가 거의 일정하여야 한다: 픽셀 값이 이산적이어서 완벽하게 상수가 될 수는 없다, 좀 더 정확히는 출력 영상의 히스토그램 누적합이 그레이 레벨에 선형 비례해야 한다는 조건으로 결정된다.

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

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

입력 영상에서 그레이 값이 0부터 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] = {0}, cumhist[256];
    for (int k = width * height; k-- > 0; ) 
        histogram[src[k]]++ ;
    
    // make 0-st order cumulative histogram;
    cumhist[0] = histogram[0];
    for (int k = 1; k < 256; k++)
        cumhist[k] = cumhist[k - 1] + histogram[k] ;
    
    // make output;
    // cumhist[255] = number of pixels;
    double v = 255.0 / cumhist[255];
    for (int k = width * height; k-- > 0; ) { 
        int  a = int(v * cumhist[src[k]]);
        dst[k] = a < 0 ? 0 : a > 255 ? 255 : a;  //clamp to [0,255] 
    }
};

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

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