C++で2種類のコンテナを同時にソートする

C++でデータを扱っていると，2種類のコンテナを同時にソートしたい場合が時々出てくる．どういうことかというと，

int main() {
    int keys[]   = {5, 3, 1, 2, 4};
    int values[] = {1, 2, 3, 4, 5};
    aoba::sort_by_key(keys, keys+5, values);
    for(int i=0; i<5; ++i) {
        std::cout << keys[i] << " ";   // 1 2 3 4 5
    }
    std::cout << std::endl;
    for(int i=0; i<5; ++i) {
        std::cout << values[i] << " "; // 3 4 2 5 1
    }
    std::cout << std::endl;
    return 0;
}

と，片方のソート結果を元にもう片方をソートするような感じ．これはCUDAライブラリであるthrustに標準搭載されており，sort_by_keyという名前がついている．

さて，こういうときに一からソートアルゴリズムを組み直すなんてことは馬鹿らしい．じゃあどうするか．基本的には既存のstd::sortに渡すイテレータを変形することで，なんとか解決しようとする．そうなると，boostにzip_iteratorが存在することに気付けば，イテレータをboost::tupleの組で表現することによって，まとめてソートさせればいいような気がしてくる．こんな風に:

std::sort(
    boost::make_zip_iteator(
        boost::make_tuple(
            keys_first, values_first)),
    boost::make_zip_iterator(
        boost::make_tuple(
            keys_last, values_last)),
    compare_functor());

しかし，そうは問屋が降ろさない．何故なら，boost::zip_iteratorはWritableではない．これはソートアルゴリズムの要件を満たしていないため，適用させようとしてもコンパイルエラーが発生してしまう！

それじゃあどのようにして解決するかというと，一番楽な方法としてはWritableなコンセプトを満たしたzipイテレータを新しく作ってやって，それを用いてstd::sortを行うのがいい．とはいっても，C++でそれをやるのは若干手間のかかる作業になる．具体的にはこんな感じで実装した:
実装したといっても，実際には『Sorting two arrays simultaneously』の内容を若干手直しするだけなので基本的にはそんなに大変じゃないし，ただ利用するぶんには単純にsort_by_key()を呼び出すだけでよい．

しかし，なんでこんなにC++は行数が多くなってしまうのか．haskellなんか一行で済むのに．

超一様分布列をBoost.Randomを用いて実装する

ご存知の通り、C++のライブラリには高性能な擬似乱数生成エンジンBoost.Randomが存在している。Boost.Randomは擬似乱数を生成するエンジン部分と、一様分布、ガウス分布、ベルーヌーイ分布など好みの分布を生成する出力部分の2つで構成されている。

しかし、ことモンテカルロ法においては主に収束速度を早める目的として、擬似乱数ではなく低食い違い量列(Low-Discrepancy Sequence, LDS)と呼ばれる数列が用いられる場合がある。これを用いると低次元での積分計算の際に収束速度を大きく改善することができ、具体的には従来のモンテカルロ法で(1/√N)だった収束速度が約(1/N)くらいにまで落とすことができる。

さて、このような低食い違い量列をBoost.Randomで用いるためにはエンジン部分を製作すればいいだけで、実際には以下のような形で実装すれば良い。今回は簡単なHalton列を用いて実装を行った:

具体的には、対象のファンクタがUniform Random Number Generatorのコンセプトを満たすためには以下の記述が必要となる:

• X::result_typeにoperator()が返す型を指定
• operator()にエンジン部分を記述
• min(), max()にエンジンが生成する分布の最小値、最大値を指定

さらに、対象のファンクタが上述の内容に加えPseudo-Random Number Generatorのコンセプトを満たすためには:

• X()にデフォルトの状態で生成されるようなコンストラクタを指定
• X(...)にユーザ指定の状態で生成されるコンストラクタを指定
• seed(...)にエンジンの状態を変更するためのメンバ関数を指定

するだけで良い。後はBoost側が良きに計らってくれる。非常に便利だと思う。

余談
C++0xには同等のstd::randomが付属しているため、これを利用すればいいのではという声もあるとは思うが、std::uniform_real_distribution<>にこのLDSを食わせたところ、どうみてもLDSではない値が出力される結果となった。恐らくresult_typeがdoubleという一般の擬似乱数ではない型を指定しているために起こる現象だと思うが、ソースを見ていないので一先ずは従来のBoost.Randomの手法を用いることにする。

Restricted Boltzmann Machineの学習手法についての簡単なまとめ

近年の機械学習ではDeep Learningと呼ばれる分野が一世を風靡しています．コンピュータビジョンや自然言語処理，音声認識などの分野では何らかの問題を解こうとした際に，まず対象の入力データからSIFTやケプストラムといった何らかのアルゴリズムを用いて特徴ベクトルを抽出し，ごりごりと判別していくといった流れが一般的です．しかし，その特徴ベクトルを生成するという生のデータから本質となる部分を抽出するアルゴリズム自体は研究者が一生懸命考えながら作るのが普通でした．

Deep Learningの分野で最も有名な手法の一つであるDeep Belief Nets(DBN) [Hinton06]は，研究者がアルゴリズムを作るのではなく，それ自体も機械学習にやらせましょうという動機で生まれたアルゴリズムです．DBNではまるで一昔前にやたら流行ったニューラルネットワークのように各ノードを層状に配置し，それらのパラメータをRestricted Boltzmann Machine(RBM)と呼ばれるモデルを用いて学習することによって自動的に特徴ベクトルを生成します．これが既存手法を大きく超える精度を出してしまったためにDeep Learningは一躍有名となり，ニューラルネットの再来か，機械学習は人工知能の夢を見るか，みたいに言われているような状況です．最近ではGoogleがネコのニューロンを自動的に学習したとかでニュースになっていました．

さて，そのような状況にもかかわらず，そのDeep Belief Netsやその改良版であるDeep Boltzmann Machineを学習するためには，まずその基礎となっているRBMについて学ぶ必要があります．しかしそのRBMですらお世辞にもわかり易いものとは言えず，しかもその最終的な更新式やその導出過程についてきちんと書いている文献は今まで確認したところ存在していません(*1)．「チュートリアルや論文だけでは良く分からないので，実際に手を動かして確認したい．だがその計算過程は公開されていないため，結局良く分からないままで終わる」ことが多いように思われました．

ということで，RBMの計算過程について書いた個人用のpdfファイルを公開してみることにします．外部に見せるということで多少文章を詳しくしましたが，これ単体で学ぶのではなく他の資料と比較しながら見るとベターかと思います．多分数式は間違ってないと思いますが，何かありましたらご一報頂けると幸いです．

http://dl.dropbox.com/u/2048288/RestrictedBoltzmannMachine.pdf

*1 Deeplearning.netの式は比較的詳しいですが，実際に計算したところ微妙に間違えている箇所を数ヶ所発見しています

直積量子化(Product Quantization)を用いた近似最近傍探索についての簡単な解説

"aka motsu-nabe" by chatani

概要
冬の寒さも一段と厳しくなってまいりました。おでんや鍋が恋しくなる季節です。

さて、最近ようやっと一仕事が終わりまして、長ったらしい記事が書けるようになりました。ですので、今回は2011年にTPAMIで発表された、近似最近傍探索についての論文『Product quantization for nearest neighbor search』について簡単に紹介したいと思います。

この論文は2011年に発表された、最近傍探索アルゴリズムの決定打です。シンプルな理論でありながら既存手法を打ち破るほどの強力な性能を有し、速度も非常に高速、かつ省メモリなのでスマートフォンに載せ、リアルタイムで動作させることも可能です。

以前この手法はCV勉強会@関東で紹介されたらしいのですが、具体的に紹介しているページは(最近すぎるので当たり前ですが)現在のところあまり見かけていません。ただ個人的にこの手法はあと5年くらいは本線で戦えるものと思っておりますので、きちんと解説した記事を出すことはある程度有意義なものかなぁと、そういう風に考えています。

最近傍探索の課題
最近傍探索は最も古典的な手法でありながら、現在においても頻繁に利用されているアルゴリズムの一つです。近年のプロセッサの性能向上によって大規模なデータ処理が家庭のパソコンでもできるようになり、コンピュータビジョンにおいても数億件規模(10^8~9)の高次元ベクトルデータを最近傍探索を用いて処理することで、良い精度の結果を得ることが当たり前の時代になってきました。

ただそうは言っても、最近傍探索をそのまま用いることは現実的ではなく、いくつかの課題が浮上してきます。具体的には以下のとおりです:

• いかに圧縮した状態でデータを格納するか: ベクトルデータをそのままの状態で保持することはメモリ容量の圧迫へ繋がります。例えば128要素のベクトルを考えた場合、floatを使用した場合でも32×128=4096bitを1エントリで消費してしまいます。IOアクセスが生じてしまうと速度は一気に低下してしまうため、ベクトルデータはすべてメモリ上に展開しておくことが望ましいー
  しかし、今回我々が扱いたいエントリ数というのは数億とかそういったオーダーです。そのオーダーのベクトルをそのまま格納するのは現在のメモリ容量だとなかなか厳しい。ですので、情報量を保ったままベクトルをいかに圧縮できるかが、実用的な面において重要な課題となります。
• いかに高速に解を見つけられるか: 次元数の多いベクトルの距離計算は基本的に低速です。さらにそのまま実装してしまうとすべてのエントリに対して距離の計算をする必要があるため速度はO(N)となり、実用的ではありません。1エントリあたりの計算コストを削減したり、そもそも明らかに異なる(クエリより離れすぎている)エントリを計算しないことで、計算を高速に行うことも同じく重要となります。

