映像奮闘記

URL:http://mglab.blogspot.com/
RSS:http://mglab.blogspot.com/feeds/posts/default

[変更]
最終更新日2014/03/02 21:36:00

タグ

映像 

記事

Numpy.hpp - .npyファイルを手軽に扱うためのヘッダライブラリ

2014-03-02 21:36:00
Introduction ご存知の通り,C++は画像処理などの大規模計算に適した言語であり,Pythonはこれらの結果を簡単にscipyやmatplotlibを使って可視化できます.そのため,多くの研究者はコアの計算をC++で行う一方で,データ整理や実験結果の可視化にはPythonを用いるアプローチが一般的でした. その際に重要となってくるのは,Numpy標準のファイル形式である.npyファイルを読み書きするための,C++ライブラリの存在です.これまでにもlibnpyなどに代表されるいくつかのライブラリが存在していましたが,予めコンパイルする必要があるため気軽に使えない,またMSVCなどの環境で使用できないなどのいくつかの欠点がありました. Numpy.hppはこれらの問題を解決するために作られたライブラリです.Numpy.hppの特徴は次の2つです. Header-only: Numpy.hppはコンパイルする必要がありません.Numpy.hppを使用するためにあなたがすることは,ただ#include "Numpy.hpp"をあなたのソースファイルに加えるだけです.外部ライブラリに依存しない: C++03のSTL環境だけで完結しているので,MSVCを含むどの環境下でも容易に動作します. Numpy.hppは以下のURLよりダウンロードできます. https://gist.github.com/rezoo/5656056 使用方法 このセクションでは,Numpy.hppの使い方について簡単に説明します. .npyファイルの読み込み template<typename Scalar>void LoadArrayFromNumpy( const std::string& filename, std::vector<int>& shape, std::vector<Scalar>& data); .npyファイルを読み込みたい場合は,aoba::LoadArrayFromNumpy()関数を使用してください. const std::string& filename: 入力先のファイル名 std::vector<int>& shape: テンソルの階数とその長さを格納するベクトルの出力先std::vector<Scalar>& data: データ本体の出力先 std::vector<T>の型Tは指定された.npyファイルの型と同一でなければなりません.もし異なる型を指定した場合には,Numpy.hppは例外を送出します. aoba::LoadArrayFromNumpy<T>()関数の簡単なサンプルは次の通りです: #include <iostream>#include "Numpy.hpp"// Now I assume np.load("file.npy").shape == (4, 5)std::vector<int> s;std::vector<double> data;aoba::LoadArrayFromNumpy("/path/to/file.npy", s, data);std::cout << s[0] << " " << s[1] << std::endl; // 4 5std::cout << data.size() << std::endl; // 20 Numpy.hppはまた,入力先のテンソルの階数が必ず固定である場合に使えるいくつかのショートカット関数を提供しています. #include "Numpy.hpp"int rows, cols;std::vector<double> data;aoba::LoadArrayFromNumpy("/path/to/file1.npy", rows, cols, data);std::cout << rows << " " << cols << std::endl;aoba::LoadArrayFromNumpy( "/path/to/file2.npy", a0, data); // 1st-order tensoraoba::LoadArrayFromNumpy( "/path/to/file3.npy", b0, b1, b2, data); // 3rd-order tensoraoba::LoadArrayFromNumpy( "/path/to/file4.npy", c0, c1, c2, c3, data); // 4th-order tensor .npyファイルの保存 template<typename Scalar>void SaveArrayAsNumpy( const std::string& filename, bool fortran_order, int n_dims, const int shape[], const Scalar* data); .npyファイルを保存したい場合は,aoba::SaveArrayAsNumpy()関数を使用してください. const std::string& filename: .npyファイルの出力先bool fortran_order: このデータの格納順がcolumn-majorである場合は,このフラグをtrueに設定してください.int n_dims: テンソルの階数const int shape[]: テンソルの長さを表す配列のポインタconst Scalar* data: テンソルの実データを指すポインタ aoba::SaveArrayAsNumpy<T>()関数の簡単なサンプルは次の通りです. #include "Numpy.hpp"// Now I assume s.size() == 2aoba::SaveArrayAsNumpy("/path/to/file.npy", false, 2, &s[0], &data[0]); aoba::LoadArrayFromNumpy<T>()と同様に,Numpy.hppはまたいくつかのショートカット関数を提供しています. // the following functions assume that fortran_order == false.aoba::SaveArrayAsNumpy("/path/to/file.npy", rows, cols, &data[0]);aoba::SaveArrayAsNumpy( "/path/to/file2.npy", a0, &data[0]); // 1st-order tensoraoba::SaveArrayAsNumpy( "/path/to/file3.npy", b0, b1, b2 &data[0]); // 3nd-order tensoraoba::SaveArrayAsNumpy( "/path/to/file4.npy", c0, c1, c2, c3, &data[0]); // 4th-order tensor もし何か質問などありましたら,お気軽に私(rezoolab at gmail.com)に連絡ください :) Enjoy!

Numpy.hpp - Simple Header-only Library for handling .npy files

