ガウシアンフィルタの並列処理実装の性能比較(リベンジ)

こんにてぃーは。前回ポストした記事にて、ガウシアンフィルタをSSEとOpenMPで並列処理してみたところ、SSEの効果が思ったほど良くなかったり、OpenCVのガウシアンフィルタ GaussianBlur() に比べて小カーネルでの性能がよろしくありませんでした。今日はそのリベンジであります。

続きを読む

ガウシアンフィルタの並列処理実装の性能比較

こんにてぃは。今日は、画像処理における超基本ツールの一つと言えるガウシアンフィルタについて、C++でSSE(SIMD演算)とOpenMP(マルチコア演算)を使ったときの性能を比較してみました。

(【2012/04/23追記】コード内に無駄なキャストが残っていたので除去しました。)

続きを読む

Google Code Jam 2012 - Qualification Round

今年もGCJの季節がやってきました。年々参加者が増えて競争が激化しつつありますね。自分はGCJはこれで3回目なのですが、毎度多倍長演算とかに悩まされてるもんで、ここ最近はPythonの練習がてらに参加しています。研究でも細々したプログラムはPythonで作ってるのですが、この活用範囲を広げるための練習という感じです。まだまだPythonの文法や関数を調べつつやってるのでコーディングが激遅ですが、少しずつ慣れつつあります。今回は自宅でゆっくりと挑戦しましてジャスト300位でした。

A. Speaking in Tongues

各アルファベットを別のアルファベットに置き換えるような暗号があって、その暗号文が与えられる。それを復号して表示しなさい、な問題。肝心の変換ルールが記載されてないけど、それはサンプルから自分で解析しなくちゃいけない。サンプルに乗ってない‘q’と‘z’は、文中にコソッとヒントがあるのでそれを見逃さないこと。実はこれでWAくらった(笑)。親切にもPythonに文字置換関数があったのでそれ使って終了。・・・Smallしかない問題ってGCJ初ですね。

B. Dancing With the Googlers

3人の審査員がダンサー達をジャッジしスコアをつける。審査基準があるので各審査員のスコアはだいたい1点差内に収まるが、たまに2点差があってそのとき観客は驚く。各ダンサーについて各審査員が何点をつけたか覚えてないけど、各ダンサーの合計得点と驚愕回数だけは覚えてるので、それらの参考に、審査員の誰かが×点以上つけたダンサーの数を推測し、その考えられる最大人数を出力。(説明むずい)
1点差スコアの場合、合計得点から審査員のつけたスコアの組み合わせは一意に決まる。この時点で審査員の誰かが×点以上つけてるダンサーをまずカウントしておく。残るダンサーのうち、もし2点差だったら最大点が×点以上に達するダンサーの数をカウントしておいて、驚愕回数と比べて、小さいほうを先ほどのカウントに足したのが答え。

C. Recycled Numbers

問題説明メンドイので解き方だけ(爆)。[A,B]範囲内の各数 n について、n

D. Hall of Mirror

幾何。怖いねぇ(笑)。案の定、ほとんどの人がスルーしてるし。でもこれを解いたら一気にヒーローだと思って頑張ることに(笑)。よくよく読んでみたら、Smallでは部屋の四方にしか鏡がないという制約らしくて、Xさんを反転コピーみたいに複製して、その鏡の中のXさんの位置をリストアップし、鏡の中のXさんとオリジナルXさんとの各距離を計算し、D以下なら到達不可とする、でイケると発見。Largeはどうにも実装できる気がしなかったので放置。そして一発AC!。
コンテスト後に@先輩の実装見たら、まー綺麗なこと。鏡の中のXさんの位置は各グリッド中心に限定されるはずなので、オリジナルXさんの座標位置から半径Dの円内にある各グリッド中心に向かって光線=レイを飛ばして、その光線をトレース(鏡面反射のシミュレート)して、オリジナルXさんに戻って来れればOK!、っていうのを実装すれば良かったらしい。当初思ってたよりは単純なアルゴリズムだったけど、そのレイトレースをバグなしで実装する自信は僕にはありません。後ほど修業します。

むすび

とりあえず無事に通過できて良かった(年々レベルアップしてる印象があるので)。次は実装時間がマヂでシビアになるけど、Pythonで行くべきかどうか悩み中。PythonのほうがC++よりデバッグが速いんだけど、不慣れって怖いしねぇ。

