Moment-preserving Thresholding

영상의 히스토그램($h[z]$)이 bimodal로 주어지는 경우 적절한 threshold 값을 선택해서 전경과 배경을 분리할 수 있다. 전경을 대표하는 픽셀 값을 $z_f$, 배경을 대표하는 픽셀 값을 $z_b$라면 이진화 후 정규화된 히스토그램은 

$$\tilde{h}[z]=p_b{\delta }_{z,z_b}+p_f{\delta }_{z,z_f}$$

로 표현된다. $p_b$은 배경에 해당하는 픽셀 비율이고, $p_f$은 전경에 해당하는 픽셀 비율이다.

threshold 값을 어떻게 선택하면 이진화된 영상의 히스토그램이 원 영상의 히스토그램의 특성을 최대한 담게 할 수 있을까? 이에 대한 기준으로 두 히스토그램의 $n$차 moment가 같은 값을 갖도록 하자. 주어진 미지수가 $p_b$, $p_f$, $z_b$, $z_f$이 있으므로 최소한 4개의 moment가 같도록 만들어야 한다. 가장 낮은 찾수의 moment로부터 시작해서 3차까지 4개의 moments가 같다는 조건에서 아래의 식들을 얻을 수 있다.

$$\text{0-차 moment: }{m}_0\equiv \sum _{z=0}^{255}h[z]=p_b+p_f=1$$

$$\text{1-차 moment: }{m}_1\equiv \sum _{z=0}^{255}zh[z]=p_b{z}_b+p_f{z}_f$$

$$\text{2-차 moment: }{m}_2\equiv \sum _{z=0}^{255}{z}^{2}h[z]=p_b{z}_b^2+p_f{z}_f^2$$

$$\text{3-차 moment: }{m}_3\equiv \sum _{z=0}^{255}{z}^{3}h[z]=p_b{z}_b^3+p_f{z}_f^3$$

원 영상의 moment $m_0$, $m_1$, $m_2$, $m_3$을 계산해서 풀면

$$\begin{array}{c}{c}_0=\frac{m_3{m}_1-m_2^2}{m_0{m}_2-m_1^2},c_1=\frac{m_1{m}_2-m_0{m}_3}{m_0{m}_2-m_1^2}\\ z_b=\frac{1}{2}(-c_1-\sqrt{c_1^2-4c_0})\\ z_f=\frac{1}{2}(-c_1+\sqrt{c_1^2-4c_0})\\ p_b=\frac{z_f-m_1}{z_f-z_b}\\ p_f=1-p_b\end{array}$$

따라서 threshold 값

$$\sum _{z=0}^{T-1}h[z]=p_b$$

을 만족하는 $T$을 선택하면 된다.

 

Ref: W. Tsai, "Moment-preserving thresholding: a new approach," Computer Vision, Graphics, and Image Processing, vol. 29, pp. 377-393, 1985.

int MomentsPreseving_threshold(int histogram[256]) {
    int tot = 0;
    for (int i = 0; i < 256; i++)
        tot += histogram[i];
    //normalised histogram
    double hist[256];
    for (int i = 0; i < 256; i++)
        hist[i] = double(histogram[i]) / tot;
    /* moments calculation: zero moment is 1 by defintion*/
    double m0 = 1, m1 = 0, m2 = 0, m3 = 0;
    for (int i = 0; i < 256; i++ ) {
        double h = hist[i];
        m1 += i * h;
        m2 += i * i * h;
        m3 += i * i * i * h;
    }
    double det = m0 * m2 - m1 * m1;
    double c0 = (m1 * m3 - m2 * m2) / det;
    double c1 = (m2 * m1 - m3 * m0) / det;
    double zb = 0.5 * (-c1 - sqrt (c1 * c1 - 4.0 * c0));
    double zf = 0.5 * (-c1 + sqrt (c1 * c1 - 4.0 * c0));
    double pb = (zf - m1) / (zf - zb);  
    double s = 0;
    for (int i = 0; i < 256; i++) {
        s += hist[i];
        if (s > pb)
            return i; // threshold
    }
    return 0;
}