2014-03-02 16:36:00
Introduction As you know, C++ is a suitable language for computing the large-scale problems, and python is a great language in terms of plotting and visualizing these results by using several attractive libraries such as scipy and matplotlib. For that reason, it is quite a common approach for researchers to use both languages – C++ and python – for the large-scale computation and the visualization. The important thing in this case is an existence of a C++ library for reading/writing .npy files (a default file extension of numpy). Although several libraries have been proposed to handle such .npy files, they require precompiling and cannot work well in several environments (e.g. MSVC). To handle these problems, I reimplemented a new C++ library called Numpy.hpp. Comparing with the existing libraries, Numpy.hpp has two advantages: Header-only: you do not need to precompile numpy.hpp. All you have to do is just add #include "Numpy.hpp" into your source file.No external dependencies: numpy.hpp uses only STL libraries in C++03. This means that numpy.hpp can work well in any environments including MSVC. You can download Numpy.hpp from the following URL: https://gist.github.com/rezoo/5656056 Usage In this section I will briefly introduce the usage of Numpy.hpp. Load .npy files template<typename Scalar>void LoadArrayFromNumpy( const std::string& filename, std::vector<int>& shape, std::vector<Scalar>& data); If you wish to load .npy files in C++, please use aoba::LoadArrayFromNumpy() functions. const std::string& filename: the filename of the input. std::vector<int>& shape: the output of the tensor’s rank and its lengths.std::vector<Scalar>& data: the destination of the data. Please note that the type of std::vector<T> must be equivalent to the type of the specified .npy file. Numpy.hpp will send an exception if you specify the different type. A simple example of aoba::LoadArrayFromNumpy<T> function is as follows: #include <iostream>#include "Numpy.hpp"// Now I assume np.load("file.npy").shape == (4, 5)std::vector<int> s;std::vector<double> data;aoba::LoadArrayFromNumpy("/path/to/file.npy", s, data);std::cout << s[0] << " " << s[1] << std::endl; // 4 5std::cout << data.size() << std::endl; // 20 Numpy.hpp also provides several shortcuts that can be used where the order of tensor is fixed. #include "Numpy.hpp"int rows, cols;std::vector<double> data;aoba::LoadArrayFromNumpy("/path/to/file1.npy", rows, cols, data);std::cout << rows << " " << cols << std::endl;aoba::LoadArrayFromNumpy( "/path/to/file2.npy", a0, data); // 1st-order tensoraoba::LoadArrayFromNumpy( "/path/to/file3.npy", b0, b1, b2, data); // 3rd-order tensoraoba::LoadArrayFromNumpy( "/path/to/file4.npy", c0, c1, c2, c3, data); // 4th-order tensor Save .npy files template<typename Scalar>void SaveArrayAsNumpy( const std::string& filename, bool fortran_order, int n_dims, const int shape[], const Scalar* data); If you with to save .npy files, please use aoba::SaveArrayAsNumpy() functions. const std::string& filename: the destination of the .npy file.bool fortran_order: if you store the data as column-major order, please set this flag as true.int n_dims: the rank of the tensor.const int shape[]: the array that contains the tensor’s lengths.const Scalar* data: the raw data of the tensor. A simple example of aoba::SaveArrayAsNumpy<T> function is as follows: #include "Numpy.hpp"// Now I assume s.size() == 2aoba::SaveArrayAsNumpy("/path/to/file.npy", false, 2, &s[0], &data[0]); As with aoba::LoadArrayFromNumpy<T>(), Numpy.hpp also provides several shortcut functions. // the following functions assume that fortran_order == false.aoba::SaveArrayAsNumpy("/path/to/file.npy", rows, cols, &data[0]);aoba::SaveArrayAsNumpy( "/path/to/file2.npy", a0, &data[0]); // 1st-order tensoraoba::SaveArrayAsNumpy( "/path/to/file3.npy", b0, b1, b2 &data[0]); // 3nd-order tensoraoba::SaveArrayAsNumpy( "/path/to/file4.npy", c0, c1, c2, c3, &data[0]); // 4th-order tensor If you have any questions, do not hesitate to ask me (rezoolab at gmail.com) :) Enjoy!

ノンパラメトリック平均場近似と確率伝播法の提案

2013-12-31 23:06:00
2013年も残り少なくなってきました.そこで,今回は研究の締めくくりとして,今年のCVPRで発表した論文について紹介したいと思います.僕の研究は,大雑把に言うとグラフィカルモデルの最適化手法に関するものです.本研究の貢献を簡単にまとめると,次のようになります.グラフィカルモデル内の確率変数が連続的である場合,計算上の理由から,変分ベイズ法ならびに確率伝播法で扱える分布のモデルには制限があります.提案手法を使うことで,非指数型の分布族に関しても変分ベイズ法を効率的に適用できるようになりました.さらに,同様の議論を確率伝播法にも適用することで,確率伝播法をノンパラメトリックに扱うこともできるようになりました.実際の論文では対象のグラフィカルモデルを一階のマルコフ確率場に限定していますが,変分ベイズ法などの別の領域にも容易に適用できます.また,議論を簡単にするため実際の論文では確率変数の離散化を焦点に当てた形で書いていますが,その流れに沿って書いてもあまり面白く無いので,ここではどのような過程でこのような研究に取り組んだのか簡単に書いていきたいと思います.グラフィカルモデルは機械学習の多くの分野で使われてきた確率モデルです.この確率モデルを用いて対象の問題を高精度に解くためには,確率変数の周辺化が重要となってきます.周辺化の計算を直接行うことは困難であるため,大きく分けて2つの方法論で周辺分布を近似的に求めることが一般的です.一つはギブスサンプリングに代表される,大量のサンプルをばらまいてモンテカルロ的に推定解を求める方法,もう1つは平均場近似(変分ベイズ法)や確率伝播法に代表される,変分原理に基づいた近似的推論法を使う方法です.後者は得られる推定解の精度が比較的悪いものの,特定の反復方程式に従って周辺分布を更新するだけで良いため,高速に周辺分布の推定解が求められる利点を持ちます.本研究は後者の推定手法に関するものです.変分原理に基づいた周辺分布の反復計算は,グラフィカルモデル内の確率変数が連続的であるか離散的であるかで異なります.変数が連続的である場合,周辺分布の密度関数は少数のパラメータからなるパラメトリックな連続分布として表現されます.そして,平均場近似,確率伝播法の両方とも,実際の最適化には変分原理に従って導出された特定の反復方程式に従って,対象のパラメータを更新していきます.しかしながら,このような更新は更新後の分布が更新前の分布と同じパラメータで表現される必要があるため,扱える分布モデルの範囲は非常に限られていました.以上の背景から,僕は平均場近似ならびに確率伝播法をノンパラメトリックに扱えるような方法論を新しく提案すれば,良いインパクトになるのではないかと考えました.このような元々の平均場近似や確率伝播法では扱えない分布をどうにかして扱えるよう拡張できないだろうかという要望は多く,これまでにいくつかの試みが提案されています.これらの試みはすべて次の混合分布の仮定を下にしています.すなわち,平均場近似や確率伝播法で扱う近似分布やメッセージを,次の混合ガウス分布で仮定することです.もし上の仮定の下で混合ガウス分布のパラメータに関する反復方程式を導出できたのであれば,対象のグラフィカルモデルがどのような分布であったとしても平均場近似や確率伝播法を適用できる.つまり,平均場近似や確率伝播法の可能性を一気に広げることができます.しかし,このアイデアを直接適用してもあまり上手くいかないことが知られていました.なぜなら,変分原理は計算の段階で混合ガウス分布のエントロピーを計算する必要があるためです.混合ガウス分布のエントロピーは計算困難であるため,このままでは反復方程式を解析的に導出できません.このような問題に対応するための方法として,これまで大きく分けて2つのアプローチが提案されてきました.1つは解析解を出すために,混合分布のエントロピーにJensenの不等式を用いてさらなる近似を行い,解析解を求める方法です[3].しかしながら,この方法論は平均場近似や確率伝播法で導入した近似分布の他に,さらなる近似をおくことになります.もう1つはモンテカルロ法を用いて,混合ガウス分布のパラメータの反復解を力技で求める方法です[1][2].しかしながら,この方法論ではサンプル数に依存せず高速に推定解を求められる,確率伝播法の利点を失うことになります.これ以外にエントロピーの計算をどうにかして解決できないだろうかと思い,思考を巡らせました.前述の通り,そもそもの問題は混合ガウス分布のエントロピーが計算困難であることでした.ガウス分布の裾が別のガウス分布に干渉してしまうため,それぞれ独立した分布として扱うことができないのです.そうであれば,最初の出発点に混合ガウス分布を選ぶのではなく,互いに干渉しない別の分布を選べば良いのではないでしょうか?例えばヒストグラムのような,互いに干渉せず,エントロピーが容易に計算できる混合矩形分布を採用するのは?ここで,hは変数iにおける,s番目の矩形分布を表します.本研究で新しく提案したアイデアは基本的にはこれだけです.この矩形分布の位置とサイズは,重なっていなければ変数空間のどの箇所に配置しても構いません.すなわち,この表現を用いて変数空間の重要な箇所を密に離散化し,重要でない箇所を疎に離散化することもできます.さらに,変数毎に配置する矩形分布の個数は異なっていても構いません.すなわち,変数の重要性に従って,変数の離散化の度合いを変えることもできます. 次に,上の表現を用いて本来の連続的な最適化問題をパラメータαの最適化問題へ変換し,パラメータに関する反復方程式を導出しました.詳しくは言及しませんが,最終的に導出した反復方程式の形は対象の変数が離散的な場合の平均場近似,確率伝播法の形と非常に似通っています.唯一の違いはエネルギー関数に加わる補正項の存在です.この補正項は変数空間の離散化の「非一様性」を補正する役割をもちます.すなわち,提案手法によって,グラフィカルモデルの連続的な変数空間を非一様に離散化し,ノンパラメトリックに扱えるようになったのです.以上,急ぎですが本研究の解説記事を書かせていただきました.興味を持っていただけた方は,次のURLからダウンロードしていただければ幸いです.http://www.vision.is.tohoku.ac.jp/index.php/download_file/view/35/177/167/参考文献[1] A. Ihler et. al., Particle Belief Propagation, AISTATS, 2009[2] E. B. Sudderth et. al., Nonparametric Belief Propagation, CVPR, 2003[3] S. J. Gershman et. al., Nonparametric Variational Inference, ICML, 2012

