Bright Preserving Histogram Equalization with Maximum Entropy

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  결과: