Statistical Region Merging

Statistical Region Merging(SRM)은 인접한 픽셀들을 점진적으로 하나의 영역으로 통합해 나가는 bottom-up 방식의 영상 분할(image segmentation) 알고리즘이다. 알고리즘은 처음에 모든 픽셀을 각각 하나의 독립된 영역으로 간주한 뒤, 서로 인접한 두 영역의 밝기 값이 통계적으로 충분히 유사하다고 판단되면 하나의 영역으로 합병한다. 이러한 과정을 더 이상 합병할 수 있는 영역이 없을 때까지 반복하여 최종적인 영상 분할 결과를 얻는다.

두 영역 \(R_1\)과 \(R_2\)의 평균 픽셀 값을 각각 \(\bar{I}(R_1)\), \(\bar{I}(R_2)\)라 하면, 두 영역은

\[\left|\bar{I}(R_1)-\bar{I}(R_2)\right|\le g\sqrt{\frac{\ln(2/\delta)}{2Q} \left(\frac{1}{|R_1|}+\frac{1}{|R_2|}\right)} \]

를 만족할 때 하나의 영역으로 합병된다.

여기서

• \(g\)는 영상의 밝기 범위(dynamic range)를 나타내며, 8비트 grayscale 영상에서는 \(g=256\)이다.
• \(|R_i|\)는 영역 \(R_i\)에 포함된 픽셀의 개수이다.
• \(\delta\)는 잘못된 영역 합병이 발생할 확률을 제한하기 위한 작은 양수이다. 실제 구현에서는 일반적으로
  \[\delta=\frac{1}{(6\times\text{width}\times\text{height})^2}\]와 같이 설정한다.
• \(Q\)는 사용자가 설정하는 복잡도(complexity) 파라미터로, 영상을 어느 정도의 세밀함으로 분할할 것인지를 결정하는 역할을 한다.

합병 조건의 우변은 두 영역의 평균 밝기 차이에 대해 허용되는 최대 오차를 의미한다. 식에서 알 수 있듯이 허용 오차는 영역의 크기에 따라 달라진다. 즉, \(|R_1|\)과 \(|R_2|\)가 커질수록

\[\frac{1}{|R_1|}+\frac{1}{|R_2|}\]

의 값이 작아지므로 허용되는 평균 밝기 차이도 감소한다. 따라서 작은 영역은 평균 밝기 값에 다소 차이가 있어도 합병될 수 있지만, 큰 영역은 평균 밝기 값이 거의 같은 경우에만 합병된다. 이는 작은 영역에서는 평균 밝기의 통계적 불확실성이 상대적으로 큰 반면, 많은 픽셀을 포함하는 큰 영역에서는 평균값의 신뢰도가 높아지기 때문이다.

또한 \(Q\)는 분할 결과의 세밀함을 결정하는 가장 중요한 매개변수이다. 식에서 허용 오차는 \(1/\sqrt{Q}\)에 비례하므로, \(Q\)가 작을수록 합병 기준이 완화되어 서로 다른 영역도 쉽게 합쳐진다. 그 결과 최종 영역의 개수가 적어지는 undersegmentation이 발생하기 쉽다. 반대로 \(Q\)가 클수록 합병 조건이 엄격해져 평균 밝기 값이 조금만 달라도 합병되지 않으므로 많은 작은 영역이 생성되는 oversegmentation이 나타난다. 따라서 실제 응용에서는 영상의 특성과 원하는 분할 수준에 맞추어 적절한 \(Q\) 값을 선택하는 것이 중요하다.

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