ビットコインの知られざる技術的魅力とその可能性

2013-11-30 19:33:00
概要ご存知の通り,ビットコイン(Bitcoin)と呼ばれる謎の仮想通貨が今世界中で大きな話題を呼んでいる.日本においてもその流れは例外でないものの,その話題の殆どがネガティブなイメージ(リンデンドルやチューリップの球根など)とともに報道されることが多い.しかしながら,これらの報道によってビットコインをただの投機先の一つと断定することは早計である.ビットコインがもつ本当の魅力は値段ではない.本当の魅力はビットコインの背後に隠れた優れたアイデアと,その可能性にある.本記事では,一般に語られるビットコインの値動きではなく,これまであまり知られていない,ビットコインが持つ別の可能性について簡単に紹介する.なお,本記事は読者が,ビットコインやネットワークに関する初歩的な知識を持つものとして話を進める.ビットコインについての概要はビットコインとは?,もしくはbitcoins.comが詳しい.また余談ではあるが,はじめてビットコインの概要を知った方がもつ典型的な疑問や疑惑のほとんどは,Bitcoin wikiのMythsの項で詳しく言及されている.興味がある方はそちらを参照されたい.イントロダクション技術的な詳細を省くと,一般のユーザにとって,ビットコインは円やドルなどの従来の法定通貨には無い4つの利点があることが知られている.世界中の何処からでも送金できる.国境による制約を受けない.通貨の偽造が非常に困難であることが数学的に保証されている. 銀行口座は決して凍結されない.誰に対しても非常に安い手数料で送金できる.しかしながら,これらの利点は日本にいる一般的な消費者にとってあまり魅力的でない.偽札が氾濫し投資行動も制限されている現在の中国では,決して偽造できず,世界中に送金できるビットコインは大きな利点がある.その一方,日本の大多数の消費者は海外送金には興味がなく,日本円の価値も揺るぎないものと信じられている.中国と比較して,日本でビットコインが相対的にあまり話題にならない理由の一つはそれらの違いにあると考えられる.だが,現在一般に知られる上の利点は,ビットコインがもつ多くの性質のごく一部にすぎない.事実,ビットコインのプロトコルは送金処理だけでなく,借用証やエスクロー,株式取引にも応用できることが分かっている.特に,ビットコインのプロトコルを利用して,株取引や物品の売買などをすべてビットコインのネットワーク上で実現しようと試みる,BitcoinX (Colored Coins)と呼ばれる取り組みがある.以降はこれらの機能に焦点を当てて説明する.ビットコインネットワークビットコインは中央管理が存在しない,P2Pのデジタル通貨である.ビットコインを使うことで,ユーザは自分が持つウォレットのビットコインを,他のユーザのウォレットへ送金できる.「アドレス」「送金」などの用語から,ビットコインが行っている処理はまるでAliceがBobにメールを送信しているかのようにビットコインのデータを直接送っているように見えるが,その実体は大きく異なる.ビットコインの実体は巨大なP2Pの分散型データベースである.例えば,AliceがCarolから貰った1 BTCを,Bobに送金する場合を考える.この場合ビットコインが行っている処理は,AliceがBobに対して1 BTCのデータを何かの形で送信しているのではない.実際は,Aliceがビットコインネットワークの全体に対して,「AliceがCarolから貰った1 BTCのうち,1 BTCをBobに送金する」という内容のトランザクションを宣言している.次に,ビットコインの参加者はこの宣言を検証して矛盾がないことを確認し,分散型データベースに対象のトランザクションを記入する.二重払いなどの矛盾がある場合は,単にそのトランザクションは無視される.分散型データベースに記入されているのはAliceやBobの口座残高ではなく,「XがYから貰ったM BTCのうち,N BTCをZに送金する」といったトランザクションの集合である.興味深いことに,このトランザクションの内容は表現力を持ったスクリプト言語で記述されている.これは,トランザクションに書き込める内容は送金処理といった単純な取引だけでなく,より複雑な取引にも対応できることを意味する.Contractsではその具体例としてデポジット,エスクロー取引,保証契約にも応用できると提案している.トランザクションの内容は誰もが見られる形で保存されるため,契約書の内容は決して偽造されない.そしてこのスクリプト言語を別の形で利用することで,ビットコインを実際の資産売買や権利譲渡にも応用しようとする取り組みがBitcoinX,またはColored Coin(色をつけたコイン)と呼ばれるプロジェクトである.BitcoinX (Colored Coins)BitcoinX (またはColored Coins)はトランザクションのスクリプトを利用して,特定のビットコインにビットコイン以外の価値―例えば債権やコモディティなど―を結びつけるプロジェクトである.この別の価値を結びつける行為のことを,発案者は「色をつける(Colored)」と表現している.BitcoinXを使うことで,ユーザはビットコインの送金と同じ感覚で,安全かつ確実に資産の受け渡しや売買を行える.クリスマスコンサートのチケットを例に説明する.Aliceは有名なクラシックコンサートの主催者で,チケットをビットコインでも発売したいと考えている.この場合,AliceはまずBitcoinXを使って特定のビットコインに「色をつける」.例えば,あるコインには「2013年12月24日 六本木クリスマスコンサート A-20席」,もう一方のコインには「2013年12月24日 六本木クリスマスコンサート S-3席」といった色である.色をつけるビットコインの額自体は少額で構わない.これは原価の安いただの紙切れが,コンサートの入場権という価値を持つことで高い値段が付けられる現象と似ている.ビットコインに色をつける行為自体は誰でもできるものの,誰が色を付けたのかという署名も一緒に記入される.この署名の偽造が困難であることは数学的に保証されているので,攻撃者によるチケットの偽造も同様に困難である.次に,Aliceはこの色のついたビットコインを販売して,多数の顧客に送信する.顧客が入場の際に行う操作は,単にこの色がついたコインを携帯のウォレットに入れ,入場口の端末にかざすだけでよい.なぜなら,分散型データベースにはこの色のついたコインが最終的にどのウォレットに入っているのか,誰もが見られる形で記録されているためである.ここで,顧客の一人であるBobは急な仕事でコンサートへ行けなくなったため,手持ちのチケット(色のついたコイン)を売却したい場合を考える.このとき,Bobはチケットを購入したい第三者を見つけ,スクリプトを利用したエスクロー取引で安全にビットコインとチケットを交換する.このとき実際に行われるのはビットコインネットワークを使用したビットコインの交換だけであるので,手数料は非常に安く抑えられる(詳細はContractsを参照されたい).このアイデアは単純だが幅広い応用が考えられる.例えば債権の売買や,オークション,株式の発行などである.これまで,ユーザはこのような取引を安全に行うために,国や企業などの仲介者に対して多額の手数料を払う必要があった.BitcoinXを使用した場合,ユーザはこのような手数料を払う必要はない.かつ,これらの取引は誰もが見られる形で公開されるため,偽造などの不正は困難である.さらにこのアイデアを発展させて,特定のビットコインと実際の通貨を結びつけることができる.ここで,ある企業が特定のコインに「100円」という色をつけたとする.このコインを実際にその企業の店頭に持ち込めば,100円と交換してもらうことができる.すなわち,このコインは100円の価値を持つ(この色をつけるBTCの額は,市場で取引される日本円の価格と対応していないことに注意されたい.すなわち,色をつけるBTCの額は,コンサートチケットと同様に少額で構わない).結果として,日本円はビットコインと同じような形で流通させることができるようになり,かつその流通にかかる手数料は非常に安く抑えられる.もちろん,ドル,ユーロ,ウォンなどの他の通貨も,同様の手続きでビットコインネットワークに流通できる.もはやビットコインは仮想通貨をやりとりするためだけの存在ではなくなった,ビットコインは仲介者が存在せずとも,安全かつ確実にあらゆる取引を遂行するためのプラットフォームにも成り得るのだ.まとめビットコインはよくバブルと共に語られることが多いが,ビットコインの本当の魅力は値段でなく,その背後にある技術と可能性である.その可能性の一部として,本記事ではビットコインを介したエクスロー取引や,資産取引の具体例について説明した.ビットコインを使うことで,ユーザはあらゆる取引を,仲介手数料を殆ど取られることなく安全かつ確実に遂行できる.現在,BitcoinXは実際のBitcoinクライアントには統合されていない,開発途中のプロジェクトである.この開発を積極的に進め,将来的にはモバイルなどの端末機器でも使えるようにすることが求められる.また,課税問題やマネーロンダリングなどといった,解決すべき法的課題も多い.しかしながら,将来的にはこれらの懸念は段階的に解決されるであろうと考えられる.ビットコインはあらゆるトランザクションが公開されるプラットフォームであるため,将来的にはむしろ,従来よりも透明性の高い資産取引とその課税ができるようになるだろう.ビットコインの未来は誰にも分からない.しかし,匿名のエンジニア―Satoshi Nakamoto―が発表したこのアイデアが画期的であることは疑いようがない.たとえビットコインの価格が崩壊しその価値を失ったとしても,この先進的な技術を他の分野へ応用することは重要であると考えている.

