Graph-based Segmentation

Graph-Based Image Segmentation은 Felzenszwalb와 Huttenlocher가 2004년에 제안한 영상 분할 알고리즘이다. 영상을 그래프로 표현하여, 각 픽셀을 vertex, 인접한 픽셀 사이의 연결을 edge로 간주한다. 간선의 가중치는 두 픽셀의 색상 차이(또는 밝기 차이)로 정의되며, 색상이 비슷할수록 작은 값을 갖는다.

알고리즘은 먼저 모든 간선을 가중치의 오름차순으로 정렬한 후, Kruskal의 MST 알고리즘과 유사한 방식으로 영역을 병합한다. 처음에는 각 픽셀이 하나의 독립된 영역이며, 가중치가 작은 간선부터 차례로 검사한다. 서로 다른 두 영역을 연결하는 간선의 가중치가 두 영역의 내부 차이를 나타내는 threshold 보다 작으면 두 영역을 하나로 합친다. 이 임계값은

$$T(R)=\operatorname{Int}(R)+\frac{c}{|R|}$$로 정의된다. 여기서 $\operatorname{Int}(R)$은 영역 내부에서 선택된 최대 간선 가중치(내부 변화량), $|R|$은 영역의 픽셀 수, $c$는 사용자가 지정하는 상수이다. 작은 영역은 $\tfrac{c}{|R|}$ 항이 커서 쉽게 병합되고, 큰 영역은 더 엄격한 조건에서만 병합되므로 노이즈에 강하면서도 의미 있는 경계를 잘 유지할 수 있다.

모든 간선을 처리한 뒤에는 사용자가 지정한 최소 크기보다 작은 영역을 인접 영역과 추가로 병합하여 작은 잡음 영역을 제거한다. 이 과정은 Union-Find(Disjoint Set) 자료구조를 이용하므로 영역의 병합과 검색을 매우 효율적으로 수행할 수 있으며, 전체 시간 복잡도는 간선 정렬을 포함하여

$$O(E\log E)$$이다. 4- 또는 8-연결 격자 영상에서는 edge 수 $E$가 픽셀 수 $N$에 비례하므로 전체 복잡도는 사실상 $O(N\log N)$에 해당한다.

이 코드에서는 RGB 색 공간에서 두 픽셀의 유클리드 거리

$$w=\sqrt{(R_1-R_2)^2+(G_1-G_2)^2+(B_1-B_2)^2}$$를 간선의 가중치로 사용하고, 분할 전에 $\sigma=0.5$의 Gaussian smoothing을 적용하여 잡음의 영향을 줄인 후 그래프 분할을 수행한다. 마지막으로 최소 크기보다 작은 영역을 병합하고, 각 영역을 대표(root) 픽셀의 색으로 채워 분할 결과를 생성한다.

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