Geometry & Recognition

Ref: Efficient Graph-Based Image Segmentation, Pedro F. Felzenszwalb and Daniel P. Huttenlocher, International Journal of Computer Vision, 59(2) September 2004.

struct edge {
    float w;	//weight=color_distance
    int a, b;   //vertices=pixel positions;
    edge() {}
    edge(int _a, int _b, float _w): a(_a), b(_b), w(_w) {}
};
bool operator<(const edge &a, const edge &b) {
  return a.w < b.w;
}
#define SQR(x) ((x)*(x)) 
#define DIFF(p,q) \
    (sqrt(SQR(red[(p)]-red[(q)])+SQR(green[(p)]-green[(q)])+SQR(blue[(p)]-blue[(q)])))
int segment_image(CRaster& raster, float c, int min_size, CRaster& out) {
    if (raster.GetBPP() != 24) return 0;
    const CSize sz = raster.GetSize();
    const int width = sz.cx, height = sz.cy;
    std::vector<float> red(width * height);
    std::vector<float> green(width * height);
    std::vector<float> blue(width * height);
    for (int y = 0, curr = 0; y <height; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < width; x++, curr++) 
            blue[curr] = *p++, green[curr] = *p++, red[curr] = *p++;
    }
    // gaussian smoothing;
    const float sigma = 0.5F;
    smooth(blue,  width, height, sigma);
    smooth(green, width, height, sigma);
    smooth(red,   width, height, sigma);

    std::vector<edge> edges;
    edges.reserve(width * height * 4);
    for (int y = 0, curr = 0; y < height; y++) {
        for (int x = 0; x < width; x++, curr++) {
            if (x < width-1)
                edges.push_back(edge(curr, curr+1, DIFF(curr, curr+1)));
            if (y < height-1) 
                edges.push_back(edge(curr, curr+width, DIFF(curr, curr+width)));
            if ((x < width-1) && (y < height-1))
                edges.push_back(edge(curr, curr+width+1, DIFF(curr, curr+width+1)));
            if ((x < width-1) && (y > 0)) 
                edges.push_back(edge(curr, curr-width+1, DIFF(curr, curr-width+1)));
        }
    }

    if (c < 0) c = 1500; 
    if (min_size < 0) min_size = 10;
    universe *u = segment_graph(width*height, edges, c);
    // join small size regions;
    for (int i = edges.size(); i--> 0;) {
        edge &e = edges[i];
        int a = u->find(e.a), b = u->find(e.b);
        if ((a != b) && ((u->size(a) < min_size) || (u->size(b) < min_size)))
            u->join(a, b);
    }
    int num_rgns = u->count();

    // paint segmented rgns with root pixel color;
    out = raster;
    for (int y = height-1; y-->0;) 
        for (int x = width-1; x-->0;) {
            int a = u->find(y * width + x);
            out.SetPixel0(x, y, RGB(blue[a],green[a],red[a]));
        } 
    delete u;
    return num_rgns;
}
// Segment a graph
// c: constant for treshold function.
universe *segment_graph(int num_vertices, std::vector<edge>& edges, float c) { 
    // sort by weight
    std::sort(edges.begin(), edges.end());
    // disjoint-set
    universe *u = new universe(num_vertices);
    // init thresholds = {c};
    std::vector<float> threshold(num_vertices, c);

    for (int i = 0; i < edges.size(); i++) {
        edge &e = edges[i];
        int a = u->find(e.a), b = u->find(e.b);
        if (a != b)
            if ((e.w <= threshold[a]) && (e.w <= threshold[b])) {
                u->join(a, b);
                a = u->find(a);
                threshold[a] = e.w + c / u->size(a);
            }
    }
    return u;
};

https://kipl.tistory.com/460

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

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

를 넘지 않아야 한다. $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;
}

이미지를 균일한 컬러 영역으로 분할하는 것은 압축이나 segmentation 등에 매우 유용하다. 4개의 child node 가운데 parent node와 차이가 많이 나지 않는 것은  부분적으로 parent node에 merge 시키는 방법도 있다.