これらの課題を克服するため、近年では近似最近傍探索と呼ばれる手法が盛んに研究されております。その中でも有名なアルゴリズムはLSH(Locality-Sensitive Hashing)であり、ハッシュ関数によって抽出される近傍エントリのみを計算することで、高速に最近傍探索を行います。ただこの手法は計算用のハッシュを別に格納する構造であるため、全体のサイズは減少するのではなく増大してしまいます。

また、Preferred Researchで紹介されているMinHashは、ベクトルの要素数が定まっておらず、成分が{0,1}の2値しかとらないような特殊なベクトルを計算する際においては有効な手法でありますが、そうでない場合においてはいったんベクトルをバイナリ表現に落とし込まなければならないので一般的に精度は低下してしまいます。今回紹介する直積量子化は、対象のベクトルが高次元であり、かつユークリッド距離による最近傍探索を行いたい場合において威力を発揮する手法です。

スカラ量子化
直積量子化について説明を行う前に、まず量子化とは何かについて説明していきます。物理をかじったことのある方ですと『量子化 = 物理量を演算子へ変換する』ものと思っちゃいますけど、そっちの量子化ではなくある数を特定のビット列で表現することの量子化です。

分割された領域に、特定のビットを割り当てる

例えば、0〜1の間をとる数を2ビットで表現しろという問題を考えた場合、普通は上図のように0〜1を4分割し、それぞれの数値を割り当てることで表現を行います。このように、少ないビット数で対象の数値を表現することを一般に『スカラ量子化(Scalar Quantization)』といいます。スカラ量子化によってfloat型(32bit)の数値は、精度を犠牲にすることによって僅か2bitで表現することができるようになりました。

ベクトル量子化
ところで、このスカラ量子化は対象のデータが数値であった場合には効率的ですが、ベクトルの場合だと効率的とは言えません。

スカラ量子化では、ベクトルを効率良く表現することができない

例えば上の図について考えてみましょう。この点群を1エントリ辺り1ビットで表現しようとした場合、スカラ量子化によって各成分を量子化してしまっては最低でも2ビットが必要となりますが、物事を空間的に考えることで発想を転換します。つまり、それぞれのクラスタを代表する点を予めコードブックとして保持しておき、各点は一番近い代表点のインデックスを保持しておくのです。これによってそれぞれの点は僅か1ビットで表現することができるようになりました。このような手法のことを一般に『ベクトル量子化(Vector Quantization)』といいます(*1)。

ベクトル量子化によって、ベクトルを2つの代表点によって表せるようになった

*1: 本手法のクラスタリングにはk-means algorithmを使用しています。k-meansに関しては『クラスタリングの定番アルゴリズム「K-means法」をビジュアライズしてみた』が視覚的に分り易いのでそちらを参照してください。

直積量子化
ベクトル量子化はスカラ量子化よりも少ないビット数で効率良くベクトルを表現できるのですが、同時に欠点もあります。次元が増大するにつれて使用するコードブックが指数的に増えてしまうのです。

これは次元の呪いに起因するものです。数次元〜数十次元くらいだったらそんなに問題とはならないんですけど、数百とか数千次元になってくるととんでもない量のコードブックが必要となり現実的でない。『直積量子化(Product Quantization)』はこの次元の呪いを解決するために用いられる、スカラ量子化とベクトル量子化の中間を行く手法です。

高次元のベクトルを直積表現によって分解する

具体的な手法は下記の通りです。まず、N次元のベクトルをM分割し、N/M次元のサブベクトルを生成します。そして、そのサブベクトルをそれぞれベクトル量子化することによって、インデックスを(i1, i2, ..., iM)のタプルで表現します。そうすれば次元の呪いは発生せず、データ表現も比較的低ビットで抑えられます。

例として128次元のベクトルを考えてみましょう。まず、このベクトルを16分割することで8次元のサブベクトルへ落とし込みます。そして、そのサブベクトルをベクトル量子化によって8bit表現に持ってきた場合、対象のベクトルは8×16=128bitで表現でき、かつ使用するコードブックサイズも16×256=4096と少数で済みます。float型のサイズで考えると30倍以上の圧縮効率です。このように、少数のコードブックサイズで効率よくベクトルを表現できることが直積量子化の利点です。

距離計算
さて、直積量子化によって各ベクトルのエントリを圧縮した状態で格納したとしましょう。あるクエリベクトルが与えられたとして、そのユークリッド距離が最小となるエントリを見つけるにはどうすればいいでしょうか？

幸運なことに、ユークリッド距離の計算は各要素毎に独立しています。つまり、ただ単にクエリベクトルを各サブベクトルに分け、それぞれの領域で距離を計算し、その和を計算するだけで良いのです。よって、直積量子化を用いて最近傍探索を近似的に行うためには:

1. クエリベクトルをM分割し、N/M次元のサブベクトルに分解する
2. クエリのサブベクトルと、それぞれのコードブックに記された代表ベクトル間の距離を予め計算しておく
3. それぞれのエントリがどのインデックスに位置しているのか確認し、2で計算した距離を参照して足し算を行う
4. 最小距離をもつエントリを候補として提出する

ことで実現できます。これが直積量子化を用いた、近似最近傍探索の核となるアルゴリズムです。

詳細について
この距離計算は非対称性からある程度のバイアスがかかっているため、実際の距離を正確に見積もるためにはこのバイアスを取り除く必要があったりします。また、具体的な距離計算がコードブックだけで済むので省力化されているとはいえ、計算量は依然としてO(N)であり、そこらへんも大規模なデータ数では足を引っ張る原因となります。

んじゃあこれをどうやって解決するのかについては、現論文

http://lear.inrialpes.fr/pubs/2011/JDS11/jegou_searching_with_quantization.pdf

を参照してください。上記2つの懸念についての解決案が記述されており、上の知識があれば多分すらすらと読めるはずです。
本当はこれも解説しなきゃなんないんですけど、書いている内に疲れたちゃったのでとりあえずコアなとこだけということで勘弁してください。気が向いたら更新します。

OpenCVアルゴリズムを高速に記述できる補助ライブラリcvutilを公開します

"The Dunny Collective - Brisbane River" by KayVee.INC

概要
OpenCVを用いてアルゴリズムを組む際、ある画像からある画像へ変換処理を行ったり、関数fを用いて新規に画像を生成したり、画像の結果を1つの値へと集約するなどの定型的な処理は非常に多く行われておりますが、並列処理を行いかつ速度を考慮したプログラムを組もうとするとある程度のコード量を必要とするためそれなりの手間となります。また画像処理によくあるループ文のネストが大量に発生するため、普通に組んでしまうと一見何をやっているのか非常に分かりづらいコードとなってしまいます。

cvutilはこのような冗長な作業を定型化し、複数のアルゴリズム関数とC++0xに新しく搭載されたラムダ式を使用して高速にアルゴリズムを記述することができます。このようにして作られたアルゴリズムは自動的にOpenMPによって並列化されるため、研究用途などの速度をある程度担保しつつもメンテナンスのし易いコードを記述するための手助けとなります。

インストール
cvutilはgithubにて公開しております。ダウンロード先は下記のURLを参照してください。

https://github.com/rezoo/cvutil

cvutilはヘッダオンリーのライブラリです。バージョンがOpenCV 2.x系統であれば、対象のincludeディレクトリにcvutilをコピーするだけで使用できます。
C++0xのラムダ式を用いることで端的に記述することができますが、関数オブジェクトを使用すればある程度古いコンパイラでも正常に動作するでしょう(そこまでして使う利点があるのかどうかは分かりませんが)。

ちなみに、元々は自分用に作られた適当なライブラリなので、バージョンアップに伴ってある程度破壊的な変更を行うこともありますのでご了承ください。その際はバージョンを0.1上げることによって対応します。

アルゴリズム一覧
ここではcvutilに搭載されているアルゴリズム一覧の一部を紹介いたします。これらのアルゴリズムはSTLのそれと似ているため、C++を触ったことのある方は名前から直感的に機能を類推できるでしょう。

transform_image
transform_image関数は対象の画像を、任意の関数を通して出力先へ変換します。cvutilの中では最も汎用性の高い関数となります。

cvutil::transform_image(src, dst, [](uchar x) {
    return cv::saturate_cast<uchar>(x*2); });

この関数は複数の画像を高度に合成するような場合においても利用できます。アリティ(渡される引数の数)によって、最も適したアルゴリズムが自動的に呼び出されます。

cvutil::transform_image(src1, src2, dst, [](uchar x, uchar y) {
    return (x + y)/2; });

generate_image
generate_image関数は任意の関数を用いて、出力先に画像を新しく生成します。generate_imageをそのまま用いることは通常ありませんが、x, y座標を元に画像を生成する関数generate_image_xyや、画像中央に原点がある右手系の座標系u, v座標を元とした関数generate_image_uvは比較的使い勝手が良いと思います。

int width = 512;
int height = 512;
cv::Mat_<cv::Vec3b> dst2(height, width);
cvutil::generate_image_xy(dst2, [=](int x, int y) {
    return cv::Vec3b(0, x*256/width, y*256/height); }); // BGR

reduce_image
reduce_image関数は対象の画像を、任意の関数を通して1つの結果へと集約し、その値を返します。

cv::Mat_<uchar> src3(3, 3);
src3(0, 0) = 0; src3(0, 1) = 1; src3(0, 2) = 2;
src3(1, 0) = 3; src3(1, 1) = 4; src3(1, 2) = 5;
src3(2, 0) = 6; src3(2, 1) = 7; src3(2, 2) = 8;
int result = cvutil::reduce_image(src3, (int)0, std::plus<int>()); // 36

上のサンプルは全要素に対しての足し算を行っています。transfromの並列化は兎も角、reduceの並列化は一般的に面倒なのですが、reduce_imageはそこらへんの処理もきちんと並列化してくれます。