条件付き確率場の推論と学習

2013-04-26 00:03:00
条件付き確率場の推論と学習 from Masaki Saito  前回の研究室ゼミで話した,条件付き確率場についての簡単な紹介スライドを公開します.内容は基本的なPairwise CRFについての話と,それが実際にどのような問題に対して応用がなされているか,またもっとも使われる最適化手法の一つである,平均場近似と確率伝搬法についての簡単な解説と,CRFの初歩的な学習法といった流れです.補足: CRFの学習は自然言語処理など他の分野でもよく使われている学習法なのでおそらく知っている方は多いと思うのですが,普通はダイレクトに対数尤度を最大化することでパラメータを求めています.しかし,今回はKL距離という確率分布間の近さを図る指標を用いて対数尤度を拡張し,パラメータを求めています.なぜわざわざそうするのかというと,だいたいの最適化手法がKL距離の最小化という議論に帰着できて,すっきりするからです.KL距離を使えば平均場近似,確率伝搬法,最尤推定,あと時間の都合上話せなかったTRWやGeneralized Belief Propagationなどの最適化手法がすべて統一的に理解できます.これによって覚えることが少なくなって楽できるだけでなく,このKL距離を使った新しい最適化手法を提案することもできます.

近況報告

2013-04-25 02:48:00
お久しぶりです.細々とした近況はTwitterに書いているのですが,あくまでTwitterはつぶやきであって,記事でない.定期的に自分の中でまとめた内容をブログに吐き出す必要がありますね.修士論文を提出し,晴れて今年度から博士後期課程の1年生です.不安が全くないと言えば嘘になりますが,それでも毎日好きなことを学び,研究できる環境は非常に楽しい.みなさんもぜひ博士課程に進学しましょう.唯一不満に思うところは頻繁に書類を書く業務が入るあたりですが,仕方のないことです.後できちんとした内容を書きますが,2013年のCVPR(コンピュータビジョンのトップカンファレンス)に採択されました.去年を含めると今回で2回目の採択ですが,前回と異なるのが,今回はポスターではなくオーラルセッションでの採択というところです.オーラルセッションは非常に狭き門(採択率 60/1870 = 3.2%)でありますので,そのような研究成果を得られたことは非常に嬉しいですし,来年もぜひ通しておきたいところです.博士後期課程になりましたので,これからは本名メインで活動していきます.プロフィール画像も更新しました.前に使ってた画像は高校時代からのものでしたから,大体7年位使っていたんでしょうか.というかこのブログ,今年で7年目になるのですね.自身の宣伝もかねて,論文やコードなどはこれから作る自己紹介のページにて積極的に公開していきます.論文についてはできるだけ読んでもらえるよう,簡単な解説も記事にしていきたいと思っております.

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