#include "QuadTree.h"
#include <queue>
#include <math.h>
QuadNode::QuadNode() {
    for (int i = 0; i < 4; i++) children[i] = NULL;
    rect = CRect(0, 0, 0, 0);
    hasChildren = false;
    mean = -1;
    stdev = 0;
}
QuadNode::QuadNode(CRect& rc) {
    for (int i = 0; i < 4; i++) children[i] = NULL;
    rect = rc;
    hasChildren = false;
    mean = -1;
    stdev = 0;
}
void QuadNode::Subdivide() {
    hasChildren = true;
    children[0] = new QuadNode(CRect(rect.TopLeft(), rect.CenterPoint()));
    children[1] = new QuadNode(CRect(rect.left, rect.CenterPoint().y, rect.CenterPoint().x, rect.bottom));
    children[2] = new QuadNode(CRect(rect.CenterPoint(), rect.BottomRight()));
    children[3] = new QuadNode(CRect(rect.CenterPoint().x, rect.top, rect.right, rect.CenterPoint().y));
}
// quadtree에 할당된 메모리를 정리함; root node는 별도로 정리.
void QuadNode::Clear() {
    if (this->hasChildren)
        for (int i = 0; i < 4; i++)
            if (children[i]) {
                children[i]->Clear();
                delete children[i];
            }
}
QuadTree::QuadTree(CRect rc, CRaster& image) {
    rect = rc & image.GetRect(); // rect는 이미지 안에.
    raster = image;
}
void QuadTree::Build() {
    root = new QuadNode();
    root->rect = this->rect;
    std::queue<QuadNode *> TreeQueue;
    TreeQueue.push(root);
    while (!TreeQueue.empty()) {
        QuadNode* node = TreeQueue.front();
        TreeQueue.pop();
        CRect rc = node->rect;
        if ((rc.Width() > 2 && rc.Height() > 2) && !IsHomogenous(node)) {
            node->Subdivide();
            for (int i = 0; i < 4; i++) 
                TreeQueue.push(node->children[i]);
        }
    }
}
// subdivision 조건: 편차, 최대-최소 차이등의 criteria를 사용할 수 있다.
#define THRESH 15
bool QuadTree::IsHomogenous(QuadNode* node) {
    CRect rc = node->rect;
    double s = 0, ss = 0;
    int minval = 255, maxval = 0;
    for (int y = rc.top; y < rc.bottom; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y) + rc.left;
        for (int x = rc.left; x < rc.right; x++) {
            int a = *p++; 
            s += a; ss += a * a;
            maxval = max(maxval, a); minval = min(minval, a);
        }
    }
    int n = rc.Width() * rc.Height();
    node->mean = s / n;
    node->stdev = sqrt((ss - s * s / n) / n);
    return (node->stdev >= THRESH) ? false: true;
    //return (maxval - minval) > THRESH ? false: true;
};
// leaf node를 찾아 list에 넣음;
void QuadTree::FindLeaf(QuadNode* node) {
    if (node == NULL) return;
    if (!node->hasChildren) {
        LeafList.push_back(node);
        return;
    }
    for (int i = 0; i < 4; i++) 
        FindLeaf(node->children[i]);
}

사용 예:

void test_main(CRaster& raster, CRaster& out) {
    ASSERT(raster.GetBPP() == 8);  // gray input;
    QuadTree qtree(raster.GetRect(), raster);
    qtree.Build();
    qtree.FindLeaf(qtree.root);
    out.SetDimensions(raster.GetSize(), 24); //RGB-color;
    // segment image;
    for (int i = 0; i < qtree.LeafList.size(); i++) {
        QuadNode *node = qtree.LeafList[i];
        int a = int(node->mean + .5);
        FillRegion(out, node->rect, RGB(a, a, a));  //gray average;
        // draw rectangle with red color;
        DrawRectangle(out, node->rect, RGB(0, 0, 0xFF));
    }
    qtree.root->Clear();
    delete qtree.root;
}