transform_reduce_image
transform_reduce関数は1枚または複数の画像を対象の関数で1つの値へ変換した後に、reduce関数を適用します。2種類の結果を合併した集約結果を返す場合などに役立つと思います。

minmax_element_image
minmax_element_image関数は画像の最小値と最大値を返します。

std::pair<uchar, uchar> minmax_result =
    cvutil::minmax_element_image(src); // (min, max)

integrate_image
integrate_image関数は画像を2次元に積分します。具体的には、inclusive_scan関数を2次元的に適用します。バイナリー関数にstd::plusを適用した場合、この関数は積分画像の生成と等価です。また、画像の型が対象の二項演算子に関してアーベル群を成していた場合、任意の矩形領域に関してのreduceはintegrate_imageによって生成された画像を用いて、定数時間による計算が可能となります(参考: 並列環境におけるreduceとscanアルゴリズム)。

OpenCVにも積分画像を生成する関数は存在しているのですが、この関数は並列化を行いかつ任意の二項演算子を適用できる点が異なります。

cv::Mat_<uchar> src3(3, 3), dst3(3, 3);
src3(0, 0) = 0; src3(0, 1) = 1; src3(0, 2) = 2;
src3(1, 0) = 3; src3(1, 1) = 4; src3(1, 2) = 5;
src3(2, 0) = 6; src3(2, 1) = 7; src3(2, 2) = 8;
cvutil::integrate_image(src3, dst3, std::plus<uchar>());
// [0, 1, 2;     [0,  1,  3;
//  3, 4, 5;  =>  3,  8,  15;
//  6, 7, 8]      9,  21, 36]

cvutilはその他にも色々と関数がありますが、多分名前からどんなことやってるのか大体は理解できると思いますので、とりあえず代表的な関数をいくつか紹介いたしました。

License
自分用に作っているだけなのですが、一応ライセンスはMIT Licenseを適用します。誰でも自由に無制限に使ってかまいません。

とりあえずはそんな感じです。もし何か使っていく上で質問や疑問、あるいはご指摘等ございましたらお気軽にコメントまたはメールをお願いします。

追記
お気づきの通り、本来の文脈ならば、アルゴリズムに渡す引数はcv::Mat_<T>ではなく何かしらの抽象構造Viewを渡すべきなのです。C++のSTL Algorithmはコンテナでなくイテレータという抽象に依存することによって、任意のデータ構造による適用が可能となり、さらにtransformed_iteratorなど、変形を遅延させるような構造を取ることができる。つまり、cvutilにおいても本来はそのような『抽象に依存する』データ構造を起点にアルゴリズムを組んでいくべきなのです。

// こんな感じで書けたら最高なんだけど…
cvutil::transform_image(
    cvutil::make_resized_view(src, dst.size()),
    dst,
    any_functor);

ただこのアナロジーを画像に適用しましょうとなるとなかなか難しく、現状としてはboost.GILしか存在していない。そういった意味で確かにboost.GILは魅力的なライブラリなのですが、GILに適用するalgorithmが殆ど存在していませんし、なによりコンパイル時間がboooooooooooostしてしまう。ぱぱっと書いてぱぱっと書き捨てられるような構造がむしろ研究にとっては重要なわけで、そんな面倒なことするよりも多少汚いOpenCVを採用したほうが使い勝手としては向上するよねという、まぁそういった言い訳です。

OpenCVの画素値アクセスについて

"Fire King" by A Continuous Lean

何度も何度も繰り返されている議題ですが、今回は画像処理に良くある『全画素を変換する』という処理のパフォーマンスについて考察します。
あなたはOpenCVのcv::Mat_<T>を使って、全ての画素値をある関数funcで変換するという処理を行いたいとします。それで、たぶんなんですけど、普通に何も考えず組むとしたらこんな感じのコードになると思います:

for(int y=0; y<src.rows; ++y) {
    for(int x=0; x<src.cols; ++x) {
        dst(y, x) = func(src(y, x));
    }
}

なんの変哲もない普通のループです。このコードは普通に動きます。おわり。ちゃんちゃん。

…そういきたいのですが、パフォーマンスを気にする方々はここでいろいろとした疑問が沸き起こってくるわけです。それじゃあどこが悪いのか。1つずつ疑問点を上げ、その上で上記のコードがある程度の合理性を持っているコードであることを示していきます。

疑問1: 関数呼び出しを行っているからその分遅い？
src(y, x)という処理は結局cv::Mat<T>::operator()という関数を呼び出していることに相当します。なので、関数呼び出しを数十万のピクセルに対して行っているため、関数呼び出しにかかるコストが無視できないレベルで発生していくのではないか？？

回答: 遅くならない
OpenCVのcv::Mat_<T>::operator()はインライン展開されます。なので、関数の呼び出しに対して発生する速度的なコストは発生しません。安心して使いましょう。

疑問2: 掛け算を余計に行なっているために遅い？
結局cv::Mat_<T>::operator()の行っていることはコードに直すとこんな感じです。

src(y, x) == ((T*)(src.data + src.step.p[0]*y))[x]

でも、この項をよく見てみると、xのループに対してyの値は固定されているため、第二項の計算は正直不要でしょう。なのに掛け算の計算でコストが余分にかかってしまうため、速度が遅くなるのでは？？

回答: 遅くなるけどそんなに問題じゃない
確かに遅くなりますが、そこまで問題となる箇所ではありません。時代は3GHzに突入しているのです。整数の掛け算にかかるコストなんてたかが知れています。

そんなところよりも、別の箇所を改善したほうが大抵の場合劇的に向上します。仮に、funcの中にstd::expやstd::logが含まれていたとすると…

dst(y, x) = std::exp(-(src(y, x) - mu)*(src(y, x) - mu)); // ぐぎゃー
dst(y, x) = -src(y, x)*std::log(src(y, x)); // おもいー

…それだけで大分重くなってしまいます。sqrt, sin, cosはまだ我慢できるにしても(それでも掛け算よりは結構重い)、指数、対数関数は普通に使うとかなり重いんです。CV系統などの精度を求めないような計算でしたら、スプライン関数などの近似式に置き換えることを検討しましょう。

at<T>(), operator()は範囲チェックを行ってくれる
昔のOpenCVにはat<T>()と同じような処理をしてくれるCV_IMAGE_ELEMマクロなんて言うものもありましたが、これはもう時代遅れの代物です。なぜなら、cv::Matのat<T>(), cv::Mat_<T>のoperator()は与えられたx, y座標が範囲内のものであるかどうか、きちんとチェックを行ってくれるからです。範囲外である場合には警告を送ってくれます。そういった意味でもCV_IMAGE_ELEMよりatやoperator()を利用していくべきでしょう。

『余計な範囲チェックを行っているために遅くならないのか？』という疑問もあるかもしれませんが、この範囲チェックはReleaseモードで消去されます。なので実際の動作には全く影響されないのです。

cv::Mat::at<T>(), cv::Mat_<T>::operator()は早い
at<T>(), operator()はきちんと考えられているメンバ関数です。ポインタを用いてがりがり書いたコードよりは少し遅いかもしれませんが、この遅さは全体のアルゴリズムに影響するような遅さではありません。ほっと胸を撫で下ろし、安心して使いましょう。

でもやっぱり速度が気になる
とはいっても、きちんと書かれていないという神経質な方もいるかもしれません。at<T>()なんて負けた気がして使いたくない……その場合のコードは大体以下のようになります。

for(int y=0; y<src.rows; ++y) {
    uchar* src_ptr = src.ptr<uchar>(y);
    uchar* dst_ptr = dst.ptr<uchar>(y);
    for(int x=0; x<src.cols; ++x) {
        dst_ptr[x] = func(src_ptr[x]);
    }
}

これだけじゃあ高速化されたといっても微々たるものなので、せっかくですからOpenMPによる並列化を行ってみましょう。さらにassertでサイズが合っているかチェックを行うようにして…そうすると最終的なループはこんな感じになります。

assert(src.size == dst.size);
#pragma omp parallel for
for(int y=0; y<src.rows; ++y) {
    uchar* src_ptr = src.ptr<uchar>(y);
    uchar* dst_ptr = dst.ptr<uchar>(y);
    for(int x=0; x<src.cols; ++x) {
        dst_ptr[x] = func(src_ptr[x]);
    }
}

…大分行数が増えてしまいました。ただ自分がやりたいことは全画素に対してfuncを適用させたい、ただそれだけなのに、ここまで行数を書く必要があるのです。しかもぱっと見で何をしたいのか良く分からない。面倒ですし分かりにくい。

ここで、C++を少しかじったことのある人でしたら、ユニタリー関数を適用させる表式のアルゴリズム関数を適用すればいいんじゃね？といった考えに至るのは自然な発想でしょう。こんなふうに:

template<typename SrcType,
         typename DstType,
         typename UnaryFunction>
void transform_image(const cv::Mat_<SrcType>& src,
                     cv::Mat_<DstType>& dst,
                     UnaryFunction unary_op) {
    assert(src.size == dst.size);
    #pragma omp parallel for
    for(int y=0; y<src.rows; ++y) {
        // equivalent to ptr<T>(y)
        const SrcType* src_ptr = src[y];
        DstType* dst_ptr = dst[y];
        for(int x=0; x<src.cols; ++x) {
            dst_ptr[x] = unary_op(src_ptr[x]);
        }
    }
}

上の形の式をヘッダファイルあたりに定義してやれば、後は以下のようなラムダ式を利用して書くだけで自動的に最適なループ処理を行ってくれます:

transform_image(src, dst, [](uchar x) -> uchar {
    return cv::saturate_cast<uchar>(x*2);
});