2012-10-20 21:19:00
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を用いて実装する

2012-08-21 23:58:00
ご存知の通り、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の学習手法についての簡単なまとめ

2012-08-18 22:58:00
近年の機械学習では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)を用いた近似最近傍探索についての簡単な解説

2011-11-24 20:11:00
"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倍以上の圧縮効率です。このように、少数のコードブックサイズで効率よくベクトルを表現できることが直積量子化の利点です。距離計算さて、直積量子化によって各ベクトルのエントリを圧縮した状態で格納したとしましょう。あるクエリベクトルが与えられたとして、そのユークリッド距離が最小となるエントリを見つけるにはどうすればいいでしょうか?幸運なことに、ユークリッド距離の計算は各要素毎に独立しています。つまり、ただ単にクエリベクトルを各サブベクトルに分け、それぞれの領域で距離を計算し、その和を計算するだけで良いのです。よって、直積量子化を用いて最近傍探索を近似的に行うためには:クエリベクトルをM分割し、N/M次元のサブベクトルに分解するクエリのサブベクトルと、それぞれのコードブックに記された代表ベクトル間の距離を予め計算しておくそれぞれのエントリがどのインデックスに位置しているのか確認し、2で計算した距離を参照して足し算を行う最小距離をもつエントリを候補として提出することで実現できます。これが直積量子化を用いた、近似最近傍探索の核となるアルゴリズムです。詳細についてこの距離計算は非対称性からある程度のバイアスがかかっているため、実際の距離を正確に見積もるためにはこのバイアスを取り除く必要があったりします。また、具体的な距離計算がコードブックだけで済むので省力化されているとはいえ、計算量は依然としてO(N)であり、そこらへんも大規模なデータ数では足を引っ張る原因となります。んじゃあこれをどうやって解決するのかについては、現論文http://lear.inrialpes.fr/pubs/2011/JDS11/jegou_searching_with_quantization.pdfを参照してください。上記2つの懸念についての解決案が記述されており、上の知識があれば多分すらすらと読めるはずです。本当はこれも解説しなきゃなんないんですけど、書いている内に疲れたちゃったのでとりあえずコアなとこだけということで勘弁してください。気が向いたら更新します。

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

2011-09-24 01:44:00
"The Dunny Collective - Brisbane River" by KayVee.INC概要OpenCVを用いてアルゴリズムを組む際、ある画像からある画像へ変換処理を行ったり、関数fを用いて新規に画像を生成したり、画像の結果を1つの値へと集約するなどの定型的な処理は非常に多く行われておりますが、並列処理を行いかつ速度を考慮したプログラムを組もうとするとある程度のコード量を必要とするためそれなりの手間となります。また画像処理によくあるループ文のネストが大量に発生するため、普通に組んでしまうと一見何をやっているのか非常に分かりづらいコードとなってしまいます。cvutilはこのような冗長な作業を定型化し、複数のアルゴリズム関数とC++0xに新しく搭載されたラムダ式を使用して高速にアルゴリズムを記述することができます。このようにして作られたアルゴリズムは自動的にOpenMPによって並列化されるため、研究用途などの速度をある程度担保しつつもメンテナンスのし易いコードを記述するための手助けとなります。インストールcvutilはgithubにて公開しております。ダウンロード先は下記のURLを参照してください。https://github.com/rezoo/cvutilcvutilはヘッダオンリーのライブラリです。バージョンがOpenCV 2.x系統であれば、対象のincludeディレクトリにcvutilをコピーするだけで使用できます。C++0xのラムダ式を用いることで端的に記述することができますが、関数オブジェクトを使用すればある程度古いコンパイラでも正常に動作するでしょう(そこまでして使う利点があるのかどうかは分かりませんが)。ちなみに、元々は自分用に作られた適当なライブラリなので、バージョンアップに伴ってある程度破壊的な変更を行うこともありますのでご了承ください。その際はバージョンを0.1上げることによって対応します。アルゴリズム一覧ここではcvutilに搭載されているアルゴリズム一覧の一部を紹介いたします。これらのアルゴリズムはSTLのそれと似ているため、C++を触ったことのある方は名前から直感的に機能を類推できるでしょう。transform_imagetransform_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_imagegenerate_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); }); // BGRreduce_imagereduce_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_imagetransform_reduce関数は1枚または複数の画像を対象の関数で1つの値へ変換した後に、reduce関数を適用します。2種類の結果を合併した集約結果を返す場合などに役立つと思います。minmax_element_imageminmax_element_image関数は画像の最小値と最大値を返します。std::pair<uchar, uchar> minmax_result = cvutil::minmax_element_image(src); // (min, max)integrate_imageintegrate_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の画素値アクセスについて

2011-09-16 14:45:00
"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 forfor(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側の意図しているところだろうと判断し、上の式を書き換えました。

boost.gilを利用して透過画像を作る

2011-06-26 01:20:00
『対象の画像から、特定の色の箇所だけを透明にした透過画像を作る』こういった課題に一ヶ月近くかかって出来なかったとかいう話を聞いて、流石にそれはかかり過ぎだろうと脊髄反射で言ってしまったわけですけど、それじゃあ自分だったらどうだろうとか後で思ってちょっとやってみることにしました。好き勝手言って出来なかったとか恥ずかしすぎですから、やっぱりそこらへんは実際にやってみる責任みたいなものは感じているわけです。というわけでboost.gilを使って書いてみた結果が以下。#include <boost/gil/gil_all.hpp>#include <boost/gil/extension/io/png_io.hpp>using namespace boost::gil;template<typename Pixel>struct transparent_functor {public: explicit transparent_functor(Pixel bg_pixel): bg_color_(bg_pixel) {}; Pixel operator()(const Pixel& src) { if(src == bg_color_) { Pixel dst = src; get_color(dst, alpha_t()) = 0; return dst; } else { return src; } }private: Pixel bg_color_;};int main() { const char* filename = "src.png"; const rgba8_pixel_t background_color(0, 255, 0, 255); rgba8_image_t src_image; png_read_image(filename, src_image); rgba8_image_t dst_image(src_image.dimensions()); transform_pixels( const_view(src_image), view(dst_image), transparent_functor<rgba8_pixel_t>(background_color)); png_write_view("dst.png", view(dst_image)); return 0;}結果としてはこんな感じです。上が元の画像で、下が結果の画像。正常に抜き出されていることが分かります:余談プログラムより文章を紡ぎ出すのに時間がかかって、改めて自分の文章力の無さにうんざりしていますというか最近暑い!自分は寝たいのに梅雨ですごい蒸し暑いから寝付けず気になってこんな不毛なことやっちゃってるわけです。今度こそ寝ます!!!

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