「言語技術」が日本のサッカーを変える

「言語技術」が日本のサッカーを変える (光文社新書)

「言語技術」が日本のサッカーを変える (光文社新書)

久々のブログ(広告ウザいね><)。今週読んだ本がめっちゃ面白かったので紹介する。日本サッカー協会の強化委員長の方が書かれた本で、日本のサッカーをより強くするためには言語技術のトレーニングが大事!と説いた本。正直言うとこの本読んで日本サッカー協会を見直した(笑)。↓の記事が概略つかむのに良さげ。

この本のいう『言語技術』ってのはいわゆる論理的な情報伝達のことで、このスキルが今後の日本サッカーの成長に欠かせない!と著者は主張する。論理的な情報伝達というのは、英文ライティングとかでよく出てくるPREP法のように、「最初にまず主張を述べて、続けてその理由・根拠を述べる」という発言スタイルや、5W1Hを常にクリアにするといったこと。まるで研究発表とかTOEFLのハウツー本みたいなトピックなのに、それがサッカーと結びついてるのが興味深く、購入後あっという間に読了してしまった。
サッカーの強化がテーマであるはずが「あーこれ、まさに研究と同じだよね」っていう場面が多々あって、全然異分野なのになぜか親近感が湧いてしまう。例えば「学生が質疑応答中に固まる(with 先生に助けを乞う眼差し?)」とか「論文やプレゼンの構成が起承転結型でイライラ」とか「厳しい質問するとふてくされる」とか。それと似たような事が日本のサッカー界でも起こってると。一方で欧米では、小さい頃からこっち方面の教育が徹底されており、選手が小さいうちから自分の意見を持ってるし、しかもそれを監督であろうがぶつけてくると。これが相互理解の訓練となり『creativity』を発揮できる土壌につながってるんじゃないかと。確かに。
また2012年に読むからこそ感じとれることもあった。執筆時期はオシム監督時代のようで、その後の南アフリカW杯の結果や日本人選手の欧州での活躍を知っていることで、より一層楽しめたと思う。また本の中で何度も言及されるJヴィレッジについても今読むと複雑な心境だ。Jヴィレッジは福島県にあるスポーツ施設で、この本の思想を実現したサッカースクールもそこにあるそうだ。福島県?と気になって調べてみたらこんなの出てきた。

実はJヴィレッジは原発20km圏内にある東京電力のサポートでできた施設だそうだ。今では原発作業員の宿舎として使われてるらしい。原発事故によって夢を邪魔されたサッカー少年はさぞ悔しかったろうと想像する。
代表選観戦する程度のサッカーファンのオレが楽しめたので、サッカーに興味ある理系大学院生や先生方はもっと面白く読めると思う。オススメ。

Google Code Jam Japan 2011 - Final Round

GCJの日本大会でごわす。200位以上でGoogleTシャツがもらえるとかでマジ狙いしてたのですが、残念ながら211位で無理でした。終盤、Small一つを狙う作戦に切り替えたのですが、切り替えが遅すぎてデバッグが間に合わず(笑)。最近の若い子たちのスピードに、おっちゃんついていけてません。

以下各問題について感じたことなど。

続きを読む

Color Coherence Vector の実装