この表式ですと、どんな変数がなんの変数に影響を及ぼすのか、全画素に対してどのような処理を行うのか一目瞭然です。また、引数に渡すunary_opはコンパイラの最適化処理によって自動的にインライン展開されるため速度的なオーバーヘッドはありませんし、OpenMPによる並列化も行ってくれます。また任意の関数オブジェクトを作用させることができるので、非常に使い回しの効く関数となっています。

画素云々について言及している方はこれまで数多く見かけましたが、このような画像に特化したアルゴリズム関数について言及している方は、何故か誰もいないというのが現状です。なのでまぁ自分用にざっくり書いているコードなのでいろいろと不備はありますけど、上記のようなアルゴリズムを詰め合わせた補助ライブラリを近いうちに上げていきたいと思います。

追記(2011/09/24): 記事のとおり、補助ライブラリであるcvutilをgithubに上げました。詳細は別記事『OpenCVアルゴリズムを高速に記述できる補助ライブラリcvutilを公開します』を参照してください。

補足
『イテレータを使って標準のアルゴリズムを適用すればいいだけなのでは？』と思う方もいるかもしれません。つまりこんな感じですね:

std::transform(src.begin(), src.end(), dst.begin(), [](uchar x) -> uchar {
    return cv::saturate_cast<uchar>(x*2);
});

確かにこの式でも動くんですけど、速度は遅いです。なぜならこのイテレータの指しているデータは、メモリ上に連続していないからです。言い換えると、任意の自然数Nに対して『&(*(iter + N)) == &(*iter) + N』であることが保障されていないため、きちんと対象のデータを指し示すように補正しなければならない。その補正処理に余分な負担がかかって遅くなってしまう。ただ、もしメモリ上に連続であるならば直接ポインタをイテレータとして入れてやることで、(きちんと動作する保障はないですが)高速に動作してくれるでしょう。並列化はされていませんけれど。

また、上記のtransform_imageは並列化は行われておりますが、SIMDによる高速化は行われておりません。これをジェネリックに行うというのは現状としてなかなか厳しいのですが、Boost.SIMDの登場によってある程度の希望の光は見えてくるのではないかと思っています。

追記:

• (2011/09/16)コードをコンパイルが通るように若干の修正。ucharを暗黙的に使うのではなく明示的に指定するようにした
• (2011/09/17)記事を書いていて初めて知ったんですけど、ptr<T>(y)はcv::Mat_<T>だとoperator[]でオーバーロードされていたんですね。恐らくこれを用いることがOpenCV側の意図しているところだろうと判断し、上の式を書き換えました。

Thrustから見る、reduceを用いたアルゴリズム実装

"The Great Day of His Wrath" by John Martin

はじめに

世はまさに大並列時代。火を手にした人類が飛躍的な進化を遂げたように、並列化のパラダイムは、今まで到底不可能と思われていた速さで計算を行うことができるようになりました。

ところが、どんな手法にも欠点は存在するもので、実際に実装しようとすると非常に難しい。何故かというと、並列には並列化特有の問題が存在しており、愚直に実装してしまうとCPUより早くなってしまうどころか遅くなってしまうことだってあり得るのです。これを回避するにはGPUの内部構造についてきちんと理解をした上で、実装したいアルゴリズムそれぞれの場合に特化したコーディングを行う必要がある。

しかしよくよく考えるとおかしな話です。私たちが実装したいのはあくまで手法であり、ハードウェアではありません。なぜこのような詳細について把握する必要があるのでしょう。これはきちんとした抽象化がなされていない結果生み出された、言ってみれば妖です。この妖によって余計な手間が生み出され、コーディング量は増大し、余計なバグが生成され、絶望した開発者はうつになり、自殺者は増え、世界は暗雲に包まれるー

そのような混沌の中、ある一筋の光、Thrustが現れました。ThrustはCUDAを抽象化し、C++のパラダイムを用いることで本質にのみ注力できるよう開発されたライブラリです。ユーザはただIteratorに対してアルゴリズムを適用するかの如くdevice_vector<T>と各種アルゴリズムを駆使することで、流れるようなプログラムを組むことが可能となります。このような抽象を用いることで、ユーザは詳細な内部構造について頭を悩ませる必要が無くなりました。誰が書いているのか分からないから使いたくない？ThrustはNVIDIAに所属している開発者Jared HoberockとNathan Bell氏によって書かれています。2人とも並列化のプロフェッショナルなので、普通の人が愚直に実装するよりも大抵の場合効率的でしょう。

さて、Thrustには並列化を行うための多種多様なアルゴリズムが存在しています。そして、Thrustはオープンソースです。これはつまり、このソースコードを読めば並列化のプロがどのようにして各種アルゴリズムの並列化を実装しているのかが分かるということです。

Reductions

ということで、最近はThrustのソースコードをちょくちょく読んでいるところです。まだ読み始めた最中なのですが、特に感動したのはReductionsの項でした。

数々あるSTLのアルゴリズムを普通に実装しようとすると、一つ一つのアルゴリズムすべてに対して、コアレッシングなどのGPGPU特有の問題を回避したコードを記述しなければならず、非常に面倒です。Thrustはこの問題をreduce関数によって解決しました。この関数は並列化が可能なので、これを用いて書かれた関数は自動的に並列化されます。つまり、Thrustはいわば『全から一を返す』アルゴリズムすべてをreduceを用いて記述することによって、冗長な記述を回避しているのです。

このreduceを用いた並列化の実装について簡単に説明していくことが、今回の記事の目的です。reduceを初めて聞いた方は前回の記事『並列環境におけるreduceとscanアルゴリズム』を参照してください。

ここで、本来のThrustはC++によって記述しているのですが、そのままC++で書くと非常に冗長で解り辛くなってしまうので、今回は本質のみを抽出するためHaskellを用いて説明します。あまりHaskellに慣れていないので不適切な書き方等あるかもしれませんが、そのへんはご指摘頂けると幸いです。また、今回のreduceにはfoldr関数を代わりに用いています。分かると思いますが一応。

実装

• max_element, min_element, minmax_element

その名の通りリストの中の最大値、最小値、最大と最小値のタプルを返します。

*Main> max_element [1,3,2,4,3,1]
4
*Main> min_element [1,3,2,4,3,1]
1
*Main> minmax_element [1,3,2,4,3,1]
(1,4)

ここらへんはまぁ普通に分かります:

max_element :: Ord a => [a] -> a
max_element list = foldr max (head list) list

min_element :: Ord a => [a] -> a
min_element list = foldr min (head list) list

minmax_element :: Ord a => [a] -> (a, a)
minmax_element list = foldr minmax_func (first, first) $ zip list list where
    first = head list
    minmax_func (min1, max1) (min2, max2) = (min min1 min2, max max1 max2)

• all_of, any_of, none_of

all_of, any_of, none_of関数は対象のリストと真偽値を判定する関数を引数に取り、

• すべて当てはまっている
• 一つ以上当てはまっている
• すべて当てはまっていない

場合にはTrueを返す関数です。pythonではall(), any(), C++ではまんまのやつがありますね。

*Main> all_of (\x -> x > 3) [1,3,2,4,3,1]
False
*Main> any_of (\x -> x > 3) [1,3,2,4,3,1]
True
*Main> none_of (\x -> x > 3) [1,3,2,4,3,1]
False

これも実装的にはらくちん:

all_of :: (a -> Bool) -> [a] -> Bool
all_of binary_pred list = foldr (&&) True $ map binary_pred list

any_of :: (a -> Bool) -> [a] -> Bool
any_of binary_pred list = foldr (||) False $ map binary_pred list

none_of :: (a -> Bool) -> [a] -> Bool
none_of binary_pred list = all_of (not . binary_pred) list

関数合成を用いてエレガントに解決。

• inner_product, transform_reduce

inner_productは2つのリストに二項演算子を作用させてから、その結果からなるリストをreduceする関数です。transform_reduceはリストをいったんユニタリー関数を用いて変形させてから、reduceを実行する関数です。

*Main> inner_product (*) (+) 0 [0,1,2] [3,4,5]
14 -- = 0*3 + 1*4 + 2*5

*Main> transform_reduce (\x -> x*2) (+) 0 [1,2,3]
12 -- = 1*2 + 2*2 + 3*2

C++だといろいろ冗長に書かないといけないですけど、Haskellだったら余裕です:

inner_product :: (a -> b -> c) -> (c -> c -> c) -> c -> [a] -> [b] -> c
inner_product binary_op1 binary_op2 init list1 list2 =
    foldr binary_op2 init $ map operate $ zip list1 list2 where
        operate (x, y) = binary_op1 x y

transform_reduce :: (a -> b) -> (b -> b -> b) -> b -> [a] -> b
transform_reduce unary_op binary_op init list =
    foldr binary_op init $ map unary_op list

• count_if, count

count_ifは条件に合致している要素数を返す関数で、countは指定した値と同じ要素の数を返す関数です。

*Main> count_if (\x -> x < 3) [1,2,3,4,5]
2
*Main> count 3 [1,3,2,2,3]
2

count_if :: (a -> Bool) -> [a] -> Int
count_if binary_pred list = foldr (+) 0 count_list where
    count_list = map transform_func list where
        transform_func x = if binary_pred x then 1 else 0

count :: Eq a => a -> [a] -> Int
count value list = count_if (== value) list

ThrustではTrueだったら1, Falseだったら0となるリストを作成して、それを足していく手法をとっていました。
countはcount_ifが思いつけば楽勝でしょう。

• find_if, find_if_not

find_ifはリストの中で条件に合致している『一番最初の要素』を示すイテレータを返します。find_if_notは逆で、条件に合致していない『一番最初の要素』を示すイテレータを返します。今回は単純に要素のラベル番号を返すようにしました。なかったら最後のラベル+1が帰ってきます。

*Main> find_if (\x -> x > 3) [1,3,2,5,4]
3
*Main> find_if (\x -> x > 7) [1,3,2,5,4]
5
*Main> find_if_not (\x -> x > 3) [1,3,2,5,4]
0

これは微妙に入り組んでいます:

find_if :: (a -> Bool) -> [a] -> Int
find_if binary_pred list = snd result where
    result = foldr find (False, length list) tuple_list where
        find (b1, i1) (b2, i2)
            | b1 && b2  = (True, min i1 i2)
            | b1        = (b1, i1)
            | b2        = (b2, i2)
            | otherwise = (False, max i1 i2)
        tuple_list = zip (map binary_pred list) [0..]

find_if_not :: (a -> Bool) -> [a] -> Int
find_if_not binary_pred list = find_if (not . binary_pred) list

ただ、find_ifはちょっと本家と違います。インデックスが小さい項が常に左の引数となるため省略してるんだとは思いますけど、対称性を満たしていないのがちょっと気持ち悪いので、ここでは交換則を満たすように関数を変形しています。

実際の実装と課題
ここではhaskellを用いて実装しましたが、もちろん実際の実装は異なっており、zipはzip_iterator, mapはtransform_iterator, タプルはtupleを用いて実装しています。これらのイテレータの評価は遅延的に行われるので、複数のメモリアクセスが発生せず、アルゴリズムの高速化に貢献します。

これらのfancy iteratorを使用した戦略は非常に賢いとは思いますが、同時に欠点もあります。複数のイテレータを大量に組み合わせるために、通常のSTL algorithmのようにbeginとendを使用してしまうと、両方に対してfancy iteratorを組み合わせなければならないため、2倍の手間がかかってしまうのです。

これらの解決案としてはboost::rangeのようなRange構造を実装することです。これをパイプ演算子で繋げていくような形をとっていけば、明快で分かりやすく、速度も十分保証されたプログラムを記述できることでしょう。

Wiresharkを用いたHTMLの初歩的なスクレイピング

前回はぽんとスクリプトを上げて「これで画像がダウンロードできますよ」とやったわけだけど、もちろん実際にGoogle Art Projectの資料が上げられているわけでもないので、きちんとHTMLやJavascriptのソースを見ながら読解していくことになる。

簡単な動作機構だったら単にHTMLのソース見るだけで良いのだけど、Google Art Projectのような複雑な機構だとそうもいかない。行数が多いと面倒だし大抵の場合難読化されてたりで一筋縄ではいかないからだ。そこで、今回はパケット解析というもうひとつの武器を利用することで内部構造を別の側面から理解することにした。思考の流れを書き綴っていくので、HTMLのスクレイピングをやってみたい方にとって何か参考になれば幸いです。

目標: 高画質な絵画をGoogle Art Projectを利用して手に入れる

まずGoogle Art Projectのサイトへアクセス。絵画ごとに固有のURLが割り当てられており、このURLを解析すれば絵画を入手できると判断。見るからにソースコードからの解析は面倒くさそうだったので、パケット解析ソフトウェアWiresharkを用いて解析を行った。

本質のみに注力するためhttpプロトコルのみキャプチャ。適当なURLにアクセスし、ズームを1回行い、パケットキャプチャを中断。解析を開始。

画像本体らしきURLにアクセスしているパケットを発見。複数のアクセスがあるってことはたぶん画像は複数に分割されていてそれを読み込んでいるんだろう。調べると/unique-id=xi-yj-zkの形でアクセスしていたので、対応するxy座標とズーム値zによって対応する画素値にアクセスすることができると推測。試しにURLにアクセスしてみると、問題なくダウンロードができた。

また、x,yの値を変更してのアクセスを行った結果も問題なくアクセスできた。ただ、どうやら範囲外のxy値にアクセスすると404が返されるらしい。当たり前といえば当たり前。アクセス禁止などのペナルティはないっぽいので、(ちょっと礼儀悪いけど)404が返されたらそれが範囲外と解釈することにする。

一旦まとめてみる。画像は分割されており、これらのアクセスはなんかよく分からないunique-idという、おそらく絵画について表しているパラメータとx,y,z座標の2つだけでダウンロードができる。ついでに、リファラー, cookie, 認証などの面倒な作業をせずとも直接のダウンロードが許されていることも分かる。そこらへんはGoogleさん寛容なんだなぁと思い、まぁAPIを提供してるくらいだからなぁ、とかいろいろ考える。

他の画像に対してのホストを確認してみると、どうやら[lh3〜lh6]までの複数のホストへとアクセスしているらしい。このホストアクセスがランダムか、それとも固有のホストを呼び出さなければならないのか確認してみたかったので、試しにホストだけを変更して、その他のパラメータは固定したままでのアクセスを試みる。結果としては問題なくダウンロードできた。たぶん負荷を分散するという単純な理由から、ランダムにJavascript側で割り振っているんだろう。ここで初めてHTMLのソースを閲覧。"cmn/js/content.js"のソースを開き(難読化はされてないようだ)"ggpht"を検索。

リストに入れているところからすると、画像へのアクセスはたぶんこのリストの中だけに限られているんじゃないかなと思う。試しにlh7へアクセス。404が正常に返された。的中。この指針で間違いないので次へ進む。

次に気になるのはunique-idだ。一体何処から仕入れているのか知らないと全く手のつけようがない。考えられるのは2つで、

• アクセス毎に固有のunique-idを発行し、一定期間のみアクセスできるようにしている
• HTMLかどこかに参照しやすいように書かれており、それを抜き出してアクセスする

という発想だ。自分は前者を疑っていたので、それっぽいパケットがあるかどうか監視。結果としてはまったく存在していなかった。というと後者かなと思い、対象のHTMLからunique-id(今回の場合だと"vtzd432xqzOpi_3X8JAJ_buly36QKI0n9zGeeWNgBcwsYJTsAKK5mbM9Pgg")を検索。

ビンゴ！これで駒はすべて揃った。まとめに入る。対象の絵画をダウンロードするためには、

1. 取得したい絵画のアドレスへアクセス
2. <div id="microscope" data-microId="unique-id">となっている項を取得、unique-idを得る
3. http://lh3.ggpht.com/ - http://lh6.ggpht.com/のいづれかにアクセス。
   http://lh3.ggpht.com/unique-id=xi-yj-zk のような形で画像を取得
4. 取得した複数の画像を繋ぎ合わせる

後は非常に簡単でただこの単純なアルゴリズムに従って書き下していくだけ。その産物がこいつ。こんなに早く仕上がるってことは、たぶんGoogleさんは本格的にダウンロードを禁止したいわけじゃないんだろう。

スクレイピングとしては楽な部類で、普通だともうちょっと複雑な認証をかましてくるので自動化しにくい。Captcha入ると流石に無理だけど、おそらく現在の文字認識技術をきちんと応用すれば解けるんじゃないかと思っている。

余談

• こういう解析は楽しい。どんな感じの楽しさかっていうと、与えられたパズルを解いていく感覚にちょっと近い
• 自分はだいたい1時間くらいで解いたけど、早い人はたぶんもっと早いんだろうなって思う。まだまだ精進が足りない

Google Art Projectを利用して高画質の絵画を手に入れる

"The Starry Night" by Vincent van Gogh

2011年2月2日、GoogleさんはGoogle Art Projectというサイトをリリースしました。どういうサイトかっていうとMoMAニューヨーク近代美術館やTate Britainなど、世界中の美術館にある有名な絵画を渡り歩くことができるという趣旨のサイトで、誰もが一度は見たことがある絵画がー家に居ながらー非常にリアリスティックな解像度で見ることができる。恐ろしい時代になったもんです。

ところがこのサイト、高解像度で見られるのは非常にいいことなのですが、肝心のダウンロードができない。せっかく高解像度なんですから保存できてもいいはずなのに、だけどできない。Google Mapなんかは常に最新の情報に更新しなければならないので、ダウンロードに制限をかけているのは合理性があるんですけど、絵画とか全く更新する必要がないでしょう。著作権も切れてるのに。よく分からないです。

ということで、なければ自分でなんとかしようという精神で、1時間くらいでちゃっちゃとPythonでスクリプト組んでダウンロードするようにしてみました。即興で組んだのでいろいろ粗い面はありますけど、そこら辺は勘弁してください。

#!/usr/bin/env python
#-*- coding: utf-8 -*-

import sys
import re
import urllib2
import random
from PIL import Image

def usage():
    print "%s [URL]" % sys.argv[0]
    return 1

def get_id(url):
    response = urllib2.urlopen(url)
    if response.code != 200:
        raise TypeError
    data = response.read()
    return re.search(
        r'<div id="microscope" data-microId="(.+)">', data).group(1)

def get_filename(x, y):
    return "%i_%i.jpg" % (x, y)

def download_image(target_id, x, y, z):
    urllist = [
        "http://lh3.ggpht.com/",
        "http://lh4.ggpht.com/",
        "http://lh5.ggpht.com/",
        "http://lh6.ggpht.com/"
    ]
    base_url = random.choice(urllist)
    url = "%s%s=x%i-y%i-z%i" % (base_url, target_id, x, y, z)

    try:
        response = urllib2.urlopen(url)
    except urllib2.HTTPError:
        return False
    with file(get_filename(x, y), "wb") as fp:
        fp.write(response.read())
    return True

def download_and_get_max_x(target_id, x, z):
    result = download_image(target_id, x, 0, z)
    if result:
        return download_and_get_max_x(target_id, x+1, z)
    else:
        return x

def download_and_get_max_y(target_id, y, z):
    result = download_image(target_id, 0, y, z)
    if result:
        return download_and_get_max_y(target_id, y+1, z)
    else:
        return y

def combine_images(name, max_x, max_y):
    x_size = 512
    y_size = 512

    result = Image.new("RGB", (max_x*x_size, max_y*y_size))
    for y in range(max_y):
        for x in range(max_x):
            image = Image.open(get_filename(x, y))
            result.paste(image, (x*x_size, y*y_size))
    result.save(name)

def main():
    if len(sys.argv) == 1:
        return usage()
    depth = 3
    target_addr = sys.argv[1]
    target_id = get_id(target_addr)
    max_x = download_and_get_max_x(target_id, 0, depth)
    max_y = download_and_get_max_y(target_id, 1, depth)
    for y in range(1, max_y):
        for x in range(1, max_x):
            download_image(target_id, x, y, depth)
    combine_images("result.png", max_x, max_y)

if __name__ == "__main__":
    sys.exit(main())

使い方は簡単で、

$ python artproject.py http://www.googleartproject.com/museums/tate/mariana-248

のようにURLを指定してあげると勝手に取得して結合してくれます。こんな感じで:

"Mariana" by John Everett Millais

あぁそれにしてもお美しい……巷ではオフィーリアちゃんが有名だけど、やっぱりマリアナさんのほうが綺麗だよ。絶対。

余談

• もともとの動機としては高解像度な絵画の壁紙が欲しくて、いろいろと探し回ったわけですが、どうもGoogle Art Projectしか高解像度の絵画が載っていないので仕方なくといった次第です。
• 利用はどうやらGoogle's Terms of Serviceに準じているようです。ざっと見した感じだとダウンロードを禁止する条項はないですし著作権も切れているので全く問題ないように思えるんですが、もし警告とかきたらチキンですから遠慮無く取り下げます。ご了承ください。
• このネタを題材に時間があればHTMLのスクレイピング、パケット解析についての入門記事を書いてみたいです。
• なんだかんだ言ってちゃんと実物見たほうが綺麗。ミレイの作品は以前東京に行ったときにBunkamuraで見ました。Bunkamuraはそういう催し物よくやっていて、そこらへんが東京羨ましかったり

並列環境におけるreduceとscanアルゴリズム

"fireking" by cmyk1219

1. 概要
reduceとscanはGPGPUなどの並列環境下において重要な役割を果たす関数だけど、あまりこれらの関数、特にscanについて言及している記事はあまり見かけない。なので今回はreduceとscanについて調べた結果をまとめていきたいと思う。

2. reduce, scan関数の概要
2.1. reduce
scanについてはreduceが分かればより理解しやすくなるので、ひとまずはreduceから始めることにする。

reduceとは、端的に言うと、対象のリストを与えられた二項演算子で纏め上げる関数だ。とは言っても、いきなりそんなことを言われてもよく分からない。まずは例を見ていこう:

>>> import operator
>>> reduce(operator.add, [1, 2, 3, 4, 5])
15
>>> reduce(operator.sub, [1, 2, 3, 4, 5])
-13
>>> reduce(operator.mul, [1, 2, 3, 4, 5])
120

reduceにおいて重要なのは二項演算子であり、上記の場合、reduceは以下のような演算を行う。

(((1 + 2) + 3) + 4) + 5 = 15
(((1 - 2) - 3) - 4) - 5 = -13
(((1 * 2) * 3) * 4) * 5 = 120

これにより、最初の演算はsum関数、3番目の関数はprod関数と等価であることが分かる。
リストの型は何も数値に限られている訳ではない:

>>> reduce(operator.and_, [True, False, True])
False
>>> reduce(operator.or_, [True, False, True])
True

これより、上記2つの演算はall, any関数と機能的に等価であることが分かるのだけれど、all, any関数のほうが効率的なので、実際はちゃんとall, any関数を使うべきだ。

他にもreduceを用いて様々な関数、アルゴリズムを記述できるが、あまりにも多いので今回は後回しにして、scanを重点的に解説していく。

2.2. scan
scan関数はreduceと比べて文献が少ない。言及している文章もあまり見かけないが、並列処理においてscanは重要な関数だ。残念ながらpythonで実装されていないけれど、仮に実装されていたとしたら以下のような働きをする:

>>> scan(operator.add, [1, 2, 3, 4, 5])
[1, 3, 6, 10, 15]
>>> scan(operator.sub, [1, 2, 3, 4, 5])
[1, -1, -4, -8, -13]
>>> scan(operator.mul, [1, 2, 3, 4, 5])
[1, 2, 6, 24, 120]

きちんと書き下せばどのような働きをしているのか分かりやすい:

[1, 1+2, (1+2)+3, ((1+2)+3)+4, (((1+2)+3)+4)+5]

結局のところ、scanは先頭からreduceしていった結果をリストにまとめ、それを返す関数となる。scanは別名prefix sumとも呼ばれている。
ちなみに、haskellではscanl1(scanl)という名前で実装されている。

Prelude> :type scanl1
scanl1 :: (a -> a -> a) -> [a] -> [a]
Prelude> scanl1 (+) [1,2,3,4,5]
[1,3,6,10,15]

また、CUDAライブラリであるthrustにはinclusive_scan(exclusive_scan)という名前で実装されている。

#include <thrust/scan.h>
int data[6] = {1, 0, 2, 2, 1, 3};
thrust::inclusive_scan(data, data + 6, data); // in-place scan
// data is now {1, 1, 3, 5, 6, 9}

3. 並列による制限
何故並列下においてreduce, scanが重要になってくるのか、それは、両方ともメニーコアな環境において効率的に動作するアルゴリズムが複数提唱されているからだ。reduce, scanは並列と相性がいい。つまり、既存のアルゴリズムをどうにかしてreduce, scanに落とし込むことができれば、並列化の恩恵を受けれられると。

ただし、並列環境においては、演算にいくつかの制限が生じてくる。具体的には、集合Gと演算子"・"の組において、

• 演算子"・"は交換則"a・b = b・a"を満たしていなければならない
• 演算子"・"は結合則"(a・b)・c = a・(b・c)"を満たしていなければならない

後者に関してはあたり前の話で、これが成り立っていなければ並列に計算を行うことができない。前者については多少曖昧で、実際には満たさなくてもよいアルゴリズムも作ることはできるけど、コアレッシングなどの問題を考えると満たしていることが強く求められる。実際、thrustのreduceは交換則を満たしていることを前提としている。

reduceの場合はこの2つの制限で十分だけれど、scanの場合は上記2つの制限に加えてさらに2つの制限を満たしておくと都合がいい。具体的には、集合Gと演算子"・"の組において、

• ただ一つの単位元"e"が存在する
• 演算子"・"は任意の元"a"に対して"a・a^{-1 =} a^-1・a = e"を満たす、ただ一つの逆元"a^-1"が存在する

これら4つの制限を満たしている集合Gと演算子"・"の組を『アーベル群』という。要するに、scan関数は、対象の型から成る集合と二項演算子の組がアーベル群を成していれば都合がいいということになる。

具体例を出そう。整数と和の組(Integer, (+))は明らかに上記4つの制限を満たす。よって、この組は並列にreduce, scanオペレーションを適用でき、かつ都合が良い。

1 + 2 == 2 + 1
(1 + 2) + 3 == 1 + (2 + 3)
1 + (-1) == (-1) + 1 == 0

真偽値と論理積の組(Bool, (&&))は上記2つの制限を満たすため並列にreduce, scanオペレーションを実行可能だが、後半2つの制限は満たしていないため都合は良くない。

3次元の回転行列と積の組(Matrix, (*))は交換則を満たしていないため―少なくともthrustでは―並列にreduce, scanオペレーションを実行することはできない。

3.1. 都合が良いとは？
とは言っても、一体アーベル群であればどこがどう都合がいいのか？それは、リスト中において局所的にreduceを行いたい場合に役に立つ。

リスト[a₀, ... , a_n]にscanオペレーションを実行して

を手に入れたとする。この場合、i < jにおいてb_j+b_i^-1を計算すると、

となる。つまり、i+1からjまでの要素についてのreduceを僅かO(1)で求めることができる。

これには大きな利点がある。例えば、大規模な時系列データを、その周りのデータから平均化(平滑化)したいとする。各要素ごとに計算を行うのは明らかに非効率であるので、まずscanオペレーションを(+)に対して適用し、その後で差分を計算し、正規化を行ってやればよい。結局、scan (+)のやっていることはf(x)の不定積分F(x)を求めることと本質的に等価となる。この考え方はCVでは、積分画像として様々な箇所で用いられる。

そしてさらに、この操作はアーベル群であれば何でも良い……とは言っても、アーベル群であることは結構厳しい制限となる。

例を出すと、リスト中の任意の範囲ーiからjの成分がすべてある条件式に適合しているかどうか確かめたいとする。
関数型の考え方だと、まずはmap関数を用いて(C++でならtransform_iteratorを用いて)リストの成分をBoolに変え、その後でscanを利用することを思いつくかもしれない。しかし、残念ながら(Bool, (&&))はアーベル群でないので、そのまま愚直に適用することができない。

Prelude> scanl1 (&&) $ map (<3) [1,2,3,2,1]
[True,True,False,False,False]
-- 最後の2要素については条件を満たしているのに、Falseで埋もれてしまって分からない！

この場合は、単純にアーベル群である(Integer, (+))をscanに適用してやればよい。

Prelude> scanl1 (+) $ map (\x->if x<3 then 1 else 0) [1,2,3,2,1]
[1,2,2,3,4]
-- きちんと第3要素のみ条件から外れていることが分かる

Pythonっぽく書くとたぶんこんな感じ:

>>> import operator
>>> f = lambda x: 1 if x < 3 else 0
>>> scan(operator.add, [f(x) for x in [1,2,3,2,1]])

4. より詳しく知りたい方は
scanアルゴリズムのGPGPU実装に関する詳細については、現在のところ『Parallel Prefix Sum (Scan) with CUDA(Mark Harris, NVIDIA)』が一番分かりやすい。より詳しく知りたい方はそちらをどうぞ。