2011-06-02 00:40:00
"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] -> amax_element list = foldr max (head list) listmin_element :: Ord a => [a] -> amin_element list = foldr min (head list) listminmax_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_ofall_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] -> Boolall_of binary_pred list = foldr (&&) True $ map binary_pred listany_of :: (a -> Bool) -> [a] -> Boolany_of binary_pred list = foldr (||) False $ map binary_pred listnone_of :: (a -> Bool) -> [a] -> Boolnone_of binary_pred list = all_of (not . binary_pred) list関数合成を用いてエレガントに解決。inner_product, transform_reduceinner_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*2C++だといろいろ冗長に書かないといけないですけど、Haskellだったら余裕です:inner_product :: (a -> b -> c) -> (c -> c -> c) -> c -> [a] -> [b] -> cinner_product binary_op1 binary_op2 init list1 list2 = foldr binary_op2 init $ map operate $ zip list1 list2 where operate (x, y) = binary_op1 x ytransform_reduce :: (a -> b) -> (b -> b -> b) -> b -> [a] -> btransform_reduce unary_op binary_op init list = foldr binary_op init $ map unary_op listcount_if, countcount_ifは条件に合致している要素数を返す関数で、countは指定した値と同じ要素の数を返す関数です。*Main> count_if (\x -> x < 3) [1,2,3,4,5]2*Main> count 3 [1,3,2,2,3]2count_if :: (a -> Bool) -> [a] -> Intcount_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 0count :: Eq a => a -> [a] -> Intcount value list = count_if (== value) listThrustではTrueだったら1, Falseだったら0となるリストを作成して、それを足していく手法をとっていました。countはcount_ifが思いつけば楽勝でしょう。find_if, find_if_notfind_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] -> Intfind_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] -> Intfind_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の初歩的なスクレイピング

2011-05-18 01:55:00
前回はぽんとスクリプトを上げて「これで画像がダウンロードできますよ」とやったわけだけど、もちろん実際に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")を検索。ビンゴ!これで駒はすべて揃った。まとめに入る。対象の絵画をダウンロードするためには、取得したい絵画のアドレスへアクセス<div id="microscope" data-microId="unique-id">となっている項を取得、unique-idを得るhttp://lh3.ggpht.com/ - http://lh6.ggpht.com/のいづれかにアクセス。http://lh3.ggpht.com/unique-id=xi-yj-zk のような形で画像を取得取得した複数の画像を繋ぎ合わせる後は非常に簡単でただこの単純なアルゴリズムに従って書き下していくだけ。その産物がこいつ。こんなに早く仕上がるってことは、たぶんGoogleさんは本格的にダウンロードを禁止したいわけじゃないんだろう。スクレイピングとしては楽な部類で、普通だともうちょっと複雑な認証をかましてくるので自動化しにくい。Captcha入ると流石に無理だけど、おそらく現在の文字認識技術をきちんと応用すれば解けるんじゃないかと思っている。余談こういう解析は楽しい。どんな感じの楽しさかっていうと、与えられたパズルを解いていく感覚にちょっと近い自分はだいたい1時間くらいで解いたけど、早い人はたぶんもっと早いんだろうなって思う。まだまだ精進が足りない

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

2011-05-16 01:20:00
"The Starry Night" by Vincent van Gogh2011年2月2日、GoogleさんはGoogle Art Projectというサイトをリリースしました。どういうサイトかっていうとMoMAニューヨーク近代美術館やTate Britainなど、世界中の美術館にある有名な絵画を渡り歩くことができるという趣旨のサイトで、誰もが一度は見たことがある絵画がー家に居ながらー非常にリアリスティックな解像度で見ることができる。恐ろしい時代になったもんです。ところがこのサイト、高解像度で見られるのは非常にいいことなのですが、肝心のダウンロードができない。せっかく高解像度なんですから保存できてもいいはずなのに、だけどできない。Google Mapなんかは常に最新の情報に更新しなければならないので、ダウンロードに制限をかけているのは合理性があるんですけど、絵画とか全く更新する必要がないでしょう。著作権も切れてるのに。よく分からないです。ということで、なければ自分でなんとかしようという精神で、1時間くらいでちゃっちゃとPythonでスクリプト組んでダウンロードするようにしてみました。即興で組んだのでいろいろ粗い面はありますけど、そこら辺は勘弁してください。#!/usr/bin/env python#-*- coding: utf-8 -*-import sysimport reimport urllib2import randomfrom PIL import Imagedef usage(): print "%s [URL]" % sys.argv[0] return 1def 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 Truedef 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 xdef 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 ydef 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アルゴリズム

