Bright Preserving Histogram Equalization with Maximum Entropy

Histogram Equalization(HE)은 영상의 누적 히스토그램을 이용하여 픽셀값을 재분배함으로써 명암 대비를 향상시키는 대표적인 영상처리 기법이다. 이상적인 연속적 상황에서는 변환된 픽셀값이 전체 명암 범위에서 균일하게 분포하도록 하지만, 실제 디지털 영상에서는 픽셀값의 이산성과 중복 때문에 균일분포에 가깝게 재분배한다. 픽셀값이 일정한 구간에서 균일하게 분포하는 경우, 별도의 제약조건이 없다면 그 분포는 최대 엔트로피를 갖는다. 따라서 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);
}

std::vector<int> hist_spec(const std::vector<double> &chist1,
                           const std::vector<double> &chist2) 
{
    std::vector<int> lut(chist1.size());                           
    for (int i = 0; i < chist1.size(); i++) {
        int j = chist1.size() - 1;
        while (chist2[j] > chist1[i]) j--;
        lut[i] = j < 0 ? 0 : j
    }
    return lut;
}

std::vector<BYTE> BPHEME(const std::vector<BYTE> &src) {
    std::vector<double> hist = make_hist(src);
    normalize_hist(hist);
    // cumulative histogram;
    std::vector<double> chist = make_cumulative_hist0(hist);
    // gray-mean; note, pixel range should be changed [0,255] -> [0,1]
    double mu = hist_mean(hist);
    // determine lambda;
    double lambda, th = 1.e-15;
    FindRoot(mu, th, lambda);
    // dst-cdf;
    std::vector<double> chist2(256);
    for (int i = 0; i < 256; i++) 
        chist2[i] = cdf(double(i) / 255., mu, lambda);
    // histogram-specification;
    std::vector<int> lut = hist_spec(chist, chist2);
    // make dst image;
    std::vector<BYTE> dst(src.size());
    for (int i = src.size(); i-->0;) dst[i] = lut[src[i]];
    return dst;
};

// Root finding procedure ::

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_src) w.r.t. s;
double DF(double s, double mu) {   
    if (s == 0.) return -1.0 / 12.0;
    double ev = exp(s), ss = s * s;
    return -(ev * (ev - ss - 2.) + 1.) / (ss * (ev - 1.) * (ev - 1.));
}
// Find root of F(s,mu_src) based on Newton-Raphson method.
int FindRoot(double mu_src, double threshold, double& s) {
    int iter = 0, max_iter = 1000;
    s = 0. ; //initial guess; problem specific!
    for (iter = 0; iter < max_iter; ++iter) {
        double sold = s;
        // ASSERT(DF)
        s = s - F(s, mu_src ) / DF(s, mu_src) ;
        if (fabs(s - sold) < threshold)
            break;
    } ;
    TRACE("terminating iters=%d\n", iter);
    return (iter < max_iter) ;
};

histogram equaliztion 결과;

BPHEME  결과: