Geometry & Recognition

Statistical region merging은 이미지의 픽셀을 일정한 기준에 따라 더 큰 영역으로 합병하는 bottom-up 방식의 과정이다. 두 영역 $R_1$과 $R_2$가 하나의 영역으로 합병이 되기 위해서는 두 영역의 평균 픽셀 값의 차이가 

$$g\sqrt{\frac{\ln (2/\delta )}{2Q}(\frac{1}{|R_1|}+\frac{1}{|R_2|})}$$

를 넘지 않아야 한다. $g=256$으로  gray 레벨의 갯수를 의미하고, $|R_i|$는 $R_i$ 영역에 포함된 픽셀 수를 나타낸다. $\delta$는 작은 수로 이미지의 픽셀 수의 제곱에 반비례한다. 보통 $\delta =1/(6\times \text{width}\times \text{height}{)}^2$로 선택한다. $Q$는 이미지의 통계적인 복잡성을 정량화하는 양으로 이 알고리즘에서는 외부에서 설정이 되는 값이다. 낮은 $Q$값을 선택하면 분할된 영역의 수가 작아지고(undersegmentation), 반대로 높은 $Q$ 값을 입력하면 분할된 영상에 너무 많은 영역이 나타나게 된다(oversegmentation).

Ref:

https://en.wikipedia.org/wiki/Statistical_region_merging

http://www.lix.polytechnique.fr/~nielsen/Srmjava.java

Srm::Srm(int width, int height, BYTE *image) {
    this->width  = width;
    this->height = height;
    int n        = width * height;
    this->count = new int [n];
    this->Ravg  = new float [n];
    this->Gavg  = new float [n];
    this->Bavg  = new float [n];
    this->image = image;
    // disjoint sets with n pixels;
    this->UF = new Universe(n);
    // initialize to each pixel a leaf region;
    for (int i = 0, pos = 0; i < n; i++, pos += 3) {
        count[i] = 1;
        Bavg[i] = image[pos    ];
        Gavg[i] = image[pos + 1];
        Ravg[i] = image[pos + 2];
    }
    this->Q = 32;		// adjustable.					    
    this->g = 256.0;
    this->logDelta = 2. * log(6.0 * n);
}
bool Srm::Predicate(int rgn1, int rgn2) {
    double dR = (Ravg[rgn1] - Ravg[rgn2]); dR *= dR;
    double dG = (Gavg[rgn1] - Gavg[rgn2]); dG *= dG;
    double dB = (Bavg[rgn1] - Bavg[rgn2]); dB *= dB;
    double logreg1 = min(g, count[rgn1]) * log(1.0 + count[rgn1]);
    double logreg2 = min(g, count[rgn2]) * log(1.0 + count[rgn2]);
    double factor = g * g / (2.0 * Q);
    double dev1 = factor * (logreg1 + logDelta) / count[rgn1] ;
    double dev2 = factor * (logreg2 + logDelta) / count[rgn2] ;
    double dev = dev1 + dev2;
    return ( (dR < dev) && (dG < dev) && (dB < dev) );
}
void Srm::Merge(int rgn1, int rgn2) {
    if (rgn1 == root2) return;
    int w1 = count[rgn1], w2 = count[rgn2];
    int root = UF->Union(rgn1, rgn2);
    //update the merged region;
    count[root] = w1 + w2;
    double count_sum = w1 + w2;
    Ravg[root] = (w1 * Ravg[rgn1] + w2 * Ravg[rgn2]) / count_sum;
    Gavg[root] = (w1 * Gavg[rgn1] + w2 * Gavg[rgn2]) / count_sum;
    Bavg[root] = (w1 * Bavg[rgn1] + w2 * Bavg[rgn2]) / count_sum;
}
Edge* Srm::Pairs(int nedge) {
    // 4-connectivity;
    int ymax = height - 1, xmax = width - 1;
    Edge* edgeList = new Edge[nedge];
    int cnt = 0;
    for (int y = 0; y < ymax; y++) {
        for (int x = 0; x < xmax; x++) {
            int pos = y * width + x;
            int b1 = image[3 * pos + 0];
            int g1 = image[3 * pos + 1];
            int r1 = image[3 * pos + 2];
            //right: x--x
            edgeList[cnt].r1 = pos;     //current
            edgeList[cnt].r2 = pos + 1; //right
            int bdiff = abs(b1 - image[3 * (pos + 1) + 0]);
            int gdiff = abs(g1 - image[3 * (pos + 1) + 1]);
            int rdiff = abs(r1 - image[3 * (pos + 1) + 2]);
            edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff) ;
            //below: x
            //       |
            //       x
            edgeList[cnt].r1 = pos;
            edgeList[cnt].r2 = pos + width;
            bdiff = abs(b1 - image[3 * (pos + width) + 0]);
            gdiff = abs(g1 - image[3 * (pos + width) + 1]);
            rdiff = abs(r1 - image[3 * (pos + width) + 2]);
            edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff);
        }
    }
    //x=width-1;
    for (int y = 0; y < ymax; y++) {
        int pos = y * width + (width - 1); // (x,y) = (width-1, y)
        // x
        // |
        // x
        edgeList[cnt].r1 = pos;
        edgeList[cnt].r2 = pos + width;
        int bdiff = abs((int)image[3 * pos + 0] - image[3 * (pos + width) + 0]);
        int gdiff = abs((int)image[3 * pos + 1] - image[3 * (pos + width) + 1]);
        int rdiff = abs((int)image[3 * pos + 2] - image[3 * (pos + width) + 2]);
        edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff);
    }
    //y=height-1;
    for (int x = 0; x < xmax; x++) {
        int pos = (height - 1) * width + x;      //(x,y)=(x, height-1);
        //right; x--x
        edgeList[cnt].r1 = pos;
        edgeList[cnt].r2 = pos + 1;
        int bdiff = abs((int)image[3 * pos + 0] - image[3 * (pos + 1) + 0]);
        int gdiff = abs((int)image[3 * pos + 1] - image[3 * (pos + 1) + 1]);
        int rdiff = abs((int)image[3 * pos + 2] - image[3 * (pos + 1) + 2]);
        edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff);
    }
    return edgeList;
}
int Srm::Segment() {
    // 4-connectivity 
    int nedge = 2 * (width - 1) * (height - 1) + (height - 1) + (width - 1);
    Edge* edgeList = Pairs(nedge);
    BucketSort(edgeList, nedge);
    for (int i = 0; i < nedge; i++) {
        Edge &e = edgeList[i];
        int r1 = UF->Find(e.r1);
        int r2 = UF->Find(e.r2);
        if ((r1 != r2) && (Predicate(r1, r2)))
            Merge(r1, r2);
    }
    delete [] edgeList;
    int rgn_count = 0;
    for (int node = width * height; node-- > 0;)
        if (UF->IsRoot(node)) rgn_count++;
    return rgn_count;
}
// sorting with buckets; returns an ordered edgeList;
void BucketSort(Edge* &edgeList, int n) {
    int hist[256] = {0}, chist[256];
    for (int i = 0; i < n; i++) hist[edgeList[i].diff]++;
    // cumulative histogram
    chist[0] = 0;  // Note, chist[0] ne hist[0];
    for (int i = 1; i < 256; i++)
        chist[i] = chist[i - 1] + hist[i - 1];

    Edge *ordered = new Edge [n];
    for (int i = 0; i < n; i++)
        ordered[chist[pair[i].diff]++] = pair[i];        
    delete[] edgeList;
    edgeList = ordered;
}

영상의 히스토그램($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;
}

$$threshold=\underset{0\le t<256}{\text{argmin}}[-p_f(t)\log ({\mu }_f(t))-p_b(t)\log ({\mu }_b(t))]$$

where

$$p_b(t)={\int }_0^{t}h(z)dz,p_f(t)={\int }_t^{255}h(z)dz$$

$${\mu }_b(t)=\frac{1}{p_b(t)}{\int }_0^{t}zh(z)dz,{\mu }_f(t)=\frac{1}{p_f(t)}{\int }_t^{255}zh(z)dz$$

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;
}