Minimum Cross Entropy Thresholding

$$threshold = \underset{0\le t < 256}{\text{argmin}} \left[-p_f(t) \log(\mu_f(t) )- p_b (t) \log( \mu_b (t)) \right] $$

where

$$ p_b (t) =\int_0^t h(z)dz, \quad  p_f(t) = \int_t^{255} h(z)dz$$

$$ \mu_b(t) = \frac{1}{p_b(t)} \int_0^{t} z h(z) dz,\quad \mu_f(t) = \frac{1}{p_f(t)} \int_t^{255} z h(z)dz $$

중간=Otsu, 마지막=MCE

Ref: Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum Cross Entropy Thresholding"Pattern Recognition Letters, 18(8): 771-776

double MCE_threshold(int hist[256]) {
    int chist[256], cxhist[256];
    chist[0] = hist[0]; cxhist[0] = 0;
    for (int i = 1; i < 256; i++) { 
        chist[i] = hist[i] + chist[i - 1];
        cxhist[i] = i * hist[i] + cxhist[i - 1];
    }
    int num = chist[255];
    double mean = double(cxhist[255]) / num;
    /* Initial estimate */
    double threshold = mean;
    while (1) {
        double old_thresh = threshold;
        int t = int(old_thresh + .5);
        /* background */
        int bgnum = chist[t];
        int bgsum = cxhist[t];
        double bgmean = bgnum == 0 ? 0: double(bgsum) / bgnum;
        /* foreground */
        int fgnum = num - bgnum;
        int fgsum = cxhist[255] - bgsum;
        double fgmean = fgnum == 0 ? 0: double(fgsum) / fgnum;
        threshold = (bgmean - fgmean) / (log(bgmean) - log(fgmean));
        // new thresh is a simple round of theta;
        ASSERT(threshold >= 0);
        if (fabs(threshold - old_thresh) < 0.5)
           break;
    }
    return threshold;
}