【追記:2013/03/25】αとβの求め方が論文と一致していないというご指摘を受け、該当箇所(//alphaと//betaの行)を修正いたしました。ご迷惑おかけしました。
Color Coherence Vector (CCV) というのは色に基づく画像検索の一手法です。結構昔の文献ながら最近の論文でもよく参照されてる有名どころの手法のようです。処理がシンプルなのでサクっと作りました。OpenCV2.1と【id:wosugi:20110628】で紹介したラベリング処理を使っています。

//ccv.hpp
#pragma once
#include <iostream>
#include <vector>
#include <cmath>
#include <C:/OpenCV2.1/include/opencv/cv.h>
#include "labeling.hpp" //前に公開したラベリング処理のヘッダ

using namespace std;
using namespace cv;

typedef vector<pair<int,int> > CCV;

const int NUM_CCV_COLOR=4*4*4;

const int W=240;
const int H=160;
Mat bufimg(H,W,CV_8UC3);

inline void setAt(Mat& img,int x,int y,uchar r,uchar g,uchar b)
{
    uchar* p=(img.data+y*img.step+3*x);
    p[0]=b;
    p[1]=g;
    p[2]=r;
}

//カラー画像からのCCV特徴の抽出
CCV extract(const Mat& image,int tau=25)
{
    resize(image,bufimg,Size(W,H));
    blur(bufimg,bufimg,Size(3,3)); //GaussianBlur(bufimg,bufimg,Size(3,3),0);

    const int MASK=0xC0;
    for(int j=0;j<H;j++)
    for(int i=0;i<W;i++)
    {
        uchar r,g,b;
        getAt(bufimg,i,j,r,g,b);
        setAt(bufimg,i,j,r&MASK,g&MASK,b&MASK);
    }

    static Mat_<int> label(H,W);
    int regions=doLabeling(bufimg,label);

    vector<int> count(regions,0);
    vector<uchar> r(regions,0);
    vector<uchar> g(regions,0);
    vector<uchar> b(regions,0);
    for(int j=0;j<H;j++)
    for(int i=0;i<W;i++)
    {
        int p=label(j,i);
        count[p]++;
        getAt(bufimg,i,j,r[p],g[p],b[p]);
    }

    CCV ccv(NUM_CCV_COLOR,pair<int,int>(0,0));
    for(int k=0;k<regions;k++)
    {
        int c=((r[k]>>6)<<4)+((g[k]>>6)<<2)+((b[k]>>6)<<0);
        if(tau<=count[k])
            ccv[c].first+=count[k];//alpha
        else
            ccv[c].second+=count[k];//beta
    }
    return ccv;
}

//CCV特徴ベクトル間の距離関数
double metric(const CCV& x,const CCV& y)
{
    double d=0.0;
    for(int k=0;k<NUM_CCV_COLOR;k++)
    {
        double a1=x[k].first;
        double b1=x[k].second;
        double a2=y[k].first;
        double b2=y[k].second;

        double da=fabs((a2-a1)/(a2+a1+1.0));
        double db=fabs((b2-b1)/(b2+b1+1.0));
        d+=da+db;
    }
    return d;
}

CCVは一枚のカラー画像から抽出される特徴ベクトルです。基本的には色の出現頻度を表すのですが、さらに色領域の面積で場合分けすることでより細かく特徴抽出しています。具体的には、まず画像サイズを適当に正規化後、色を64色(RGB空間を4x4x4ビンで表現)に量子化し同色な領域をまとめます。ここで色数分(つまり64)の要素を持つ配列αとβを用意します。各領域について面積が閾値τ以上であればα[c]+1、そうでなければβ[c]+1します【修正:2013/03/25】閾値τ以上であればα[c]に、そうでなければβ[c]に、その領域面積を加算します(論文では画素数をカウントしているので以前の実装は間違いです)。ただしcはその領域の色を示すインデックス(c∈{0,1,...,63})です。このαとβを合わせたものがCCVになります。
論文の記述との違いですが、正規化後のサイズは240x160=36800pixにしました(論文では38976pixとだけ書いてある)。τのデフォルト値は論文にある25にしてますが、画像サイズと相互依存なので注意が必要。また平滑化は、論文と同じ3x3の平均値フィルタを使ってますが、こちらの手元にあるデータでは3x3ガウスフィルタのほうが精度が高かったです(コメントアウトしてます)。
バグなどあったら教えてください。

参考文献

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

画像処理の基礎的な処理の一つにラベリング処理というのがあります。これは連結する画素・領域をくっつけていって、最終的に同じ性質を持つ領域ごとに番号(ラベル)をふるというものです。このラベリングには数多くのやり方があり、突き詰めていけばワンパス(画像を一回だけラスタスキャンする)で処理できることが知られています。しかしネット上には意外とワンパスのコードが少なく、ラスタスキャンではなく一部幅優先探索するものや、あるいはワンパスっぽいけどポインタ駆使しすぎ&コード量が多いものが多々あり*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++)
    {
        //隣接画素(4近傍)との連結チェック
        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必要ですが、メモリ効率を考慮して空気を読んで使ってください。あと8近傍にしたい方のために念のためそのコードをコメントアウトで残しています。バグなどがあったら教えてください。

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