ほぼワンパスなラベリング処理

プログラミング・プロコン 画像処理・パターン認識

画像処理の基礎的な処理の一つにラベリング処理というのがあります。これは連結する画素・領域をくっつけていって、最終的に同じ性質を持つ領域ごとに番号（ラベル）をふるというものです。このラベリングには数多くのやり方があり、突き詰めていけばワンパス（画像を一回だけラスタスキャンする）で処理できることが知られています。しかしネット上には意外とワンパスのコードが少なく、ラスタスキャンではなく一部幅優先探索するものや、あるいはワンパスっぽいけどポインタ駆使しすぎ＆コード量が多いものが多々あり*1、処理速度や読みやすさの点でもっと改善できるんじゃないかなぁ〜と感じました。
今回はできるだけ簡単な記述で、ほぼワンパスのラベリングのコードを書きましたので公開します。詳細は後回しして、とりあえずコードです。

//labeling.hpp
#pragma once
#include <vector>
#include <C:/OpenCV2.1/include/opencv/cv.h>

using namespace std;
using namespace cv;

//ラベルのバッファ長
//予想領域数より十分に大きく
const int BUF_LABEL=512;

inline bool isIn(int w,int h,int x,int y)
{
    return 0<=x && x<w && 0<=y && y<h;
}

inline unsigned int getAt(const Mat& img,int x,int y)
{
    unsigned char* p=(img.data+y*img.step+3*x);
    unsigned int b=p[0];
    unsigned int g=p[1];
    unsigned int r=p[2];
    return (r<<16)+(g<<8)+(b);
}

//aの属すグループの代表に向かって経路圧縮（代表を返す）
inline int compress(vector<int>& parents,int a)
{
    while(a!=parents[a])
    {
        parents[a]=parents[parents[a]];
        a=parents[a];
    }
    return a;
}

//aの属すグループとbの属すグループを併合（併合後の代表を返す）
inline int link(vector<int>& parents,int a,int b)
{
    a=compress(parents,a);
    b=compress(parents,b);
    if(a<b)
        return parents[b]=a;
    else
        return parents[a]=b;
}

//番号とびとびなラベルを0,1,2,...に貼り替え
inline int relabel(vector<int>& parents)
{
    int index=0;
    for(int k=0;k<(int)parents.size();k++)
    {
        if(k==parents[k])
            parents[k]=index++;
        else
            parents[k]=parents[parents[k]];
    }
    return index;
}

int doLabeling(const Mat& image,Mat_<int>& label)
{
    const int W=image.cols;
    const int H=image.rows;

    //関数を複数回呼び出すときメモリ確保が何度も生じるので、
    //例えばリアルタイム処理では関数外で確保するのがベター。
    vector<int> parents;
    parents.reserve(BUF_LABEL);

    int index=0;
    for(int j=0;j<H;j++)
    for(int i=0;i<W;i++)
    {
        //隣接画素（４近傍）との連結チェック
        unsigned int c=getAt(image,i,j);
        bool flagA=(isIn(W,H,i-1,j  ) && c==getAt(image,i-1,j  )); //左
        bool flagB=(isIn(W,H,i  ,j-1) && c==getAt(image,i  ,j-1)); //上
        //bool flagC=(isIn(W,H,i-1,j-1) && c==getAt(image,i-1,j-1)); //左上
        //bool flagD=(isIn(W,H,i+1,j-1) && c==getAt(image,i+1,j-1)); //右上

        //着目画素と連結画素を併合
        label(j,i)=index;
        if((flagA|flagB)==true) //if((flagA|flagB|flagC|flagD)==true)
        {
	    parents.push_back(index);
	    if(flagA) label(j,i)=link(parents,label(j,i),label(j  ,i-1));
	    if(flagB) label(j,i)=link(parents,label(j,i),label(j-1,i  ));
	    //if(flagC) label(j,i)=link(parents,label(j,i),label(j-1,i-1));
	    //if(flagD) label(j,i)=link(parents,label(j,i),label(j-1,i+1));
	    parents.pop_back();
        }
        else
	    parents.push_back(index++);
    }

    //再ラベリング
    int regions=relabel(parents);
    for(int j=0;j<H;j++)
    for(int i=0;i<W;i++)
        label(j,i)=parents[label(j,i)];
    return regions;
}

ラベリングそのものはワンパスですが、再ラベリングにさらにワンパス必要なのでこのコードは実質ツーパスです。再ラベリングが不要なら消してください。領域の統合と代表取得には「合併併合アルゴリズム（union find）」あるいは「disjoint set」で用いられている経路圧縮と呼ばれるテクニックを使いました。上記のコードが厳密なものかどうかは分かりませんが、一般にその処理オーダーはAckermann関数（爆発的に増加する関数）の逆関数だそうなので、ほぼO(WH)時間で処理できるとみなして構わないでしょう。また処理順序がラスタスキャンなのでハードウェアとの親和性も高そうです。
使用上の注意：パブリックドメインで公開しますので、適宜用途に合わせて変更して使ってもらって構いません。が、何か一言コメント頂けると嬉しいです＾＾。BUF_LABELは最悪な状況（市松模様）ではW*H必要ですが、メモリ効率を考慮して空気を読んで使ってください。あと８近傍にしたい方のために念のためそのコードをコメントアウトで残しています。バグなどがあったら教えてください。

参考文献

*1:さらに間違っているコードを公開しているのもあったり・・・