2011-04-19 01:28:00
"fireking" by cmyk12191. 概要reduceとscanはGPGPUなどの並列環境下において重要な役割を果たす関数だけど、あまりこれらの関数、特にscanについて言及している記事はあまり見かけない。なので今回はreduceとscanについて調べた結果をまとめていきたいと思う。2. reduce, scan関数の概要2.1. reducescanについては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])120reduceにおいて重要なのは二項演算子であり、上記の場合、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. scanscan関数は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 scanl1scanl1 :: (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を行いたい場合に役に立つ。リスト[a0, ... , an]にscanオペレーションを実行してを手に入れたとする。この場合、i < jにおいてbj+bi-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の出しているテキストは分かりやすい。すごいことです。

Moving

2011-03-31 19:17:00
本当は4/1に書く予定でしたが、そういえば4/1はエイプリルフールでしたので今日書きます。えーと、本日、東北大学理学部物理学科を卒業し、明日から東北大学大学院情報科学研究科に入学します。所属する研究室はコンピュータビジョンなど、主に視覚的な情報について研究している研究室ですので、たぶんこれから書く記事はそれ系統の内容が多くなると思います。このブログのタイトルに少し近づきました。元々このブログは映像を作り、その技術的な情報についてログを残せるよう立ち上げたブログでした。確かに現在は映像の「え」の字もないブログに仕上がっていますが、もちろん今でも映像は好きですし、好きなプロダクションさんの映像はいつも観に行くようにしています。そして全く役に立たなかったかと言われればそういうわけではなくて、個人的に仕事を頼まれることもありましたし、Webサイトや印刷物のデザインをする際にも役立ちましたし、ここで得た知見のお陰で様々な方とお話することもできました。物理学科に入ったのも元々は相対論や量子力学の深遠な世界について知りたいために入りました。正直に言って完全に理解できたかというと怪しいところでありますし、理学から工学という他分野に移るわけで、一見すると今までの知識が役立たなくなるように見えます。が、実はそうでもなく、統計力学(*1)や相対論で使っている手法をコンピュータビジョンに輸入している論文があったり、案外様々な場面で役立ったりしているわけで(*2)。全然違う領域に見えても絡んでくる領域というものは少なからずあるものです。ということで、映像と物理がお互いに絡み合っている、コンピュータビジョンの道に進むことを決めました。これからもどうぞよろしくお願いします。p.s. 家族は無事でした、よかったよかった*1 ここで宣伝ですが、田崎先生の『統計力学I, II』はお薦めです。この本はぜひ理学系だけでなく工学系の方々にも読んでもらいたい本です。値段も2冊併せて7,000円くらいと安い。新品のゲームソフト1本買うくらいならこれを買いましょう*2 もちろん線形代数や基本的な偏微分も数式や概念を学んだりする上で役立ちました。一番学んで良かったのはブラケット記法とヒルベルト空間、あれは素晴らしい

自分は無事です

2011-03-11 20:51:00
自分は出張で京都にいましたので、無事です。もし身内の方が見ていましたらご安心ください。ただ暫くは京都に滞在することになると思います。早く交通機関が復旧してほしいです。ただ家族がどうなっているのか全く分からない。怖いです、不安です。[文句]iPhoneではSoftbankの災害用伝言ダイヤルに『登録できない』!!!明らかにsoftbank側の怠慢、糞だ。糞すぎる

OMakeのContributionsに載ったった

2011-02-16 02:22:00
The OMake build system: User Contributionsやったね!余談翻訳をすることは勉強にもなって非常にいいですし皆さんにもぜひやって欲しいのですが、その分、なんというか、日本の中へと篭ってしまうことにも繋がっている気がします。英語に臆することなく積極的にアウトプットできるよう、自分自身の心構えをもう少しきちんとしなければならんですね。

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

2011-02-15 03:54:00
"Image dans le néant" by gelinhboost.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.

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

2010-12-18 18:25:00
『画像をぼかす』と一口に言っても、その方法にはいろんなものがあります。一番有名なのは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つ目が、『全てのデジタルカメラにおいて、ピクセルは輝度の対数の次元を取っているのだろうか?』という単純な疑問です。きちんと物理的/工学的な意味合いを持つ画像処理にするためには、格納されているピクセルデータの次元がきちんと揃えられている必要があります。自分自身も普通に考えて輝度の対数が格納されているんだろうということで、任意ではなく指数関数に制限された関数系を使いましょうと運びたかったのですが、そうするための証拠が不足していた。もしかしたら自分の予想もつかない方法で処理されている可能性があるわけで、そういった証拠が十分にないことから、単調増加関数と少し曖昧とした意味合いを含ませたのです。もう一つの理由は『今回の画像処理は物理的に正確である必要がない』という理由です。これが大きかった。そもそも今回の画像処理はどういう目的があるのでしょうか?ぼけて普通よりもリアルに見えて面白い、それで終わりです。物理的に計測をするわけではなく、ざっくり『それっぽく』見えればそれでいいという類のものです。そのような目的の下では、下手に物理的正確さを求めて処理時間を長くするよりは、計算速度の早いフェイクのほうがより適しています。だから何も指数関数に制限する理由はあまりないわけで、より早くてよりそれっぽく見えれば、そちらのほうが優れていると、そう言えるわけです(ただ、今回の計算量からすると最も大きなオーダーとなっているのは畳み込み積分の箇所で、単調増加関数自体のコストはそれに比べれば微々たるものでしょう)。

OMakeマニュアル 日本語訳 PDFバージョンをリリースしました

2010-10-27 02:03:00
"Aint they cute?" by Nicki Varkevisser題名の通りで、以前から取り組んでおりました『OMakeマニュアル 日本語訳』の全ページを一つに纏めたPDFバージョンをリリースしました。http://omake-japanese.sourceforge.jp/OMake.pdfよりダウンロードできます。電子書籍でまとめてご覧になりたい方や、職場などの環境で使いたいけれどネットが使えないのでごりごり印刷したい…といった場合に便利です。元々は渋川さんの『LaTeX経由でのPDF作成』という記事を拝見しまして、試験的にPDF版を作成してみたところ、思いのほか綺麗に製本されていることからきちんとリリースしてみようと思い立った次第です。殆ど素の状態から手を加えていないので暫定的ですが、後はきちんとスタイルを組んでいって、もうちょっと読みやすいレイアウトにしていけばいいかなって感じです。思えば個人的な暇潰しから始まったこのプロジェクトも現在はPDFにして約200ページにもなり、きちんと一つの事柄に集中して取り組めたことは、学部時代のいい経験になったと考えています。まあ、これからも暇があればちょくちょくレイアウト、文章含めて修正していきますので、どうぞ今後ともよろしくお願いします。

確率微分方程式を解くー幾何ブラウン運動、バシチェックモデルの場合

2010-10-04 01:58:00
"New York Fashion Week Fall 2007: Doo Ri" by Art Comments確率微分方程式という分野がある。どういうものかっていうと、要は確率過程を微分方程式によって表そうという試みであり、微分方程式の中に微小の確率変数が含まれている。具体的にはこんな感じ。このdWはブラウン運動を表していて、時間がたつにつれてうねうねと確率的に動く。期待値と標準偏差はそれぞれ0と√tで、時間がたつにつれて動く場所もどんどん広がっていく。この運動はランダムウオークとか酔歩とも呼ばれていて、酔っ払ったおっさんがうねうねとあちこちに歩いているように見えることからこういう名前がついたそうな。上のような式で表すことができる確率過程を一般には『伊藤過程』と言って、他にもストラトノビッチ確率解析とかいうのもあるらしいけど、詳しくやっていないので分かりません。これの何が便利かっていうと確率過程を通常の微分方程式のように扱えるところで、確率変数をまるで微分方程式を解くかのように求めることができる。とは言っても通常の微分方程式にはない規則がいくつかあって、初めはすごく奇妙に思えるかもしれない。具体的には、以下のように確率変数x(t)とy(t)を定めた場合、このy(t)もまた伊藤過程に従い、確率微分方程式dyは以下のルールに従う。普通だったら消えるはずの2乗の項が残ってしまう。それじゃあ(dx)^2はどのようにして計算するのかというと、以下の公式に従う。なんと、(dW)^2がdtに変化してしまう!これを『伊藤の公式』という。なんでそうなるのか?それを知りたい方は専門書を漁ればいくらでも載っているので、そこらへんを参照すればいいと思う。ここが通常の微分方程式と違うところで、逆を言えばそこさえ守ればきちんと積分によって求めることができる。専門書だとこの箇所がすごく分かり辛い表現で説明されているけど、詳細をはしょればこんな感じなはず。それじゃあ、まず手始めに幾何ブラウン運動と呼ばれるモデルを求めてみる。これは株価の変動を表す最も簡単なモデルなんだけど、何故これが株価の変動を表しているのかっていうのは、とりあえず解いてみれば判明するはず。さて、普通の考えだったら『変数分離法を用いてxを分離させて、そのまま積分すればよくね?』って発想になるんだけど、何度も言っているようにこれは確率微分方程式なので、まずは伊藤の補題を用いて変数変換を行い、xを分離させてみることにする。y = ln(x)とおくと、伊藤の補題よりとなる。(-1/2 σ^2)の項が出てくるところがミソで、普通に変数変換してしまうと出現しない。これで定数部分だけに分離できたので、そのまま愚直に積分を行うと、となる。これは対数正規分布の形となっているため、株価の変動を表す簡単なモデルとなっていることは容易に理解できると思う。それでは、次にVasicekモデルと呼ばれるモデルについて求めてみる。これは主に国債金利の利子率の時間的変動を表すモデルで、いったいこいつの何が重要なんじゃいっていうと、金融商品のリスクを推定する際に使われる。金融商品の価格は利子率によって変動するので、利子率の変動を推定することは金融商品のリスクを推定することに繋がる(*1)。つまり、今の自分のポートフォリオがどれくらいのリスクを持っているのか判断できる。そもそもこの記事を書くきっかけになったのはこいつが原因で、wikipediaさんにはただ『積分するとこれになりますよ』って書いてあるだけで、どうやって積分するのかについては書いていなかった。しょうがないので自分で手を動かして、ようやっと理解できたので公開してみようと。さて、まずはu=b-rとおいて変数変換すると、となる。ただ普通にChain Ruleやってるように見えるけど、たまたま結果が同じだっただけで、どんな変数変換でもきちんと伊藤の補題を使わなきゃいけない。そんで、物理やってる人はピンとくるかもしれないけど、これはランジュバン方程式の形となっているので、積分には同様の手法を利用できる。まぁ自分はピンとこなかったんですけど。それはさておき、v=e^(at)uとおいてさらに変数変換してやると、となり、いい感じに余計な項が消えてブラウン運動だけになった。よってこれを積分してやるとあとはごちゃごちゃ整理してやるとこれでwikipediaさんにあるのとおんなじ形になった。めでたしめでたし。なんかブログに書くと導出するのは楽ちんなように思えるけど、結構悩んだし大変だった。どんなもんでも答えを見るとすぐ分かるし思いつきそうなんだけど、実際になにもない地点から思いつくかどうかっていうのはまた別な話で、それを考えると新しい学問を生み出すっていうのは凄まじいほどの労力が必要だ。実際、未だにどういう発想で量子力学を思いついたのか本気で理解できない。あんなのヒントがいくらあっても無理だ。そういうことを考えた一日でありましたとさ。*1: もちろんこの説明は適当でないが、詳しく説明してもあんまり意味がないのでこう書くことにする追記(2010/10/04)1/2の項が抜けていたので修正。

マルコフ連鎖モンテカルロ法を実装してみよう

2010-09-28 18:31:00
"Chain Bridge, Budapest" by szeke約2ヶ月ぶりということで、ご無沙汰しております。書きたいネタというのは結構あって書こう書こうとは思っているんですが、なにより書くとなるとあれもこれも伝えなきゃいけないみたいな思いになって、結局分量の多さから諦めてしまうというのが結構続いています。もう少し気を張らずに更新していきたいものです。さて、最近の自分はマルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo, MCMC)についておりました。何にも知らない方ににとってはよく分からないかもしれませんが、この手法はマルコフ連鎖が持つ簡明さとモンテカルロ法が持つ実用性が合わさった凄まじい手法なんです。そしてなによりとてもエレガント。とりあえず知らない人のためにてきとーに解説しますと、与えられた適当な関数から確率変数をサンプリングするための公式です。べつにこれじゃなくても確率変数のサンプリングにはいろんなアルゴリズムがあるわけですが、これを使うと関数の計算コストがやけに高い場合だとか、ピーク部分を集中的にサンプリングしてくれたりとかで結構便利で、金融や工学、物理の分野でも非常に役立っております。計算物理をやる人にとっては必須に近い手法といっても過言ではありません。そんなMCMCですが、実際に提案されている手法には様々なものがあります。不変分布から遷移確率を逆算するため、自由度が増えるから遷移確率は一意じゃなくなってしまうんです。今回は最も有名なMetropolis-Hastingsアルゴリズム(*1)を用いて、任意の関数に従って乱数のサンプリングを行うコードを練習がてらpythonで書いてみました。#!/usr/bin/env python#-*- coding:utf-8 -*-import sysimport mathimport randomdef main(): f = lambda x: math.exp(-x*x / 10.0) for x, fx in metropolis(f): print x, fxdef metropolis(func, start=0.0, delta=2.0, burn_in=100, seed=1): metro_iter = _metropolis(func, start, delta, seed) for i, value in enumerate(metro_iter): if i > burn_in: yield value break for value in metro_iter: yield valuedef _metropolis(func, start, delta, seed): random.seed(seed) initial_x = start initial_fx = func(initial_x) proposed_x, proposed_fx = 0.0, 0.0 while True: proposed_x = initial_x + random.uniform(-delta, delta) proposed_fx = func(proposed_x) ratio = proposed_fx / initial_fx if ratio > 1.0 or ratio > random.random(): initial_x, initial_fx = proposed_x, proposed_fx yield proposed_x, proposed_fx else: yield initial_x, initial_fxif __name__ == '__main__': sys.exit(main())うーん・・・きれいだ(*2)・・・こいつのすごいところは、アルゴリズム自体は非常にシンプルで実装しやすいというところにあります。行数も上を見ると分かるのですがとっても少ないし、コア部分はたった25行くらいで構成されています。なのに使っている学問は結構複雑で、定常分布から遷移行列を導出するために詳細釣り合いの条件からMetropolis的制約を用いて導出、サンプリングを行っております。このアルゴリズム自体は結構知ってる人が多いかもしれませんが、実際に理論方面で理解している人はどうなんでしょう。んで、正規分布のガウス関数からサンプリングを行った結果がこれ。 赤線が元の確率分布。きちんとサンプリングが行われているようです。今回は正規分布で計算を行いましたが、非負の関数であればなんでも構いません。あとはこれを研究テーマへと落とし込むことができれば、とりあえず一区切りつけるかなって感じです。*1: 提案確率関数が対称であるという制約を持つのがMetropolis。それを非対称の場合にも拡張したのがHastings。個人的にはMetropolis-HastingよりMetropolisのほうが綺麗で好きです。というか提案確率密度関数が非対称じゃないといけないようなシチュエーションがよく分からない*2: yieldとジェネレータって便利ですよね。計算速度も効率的だし、慣れてくるとすべて遅延評価で行わせたいくらい