というより、基本NVIDIAの出しているテキストは分かりやすい。すごいことです。

boost.GILで組む、ジェネリックな画像処理のアルゴリズム

"Image dans le néant" by gelinh

boost.GILは凄い！開発者の頭の良さがビシビシと伝わってくる！
ということで、今回はGILに関しての紹介記事を書こうと思います。

概要
あなたは画像処理のエキスパート。顧客の依頼で、8bitのRGB画像を処理するアルゴリズムを記述していたとします。ところが対象となるデバイスの仕様を調べていた際に、実はRGBA画像にも対応させなければいけないことが分かりました。面倒だと思いながら書いていた矢先、さらにBGRやABGR,さらには16や24bitにも対応したアルゴリズムを記述しなければならないことが判明しました。なんということでしょう…これらの画像すべてに対してアルゴリズムを書くなんて、とてもじゃないですがやってられない。やめてくれ！って感じです。

boostに付いてくるGILを用いることで、画像に対する操作をよりジェネリックに行うことが出来ます。具体的に言うと、あなたはGILの作法に従ってアルゴリズムを書き下すだけで、メモリに展開されている画像データの色空間、色深度、レイアウト、配置順に『関わらず』、任意のデータ配置に対して適用できるアルゴリズムを作ることができるのです。最初に提示した問題が一気に解決しちゃいました。

長所と短所
長所

• 画像処理に対してジェネリックなアルゴリズムを記述でき、少ない記述量で全てをカバーできる。処理速度も普通に書き下した場合と同等(静的なポリモーフィズムしか使っていないので、コンパイル時に必要な計算はすべて終わっている)。
• 様々なアダプタが存在している(e.g. subimage_view, rotated180_view)ので、ばんばん遅延評価して一気に処理を行うことが可能。
• 普通の人が普通に書くよりも、GILのアルゴリズムに従って書いたほうが大抵の場合早い
• ソースコードが比較的読みやすい。コード量も短くなる。
• 誰が作ったかわかんないもの利用したくねーという方も安心のAdobe製。画像処理の専門家が作ってるのでとても合理的。

短所

• テンプレートを多用しており『非常に難解』。
• 全てを抽象化しているために『非常に難解』。
• C++と画像処理の両方に秀でていなければ記述できないので、使っている人が殆どいない。
• 英語のサイトですら文献が殆どない

正直なところ、短所の部分がかなり厳しい。かなりC++をやりこんでいないと使えないということは、それだけでユーザ層を大分狭めます。仕方がないことですけど。

Viewの概要
GILでは、画像に対する操作はすべて『View』を起点にして行われます。Viewは以下のような構造を持ち、ジェネリックな操作を行えるよう工夫してあります。


View
|-- Locator
|    `-- Pixel
|         |-- Channel
|         `-- Layout
|              |-- Color Space
|              `-- Channel Mapping
`-- dimensions

(上図はMemory basedで、画素のメモリ配置がInterleavedな場合のView構造)

以下のリストで、それぞれの名称がどのような役割を持っているのか説明します。

• View: (x, y)座標における画素値の他、特定のX、Y方向におけるイテレータや全ピクセルにわたって捜査できるxyイテレータなど、「様々な手法によって画素にアクセスできる手段」を提供する、一種の抽象的なクラスを表しています。画像の幅、高さを表すdimensionsと、後述するLocatorを持ちます。(注: 画像のデータを実際に保持している『わけではありません』。Viewはあくまで、対象の画素値へとアクセスするための『手段』を提供しているだけに過ぎないのです)
• Locator: 画像の2次元的な『位置』を表します。ランダムアクセスイテレータの2次元版のようなものと考えれば理解が早いと思います。具体的に言いますと、 *loc あるいは *(loc + point2(dx, dy)) などで現在の画素値や(dx, dy)だけ離れた位置の画素値を取得することは出来ますが、Locatorの現在地(xy座標)や画像の大きさ(width, height)を取得することはできません。
• Pixel: Locatorによって返される画素値を表します。後述するChannelとLayoutを持ちます。
• Channel: 画素の色深度を表します(8bit, 16bit, etc…)。各型の取りうる最大値、最小値も取得できます。
• Layout: 色空間のレイアウトを表します。ちょっと分かりにくいので、下の2つを読んだほうが早いでしょう。
• Color Space: 色空間を表します(e.g. RGB, CMYK, HSV, Gray, etc…)。
• Channel Mapping: 実メモリ上での配置順を表します(e.g. RGB, BGR, etc…)。ただし、この情報はInterleave画像のみ用います(Planar画像だと配置順は関係ないため)。

アルゴリズムを用いて書いてみる
それでは実際にGILに付属しているアルゴリズムを用いて実際に書き下してみます。以下の2つのアルゴリズムを用いてみましょう。

boost::gil::transform_pixel(const View& src, const View& dst, F functor)
対象のピクセルを受け取って別のピクセルを返す関数を、全ピクセルに対して適用します。この関数は、画素値と出力値が1対1対応である場合に用います(e.g. 明度、コントラスト, レベル補正, トーンカーブ, etc...)。

boost::gil::transform_pixel_positions(const View& src, const View& dst, F functor)
周囲のピクセルの値を元に結果となるピクセルを返す関数を、全ピクセルに対して適用します。transform_pixelとは、受け取る引数がロケータであるという点が異なっており、周囲の画素値を元に実際の出力値を求める場合に有効です(e.g. Blur, Sovel, Laplacian, etc...)。

実際にやってみれば、どんな感じか分かるでしょう:

#include <boost/gil/gil_all.hpp>
#include <boost/gil/extension/io/jpeg_io.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace boost::gil;

struct change_brightness {
public:
    change_brightness(): weight_(1.0) {};
    explicit change_brightness(float weight): weight_(weight) {};

    template<typename Pixel>
    Pixel operator()(Pixel src) {
        typedef typename channel_type<Pixel>::type channel_type;
        Pixel dst;
        for(int c=0; c<num_channels<Pixel>::value; ++c) {
             dst[c] = cv::saturate_cast<channel_type>(
                     static_cast<float>(src[c] * weight_));
        }
        return dst;
    }
private:
    float weight_;
};

struct change_brightness_locator {
public:
    change_brightness_locator(): weight_(1.0) {};
    explicit change_brightness_locator(float weight): weight_(weight) {};

    template<typename Locator>
    typename Locator::value_type operator()(Locator loc) {
        typedef typename Locator::value_type pixel_type;
        typedef typename channel_type<Locator>::type channel_type;
        pixel_type src = loc(0, 0);
        pixel_type dst;
        for(int c=0; c<num_channels<Locator>::value; ++c) {
            dst[c] = cv::saturate_cast<channel_type>(
                    static_cast<float>(src[c]) * weight_);
        }
        return dst;
    }
private:
    float weight_;
};

template<typename Locator>
struct x_gradient_functor {
public:
    explicit x_gradient_functor(Locator loc)
    : left_(loc.cache_location(-1, 0)),
      right_(loc.cache_location(1, 0)) {};

    typename Locator::value_type operator()(Locator loc) {
        typedef typename Locator::value_type pixel_type;
        typedef typename channel_type<Locator>::type channel_type;
        pixel_type src_left = loc[left_];
        pixel_type src_right = loc[right_];
        pixel_type dst;
        for(int c=0; c<num_channels<Locator>::value; ++c)
            dst[c] = (src_left[c] - src_right[c])/2;
        return dst;
    }

private:
    typename Locator::cached_location_t left_;
    typename Locator::cached_location_t right_;
};

template <typename SrcView, typename DstView>
void x_gradient(const SrcView& src, const DstView& dst) {
    assert(src.dimensions() == dst.dimensions());
    transform_pixel_positions(
            subimage_view(src, 1, 0, src.width()-2, src.height()),
            subimage_view(dst, 1, 0, src.width()-2, src.height()),
            x_gradient_functor<typename SrcView::locator>(src.xy_at(0, 0)));
}

int main() {
    const char* filename = "lena.jpg";
    rgb8_image_t src_image;

    jpeg_read_image(filename, src_image);
    rgb8_image_t dst_image(src_image.dimensions());

    transform_pixels(
            const_view(src_image),
            view(dst_image),
            change_brightness(1.5));
    jpeg_write_view("result_normal.jpg", view(dst_image));

    transform_pixel_positions(
            const_view(src_image),
            view(dst_image),
            change_brightness_locator(2.0));
    jpeg_write_view("result_locator.jpg", view(dst_image));

    x_gradient(const_view(src_image), view(dst_image));
    jpeg_write_view("result_gradient.jpg", view(dst_image));

    return 0;
}

また、画像の出力例を下記に示します:

左から、元画像, change_brightness, change_brightness_locator, x_gradientの結果

change_brightness_locator()については、for_each_pixel用のアルゴリズムをlocatorを用いて適用させた場合にどうなるか確認するためだけに実装しましたので、実務的ではありません。また、値を型の持つ最大、最小値に丸めるためにOpenCVのcv::saturate_castを使用しています。

main()についてはどんなことをやってるのか大体感覚的に理解できるはずですので、とりあえず3つのファンクタに関して説明を加えていきます。

typedef typename channel_type<Pixel>::type channel_type;

channel_type<T>::typeで実際に使われているチャネルの型を取得できます。Tにはピクセル、ロケータ、ビューのどれを指定しても構いません。今回の場合ですと8bitのチャネルを指定しているので、型にはunsigned charが入ります。

num_channels<Pixel>::value

ピクセルを構成しているチャネルの数を取得できます。この場合ですと色空間がRGBなので値には3が入ります。この値はコンパイル時に決定されるので、最適化を用いることでループ展開が行われ、速度的なオーバーヘッドはありません。

typedef typename Locator::value_type pixel_type;

ロケータが返すピクセルの型を取得できます。この場合ですとrgb8_pixel_tが入ります。

pixel_type src = loc(0, 0);

()演算子によって現在のピクセル値を取得しています。他にも*演算子でイテレータっぽくアクセスすることもできます。

explicit x_gradient_functor(Locator loc)
    : left_(loc.cache_location(-1, 0)),
      right_(loc.cache_location(1, 0)) {};

ロケータの変動値をキャッシュしておくことで、高速化を図ります。Locator::cached_location_tによってロケータのキャッシュ型を取得します。キャッシュを使ったアクセスには[]演算子を用います。
仰々しく言っていますが、要はstd::ptrdiff_tみたいなものです。移動する値は同じなんだから記憶しておきましょうということです。

subimage_view(src, 1, 0, src.width()-2, src.height())

subimage_view()によって画像の領域を切り出しています。これはアダプタで、実際の評価は遅延的に行われます。

まとめ
大体こんな感じでアルゴリズムを組んでいくのがGIL流です。遅延評価とアルゴリズムを武器に、ジェネリックなアルゴリズムを作っていきましょう。

余談
C++に関しての記事を投稿するのは正直に言って『少々怖い』です。大筋では合っていると思いますが、完全に合っているという自信があまりないからです。ですが、GILに関しては本当に、本当に参考となる文献が少ない。サンプルコードも多くない。こういった中でこの記事が少しでも参考程度になっていただければ幸いです。

それと、ざっくり組んでしまったのでこうしたほうがいいよとか、間違い等あれば遠慮なく言ってもらえると助かります。

(2/16) 追記: ソースコード若干の修正、説明や画像例の追加。
(3/31) typo: e.q. -> e.g.

画像をぼかす: いろんなボケ、リアルなボケ

『画像をぼかす』と一口に言っても、その方法にはいろんなものがあります。一番有名なのはPhotoshopの『ガウスぼかし』に見られるようなGaussian Blurですが、残念ながらこれはカメラに見られるようなリアルなぼかしと比べるとどうしても見劣りします。原因はいくつか考えられて、一つはGaussian Blurに使われているガウス関数の勾配が穏やかなため、現実のレンズのシチュエーションを考慮していないという点があります。

それじゃあ勾配を急にした関数ー例えばフェルミ分布関数のようなーをカーネルに採用すれば良いのかというとそういうわけではありません。有名なLena画像で試してみましょう。

左が元画像、中央がカーネルにガウス関数を用いた画像、右がフェルミ分布関数を用いた画像。
確かに中央よりは綺麗だけど、それでもまだリアルとは言えない

それじゃあ他に何が違うのでしょう？端的に言うと、黒も白も等価に扱って平均してしまうのがマズい。

画像をぼかすという行為は要は畳込み積分で、周囲のピクセル値を特定のウェイトで平均してやってることになります。で、普通にブラーしてしまうとピクセルの輝度に関わらず常に一定の割合で平均してしまいますので、なんか濁ったような画像が生成されてしまうと。

通常のブラー。カーネルにはよりリアルに見せるため正六角形を用いている

それではどのようなアプローチを用いればいいのでしょうか？レンズブラーを観察してみると、明るい箇所が強調されてぼかされていることが分かります。言い換えると、輝度の高いピクセルが優先的に混合されている、つまり予めRGBの全ピクセル値を適当な単調増加関数で写像してから畳込み積分を行って、その後対応する逆関数で戻してやればいい。要はhirax.netとやっていることは同じで、expしてからボカしてlogすると、こういうわけです。

結果の画像。比較的リアルにボカされていることが分かる

並べてみれば、違いがはっきり(クリックで拡大)

Lenaさんだけでは少し分かりづらいので、適当な背景写真に適用してみましょう。

森林画像

まず、通常の手法でぼかしてしまうとどうなるでしょう？

森林画像(通常のブラー)

お世辞にも綺麗なボケとは言えません。それじゃあ今回の手法でぼかしてみます。

森林画像(リアルブラー)

なんということでしょう…リアルなボケ画像ができました。

比較画像(クリックで拡大)

今回組んだOpenCVのソースコードを下記に示します。2.2でコンパイルしましたが、多分2.0以上だったら普通に通ると思います。

#include <algorithm>
#include <cmath>
#include <string>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/core/operations.hpp>

struct exp_each_element
{
public:
    exp_each_element(float factor): factor_(factor/256.0){};
    cv::Vec3f operator()(cv::Vec3b &color) {
        return cv::Vec3f(
                std::exp(static_cast<float>(color[0])*factor_),
                std::exp(static_cast<float>(color[1])*factor_),
                std::exp(static_cast<float>(color[2])*factor_));
    }
private:
    const float factor_;
};

struct log_each_element
{
public:
    log_each_element(float factor): factor_(factor/256.0){};
    cv::Vec3b operator()(cv::Vec3f &color) {
        return cv::Vec3b(
                cv::saturate_cast<uchar>(std::log(color[0])/factor_),
                cv::saturate_cast<uchar>(std::log(color[1])/factor_),
                cv::saturate_cast<uchar>(std::log(color[2])/factor_));
    }
private:
    const float factor_;
};

cv::Mat make_fermi_kernel(float radius, float range) {
    const int size = static_cast<int>(2.0*range + radius)*2 + 5;
    const float mu = static_cast<float>(size/2);
    const float range_inversed = 1.0 / range;
    cv::Mat dst(size, size, CV_32F);
    for(int i=0; i<size; ++i) {
        for(int j=0; j<size; ++j) {
            const float u = static_cast<float>(i - mu);
            const float v = static_cast<float>(j - mu);
            const float r = std::sqrt(u*u + v*v);
            dst.at<float>(i, j) =
                    1.0/(std::exp(range_inversed*(r - radius)) + 1.0);
        }
    }
    cv::normalize(dst, dst, 1.0, 0, CV_L1);
    return dst;
}

cv::Mat make_regular_polygon_kernel(int number, float radius) {
    const int size = static_cast<int>(2.0*radius) + 1;
    const float center = static_cast<float>(size)/2.0;
    std::vector<cv::Point> dst_vector;
    for(int i=0; i<number; ++i) {
        const float theta = 2.0*M_PI*i/number;
        const cv::Point point = cv::Point(
                center + radius*std::cos(theta),
                center + radius*std::sin(theta));
        dst_vector.push_back(point);
    }
    cv::Mat dst(size, size, CV_32F);
    const cv::Scalar color(1.0, 1.0, 1.0, 1.0);
    cv::fillConvexPoly(dst, &dst_vector[0], number, color);
    cv::normalize(dst, dst, 1.0, 0, CV_L1);
    return dst;
}

int main(int argc, char** argv) {
    const std::string source = argc >= 2 ? argv[1] : "lena_std.tif";
    const std::string destination = argc >= 3 ? argv[2] : "out.jpg";
    const float factor = 6.0;

    cv::Mat image = cv::imread(source);
    CV_Assert(image.type() == CV_8UC3);
    cv::Mat kernel = make_regular_polygon_kernel(6, 8.0);

    cv::Mat exp_image(image.size(), CV_32FC3);
    std::transform(image.begin<cv::Vec3b>(), image.end<cv::Vec3b>(),
            exp_image.begin<cv::Vec3f>(), exp_each_element(factor));

    cv::filter2D(exp_image, exp_image, -1, kernel);

    std::transform(exp_image.begin<cv::Vec3f>(), exp_image.end<cv::Vec3f>(),
            image.begin<cv::Vec3b>(), log_each_element(factor));

    cv::imwrite(destination, image);
    return 0;
}

余談

• 実質80行くらいで完成したので嬉しいです。可読性も多分結構読みやすい部類だと思う
• C++のコードが何やってるか分からない方は、まず関数オブジェクトとSTLについて調べてみましょう。あれはいいもんです

追記(11/03/23)

どんなデバイスであるか、たとえば、それが入力デバイスであるのか、出力でバイスであるのか、はたまた、計算処理デバイスであるのか次第で「輝度値の単位系」をどのように保持（処理）するべきかは違うことでしょう。
大切なことは、取り扱う「値」がどんな単位であるのかを見失わないことなのか、と思います。 

hirax.net::「コンピュータ画像処理」と「輝度値の単位系」 

自分がこのような『任意の単調増加関数で良い』という具合で筆を進めたのには2つ理由があります。

まず1つ目が、『全てのデジタルカメラにおいて、ピクセルは輝度の対数の次元を取っているのだろうか？』という単純な疑問です。

きちんと物理的/工学的な意味合いを持つ画像処理にするためには、格納されているピクセルデータの次元がきちんと揃えられている必要があります。自分自身も普通に考えて輝度の対数が格納されているんだろうということで、任意ではなく指数関数に制限された関数系を使いましょうと運びたかったのですが、そうするための証拠が不足していた。

もしかしたら自分の予想もつかない方法で処理されている可能性があるわけで、そういった証拠が十分にないことから、単調増加関数と少し曖昧とした意味合いを含ませたのです。

もう一つの理由は『今回の画像処理は物理的に正確である必要がない』という理由です。これが大きかった。

そもそも今回の画像処理はどういう目的があるのでしょうか？ぼけて普通よりもリアルに見えて面白い、それで終わりです。物理的に計測をするわけではなく、ざっくり『それっぽく』見えればそれでいいという類のものです。そのような目的の下では、下手に物理的正確さを求めて処理時間を長くするよりは、計算速度の早いフェイクのほうがより適しています。

だから何も指数関数に制限する理由はあまりないわけで、より早くてよりそれっぽく見えれば、そちらのほうが優れていると、そう言えるわけです(ただ、今回の計算量からすると最も大きなオーダーとなっているのは畳み込み積分の箇所で、単調増加関数自体のコストはそれに比べれば微々たるものでしょう)。