tag:blogger.com,1999:blog-10622833001158220262024-03-05T20:27:02.334+09:00映像奮闘記ちょっと映像関係でてきましたrezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.comBlogger169125tag:blogger.com,1999:blog-1062283300115822026.post-42228275189963882962014-03-02T21:36:00.001+09:002014-03-02T21:36:42.499+09:00Numpy.hpp - .npyファイルを手軽に扱うためのヘッダライブラリ<h2 id="introduction">Introduction</h2>
<p>ご存知の通り,C++は画像処理などの大規模計算に適した言語であり,Pythonはこれらの結果を簡単にscipyやmatplotlibを使って可視化できます.そのため,多くの研究者はコアの計算をC++で行う一方で,データ整理や実験結果の可視化にはPythonを用いるアプローチが一般的でした.</p>
<p>その際に重要となってくるのは,Numpy標準のファイル形式である<code>.npy</code>ファイルを読み書きするための,C++ライブラリの存在です.これまでにもlibnpyなどに代表されるいくつかのライブラリが存在していましたが,予めコンパイルする必要があるため気軽に使えない,またMSVCなどの環境で使用できないなどのいくつかの欠点がありました.</p>
<p><code>Numpy.hpp</code>はこれらの問題を解決するために作られたライブラリです.<code>Numpy.hpp</code>の特徴は次の2つです.</p>
<ul>
<li><strong>Header-only</strong>: <code>Numpy.hpp</code>はコンパイルする必要がありません.<code>Numpy.hpp</code>を使用するためにあなたがすることは,ただ<code>#include "Numpy.hpp"</code>をあなたのソースファイルに加えるだけです.</li>
<li><strong>外部ライブラリに依存しない</strong>: C++03のSTL環境だけで完結しているので,MSVCを含むどの環境下でも容易に動作します.</li>
</ul>
<p><code>Numpy.hpp</code>は以下のURLよりダウンロードできます. <br>
<a href="https://gist.github.com/rezoo/5656056">https://gist.github.com/rezoo/5656056</a></p><div class="se-section-delimiter"></div>
<h2 id="使用方法">使用方法</h2>
<p>このセクションでは,<code>Numpy.hpp</code>の使い方について簡単に説明します.</p><div class="se-section-delimiter"></div>
<h3 id="npyファイルの読み込み"><code>.npy</code>ファイルの読み込み</h3><div class="se-section-delimiter"></div>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="kwd">template</span><span class="pun"><</span><span class="kwd">typename</span><span class="pln"> </span><span class="typ">Scalar</span><span class="pun">></span><span class="pln">
</span><span class="kwd">void</span><span class="pln"> </span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="kwd">const</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">string</span><span class="pun">&</span><span class="pln"> filename</span><span class="pun">,</span><span class="pln"> std</span><span class="pun">::</span><span class="typ">vector</span><span class="str"><int></span><span class="pun">&</span><span class="pln"> shape</span><span class="pun">,</span><span class="pln">
std</span><span class="pun">::</span><span class="typ">vector</span><span class="pun"><</span><span class="typ">Scalar</span><span class="pun">>&</span><span class="pln"> data</span><span class="pun">);</span></code></pre>
<p><code>.npy</code>ファイルを読み込みたい場合は,<code>aoba::LoadArrayFromNumpy()</code>関数を使用してください.</p>
<ol>
<li><code>const std::string& filename</code>: 入力先のファイル名 </li>
<li><code>std::vector<int>& shape</code>: テンソルの階数とその長さを格納するベクトルの出力先</li>
<li><code>std::vector<Scalar>& data</code>: データ本体の出力先</li>
</ol>
<p><code>std::vector<T></code>の型<code>T</code>は指定された<code>.npy</code>ファイルの型と同一でなければなりません.もし異なる型を指定した場合には,<code>Numpy.hpp</code>は例外を送出します.</p>
<p><code>aoba::LoadArrayFromNumpy<T>()</code>関数の簡単なサンプルは次の通りです:</p><div class="se-section-delimiter"></div>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="com">#include</span><span class="pln"> </span><span class="str"><iostream></span><span class="pln">
</span><span class="com">#include</span><span class="pln"> </span><span class="str">"Numpy.hpp"</span><span class="pln">
</span><span class="com">// Now I assume np.load("file.npy").shape == (4, 5)</span><span class="pln">
std</span><span class="pun">::</span><span class="typ">vector</span><span class="str"><int></span><span class="pln"> s</span><span class="pun">;</span><span class="pln">
std</span><span class="pun">::</span><span class="typ">vector</span><span class="str"><double></span><span class="pln"> data</span><span class="pun">;</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="str">"/path/to/file.npy"</span><span class="pun">,</span><span class="pln"> s</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln">
std</span><span class="pun">::</span><span class="pln">cout </span><span class="pun"><<</span><span class="pln"> s</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> </span><span class="str">" "</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> s</span><span class="pun">[</span><span class="lit">1</span><span class="pun">]</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">endl</span><span class="pun">;</span><span class="pln"> </span><span class="com">// 4 5</span><span class="pln">
std</span><span class="pun">::</span><span class="pln">cout </span><span class="pun"><<</span><span class="pln"> data</span><span class="pun">.</span><span class="pln">size</span><span class="pun">()</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">endl</span><span class="pun">;</span><span class="pln"> </span><span class="com">// 20</span></code></pre>
<p><code>Numpy.hpp</code>はまた,入力先のテンソルの階数が必ず固定である場合に使えるいくつかのショートカット関数を提供しています.</p><div class="se-section-delimiter"></div>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="com">#include</span><span class="pln"> </span><span class="str">"Numpy.hpp"</span><span class="pln">
</span><span class="typ">int</span><span class="pln"> rows</span><span class="pun">,</span><span class="pln"> cols</span><span class="pun">;</span><span class="pln">
std</span><span class="pun">::</span><span class="typ">vector</span><span class="str"><double></span><span class="pln"> data</span><span class="pun">;</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="str">"/path/to/file1.npy"</span><span class="pun">,</span><span class="pln"> rows</span><span class="pun">,</span><span class="pln"> cols</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln">
std</span><span class="pun">::</span><span class="pln">cout </span><span class="pun"><<</span><span class="pln"> rows </span><span class="pun"><<</span><span class="pln"> </span><span class="str">" "</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> cols </span><span class="pun"><<</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">endl</span><span class="pun">;</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file2.npy"</span><span class="pun">,</span><span class="pln"> a0</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln"> </span><span class="com">// 1st-order tensor</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file3.npy"</span><span class="pun">,</span><span class="pln"> b0</span><span class="pun">,</span><span class="pln"> b1</span><span class="pun">,</span><span class="pln"> b2</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln"> </span><span class="com">// 3rd-order tensor</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file4.npy"</span><span class="pun">,</span><span class="pln"> c0</span><span class="pun">,</span><span class="pln"> c1</span><span class="pun">,</span><span class="pln"> c2</span><span class="pun">,</span><span class="pln"> c3</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln"> </span><span class="com">// 4th-order tensor</span></code></pre>
<h3 id="npyファイルの保存"><code>.npy</code>ファイルの保存</h3><div class="se-section-delimiter"></div>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="kwd">template</span><span class="pun"><</span><span class="kwd">typename</span><span class="pln"> </span><span class="typ">Scalar</span><span class="pun">></span><span class="pln">
</span><span class="kwd">void</span><span class="pln"> </span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="kwd">const</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">string</span><span class="pun">&</span><span class="pln"> filename</span><span class="pun">,</span><span class="pln"> </span><span class="kwd">bool</span><span class="pln"> fortran_order</span><span class="pun">,</span><span class="pln">
</span><span class="typ">int</span><span class="pln"> n_dims</span><span class="pun">,</span><span class="pln"> </span><span class="kwd">const</span><span class="pln"> </span><span class="typ">int</span><span class="pln"> shape</span><span class="pun">[],</span><span class="pln"> </span><span class="kwd">const</span><span class="pln"> </span><span class="typ">Scalar</span><span class="pun">*</span><span class="pln"> data</span><span class="pun">);</span></code></pre>
<p><code>.npy</code>ファイルを保存したい場合は,<code>aoba::SaveArrayAsNumpy()</code>関数を使用してください. </p>
<ol>
<li><code>const std::string& filename</code>: <code>.npy</code>ファイルの出力先</li>
<li><code>bool fortran_order</code>: このデータの格納順がcolumn-majorである場合は,このフラグを<code>true</code>に設定してください.</li>
<li><code>int n_dims</code>: テンソルの階数</li>
<li><code>const int shape[]</code>: テンソルの長さを表す配列のポインタ</li>
<li><code>const Scalar* data</code>: テンソルの実データを指すポインタ</li>
</ol>
<p><code>aoba::SaveArrayAsNumpy<T>()</code>関数の簡単なサンプルは次の通りです.</p><div class="se-section-delimiter"></div>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="com">#include</span><span class="pln"> </span><span class="str">"Numpy.hpp"</span><span class="pln">
</span><span class="com">// Now I assume s.size() == 2</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="str">"/path/to/file.npy"</span><span class="pun">,</span><span class="pln"> </span><span class="kwd">false</span><span class="pun">,</span><span class="pln"> </span><span class="lit">2</span><span class="pun">,</span><span class="pln"> </span><span class="pun">&</span><span class="pln">s</span><span class="pun">[</span><span class="lit">0</span><span class="pun">],</span><span class="pln"> </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span></code></pre>
<p><code>aoba::LoadArrayFromNumpy<T>()</code>と同様に,<code>Numpy.hpp</code>はまたいくつかのショートカット関数を提供しています.</p><div class="se-section-delimiter"></div>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="com">// the following functions assume that fortran_order == false.</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="str">"/path/to/file.npy"</span><span class="pun">,</span><span class="pln"> rows</span><span class="pun">,</span><span class="pln"> cols</span><span class="pun">,</span><span class="pln"> </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file2.npy"</span><span class="pun">,</span><span class="pln"> a0</span><span class="pun">,</span><span class="pln"> </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span><span class="pln"> </span><span class="com">// 1st-order tensor</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file3.npy"</span><span class="pun">,</span><span class="pln"> b0</span><span class="pun">,</span><span class="pln"> b1</span><span class="pun">,</span><span class="pln"> b2 </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span><span class="pln"> </span><span class="com">// 3nd-order tensor</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file4.npy"</span><span class="pun">,</span><span class="pln"> c0</span><span class="pun">,</span><span class="pln"> c1</span><span class="pun">,</span><span class="pln"> c2</span><span class="pun">,</span><span class="pln"> c3</span><span class="pun">,</span><span class="pln"> </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span><span class="pln"> </span><span class="com">// 4th-order tensor</span></code></pre>
<p>もし何か質問などありましたら,お気軽に私(rezoolab at gmail.com)に連絡ください :) Enjoy!</p>rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com1tag:blogger.com,1999:blog-1062283300115822026.post-61118450189733232602014-03-02T16:36:00.001+09:002014-03-02T22:24:29.019+09:00Numpy.hpp - Simple Header-only Library for handling .npy files<h2 id="introduction">Introduction</h2>
<p>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.</p>
<p>The important thing in this case is an existence of a C++ library for reading/writing <code>.npy</code> 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).</p>
<p>To handle these problems, I reimplemented a new C++ library called <code>Numpy.hpp</code>. Comparing with the existing libraries, <code>Numpy.hpp</code> has two advantages:</p>
<ul>
<li><strong>Header-only</strong>: you do not need to precompile numpy.hpp. All you have to do is just add <code>#include "Numpy.hpp"</code> into your source file.</li>
<li><strong>No external dependencies</strong>: numpy.hpp uses only STL libraries in C++03. This means that <code>numpy.hpp</code> can work well in any environments including MSVC.</li>
</ul>
<p>You can download <code>Numpy.hpp</code> from the following URL: <br>
<a href="https://gist.github.com/rezoo/5656056">https://gist.github.com/rezoo/5656056</a></p>
<h2 id="usage">Usage</h2>
<p>In this section I will briefly introduce the usage of <code>Numpy.hpp</code>.</p>
<h3 id="load-npy-files">Load <code>.npy</code> files</h3>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="kwd">template</span><span class="pun"><</span><span class="kwd">typename</span><span class="pln"> </span><span class="typ">Scalar</span><span class="pun">></span><span class="pln">
</span><span class="kwd">void</span><span class="pln"> </span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="kwd">const</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">string</span><span class="pun">&</span><span class="pln"> filename</span><span class="pun">,</span><span class="pln"> std</span><span class="pun">::</span><span class="typ">vector</span><span class="str"><int></span><span class="pun">&</span><span class="pln"> shape</span><span class="pun">,</span><span class="pln">
std</span><span class="pun">::</span><span class="typ">vector</span><span class="pun"><</span><span class="typ">Scalar</span><span class="pun">>&</span><span class="pln"> data</span><span class="pun">);</span></code></pre>
<p>If you wish to load <code>.npy</code> files in C++, please use <code>aoba::LoadArrayFromNumpy()</code> functions.</p>
<ol>
<li><code>const std::string& filename</code>: the filename of the input. </li>
<li><code>std::vector<int>& shape</code>: the output of the tensor’s rank and its lengths.</li>
<li><code>std::vector<Scalar>& data</code>: the destination of the data.</li>
</ol>
<p>Please note that the type of <code>std::vector<T></code> must be equivalent to the type of the specified <code>.npy</code> file. <code>Numpy.hpp</code> will send an exception if you specify the different type.</p>
<p>A simple example of <code>aoba::LoadArrayFromNumpy<T></code> function is as follows:</p><div class="se-section-delimiter"></div>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="com">#include</span><span class="pln"> </span><span class="str"><iostream></span><span class="pln">
</span><span class="com">#include</span><span class="pln"> </span><span class="str">"Numpy.hpp"</span><span class="pln">
</span><span class="com">// Now I assume np.load("file.npy").shape == (4, 5)</span><span class="pln">
std</span><span class="pun">::</span><span class="typ">vector</span><span class="str"><int></span><span class="pln"> s</span><span class="pun">;</span><span class="pln">
std</span><span class="pun">::</span><span class="typ">vector</span><span class="str"><double></span><span class="pln"> data</span><span class="pun">;</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="str">"/path/to/file.npy"</span><span class="pun">,</span><span class="pln"> s</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln">
std</span><span class="pun">::</span><span class="pln">cout </span><span class="pun"><<</span><span class="pln"> s</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> </span><span class="str">" "</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> s</span><span class="pun">[</span><span class="lit">1</span><span class="pun">]</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">endl</span><span class="pun">;</span><span class="pln"> </span><span class="com">// 4 5</span><span class="pln">
std</span><span class="pun">::</span><span class="pln">cout </span><span class="pun"><<</span><span class="pln"> data</span><span class="pun">.</span><span class="pln">size</span><span class="pun">()</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">endl</span><span class="pun">;</span><span class="pln"> </span><span class="com">// 20</span></code></pre>
<p><code>Numpy.hpp</code> also provides several shortcuts that can be used where the order of tensor is fixed.</p>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="com">#include</span><span class="pln"> </span><span class="str">"Numpy.hpp"</span><span class="pln">
</span><span class="typ">int</span><span class="pln"> rows</span><span class="pun">,</span><span class="pln"> cols</span><span class="pun">;</span><span class="pln">
std</span><span class="pun">::</span><span class="typ">vector</span><span class="str"><double></span><span class="pln"> data</span><span class="pun">;</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="str">"/path/to/file1.npy"</span><span class="pun">,</span><span class="pln"> rows</span><span class="pun">,</span><span class="pln"> cols</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln">
std</span><span class="pun">::</span><span class="pln">cout </span><span class="pun"><<</span><span class="pln"> rows </span><span class="pun"><<</span><span class="pln"> </span><span class="str">" "</span><span class="pln"> </span><span class="pun"><<</span><span class="pln"> cols </span><span class="pun"><<</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">endl</span><span class="pun">;</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file2.npy"</span><span class="pun">,</span><span class="pln"> a0</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln"> </span><span class="com">// 1st-order tensor</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file3.npy"</span><span class="pun">,</span><span class="pln"> b0</span><span class="pun">,</span><span class="pln"> b1</span><span class="pun">,</span><span class="pln"> b2</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln"> </span><span class="com">// 3rd-order tensor</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">LoadArrayFromNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file4.npy"</span><span class="pun">,</span><span class="pln"> c0</span><span class="pun">,</span><span class="pln"> c1</span><span class="pun">,</span><span class="pln"> c2</span><span class="pun">,</span><span class="pln"> c3</span><span class="pun">,</span><span class="pln"> data</span><span class="pun">);</span><span class="pln"> </span><span class="com">// 4th-order tensor</span></code></pre>
<h3 id="save-npy-files">Save <code>.npy</code> files</h3>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="kwd">template</span><span class="pun"><</span><span class="kwd">typename</span><span class="pln"> </span><span class="typ">Scalar</span><span class="pun">></span><span class="pln">
</span><span class="kwd">void</span><span class="pln"> </span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="kwd">const</span><span class="pln"> std</span><span class="pun">::</span><span class="pln">string</span><span class="pun">&</span><span class="pln"> filename</span><span class="pun">,</span><span class="pln"> </span><span class="kwd">bool</span><span class="pln"> fortran_order</span><span class="pun">,</span><span class="pln">
</span><span class="typ">int</span><span class="pln"> n_dims</span><span class="pun">,</span><span class="pln"> </span><span class="kwd">const</span><span class="pln"> </span><span class="typ">int</span><span class="pln"> shape</span><span class="pun">[],</span><span class="pln"> </span><span class="kwd">const</span><span class="pln"> </span><span class="typ">Scalar</span><span class="pun">*</span><span class="pln"> data</span><span class="pun">);</span></code></pre>
<p>If you with to save <code>.npy</code> files, please use <code>aoba::SaveArrayAsNumpy()</code> functions. </p>
<ol>
<li><code>const std::string& filename</code>: the destination of the <code>.npy</code> file.</li>
<li><code>bool fortran_order</code>: if you store the data as column-major order, please set this flag as <code>true</code>.</li>
<li><code>int n_dims</code>: the rank of the tensor.</li>
<li><code>const int shape[]</code>: the array that contains the tensor’s lengths.</li>
<li><code>const Scalar* data</code>: the raw data of the tensor.</li>
</ol>
<p>A simple example of <code>aoba::SaveArrayAsNumpy<T></code> function is as follows:</p>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="com">#include</span><span class="pln"> </span><span class="str">"Numpy.hpp"</span><span class="pln">
</span><span class="com">// Now I assume s.size() == 2</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="str">"/path/to/file.npy"</span><span class="pun">,</span><span class="pln"> </span><span class="kwd">false</span><span class="pun">,</span><span class="pln"> </span><span class="lit">2</span><span class="pun">,</span><span class="pln"> </span><span class="pun">&</span><span class="pln">s</span><span class="pun">[</span><span class="lit">0</span><span class="pun">],</span><span class="pln"> </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span></code></pre>
<p>As with <code>aoba::LoadArrayFromNumpy<T>()</code>, <code>Numpy.hpp</code> also provides several shortcut functions.</p>
<pre class="prettyprint prettyprinted" style=""><code class="language-cpp"><span class="com">// the following functions assume that fortran_order == false.</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="str">"/path/to/file.npy"</span><span class="pun">,</span><span class="pln"> rows</span><span class="pun">,</span><span class="pln"> cols</span><span class="pun">,</span><span class="pln"> </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file2.npy"</span><span class="pun">,</span><span class="pln"> a0</span><span class="pun">,</span><span class="pln"> </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span><span class="pln"> </span><span class="com">// 1st-order tensor</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file3.npy"</span><span class="pun">,</span><span class="pln"> b0</span><span class="pun">,</span><span class="pln"> b1</span><span class="pun">,</span><span class="pln"> b2 </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span><span class="pln"> </span><span class="com">// 3nd-order tensor</span><span class="pln">
aoba</span><span class="pun">::</span><span class="typ">SaveArrayAsNumpy</span><span class="pun">(</span><span class="pln">
</span><span class="str">"/path/to/file4.npy"</span><span class="pun">,</span><span class="pln"> c0</span><span class="pun">,</span><span class="pln"> c1</span><span class="pun">,</span><span class="pln"> c2</span><span class="pun">,</span><span class="pln"> c3</span><span class="pun">,</span><span class="pln"> </span><span class="pun">&</span><span class="pln">data</span><span class="pun">[</span><span class="lit">0</span><span class="pun">]);</span><span class="pln"> </span><span class="com">// 4th-order tensor</span></code></pre>
<p>If you have any questions, do not hesitate to ask me (rezoolab at gmail.com) :) Enjoy!</p>rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-82486029241828980052013-12-31T23:06:00.000+09:002014-01-01T08:24:17.245+09:00ノンパラメトリック平均場近似と確率伝播法の提案2013年も残り少なくなってきました.そこで,今回は研究の締めくくりとして,今年のCVPRで発表した論文について紹介したいと思います.<br />
<br />
僕の研究は,大雑把に言うとグラフィカルモデルの最適化手法に関するものです.本研究の貢献を簡単にまとめると,次のようになります.<br />
<br />
<ul>
<li>グラフィカルモデル内の確率変数が連続的である場合,計算上の理由から,変分ベイズ法ならびに確率伝播法で扱える分布のモデルには制限があります.提案手法を使うことで,非指数型の分布族に関しても変分ベイズ法を効率的に適用できるようになりました.</li>
<li>さらに,同様の議論を確率伝播法にも適用することで,確率伝播法をノンパラメトリックに扱うこともできるようになりました.</li>
</ul>
<br />
<div>
実際の論文では対象のグラフィカルモデルを一階のマルコフ確率場に限定していますが,変分ベイズ法などの別の領域にも容易に適用できます.また,議論を簡単にするため実際の論文では確率変数の離散化を焦点に当てた形で書いていますが,その流れに沿って書いてもあまり面白く無いので,ここではどのような過程でこのような研究に取り組んだのか簡単に書いていきたいと思います.</div>
<div>
<div>
<br /></div>
<div>
グラフィカルモデルは機械学習の多くの分野で使われてきた確率モデルです.この確率モデルを用いて対象の問題を高精度に解くためには,確率変数の周辺化が重要となってきます.周辺化の計算を直接行うことは困難であるため,大きく分けて2つの方法論で周辺分布を近似的に求めることが一般的です.一つはギブスサンプリングに代表される,大量のサンプルをばらまいてモンテカルロ的に推定解を求める方法,もう1つは平均場近似(変分ベイズ法)や確率伝播法に代表される,変分原理に基づいた近似的推論法を使う方法です.後者は得られる推定解の精度が比較的悪いものの,特定の反復方程式に従って周辺分布を更新するだけで良いため,高速に周辺分布の推定解が求められる利点を持ちます.本研究は後者の推定手法に関するものです.</div>
<div>
<br /></div>
<div>
変分原理に基づいた周辺分布の反復計算は,グラフィカルモデル内の確率変数が連続的であるか離散的であるかで異なります.変数が連続的である場合,周辺分布の密度関数は少数のパラメータからなるパラメトリックな連続分布として表現されます.そして,平均場近似,確率伝播法の両方とも,実際の最適化には変分原理に従って導出された特定の反復方程式に従って,対象のパラメータを更新していきます.しかしながら,このような更新は更新後の分布が更新前の分布と同じパラメータで表現される必要があるため,扱える分布モデルの範囲は非常に限られていました.以上の背景から,僕は平均場近似ならびに確率伝播法をノンパラメトリックに扱えるような方法論を新しく提案すれば,良いインパクトになるのではないかと考えました.</div>
<div>
<br /></div>
<div>
このような元々の平均場近似や確率伝播法では扱えない分布をどうにかして扱えるよう拡張できないだろうかという要望は多く,これまでにいくつかの試みが提案されています.これらの試みはすべて次の混合分布の仮定を下にしています.すなわち,平均場近似や確率伝播法で扱う近似分布やメッセージを,次の混合ガウス分布で仮定することです.</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJRs7fd-FHrrk-9yuuOCe7SzDS0hhkYOt0nxPeyHs9l01Ec7Sk0ktJVWUA30PK9108HFbba4u-mQSLvDN1nINkCPOVJSJE3CLMDbXeGlSeENikwb16MCmwyc3t1COJOMcxF26Wj8DJTDg/s1600/image3.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJRs7fd-FHrrk-9yuuOCe7SzDS0hhkYOt0nxPeyHs9l01Ec7Sk0ktJVWUA30PK9108HFbba4u-mQSLvDN1nINkCPOVJSJE3CLMDbXeGlSeENikwb16MCmwyc3t1COJOMcxF26Wj8DJTDg/s1600/image3.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div>
もし上の仮定の下で混合ガウス分布のパラメータに関する反復方程式を導出できたのであれば,対象のグラフィカルモデルがどのような分布であったとしても平均場近似や確率伝播法を適用できる.つまり,平均場近似や確率伝播法の可能性を一気に広げることができます.</div>
<div>
<br /></div>
<div>
しかし,このアイデアを直接適用してもあまり上手くいかないことが知られていました.なぜなら,変分原理は計算の段階で混合ガウス分布のエントロピーを計算する必要があるためです.混合ガウス分布のエントロピーは計算困難であるため,このままでは反復方程式を解析的に導出できません.</div>
</div>
<div>
<br /></div>
<div>
<div>
このような問題に対応するための方法として,これまで大きく分けて2つのアプローチが提案されてきました.1つは解析解を出すために,混合分布のエントロピーにJensenの不等式を用いてさらなる近似を行い,解析解を求める方法です[3].しかしながら,この方法論は平均場近似や確率伝播法で導入した近似分布の他に,さらなる近似をおくことになります.もう1つはモンテカルロ法を用いて,混合ガウス分布のパラメータの反復解を力技で求める方法です[1][2].しかしながら,この方法論ではサンプル数に依存せず高速に推定解を求められる,確率伝播法の利点を失うことになります.</div>
<div>
<br /></div>
<div>
これ以外にエントロピーの計算をどうにかして解決できないだろうかと思い,思考を巡らせました.前述の通り,そもそもの問題は混合ガウス分布のエントロピーが計算困難であることでした.ガウス分布の裾が別のガウス分布に干渉してしまうため,それぞれ独立した分布として扱うことができないのです.そうであれば,最初の出発点に混合ガウス分布を選ぶのではなく,互いに干渉しない別の分布を選べば良いのではないでしょうか?例えばヒストグラムのような,互いに干渉せず,エントロピーが容易に計算できる混合矩形分布を採用するのは?<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhuWJJKW0RQmV9HYamifb_hOSwqw1LY1L1R5hNpPZmMkujhmLgikHN7x-JcO0PQEoFKntsk22Muv7hJggzqfpT3N1he0JguCaRuYZhpV9WnBcnnlBqkmKGUCdaOcrK7T8ezdM9LUMnW5xw/s1600/image4.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhuWJJKW0RQmV9HYamifb_hOSwqw1LY1L1R5hNpPZmMkujhmLgikHN7x-JcO0PQEoFKntsk22Muv7hJggzqfpT3N1he0JguCaRuYZhpV9WnBcnnlBqkmKGUCdaOcrK7T8ezdM9LUMnW5xw/s1600/image4.png" /></a></div>
<br /></div>
<div>
ここで,hは変数iにおける,s番目の矩形分布を表します.本研究で新しく提案したアイデアは基本的にはこれだけです.この矩形分布の位置とサイズは,重なっていなければ変数空間のどの箇所に配置しても構いません.すなわち,この表現を用いて変数空間の重要な箇所を密に離散化し,重要でない箇所を疎に離散化することもできます.さらに,変数毎に配置する矩形分布の個数は異なっていても構いません.すなわち,変数の重要性に従って,変数の離散化の度合いを変えることもできます. </div>
<div>
<br /></div>
<div>
次に,上の表現を用いて本来の連続的な最適化問題をパラメータαの最適化問題へ変換し,パラメータに関する反復方程式を導出しました.詳しくは言及しませんが,最終的に導出した反復方程式の形は対象の変数が離散的な場合の平均場近似,確率伝播法の形と非常に似通っています.唯一の違いはエネルギー関数に加わる補正項の存在です.この補正項は変数空間の離散化の「非一様性」を補正する役割をもちます.すなわち,提案手法によって,グラフィカルモデルの連続的な変数空間を非一様に離散化し,ノンパラメトリックに扱えるようになったのです.</div>
<div>
<br /></div>
<div>
以上,急ぎですが本研究の解説記事を書かせていただきました.興味を持っていただけた方は,次のURLからダウンロードしていただければ幸いです.</div>
</div>
<div>
<br /></div>
<div>
<a href="http://www.vision.is.tohoku.ac.jp/index.php/download_file/view/35/177/167/">http://www.vision.is.tohoku.ac.jp/index.php/download_file/view/35/177/167/</a><br />
<br />
参考文献<br />
[1] A. Ihler et. al., Particle Belief Propagation, AISTATS, 2009<br />
[2] E. B. Sudderth et. al., Nonparametric Belief Propagation, CVPR, 2003<br />
[3] S. J. Gershman et. al., Nonparametric Variational Inference, ICML, 2012</div>
rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-29155019512327095162013-04-26T00:03:00.002+09:002013-04-26T02:18:24.152+09:00条件付き確率場の推論と学習<div style="text-align: center;">
<iframe allowfullscreen="" frameborder="0" height="356" marginheight="0" marginwidth="0" mozallowfullscreen="" scrolling="no" src="http://www.slideshare.net/slideshow/embed_code/19715143" style="border-width: 1px 1px 0; border: 1px solid #CCC; margin-bottom: 5px;" webkitallowfullscreen="" width="427"> </iframe> </div>
<div style="margin-bottom: 5px;">
<div style="text-align: center;">
<strong> <a href="http://www.slideshare.net/rezoolab/seminar-19715143" target="_blank" title="条件付き確率場の推論と学習">条件付き確率場の推論と学習</a> </strong> from <strong><a href="http://www.slideshare.net/rezoolab" target="_blank">Masaki Saito</a></strong> <br />
<br />
<div style="text-align: left;">
<span style="text-align: -webkit-auto;"> 前回の研究室ゼミで話した,条件付き確率場についての簡単な紹介スライドを公開します.内容は基本的なPairwise CRFについての話と,それが実際にどのような問題に対して応用がなされているか,またもっとも使われる最適化手法の一つである,平均場近似と確率伝搬法についての簡単な解説と,CRFの初歩的な学習法といった流れです.</span><br />
<span style="text-align: -webkit-auto;"><br /></span>
<span style="text-align: -webkit-auto;">補足:</span><br />
<span style="text-align: -webkit-auto;"> CRFの学習は自然言語処理など他の分野でもよく使われている学習法なのでおそらく知っている方は多いと思うのですが,普通はダイレクトに対数尤度を最大化することでパラメータを求めています.しかし,今回はKL距離という確率分布間の近さを図る指標を用いて対数尤度を拡張し,パラメータを求めています.</span>なぜわざわざそうするのかというと,だいたいの最適化手法がKL距離の最小化という議論に帰着できて,すっきりするからです.KL距離を使えば平均場近似,確率伝搬法,最尤推定,あと時間の都合上話せなかったTRWやGeneralized Belief Propagationなどの最適化手法がすべて統一的に理解できます.これによって覚えることが少なくなって楽できるだけでなく,このKL距離を使った新しい最適化手法を提案することもできます.</div>
</div>
</div>
rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-8664712575071081072013-04-25T02:48:00.005+09:002013-04-25T18:52:33.549+09:00近況報告お久しぶりです.細々とした近況は<a href="http://twitter.com/rezoolab">Twitter</a>に書いているのですが,あくまでTwitterはつぶやきであって,記事でない.定期的に自分の中でまとめた内容をブログに吐き出す必要がありますね.<br />
<ol>
<li>修士論文を提出し,晴れて今年度から博士後期課程の1年生です.不安が全くないと言えば嘘になりますが,それでも毎日好きなことを学び,研究できる環境は非常に楽しい.みなさんもぜひ博士課程に進学しましょう.唯一不満に思うところは頻繁に書類を書く業務が入るあたりですが,仕方のないことです.</li>
<li>後できちんとした内容を書きますが,2013年のCVPR(コンピュータビジョンのトップカンファレンス)に採択されました.去年を含めると今回で2回目の採択ですが,前回と異なるのが,今回はポスターではなくオーラルセッションでの採択というところです.オーラルセッションは非常に狭き門(採択率 60/1870 = 3.2%)でありますので,そのような研究成果を得られたことは非常に嬉しいですし,来年もぜひ通しておきたいところです.</li>
<li>博士後期課程になりましたので,これからは本名メインで活動していきます.プロフィール画像も更新しました.前に使ってた画像は高校時代からのものでしたから,大体7年位使っていたんでしょうか.というかこのブログ,今年で7年目になるのですね.</li>
<li>自身の宣伝もかねて,論文やコードなどはこれから作る自己紹介のページにて積極的に公開していきます.論文についてはできるだけ読んでもらえるよう,簡単な解説も記事にしていきたいと思っております.</li>
</ol>
rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-68037656088118915372012-10-20T21:19:00.000+09:002013-05-29T01:08:57.964+09:00C++で2種類のコンテナを同時にソートするC++でデータを扱っていると,2種類のコンテナを同時にソートしたい場合が時々出てくる.どういうことかというと,
<br />
<pre class="prettyprint">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;
}
</pre>
と,片方のソート結果を元にもう片方をソートするような感じ.これはCUDAライブラリである<a href="http://code.google.com/p/thrust/">thrust</a>に標準搭載されており,<span style="font-family: Courier New, Courier, monospace;">sort_by_key</span>という名前がついている.<br />
<br />
さて,こういうときに一からソートアルゴリズムを組み直すなんてことは馬鹿らしい.じゃあどうするか.基本的には既存の<span style="font-family: Courier New, Courier, monospace;">std::sort</span>に渡すイテレータを変形することで,なんとか解決しようとする.そうなると,boostに<span style="font-family: Courier New, Courier, monospace;">zip_iterator</span>が存在することに気付けば,イテレータを<span style="font-family: Courier New, Courier, monospace;">boost::tuple</span>の組で表現することによって,まとめてソートさせればいいような気がしてくる.こんな風に:
<br />
<pre class="prettyprint">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());
</pre>
しかし,そうは問屋が降ろさない.何故なら,<span style="font-family: Courier New, Courier, monospace;">boost::zip_iterator</span>はWritableではない.これはソートアルゴリズムの要件を満たしていないため,適用させようとしてもコンパイルエラーが発生してしまう!<br />
<br />
それじゃあどのようにして解決するかというと,一番楽な方法としてはWritableなコンセプトを満たしたzipイテレータを新しく作ってやって,それを用いて<span style="font-family: Courier New, Courier, monospace;">std::sort</span>を行うのがいい.とはいっても,C++でそれをやるのは若干手間のかかる作業になる.具体的にはこんな感じで実装した:<br />
<script src="https://gist.github.com/rezoo/5663872.js"></script>
実装したといっても,実際には『<a href="http://www.stanford.edu/~dgleich/notebook/2006/03/sorting_two_arrays_simultaneou.html">Sorting two arrays simultaneously</a>』の内容を若干手直しするだけなので基本的にはそんなに大変じゃないし,ただ利用するぶんには単純に<span style="font-family: Courier New, Courier, monospace;">sort_by_key()</span>を呼び出すだけでよい.<br />
<br />
しかし,なんでこんなにC++は行数が多くなってしまうのか.haskellなんか一行で済むのに.rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com3tag:blogger.com,1999:blog-1062283300115822026.post-6317047433148628602012-08-21T23:58:00.000+09:002013-05-29T01:14:39.259+09:00超一様分布列をBoost.Randomを用いて実装するご存知の通り、C++のライブラリには高性能な擬似乱数生成エンジン<a href="http://www.boost.org/libs/random/">Boost.Random</a>が存在している。Boost.Randomは擬似乱数を生成する<a href="http://www.boost.org/doc/libs/1_51_0/doc/html/boost_random/reference.html#boost_random.reference.generators">エンジン部分</a>と、一様分布、ガウス分布、ベルーヌーイ分布など好みの分布を生成する<a href="http://www.boost.org/doc/libs/1_51_0/doc/html/boost_random/reference.html#boost_random.reference.distributions">出力部分</a>の2つで構成されている。<br />
<br />
しかし、ことモンテカルロ法においては主に収束速度を早める目的として、擬似乱数ではなく<a href="http://www.math.tohoku.ac.jp/akama/2006/lds.html">低食い違い量列</a>(Low-Discrepancy Sequence, LDS)と呼ばれる数列が用いられる場合がある。これを用いると低次元での積分計算の際に収束速度を大きく改善することができ、具体的には従来のモンテカルロ法で(1/√N)だった収束速度が約(1/N)くらいにまで落とすことができる。<br />
<br />
さて、このような低食い違い量列をBoost.Randomで用いるためにはエンジン部分を製作すればいいだけで、実際には以下のような形で実装すれば良い。今回は簡単な<a href="http://lucille.sourceforge.net/blog/archives/000033.html">Halton列</a>を用いて実装を行った:<br />
<br />
<script src="https://gist.github.com/rezoo/5663945.js"></script>
具体的には、対象のファンクタがUniform Random Number Generatorのコンセプトを満たすためには以下の記述が必要となる:<br />
<ul>
<li><span style="font-family: 'Courier New', Courier, monospace;">X::result_type</span>に<span style="font-family: 'Courier New', Courier, monospace;">operator()</span>が返す型を指定</li>
<li><span style="font-family: 'Courier New', Courier, monospace;">operator()</span>にエンジン部分を記述</li>
<li><span style="font-family: 'Courier New', Courier, monospace;">min(), max()</span>にエンジンが生成する分布の最小値、最大値を指定</li>
</ul>
<div>
さらに、対象のファンクタが上述の内容に加えPseudo-Random Number Generatorのコンセプトを満たすためには:<br />
<br />
<ul>
<li><span style="font-family: 'Courier New', Courier, monospace;">X()</span>にデフォルトの状態で生成されるようなコンストラクタを指定</li>
<li><span style="font-family: 'Courier New', Courier, monospace;">X(...)</span>にユーザ指定の状態で生成されるコンストラクタを指定</li>
<li><span style="font-family: 'Courier New', Courier, monospace;">seed(...)</span>にエンジンの状態を変更するためのメンバ関数を指定</li>
</ul>
</div>
するだけで良い。後はBoost側が良きに計らってくれる。非常に便利だと思う。<br />
<br />
余談<br />
C++0xには同等の<span style="font-family: 'Courier New', Courier, monospace;">std::random</span>が付属しているため、これを利用すればいいのではという声もあるとは思うが、<span style="font-family: 'Courier New', Courier, monospace;">std::uniform_real_distribution<></span>にこのLDSを食わせたところ、どうみてもLDSではない値が出力される結果となった。恐らく<span style="font-family: 'Courier New', Courier, monospace;">result_type</span>がdoubleという一般の擬似乱数ではない型を指定しているために起こる現象だと思うが、ソースを見ていないので一先ずは従来のBoost.Randomの手法を用いることにする。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com3tag:blogger.com,1999:blog-1062283300115822026.post-32505301340208955522012-08-18T22:58:00.001+09:002012-08-18T23:44:13.421+09:00Restricted Boltzmann Machineの学習手法についての簡単なまとめ<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjdDk0ofnZOMHcfdLeUJfVMzP9-MiWzAHNTQUa_BTBlfP9P56tJs1ot2xJUjomuaANmFw7pexj8jCAp1fjuW2kGq2dq5_0r-qvXofNQEeTN6MZKG9DivnWRm-3KDvj5CzChaCNdY_UfLPg/s1600/rbm.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjdDk0ofnZOMHcfdLeUJfVMzP9-MiWzAHNTQUa_BTBlfP9P56tJs1ot2xJUjomuaANmFw7pexj8jCAp1fjuW2kGq2dq5_0r-qvXofNQEeTN6MZKG9DivnWRm-3KDvj5CzChaCNdY_UfLPg/s1600/rbm.png" style="border-bottom-width: 0px; border-color: initial; border-image: initial; border-left-width: 0px; border-right-width: 0px; border-style: initial; border-top-width: 0px;" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><br /></td></tr>
</tbody></table>
近年の機械学習では<a href="http://en.wikipedia.org/wiki/Deep_learning">Deep Learning</a>と呼ばれる分野が一世を風靡しています.コンピュータビジョンや自然言語処理,音声認識などの分野では何らかの問題を解こうとした際に,まず対象の入力データからSIFTやケプストラムといった何らかのアルゴリズムを用いて特徴ベクトルを抽出し,ごりごりと判別していくといった流れが一般的です.しかし,その特徴ベクトルを生成するという生のデータから本質となる部分を抽出するアルゴリズム自体は研究者が一生懸命考えながら作るのが普通でした.<br />
<br />
Deep Learningの分野で最も有名な手法の一つである<a href="http://www.scholarpedia.org/article/Deep_belief_networks">Deep Belief Nets</a>(DBN) [Hinton06]は,研究者がアルゴリズムを作るのではなく,それ自体も機械学習にやらせましょうという動機で生まれたアルゴリズムです.DBNではまるで一昔前にやたら流行ったニューラルネットワークのように各ノードを層状に配置し,それらのパラメータを<a href="http://deeplearning.net/tutorial/rbm.html">Restricted Boltzmann Machine</a>(RBM)と呼ばれるモデルを用いて学習することによって自動的に特徴ベクトルを生成します.これが既存手法を大きく超える精度を出してしまったためにDeep Learningは一躍有名となり,ニューラルネットの再来か,機械学習は人工知能の夢を見るか,みたいに言われているような状況です.最近では<a href="http://itpro.nikkeibp.co.jp/article/NEWS/20120627/405501/">Googleがネコのニューロンを自動的に学習した</a>とかでニュースになっていました.<br />
<br />
さて,そのような状況にもかかわらず,そのDeep Belief Netsやその改良版である<a href="http://www.mit.edu/~rsalakhu/papers/dbm.pdf">Deep Boltzmann Machine</a>を学習するためには,まずその基礎となっているRBMについて学ぶ必要があります.しかしそのRBMですらお世辞にもわかり易いものとは言えず,しかもその最終的な更新式やその導出過程についてきちんと書いている文献は今まで確認したところ存在していません(*1).「チュートリアルや論文だけでは良く分からないので,実際に手を動かして確認したい.だがその計算過程は公開されていないため,結局良く分からないままで終わる」ことが多いように思われました.<br />
<br />
ということで,RBMの計算過程について書いた個人用のpdfファイルを公開してみることにします.外部に見せるということで多少文章を詳しくしましたが,これ単体で学ぶのではなく他の資料と比較しながら見るとベターかと思います.多分数式は間違ってないと思いますが,何かありましたらご一報頂けると幸いです.<br />
<br />
<a href="http://dl.dropbox.com/u/2048288/RestrictedBoltzmannMachine.pdf">http://dl.dropbox.com/u/2048288/RestrictedBoltzmannMachine.pdf</a><br />
<br />
*1 <a href="http://deeplearning.net/tutorial/rbm.html">Deeplearning.net</a>の式は比較的詳しいですが,実際に計算したところ微妙に間違えている箇所を数ヶ所発見していますrezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com6tag:blogger.com,1999:blog-1062283300115822026.post-62824618831600944952011-11-24T20:11:00.001+09:002011-11-27T00:53:40.152+09:00直積量子化(Product Quantization)を用いた近似最近傍探索についての簡単な解説<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjfFa1hv_aB56UGH-OCP1Bf_aKTlAINGZhKJvd41XTuMFXf8T6brybloLFKQvOCzIbNr_8zS50aMiQkniQyF-e12W1J9ji-cdEO3JNQwcYwuDKK9tXvCge0t35a3rAdsN1XzLi0l1kQvjE/s1600/nabe.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjfFa1hv_aB56UGH-OCP1Bf_aKTlAINGZhKJvd41XTuMFXf8T6brybloLFKQvOCzIbNr_8zS50aMiQkniQyF-e12W1J9ji-cdEO3JNQwcYwuDKK9tXvCge0t35a3rAdsN1XzLi0l1kQvjE/s1600/nabe.jpg" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">"<a href="http://www.flickr.com/photos/chatani/303581016/">aka motsu-nabe</a>" by <a href="http://www.flickr.com/photos/chatani/">chatani</a></td></tr>
</tbody></table>
<span class="Apple-style-span" style="font-size: large;">概要</span><br />
冬の寒さも一段と厳しくなってまいりました。おでんや鍋が恋しくなる季節です。<br />
<br />
さて、最近ようやっと一仕事が終わりまして、長ったらしい記事が書けるようになりました。ですので、今回は2011年に<a href="http://www.computer.org/portal/web/tpami">TPAMI</a>で発表された、近似最近傍探索についての論文『<a href="http://lear.inrialpes.fr/pubs/2011/JDS11/jegou_searching_with_quantization.pdf">Product quantization for nearest neighbor search</a>』について簡単に紹介したいと思います。<br />
<br />
この論文は2011年に発表された、最近傍探索アルゴリズムの決定打です。シンプルな理論でありながら既存手法を打ち破るほどの強力な性能を有し、速度も非常に高速、かつ省メモリなのでスマートフォンに載せ、リアルタイムで動作させることも可能です。<br />
<br />
以前この手法は<a href="http://www.slideshare.net/ren4yu/334-8419048">CV勉強会@関東</a>で紹介されたらしいのですが、具体的に紹介しているページは(最近すぎるので当たり前ですが)現在のところあまり見かけていません。ただ個人的にこの手法はあと5年くらいは本線で戦えるものと思っておりますので、きちんと解説した記事を出すことはある程度有意義なものかなぁと、そういう風に考えています。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">最近傍探索の課題</span><br />
<a href="http://ja.wikipedia.org/wiki/%E6%9C%80%E8%BF%91%E5%82%8D%E6%8E%A2%E7%B4%A2">最近傍探索</a>は最も古典的な手法でありながら、現在においても頻繁に利用されているアルゴリズムの一つです。近年のプロセッサの性能向上によって大規模なデータ処理が家庭のパソコンでもできるようになり、コンピュータビジョンにおいても数億件規模(10^8~9)の高次元ベクトルデータを最近傍探索を用いて処理することで、良い精度の結果を得ることが当たり前の時代になってきました。<br />
<br />
ただそうは言っても、最近傍探索をそのまま用いることは現実的ではなく、いくつかの課題が浮上してきます。具体的には以下のとおりです:<br />
<ul>
<li><b>いかに圧縮した状態でデータを格納するか</b>: ベクトルデータをそのままの状態で保持することはメモリ容量の圧迫へ繋がります。例えば128要素のベクトルを考えた場合、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">float</span>を使用した場合でも<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">32×128=4096bit</span>を1エントリで消費してしまいます。IOアクセスが生じてしまうと速度は一気に低下してしまうため、ベクトルデータはすべてメモリ上に展開しておくことが望ましいー<br />しかし、今回我々が扱いたいエントリ数というのは数億とかそういったオーダーです。そのオーダーのベクトルをそのまま格納するのは現在のメモリ容量だとなかなか厳しい。ですので、情報量を保ったままベクトルをいかに圧縮できるかが、実用的な面において重要な課題となります。</li>
<li><b>いかに高速に解を見つけられるか</b>: 次元数の多いベクトルの距離計算は基本的に低速です。さらにそのまま実装してしまうとすべてのエントリに対して距離の計算をする必要があるため速度はO(N)となり、実用的ではありません。1エントリあたりの計算コストを削減したり、そもそも明らかに異なる(クエリより離れすぎている)エントリを計算しないことで、計算を高速に行うことも同じく重要となります。</li>
</ul>
<div>
これらの課題を克服するため、近年では近似最近傍探索と呼ばれる手法が盛んに研究されております。その中でも有名なアルゴリズムは<a href="http://d.hatena.ne.jp/mamoruk/20091113/p1" style="font-style: italic;">LSH(Locality-Sensitive Hashing)</a>であり、ハッシュ関数によって抽出される近傍エントリのみを計算することで、高速に最近傍探索を行います。ただこの手法は計算用のハッシュを別に格納する構造であるため、全体のサイズは減少するのではなく増大してしまいます。<br />
<br />
また、<a href="http://research.preferred.jp/">Preferred Research</a>で紹介されている<i><a href="http://research.preferred.jp/2011/02/minhash/">MinHash</a></i>は、ベクトルの要素数が定まっておらず、成分が{0,1}の2値しかとらないような特殊なベクトルを計算する際においては有効な手法でありますが、そうでない場合においてはいったんベクトルをバイナリ表現に落とし込まなければならないので一般的に精度は低下してしまいます。今回紹介する直積量子化は、対象のベクトルが高次元であり、かつユークリッド距離による最近傍探索を行いたい場合において威力を発揮する手法です。</div>
<div>
<br /></div>
<span class="Apple-style-span" style="font-size: large;">スカラ量子化</span><br />
直積量子化について説明を行う前に、まず量子化とは何かについて説明していきます。物理をかじったことのある方ですと『<i>量子化 = 物理量を演算子へ変換する</i>』ものと思っちゃいますけど、そっちの量子化ではなくある数を特定のビット列で表現することの量子化です。<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgCNHFb6fZj5Gd8KZd-N4E5hwdAA91WtG75ETo1pLg2S8A3Mixpy48rZXhQseigKBkfc3o-2LSPOtf-rT-d2rBhNsBOesssgdjNH-zYrPs8xwF5mO9mCDpF86hD2HV8wmK3tYIepv6O2-Y/s1600/pic01.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgCNHFb6fZj5Gd8KZd-N4E5hwdAA91WtG75ETo1pLg2S8A3Mixpy48rZXhQseigKBkfc3o-2LSPOtf-rT-d2rBhNsBOesssgdjNH-zYrPs8xwF5mO9mCDpF86hD2HV8wmK3tYIepv6O2-Y/s1600/pic01.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">分割された領域に、特定のビットを割り当てる</td></tr>
</tbody></table>
例えば、0〜1の間をとる数を2ビットで表現しろという問題を考えた場合、普通は上図のように0〜1を4分割し、それぞれの数値を割り当てることで表現を行います。このように、少ないビット数で対象の数値を表現することを一般に『<b>スカラ量子化(Scalar Quantization)</b>』といいます。スカラ量子化によってfloat型(32bit)の数値は、精度を犠牲にすることによって僅か2bitで表現することができるようになりました。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">ベクトル量子化</span><br />
ところで、このスカラ量子化は対象のデータが数値であった場合には効率的ですが、ベクトルの場合だと効率的とは言えません。<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYgr3W8SzXAKTU1qj9dxHonEmv8egM81ftEgD9SEIiVflxdXg5ieULe9VlWD9ZICAKyWaMpadokuDxpRiEo0w8ARHvBO3b7J1EMJqeMCvUBUbVVWlD54DgdvAK68YfncBFjR7jCseSuoo/s1600/pic02.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYgr3W8SzXAKTU1qj9dxHonEmv8egM81ftEgD9SEIiVflxdXg5ieULe9VlWD9ZICAKyWaMpadokuDxpRiEo0w8ARHvBO3b7J1EMJqeMCvUBUbVVWlD54DgdvAK68YfncBFjR7jCseSuoo/s1600/pic02.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">スカラ量子化では、ベクトルを効率良く表現することができない</td></tr>
</tbody></table>
例えば上の図について考えてみましょう。この点群を1エントリ辺り1ビットで表現しようとした場合、スカラ量子化によって各成分を量子化してしまっては最低でも2ビットが必要となりますが、物事を空間的に考えることで発想を転換します。つまり、それぞれのクラスタを代表する点を予めコードブックとして保持しておき、各点は一番近い代表点のインデックスを保持しておくのです。これによってそれぞれの点は僅か1ビットで表現することができるようになりました。このような手法のことを一般に『<b>ベクトル量子化(Vector Quantization)</b>』といいます(*1)。<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjKu95kqUAO2R-FEqLqY5A1qogq10drHXxGAmboyFpaZ94FRwBBBm3Inh8MapICm4SqpNpc3132qyFGrTmJltmcOtp4jUK9yaGhnpbeU5-RikHq4Q4-Mu9i55O3KIy_V11QvlF73ktuNVs/s1600/pic03.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjKu95kqUAO2R-FEqLqY5A1qogq10drHXxGAmboyFpaZ94FRwBBBm3Inh8MapICm4SqpNpc3132qyFGrTmJltmcOtp4jUK9yaGhnpbeU5-RikHq4Q4-Mu9i55O3KIy_V11QvlF73ktuNVs/s1600/pic03.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">ベクトル量子化によって、ベクトルを2つの代表点によって表せるようになった</td></tr>
</tbody></table>
<br />
*1: 本手法のクラスタリングには<i>k-means algorithm</i>を使用しています。k-meansに関しては『<a href="http://d.hatena.ne.jp/nitoyon/20090409/kmeans_visualise">クラスタリングの定番アルゴリズム「K-means法」をビジュアライズしてみた</a>』が視覚的に分り易いのでそちらを参照してください。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">直積量子化</span><br />
ベクトル量子化はスカラ量子化よりも少ないビット数で効率良くベクトルを表現できるのですが、同時に欠点もあります。次元が増大するにつれて使用するコードブックが指数的に増えてしまうのです。<br />
<br />
これは次元の呪いに起因するものです。数次元〜数十次元くらいだったらそんなに問題とはならないんですけど、数百とか数千次元になってくるととんでもない量のコードブックが必要となり現実的でない。『<b>直積量子化(Product Quantization)</b>』はこの次元の呪いを解決するために用いられる、スカラ量子化とベクトル量子化の中間を行く手法です。<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiUMw9nhhuXNGN2GNR4YvuXIbdj09GSoWV8lVw148TCqXlu9EYPikVGjaiwH6sE4-fJZCcLjXnZ3O5eGZRwCr85IoDkECm2QNzJ6DA93RsEwnqA1fg2vW1A35Sbh6Yuoz-RkHTnqVsifMY/s1600/pic04.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiUMw9nhhuXNGN2GNR4YvuXIbdj09GSoWV8lVw148TCqXlu9EYPikVGjaiwH6sE4-fJZCcLjXnZ3O5eGZRwCr85IoDkECm2QNzJ6DA93RsEwnqA1fg2vW1A35Sbh6Yuoz-RkHTnqVsifMY/s1600/pic04.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">高次元のベクトルを直積表現によって分解する</td></tr>
</tbody></table>
<br />
具体的な手法は下記の通りです。まず、N次元のベクトルをM分割し、N/M次元のサブベクトルを生成します。そして、そのサブベクトルをそれぞれベクトル量子化することによって、インデックスを(i<span class="Apple-style-span" style="font-size: x-small;">1</span>, i<span class="Apple-style-span" style="font-size: x-small;">2</span>, ..., i<span class="Apple-style-span" style="font-size: x-small;">M</span>)のタプルで表現します。そうすれば次元の呪いは発生せず、データ表現も比較的低ビットで抑えられます。<br />
<br />
例として128次元のベクトルを考えてみましょう。まず、このベクトルを16分割することで8次元のサブベクトルへ落とし込みます。そして、そのサブベクトルをベクトル量子化によって8bit表現に持ってきた場合、対象のベクトルは8×16=128bitで表現でき、かつ使用するコードブックサイズも16×256=4096と少数で済みます。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">float</span>型のサイズで考えると30倍以上の圧縮効率です。このように、少数のコードブックサイズで効率よくベクトルを表現できることが直積量子化の利点です。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">距離計算</span><br />
さて、直積量子化によって各ベクトルのエントリを圧縮した状態で格納したとしましょう。あるクエリベクトルが与えられたとして、そのユークリッド距離が最小となるエントリを見つけるにはどうすればいいでしょうか?<br />
<br />
幸運なことに、ユークリッド距離の計算は各要素毎に独立しています。つまり、ただ単にクエリベクトルを各サブベクトルに分け、それぞれの領域で距離を計算し、その和を計算するだけで良いのです。よって、直積量子化を用いて最近傍探索を近似的に行うためには:<br />
<ol>
<li>クエリベクトルをM分割し、N/M次元のサブベクトルに分解する</li>
<li>クエリのサブベクトルと、それぞれのコードブックに記された代表ベクトル間の距離を予め計算しておく</li>
<li>それぞれのエントリがどのインデックスに位置しているのか確認し、2で計算した距離を参照して足し算を行う</li>
<li>最小距離をもつエントリを候補として提出する</li>
</ol>
<div>
ことで実現できます。これが直積量子化を用いた、近似最近傍探索の核となるアルゴリズムです。</div>
<br />
<span class="Apple-style-span" style="font-size: large;">詳細について</span><br />
この距離計算は非対称性からある程度のバイアスがかかっているため、実際の距離を正確に見積もるためにはこのバイアスを取り除く必要があったりします。また、具体的な距離計算がコードブックだけで済むので省力化されているとはいえ、計算量は依然としてO(N)であり、そこらへんも大規模なデータ数では足を引っ張る原因となります。<br />
<br />
んじゃあこれをどうやって解決するのかについては、現論文<br />
<br />
<a href="http://lear.inrialpes.fr/pubs/2011/JDS11/jegou_searching_with_quantization.pdf">http://lear.inrialpes.fr/pubs/2011/JDS11/jegou_searching_with_quantization.pdf</a><br />
<br />
を参照してください。上記2つの懸念についての解決案が記述されており、上の知識があれば多分すらすらと読めるはずです。<br />
本当はこれも解説しなきゃなんないんですけど、書いている内に疲れたちゃったのでとりあえずコアなとこだけということで勘弁してください。気が向いたら更新します。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-16969603892474936922011-09-24T01:44:00.000+09:002012-08-19T00:15:02.580+09:00OpenCVアルゴリズムを高速に記述できる補助ライブラリcvutilを公開します<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh8aKPMmYLYUJKhlzDchIf8EMQlRP9F8R7NfZrbZ9HHWOTl8sNJSSfuwYBMK9SVE7a27AFwZ4osi3ih5825R0_fW_KAKWgaHCNQR7BXeVJ1kpVdKqH3xgERFOvpJayhyke9t-7Z7xAFaTw/s1600/dunny.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh8aKPMmYLYUJKhlzDchIf8EMQlRP9F8R7NfZrbZ9HHWOTl8sNJSSfuwYBMK9SVE7a27AFwZ4osi3ih5825R0_fW_KAKWgaHCNQR7BXeVJ1kpVdKqH3xgERFOvpJayhyke9t-7Z7xAFaTw/s1600/dunny.jpg" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">"<a href="http://www.flickr.com/photos/kayveeinc/3137473002/">The Dunny Collective - Brisbane River</a>" by <a href="http://www.flickr.com/photos/kayveeinc/">KayVee.INC</a></td></tr>
</tbody></table>
<b><span class="Apple-style-span" style="font-size: large;">概要</span></b><br />
OpenCVを用いてアルゴリズムを組む際、<b>ある画像からある画像へ変換処理</b>を行ったり、<b>関数fを用いて新規に画像を生成</b>したり、<b>画像の結果を1つの値へと集約する</b>などの定型的な処理は非常に多く行われておりますが、並列処理を行いかつ速度を考慮したプログラムを組もうとするとある程度のコード量を必要とするためそれなりの手間となります。また画像処理によくあるループ文のネストが大量に発生するため、普通に組んでしまうと一見何をやっているのか非常に分かりづらいコードとなってしまいます。<br />
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cvutil</span>はこのような冗長な作業を定型化し、複数のアルゴリズム関数と<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">C++0x</span>に新しく搭載されたラムダ式を使用して高速にアルゴリズムを記述することができます。このようにして作られたアルゴリズムは自動的にOpenMPによって並列化されるため、研究用途などの速度をある程度担保しつつもメンテナンスのし易いコードを記述するための手助けとなります。<br />
<br />
<span class="Apple-style-span" style="font-size: large;"><b>インストール</b></span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cvutil</span>は<a href="https://github.com/">github</a>にて公開しております。ダウンロード先は下記のURLを参照してください。<br />
<br />
<a href="https://github.com/rezoo/cvutil">https://github.com/rezoo/cvutil</a><br />
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cvutil</span>はヘッダオンリーのライブラリです。バージョンがOpenCV 2.x系統であれば、対象の<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">include</span>ディレクトリに<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cvutil</span>をコピーするだけで使用できます。<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">C++0x</span>のラムダ式を用いることで端的に記述することができますが、関数オブジェクトを使用すればある程度古いコンパイラでも正常に動作するでしょう(そこまでして使う利点があるのかどうかは分かりませんが)。<br />
<br />
ちなみに、元々は自分用に作られた適当なライブラリなので、バージョンアップに伴ってある程度破壊的な変更を行うこともありますのでご了承ください。その際はバージョンを0.1上げることによって対応します。<br />
<br />
<span class="Apple-style-span" style="font-size: large;"><b>アルゴリズム一覧</b></span><br />
ここでは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cvutil</span>に搭載されているアルゴリズム一覧の一部を紹介いたします。これらのアルゴリズムは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">STL</span>のそれと似ているため、C++を触ったことのある方は名前から直感的に機能を類推できるでしょう。<br />
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b>transform_image</b></span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">transform_image</span>関数は対象の画像を、任意の関数を通して出力先へ変換します。cvutilの中では最も汎用性の高い関数となります。<br />
<pre class="prettyprint">cvutil::transform_image(src, dst, [](uchar x) {
return cv::saturate_cast<uchar>(x*2); });
</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgggRKz8iMo7kn2LXax38pgXFvOBFUaUzrsxIYFWVZ48uOJKIzNk0SG0AJ_g3qGielBV1Kbz7ZP0aOe8_hGo50U1dgN_IsRRMko9FpBPReBTZIl6Bx3qslM3cEHBBIFAr81eWUyv_s_Uvg/s1600/result.jpg" imageanchor="1"><img border="0" height="160" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgggRKz8iMo7kn2LXax38pgXFvOBFUaUzrsxIYFWVZ48uOJKIzNk0SG0AJ_g3qGielBV1Kbz7ZP0aOe8_hGo50U1dgN_IsRRMko9FpBPReBTZIl6Bx3qslM3cEHBBIFAr81eWUyv_s_Uvg/s320/result.jpg" width="320" /></a></div>
この関数は複数の画像を高度に合成するような場合においても利用できます。アリティ(渡される引数の数)によって、最も適したアルゴリズムが自動的に呼び出されます。<br />
<pre class="prettyprint">cvutil::transform_image(src1, src2, dst, [](uchar x, uchar y) {
return (x + y)/2; });</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhgJi6EcXPEzjLj7OFx3nH6wfSrSct6ttB2T_ta10J6pG8qAKOOOOsUhxkdEv684oyb0WDuRkyUZ2YFy4L4fn2YmnYQMoi9j_akJbxI-t39QgOCjya9M0QEZD9wRJobz8lrmAyZve9cASA/s1600/result2.jpg" imageanchor="1"><img border="0" height="132" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhgJi6EcXPEzjLj7OFx3nH6wfSrSct6ttB2T_ta10J6pG8qAKOOOOsUhxkdEv684oyb0WDuRkyUZ2YFy4L4fn2YmnYQMoi9j_akJbxI-t39QgOCjya9M0QEZD9wRJobz8lrmAyZve9cASA/s400/result2.jpg" width="400" /></a></div>
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b><br /></b></span>
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b>generate_image</b></span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">generate_image</span>関数は任意の関数を用いて、出力先に画像を新しく生成します。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">generate_image</span>をそのまま用いることは通常ありませんが、x, y座標を元に画像を生成する関数<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">generate_image_xy</span>や、画像中央に原点がある右手系の座標系u, v座標を元とした関数<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">generate_image_uv</span>は比較的使い勝手が良いと思います。<br />
<pre class="prettyprint">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
</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiQrBH2bGdV1Dnp2iAIgHx2a3R3EERrAHXatk5yo4hS6uJbdEYYQxTawUO1pKNBrpWbIQ_S2hX3UFecyTZGDmbB28OwIJ16JmqYijS_aILFN3MRHbItJRkPCOOvuY2zBK4VbdV0l0dqER0/s1600/result3.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiQrBH2bGdV1Dnp2iAIgHx2a3R3EERrAHXatk5yo4hS6uJbdEYYQxTawUO1pKNBrpWbIQ_S2hX3UFecyTZGDmbB28OwIJ16JmqYijS_aILFN3MRHbItJRkPCOOvuY2zBK4VbdV0l0dqER0/s200/result3.jpg" width="200" /></a></div>
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b>reduce_image</b></span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce_image</span>関数は対象の画像を、任意の関数を通して1つの結果へと集約し、その値を返します。<br />
<pre class="prettyprint">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</pre>
上のサンプルは全要素に対しての足し算を行っています。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">transfrom</span>の並列化は兎も角、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>の並列化は一般的に面倒なのですが、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce_image</span>はそこらへんの処理もきちんと並列化してくれます。<br />
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b>transform_reduce_image</b></span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">transform_reduce</span>関数は1枚または複数の画像を対象の関数で1つの値へ変換した後に、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>関数を適用します。2種類の結果を合併した集約結果を返す場合などに役立つと思います。<br />
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b>minmax_element_image</b></span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">minmax_element_image関数は画像の最小値と最大値を返します。</span><br />
<pre class="prettyprint">std::pair<uchar, uchar> minmax_result =
cvutil::minmax_element_image(src); // (min, max)</pre>
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b>integrate_image</b></span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">integrate_image</span>関数は画像を2次元に積分します。具体的には、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">inclusive_scan</span>関数を2次元的に適用します。バイナリー関数に<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">std::plus</span>を適用した場合、この関数は積分画像の生成と等価です。また、画像の型が対象の二項演算子に関してアーベル群を成していた場合、任意の矩形領域に関しての<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">integrate_image</span>によって生成された画像を用いて、定数時間による計算が可能となります(参考: <a href="http://mglab.blogspot.com/2011/04/reducescan.html">並列環境におけるreduceとscanアルゴリズム</a>)。<br />
<br />
OpenCVにも積分画像を生成する関数は存在しているのですが、この関数は並列化を行いかつ任意の二項演算子を適用できる点が異なります。
<br />
<pre class="prettyprint">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]</pre>
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cvutil</span>はその他にも色々と関数がありますが、多分名前からどんなことやってるのか大体は理解できると思いますので、とりあえず代表的な関数をいくつか紹介いたしました。<br />
<br />
<span class="Apple-style-span" style="font-size: large;"><b>License</b></span><br />
自分用に作っているだけなのですが、一応ライセンスは<a href="http://ja.wikipedia.org/wiki/MIT_License">MIT License</a>を適用します。誰でも自由に無制限に使ってかまいません。<br />
<br />
とりあえずはそんな感じです。もし何か使っていく上で質問や疑問、あるいはご指摘等ございましたらお気軽にコメントまたは<a href="http://www.blogger.com/profile/09706721371333005060">メール</a>をお願いします。<br />
<br />
<span class="Apple-style-span" style="font-size: large;"><b>追記</b></span><br />
お気づきの通り、本来の文脈ならば、アルゴリズムに渡す引数は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cv::Mat_<T></span>ではなく何かしらの抽象構造<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">View</span>を渡すべきなのです。C++のSTL Algorithmはコンテナでなくイテレータという抽象に依存することによって、任意のデータ構造による適用が可能となり、さらに<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">transformed_iterator</span>など、変形を遅延させるような構造を取ることができる。つまり、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cvutil</span>においても本来はそのような『<b>抽象に依存する</b>』データ構造を起点にアルゴリズムを組んでいくべきなのです。<br />
<pre class="prettyprint">// こんな感じで書けたら最高なんだけど…
cvutil::transform_image(
cvutil::make_resized_view(src, dst.size()),
dst,
any_functor);</pre>
ただこのアナロジーを画像に適用しましょうとなるとなかなか難しく、現状としては<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><a href="http://www.boost.org/doc/libs/1_47_0/libs/gil/doc/index.html">boost.GIL</a></span>しか存在していない。そういった意味で確かに<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">boost.GIL</span>は魅力的なライブラリなのですが、GILに適用するalgorithmが殆ど存在していませんし、なによりコンパイル時間が<b>boooooooooooost</b>してしまう。ぱぱっと書いてぱぱっと書き捨てられるような構造がむしろ研究にとっては重要なわけで、そんな面倒なことするよりも多少汚いOpenCVを採用したほうが使い勝手としては向上するよねという、まぁそういった言い訳です。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-29414476733209798502011-09-16T14:45:00.000+09:002011-09-24T01:47:51.902+09:00OpenCVの画素値アクセスについて<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi0-3UI9j9l621vL7thcWHxeCv2xAFgQ1fLTrdMU7pBT-ug15HaAWDRnsU3YLqk-iLdHxvIvSZ0u_60Et1HFFuCqFYBhVfSUqx_yIPXWCvXHGq78js-aV1RrxJW59APUFWRJS6DHdtH14M/s1600/fireking.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi0-3UI9j9l621vL7thcWHxeCv2xAFgQ1fLTrdMU7pBT-ug15HaAWDRnsU3YLqk-iLdHxvIvSZ0u_60Et1HFFuCqFYBhVfSUqx_yIPXWCvXHGq78js-aV1RrxJW59APUFWRJS6DHdtH14M/s1600/fireking.jpg" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">"<a href="http://www.flickr.com/photos/mkwilliams/4627921577/">Fire King</a>" by <a href="http://www.flickr.com/photos/mkwilliams/">A Continuous Lean</a></td></tr>
</tbody></table>
何度も何度も繰り返されている議題ですが、今回は画像処理に良くある『<b>全画素を変換する</b>』という処理のパフォーマンスについて考察します。<br />
あなたはOpenCVの<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cv::Mat_<T></span>を使って、全ての画素値をある関数<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">func</span>で変換するという処理を行いたいとします。それで、たぶんなんですけど、普通に何も考えず組むとしたらこんな感じのコードになると思います:<br />
<pre class="prettyprint">for(int y=0; y<src.rows; ++y) {
for(int x=0; x<src.cols; ++x) {
dst(y, x) = func(src(y, x));
}
}</pre>
なんの変哲もない普通のループです。このコードは普通に動きます。おわり。ちゃんちゃん。<br />
<br />
…そういきたいのですが、パフォーマンスを気にする方々はここでいろいろとした疑問が沸き起こってくるわけです。それじゃあどこが悪いのか。1つずつ疑問点を上げ、その上で上記のコードがある程度の合理性を持っているコードであることを示していきます。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">疑問1: 関数呼び出しを行っているからその分遅い?</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">src(y, x)</span>という処理は結局<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cv::Mat<T>::operator()</span>という関数を呼び出していることに相当します。なので、関数呼び出しを数十万のピクセルに対して行っているため、関数呼び出しにかかるコストが無視できないレベルで発生していくのではないか??<br />
<br />
<span class="Apple-style-span" style="font-size: large;">回答: 遅くならない</span><br />
OpenCVの<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cv::Mat_<T>::operator()</span>はインライン展開されます。なので、関数の呼び出しに対して発生する速度的なコストは発生しません。安心して使いましょう。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">疑問2: 掛け算を余計に行なっているために遅い?</span><br />
結局<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cv::Mat_<T>::operator()</span>の行っていることはコードに直すとこんな感じです。<br />
<pre class="prettyprint">src(y, x) == ((T*)(src.data + src.step.p[0]*y))[x]
</pre>
でも、この項をよく見てみると、xのループに対してyの値は固定されているため、第二項の計算は正直不要でしょう。なのに掛け算の計算でコストが余分にかかってしまうため、速度が遅くなるのでは??<br />
<br />
<span class="Apple-style-span" style="font-size: large;">回答: 遅くなるけどそんなに問題じゃない</span><br />
確かに遅くなりますが、そこまで問題となる箇所ではありません。時代は3GHzに突入しているのです。整数の掛け算にかかるコストなんてたかが知れています。<br />
<br />
そんなところよりも、別の箇所を改善したほうが大抵の場合劇的に向上します。仮に、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">func</span>の中に<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">std::exp</span>や<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">std::log</span>が含まれていたとすると…<br />
<pre class="prettyprint">dst(y, x) = std::exp(-(src(y, x) - mu)*(src(y, x) - mu)); // ぐぎゃー
dst(y, x) = -src(y, x)*std::log(src(y, x)); // おもいー</pre>
…それだけで大分重くなってしまいます。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">sqrt, sin, cos</span>はまだ我慢できるにしても(それでも掛け算よりは結構重い)、指数、対数関数は普通に使うとかなり重いんです。CV系統などの精度を求めないような計算でしたら、スプライン関数などの近似式に置き換えることを検討しましょう。<br />
<br />
<span class="Apple-style-span" style="font-size: large;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">at<T>(), operator()</span>は範囲チェックを行ってくれる</span><br />
昔のOpenCVには<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">at<T>()</span>と同じような処理をしてくれる<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">CV_IMAGE_ELEM</span>マクロなんて言うものもありましたが、これはもう時代遅れの代物です。なぜなら、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cv::Mat</span>の<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">at<T>()</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cv::Mat_<T></span>の<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">operator()</span>は与えられたx, y座標が範囲内のものであるかどうか、きちんとチェックを行ってくれるからです。範囲外である場合には警告を送ってくれます。そういった意味でも<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">CV_IMAGE_ELEM</span>より<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">at</span>や<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">operator()</span>を利用していくべきでしょう。<br />
<br />
『<b>余計な範囲チェックを行っているために遅くならないのか?</b>』という疑問もあるかもしれませんが、この範囲チェックは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Release</span>モードで消去されます。なので実際の動作には全く影響されないのです。<br />
<br />
<span class="Apple-style-span" style="font-size: large;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">
cv::Mat::at<T>(), cv::Mat_<T>::operator()</span>は早い</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">at<T>(), operator()</span>はきちんと考えられているメンバ関数です。ポインタを用いてがりがり書いたコードよりは少し遅いかもしれませんが、この遅さは全体のアルゴリズムに影響するような遅さではありません。ほっと胸を撫で下ろし、安心して使いましょう。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">でもやっぱり速度が気になる</span><br />
とはいっても、きちんと書かれていないという神経質な方もいるかもしれません。<b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">at<T>()</span>なんて負けた気がして使いたくない</b>……その場合のコードは大体以下のようになります。<br />
<pre class="prettyprint">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]);
}
}</pre>
これだけじゃあ高速化されたといっても微々たるものなので、せっかくですからOpenMPによる並列化を行ってみましょう。さらに<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">assert</span>でサイズが合っているかチェックを行うようにして…そうすると最終的なループはこんな感じになります。<br />
<pre class="prettyprint">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]);
}
}</pre>
…<b>大分行数が増えてしまいました。</b>ただ自分がやりたいことは全画素に対して<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">func</span>を適用させたい、ただそれだけなのに、ここまで行数を書く必要があるのです。しかもぱっと見で何をしたいのか良く分からない。面倒ですし分かりにくい。<br />
<br />
ここで、C++を少しかじったことのある人でしたら、<i>ユニタリー関数を適用させる表式のアルゴリズム関数を適用すればいいんじゃね?</i>といった考えに至るのは自然な発想でしょう。こんなふうに:<br />
<pre class="prettyprint">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]);
}
}
}
</pre>
上の形の式をヘッダファイルあたりに定義してやれば、後は以下のようなラムダ式を利用して書くだけで自動的に最適なループ処理を行ってくれます:
<br />
<pre class="prettyprint">transform_image(src, dst, [](uchar x) -> uchar {
return cv::saturate_cast<uchar>(x*2);
});</pre>
この表式ですと、どんな変数がなんの変数に影響を及ぼすのか、全画素に対してどのような処理を行うのか一目瞭然です。また、引数に渡す<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">unary_op</span>はコンパイラの最適化処理によって自動的にインライン展開されるため速度的なオーバーヘッドはありませんし、OpenMPによる並列化も行ってくれます。また任意の関数オブジェクトを作用させることができるので、非常に使い回しの効く関数となっています。<br />
<br />
画素云々について言及している方はこれまで数多く見かけましたが、このような画像に特化したアルゴリズム関数について言及している方は、何故か誰もいないというのが現状です。なのでまぁ自分用にざっくり書いているコードなのでいろいろと不備はありますけど、上記のようなアルゴリズムを詰め合わせた補助ライブラリを近いうちに上げていきたいと思います。<br />
<br />
追記(2011/09/24): 記事のとおり、補助ライブラリである<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><b>cvutil</b></span>を<a href="https://github.com/rezoo/cvutil">github</a>に上げました。詳細は別記事『<a href="http://mglab.blogspot.com/2011/09/opencvcvutil.html"><b>OpenCVアルゴリズムを高速に記述できる補助ライブラリcvutilを公開します</b></a>』を参照してください。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">補足</span><br />
『イテレータを使って標準のアルゴリズムを適用すればいいだけなのでは?』と思う方もいるかもしれません。つまりこんな感じですね:<br />
<pre class="prettyprint">std::transform(src.begin(), src.end(), dst.begin(), [](uchar x) -> uchar {
return cv::saturate_cast<uchar>(x*2);
});
</pre>
確かにこの式でも動くんですけど、速度は遅いです。なぜならこのイテレータの指しているデータは、メモリ上に連続していないからです。言い換えると、任意の自然数Nに対して『<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">&(*(iter + N)) == &(*iter) + N</span>』であることが保障されていないため、きちんと対象のデータを指し示すように補正しなければならない。その補正処理に余分な負担がかかって遅くなってしまう。ただ、もしメモリ上に連続であるならば直接ポインタをイテレータとして入れてやることで、(きちんと動作する保障はないですが)高速に動作してくれるでしょう。並列化はされていませんけれど。<br />
<br />
また、上記の<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">transform_image</span>は並列化は行われておりますが、SIMDによる高速化は行われておりません。これをジェネリックに行うというのは現状としてなかなか厳しいのですが、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Boost.SIMD</span>の登場によってある程度の希望の光は見えてくるのではないかと思っています。<br />
<br />
追記:<br />
<ul>
<li>(2011/09/16)コードをコンパイルが通るように若干の修正。ucharを暗黙的に使うのではなく明示的に指定するようにした</li>
<li>(2011/09/17)記事を書いていて初めて知ったんですけど、<b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">ptr<T>(y)</span>は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cv::Mat_<T></span>だと<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">operator[]</span>でオーバーロードされていた</b>んですね。恐らくこれを用いることがOpenCV側の意図しているところだろうと判断し、上の式を書き換えました。</li>
</ul>
rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-37696992349488209692011-06-26T01:20:00.006+09:002011-06-26T01:37:30.672+09:00boost.gilを利用して透過画像を作る<div style="text-align: center;">『<b>対象の画像から、特定の色の箇所だけを透明にした透過画像を作る</b>』</div><br />
こういった課題に一ヶ月近くかかって出来なかったとかいう話を聞いて、流石にそれはかかり過ぎだろうと脊髄反射で言ってしまったわけですけど、それじゃあ自分だったらどうだろうとか後で思ってちょっとやってみることにしました。好き勝手言って出来なかったとか恥ずかしすぎですから、やっぱりそこらへんは実際にやってみる責任みたいなものは感じているわけです。というわけでboost.gilを使って書いてみた結果が以下。<br />
<pre class="prettyprint">#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;
}
</pre>結果としてはこんな感じです。上が元の画像で、下が結果の画像。正常に抜き出されていることが分かります:<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjbkE0A_Fn2hPMvFtB7RW7Hczash8dHP8DminnTFO3jeb1z5Q7C4Ku9MLgPPcHV_Wmu8hrvtgPTeostizVNwFeNAb_vJUBlBt-jzsE2Ifm0XFyHt5DXymnegsnyZkGI4ljoERcu6D1YOUY/s1600/src.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="150" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjbkE0A_Fn2hPMvFtB7RW7Hczash8dHP8DminnTFO3jeb1z5Q7C4Ku9MLgPPcHV_Wmu8hrvtgPTeostizVNwFeNAb_vJUBlBt-jzsE2Ifm0XFyHt5DXymnegsnyZkGI4ljoERcu6D1YOUY/s200/src.png" width="200" /></a><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEikm5-RaQEfa9vHr9358jf7cIFeIDLx6ZTkYY0CQPo4IOBoAfyTk6MwzQcosZLeXD8ibhyphenhyphenYnSjGhYP8IxX5fWKZhLwlLmjxmzvVdjqE9UIWlQlfyYMOiQunSG5npeSd7izVVJeF3moHJYI/s1600/dst.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="150" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEikm5-RaQEfa9vHr9358jf7cIFeIDLx6ZTkYY0CQPo4IOBoAfyTk6MwzQcosZLeXD8ibhyphenhyphenYnSjGhYP8IxX5fWKZhLwlLmjxmzvVdjqE9UIWlQlfyYMOiQunSG5npeSd7izVVJeF3moHJYI/s200/dst.png" width="200" /></a></div><b>余談</b><br />
<ul><li>プログラムより文章を紡ぎ出すのに時間がかかって、改めて自分の文章力の無さにうんざりしています</li>
<li>というか最近暑い!自分は寝たいのに梅雨ですごい蒸し暑いから寝付けず気になってこんな不毛なことやっちゃってるわけです。今度こそ寝ます!!!</li>
</ul>rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-73780384969806914602011-06-02T00:40:00.002+09:002011-06-02T00:43:28.750+09:00Thrustから見る、reduceを用いたアルゴリズム実装<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEivya_Ubecs-1amULvVWFTTla6hjkHH0FiYcDYhrMXE1zuMmhEX7EA78Ycph7-6gSAi_RrOHotIQt0fmPgg5QwosyWt7iMi-MATcaQ5imGmOunjERpA_LNj27NjEPNVKwC1b1Z1HFWIUjo/s1600/John+Martin+-+The+Great+Day+of+His+Wrath.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEivya_Ubecs-1amULvVWFTTla6hjkHH0FiYcDYhrMXE1zuMmhEX7EA78Ycph7-6gSAi_RrOHotIQt0fmPgg5QwosyWt7iMi-MATcaQ5imGmOunjERpA_LNj27NjEPNVKwC1b1Z1HFWIUjo/s1600/John+Martin+-+The+Great+Day+of+His+Wrath.jpg" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">"The Great Day of His Wrath" by John Martin</td></tr>
</tbody></table><span class="Apple-style-span" style="font-size: large;">はじめに</span><br />
<br />
<b>世はまさに大並列時代</b>。火を手にした人類が飛躍的な進化を遂げたように、並列化のパラダイムは、今まで到底不可能と思われていた速さで計算を行うことができるようになりました。<br />
<br />
ところが、どんな手法にも欠点は存在するもので、実際に実装しようとすると非常に難しい。何故かというと、並列には並列化特有の問題が存在しており、愚直に実装してしまうとCPUより早くなってしまうどころか遅くなってしまうことだってあり得るのです。これを回避するにはGPUの内部構造についてきちんと理解をした上で、実装したいアルゴリズムそれぞれの場合に特化したコーディングを行う必要がある。<br />
<br />
しかしよくよく考えるとおかしな話です。<b>私たちが実装したいのはあくまで手法であり、ハードウェアではありません</b>。なぜこのような詳細について把握する必要があるのでしょう。これはきちんとした抽象化がなされていない結果生み出された、言ってみれば妖です。この妖によって余計な手間が生み出され、コーディング量は増大し、余計なバグが生成され、絶望した開発者はうつになり、自殺者は増え、世界は暗雲に包まれるー<br />
<br />
そのような混沌の中、ある一筋の光、<a href="http://code.google.com/p/thrust/">Thrust</a>が現れました。ThrustはCUDAを抽象化し、C++のパラダイムを用いることで本質にのみ注力できるよう開発されたライブラリです。ユーザはただ<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Iterator</span>に対してアルゴリズムを適用するかの如く<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">device_vector<T></span>と各種アルゴリズムを駆使することで、流れるようなプログラムを組むことが可能となります。このような抽象を用いることで、<b>ユーザは詳細な内部構造について頭を悩ませる必要が無くなりました</b>。誰が書いているのか分からないから使いたくない?ThrustはNVIDIAに所属している開発者<a href="http://research.nvidia.com/users/jared-hoberock">Jared Hoberock</a>と<a href="http://graphics.cs.uiuc.edu/~wnbell/">Nathan Bell</a>氏によって書かれています。2人とも並列化のプロフェッショナルなので、普通の人が愚直に実装するよりも大抵の場合効率的でしょう。<br />
<br />
さて、Thrustには並列化を行うための多種多様なアルゴリズムが存在しています。そして、Thrustはオープンソースです。これはつまり、このソースコードを読めば並列化のプロがどのようにして各種アルゴリズムの並列化を実装しているのかが分かるということです。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">Reductions</span><br />
<br />
ということで、最近はThrustのソースコードをちょくちょく読んでいるところです。まだ読み始めた最中なのですが、特に感動したのは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Reductions</span>の項でした。<br />
<br />
数々あるSTLのアルゴリズムを普通に実装しようとすると、一つ一つのアルゴリズムすべてに対して、コアレッシングなどのGPGPU特有の問題を回避したコードを記述しなければならず、非常に面倒です。Thrustはこの問題を<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce関数</span>によって解決しました。この関数は並列化が可能なので、これを用いて書かれた関数は自動的に並列化されます。つまり、Thrustはいわば『<b>全から一を返す</b>』アルゴリズムすべてを<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>を用いて記述することによって、<b>冗長な記述を回避</b>しているのです。<br />
<br />
この<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>を用いた並列化の実装について簡単に説明していくことが、今回の記事の目的です。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>を初めて聞いた方は前回の記事『<a href="http://mglab.blogspot.com/2011/04/reducescan.html">並列環境におけるreduceとscanアルゴリズム</a>』を参照してください。<br />
<br />
ここで、本来のThrustはC++によって記述しているのですが、そのままC++で書くと非常に冗長で解り辛くなってしまうので、今回は本質のみを抽出するためHaskellを用いて説明します。あまりHaskellに慣れていないので不適切な書き方等あるかもしれませんが、そのへんはご指摘頂けると幸いです。また、今回の<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>には<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">foldr</span>関数を代わりに用いています。分かると思いますが一応。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">実装</span><br />
<ul><li><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b>max_element, min_element, minmax_element</b></span></li>
</ul>その名の通りリストの中の最大値、最小値、最大と最小値のタプルを返します。<br />
<pre class="prettyprint">*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)</pre>ここらへんはまぁ普通に分かります:<br />
<pre class="prettyprint">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)
</pre><br />
<ul><li><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b>all_of, any_of, none_of</b></span></li>
</ul><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">all_of, any_of, none_of関数は</span>対象のリストと真偽値を判定する関数を引数に取り、<br />
<ul><li>すべて当てはまっている</li>
<li>一つ以上当てはまっている</li>
<li>すべて当てはまっていない</li>
</ul>場合にはTrueを返す関数です。pythonでは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">all(), any()</span>, C++ではまんまのやつがありますね。<br />
<pre class="prettyprint">*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
</pre>これも実装的にはらくちん:<br />
<pre class="prettyprint">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
</pre>関数合成を用いてエレガントに解決。<br />
<ul><li><b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;">inner_product, transform_reduce</span></b></li>
</ul><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">inner_product</span>は2つのリストに二項演算子を作用させてから、その結果からなるリストを<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>する関数です。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">transform_reduce</span>はリストをいったんユニタリー関数を用いて変形させてから、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>を実行する関数です。<br />
<pre class="prettyprint">*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
</pre>C++だといろいろ冗長に書かないといけないですけど、Haskellだったら余裕です:<br />
<pre class="prettyprint">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
</pre><ul><li><b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;">count_if, count</span></b></li>
</ul><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">count_if</span>は条件に合致している要素数を返す関数で、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">count</span>は指定した値と同じ要素の数を返す関数です。<br />
<pre class="prettyprint">*Main> count_if (\x -> x < 3) [1,2,3,4,5]
2
*Main> count 3 [1,3,2,2,3]
2</pre><pre class="prettyprint">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
</pre>Thrustでは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">True</span>だったら1, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">False</span>だったら0となるリストを作成して、それを足していく手法をとっていました。<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">count</span>は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">count_if</span>が思いつけば楽勝でしょう。<br />
<ul><li><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: large;"><b>find_if, find_if_not</b></span></li>
</ul><div class="separator" style="clear: both; text-align: center;"></div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">find_if</span>はリストの中で条件に合致している『一番最初の要素』を示すイテレータを返します。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">find_if_not</span>は逆で、条件に合致していない『一番最初の要素』を示すイテレータを返します。今回は単純に要素のラベル番号を返すようにしました。なかったら最後のラベル+1が帰ってきます。<br />
<pre class="prettyprint">*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
</pre>これは微妙に入り組んでいます:<br />
<pre class="prettyprint">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
</pre>ただ、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">find_if</span>はちょっと本家と違います。インデックスが小さい項が常に左の引数となるため省略してるんだとは思いますけど、対称性を満たしていないのがちょっと気持ち悪いので、ここでは交換則を満たすように関数を変形しています。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">実際の実装と課題</span><br />
ここではhaskellを用いて実装しましたが、もちろん実際の実装は異なっており、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">zip</span>は<a href="http://wiki.thrust.googlecode.com/hg/html/classthrust_1_1zip__iterator.html"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">zip_iterator</span></a>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">map</span>は<a href="http://wiki.thrust.googlecode.com/hg/html/classthrust_1_1transform__iterator.html"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">transform_iterator</span></a>, タプルは<a href="http://wiki.thrust.googlecode.com/hg/html/classthrust_1_1tuple.html"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">tuple</span></a>を用いて実装しています。これらのイテレータの評価は遅延的に行われるので、複数のメモリアクセスが発生せず、アルゴリズムの高速化に貢献します。<br />
<br />
これらの<a href="http://wiki.thrust.googlecode.com/hg/html/group__fancyiterator.html"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">fancy iterator</span></a>を使用した戦略は非常に賢いとは思いますが、同時に欠点もあります。複数のイテレータを大量に組み合わせるために、通常のSTL algorithmのように<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">begin</span>と<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">end</span>を使用してしまうと、両方に対して<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">fancy iterator</span>を組み合わせなければならないため、2倍の手間がかかってしまうのです。<br />
<br />
これらの解決案としては<a href="http://www.boost.org/doc/libs/1_46_1/libs/range/doc/html/index.html"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">boost::range</span></a>のようなRange構造を実装することです。これをパイプ演算子で繋げていくような形をとっていけば、明快で分かりやすく、速度も十分保証されたプログラムを記述できることでしょう。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-27750516862655624462011-05-18T01:55:00.004+09:002011-05-18T02:01:13.145+09:00Wiresharkを用いたHTMLの初歩的なスクレイピング<a href="http://mglab.blogspot.com/2011/05/google-art-project.html">前回</a>はぽんとスクリプトを上げて「これで画像がダウンロードできますよ」とやったわけだけど、もちろん実際に<a href="http://www.googleartproject.com/">Google Art Project</a>の資料が上げられているわけでもないので、きちんとHTMLやJavascriptのソースを見ながら読解していくことになる。<br />
<br />
簡単な動作機構だったら単にHTMLのソース見るだけで良いのだけど、<a href="http://www.googleartproject.com/">Google Art Project</a>のような複雑な機構だとそうもいかない。行数が多いと面倒だし大抵の場合難読化されてたりで一筋縄ではいかないからだ。そこで、今回はパケット解析というもうひとつの武器を利用することで内部構造を別の側面から理解することにした。思考の流れを書き綴っていくので、HTMLのスクレイピングをやってみたい方にとって何か参考になれば幸いです。<br />
<br />
<div style="text-align: center;"><b>目標: 高画質な絵画をGoogle Art Projectを利用して手に入れる</b></div><br />
まず<a href="http://www.googleartproject.com/">Google Art Project</a>のサイトへアクセス。絵画ごとに固有のURLが割り当てられており、このURLを解析すれば絵画を入手できると判断。見るからにソースコードからの解析は面倒くさそうだったので、パケット解析ソフトウェア<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><a href="http://ja.wikipedia.org/wiki/Wireshark">Wireshark</a></span>を用いて解析を行った。<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi2SfOXEFBlogniEHgqbL_92bfBoAxbEIQSGVDgv3zVSJusH0RWdu5sC_KOUP3G3Sdb1LLD5oRN_Ory4Z6h6LjJMIY3F0J5Fos_oXoLsudNuYOWnoim5i-DTcqisuyE-Osr161RX09BUsA/s1600/initial.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="276" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi2SfOXEFBlogniEHgqbL_92bfBoAxbEIQSGVDgv3zVSJusH0RWdu5sC_KOUP3G3Sdb1LLD5oRN_Ory4Z6h6LjJMIY3F0J5Fos_oXoLsudNuYOWnoim5i-DTcqisuyE-Osr161RX09BUsA/s400/initial.jpg" width="400" /></a></div>本質のみに注力するため<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">http</span>プロトコルのみキャプチャ。適当なURLにアクセスし、ズームを1回行い、パケットキャプチャを中断。解析を開始。<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjszwTJNX2Y6iILJvQb2-rxm2BwmXYy5KALiO9ULGkNA-wWjnOhhlGDVufcLtBTSUDqj935ym-u649LuWIDjVScnBpOHeMgcA2Ms57FkoBDJGp1L-aw8dpiIvRuZWvmPtQyoaIms-qE3LQ/s1600/wireshark.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="356" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjszwTJNX2Y6iILJvQb2-rxm2BwmXYy5KALiO9ULGkNA-wWjnOhhlGDVufcLtBTSUDqj935ym-u649LuWIDjVScnBpOHeMgcA2Ms57FkoBDJGp1L-aw8dpiIvRuZWvmPtQyoaIms-qE3LQ/s400/wireshark.png" width="400" /></a></div>画像本体らしきURLにアクセスしているパケットを発見。複数のアクセスがあるってことはたぶん画像は複数に分割されていてそれを読み込んでいるんだろう。調べると<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">/unique-id=xi-yj-zk</span>の形でアクセスしていたので、対応する<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">xy</span>座標とズーム値<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">z</span>によって対応する画素値にアクセスすることができると推測。試しにURLにアクセスしてみると、問題なくダウンロードができた。<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEik5LlISsu3kOapIVHCK4H0oEXgApf2LW_M73VmVMvjCl-BZBGhYY83yWDkiK-v2EFs7pWRCPjSUom8jZORIKOtr_4uSKOUsu4UTp9W_KtthbbuZwoehUMV0RdsUI4BaVOzE1CTN9G6u6M/s1600/unnamed.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEik5LlISsu3kOapIVHCK4H0oEXgApf2LW_M73VmVMvjCl-BZBGhYY83yWDkiK-v2EFs7pWRCPjSUom8jZORIKOtr_4uSKOUsu4UTp9W_KtthbbuZwoehUMV0RdsUI4BaVOzE1CTN9G6u6M/s200/unnamed.jpg" width="200" /></a></div>また、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">x,y</span>の値を変更してのアクセスを行った結果も問題なくアクセスできた。ただ、どうやら範囲外の<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">xy</span>値にアクセスすると404が返されるらしい。当たり前といえば当たり前。アクセス禁止などのペナルティはないっぽいので、(ちょっと礼儀悪いけど)404が返されたらそれが範囲外と解釈することにする。<br />
<br />
一旦まとめてみる。<b>画像は分割されており、これらのアクセスはなんかよく分からない<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">unique-id</span>という、おそらく絵画について表しているパラメータと<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">x,y,z</span>座標の2つだけでダウンロードができる</b>。ついでに、<b>リファラー, cookie, 認証などの面倒な作業をせずとも直接のダウンロードが許されている</b>ことも分かる。そこらへんはGoogleさん寛容なんだなぁと思い、まぁAPIを提供してるくらいだからなぁ、とかいろいろ考える。<br />
<br />
他の画像に対してのホストを確認してみると、どうやら<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">[lh3〜lh6]</span>までの複数のホストへとアクセスしているらしい。このホストアクセスがランダムか、それとも固有のホストを呼び出さなければならないのか確認してみたかったので、試しにホストだけを変更して、その他のパラメータは固定したままでのアクセスを試みる。結果としては問題なくダウンロードできた。たぶん負荷を分散するという単純な理由から、ランダムに<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Javascript</span>側で割り振っているんだろう。ここで初めてHTMLのソースを閲覧。"<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cmn/js/content.js</span>"のソースを開き(難読化はされてないようだ)"<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">ggpht</span>"を検索。<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj_6xeRdVFU8WhwSoF4OxCTGfCeuWoZOxMz9Zd05Bk6BYikZbJHAv6UFrrBuHJzubNh2FvO6OXVGzUfXMCn-tQEkUgV4P8J9pzAi4Dv-Bx-RyYa7nEF3BaGgt5b7tVTBGncVri4aALJmiI/s1600/ggpht.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="222" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj_6xeRdVFU8WhwSoF4OxCTGfCeuWoZOxMz9Zd05Bk6BYikZbJHAv6UFrrBuHJzubNh2FvO6OXVGzUfXMCn-tQEkUgV4P8J9pzAi4Dv-Bx-RyYa7nEF3BaGgt5b7tVTBGncVri4aALJmiI/s320/ggpht.png" width="320" /></a></div>リストに入れているところからすると、画像へのアクセスはたぶんこのリストの中だけに限られているんじゃないかなと思う。試しに<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">lh7</span>へアクセス。404が正常に返された。的中。この指針で間違いないので次へ進む。<br />
<br />
次に気になるのは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">unique-id</span>だ。一体何処から仕入れているのか知らないと全く手のつけようがない。考えられるのは2つで、<br />
<ul><li>アクセス毎に固有の<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">unique-id</span>を発行し、一定期間のみアクセスできるようにしている</li>
<li>HTMLかどこかに参照しやすいように書かれており、それを抜き出してアクセスする</li>
</ul>という発想だ。自分は前者を疑っていたので、それっぽいパケットがあるかどうか監視。結果としてはまったく存在していなかった。というと後者かなと思い、対象のHTMLから<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">unique-id(今回の場合だと"vtzd432xqzOpi_3X8JAJ_buly36QKI0n9zGeeWNgBcwsYJTsAKK5mbM9Pgg")</span>を検索。<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh4dcZvfzKM_TbWshWIAQ899xrgqY6Vf8LeS-euubjbENHZk_55dHscKZxJHnVUXbC1pm-ApxoJzDHRohopMdhyeJNH2bViwADgDzd_yfRcCFeD3VdO3m8E_QQ1TDdTLz1bTucbdD9xS1o/s1600/unique.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="222" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh4dcZvfzKM_TbWshWIAQ899xrgqY6Vf8LeS-euubjbENHZk_55dHscKZxJHnVUXbC1pm-ApxoJzDHRohopMdhyeJNH2bViwADgDzd_yfRcCFeD3VdO3m8E_QQ1TDdTLz1bTucbdD9xS1o/s320/unique.png" width="320" /></a></div>ビンゴ!これで駒はすべて揃った。まとめに入る。対象の絵画をダウンロードするためには、<br />
<ol><li>取得したい絵画のアドレスへアクセス</li>
<li><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><div id="microscope" data-microId="unique-id"></span>となっている項を取得、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">unique-id</span>を得る</li>
<li><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">http://lh3.ggpht.com/ - http://lh6.ggpht.com/</span>のいづれかにアクセス。<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">http://lh3.ggpht.com/unique-id=xi-yj-zk</span> のような形で画像を取得</li>
<li>取得した複数の画像を繋ぎ合わせる</li>
</ol>後は非常に簡単でただこの単純なアルゴリズムに従って書き下していくだけ。その産物が<a href="http://mglab.blogspot.com/2011/05/google-art-project.html">こいつ</a>。こんなに早く仕上がるってことは、たぶんGoogleさんは本格的にダウンロードを禁止したいわけじゃないんだろう。<br />
<br />
スクレイピングとしては楽な部類で、普通だともうちょっと複雑な認証をかましてくるので自動化しにくい。<a href="http://ja.wikipedia.org/wiki/CAPTCHA">Captcha</a>入ると流石に無理だけど、おそらく現在の文字認識技術をきちんと応用すれば解けるんじゃないかと思っている。<br />
<br />
<b>余談</b><br />
<br />
<ul><li>こういう解析は楽しい。どんな感じの楽しさかっていうと、与えられたパズルを解いていく感覚にちょっと近い</li>
<li>自分はだいたい1時間くらいで解いたけど、早い人はたぶんもっと早いんだろうなって思う。まだまだ精進が足りない</li>
</ul>rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-39047432855545024432011-05-16T01:20:00.002+09:002011-05-16T01:22:08.404+09:00Google Art Projectを利用して高画質の絵画を手に入れる<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhQDoCyMony9zLzMgr64TcJ6CXlaPwwM1Q2-r_l7txVSuHsnkZ_940C8JN2iNxCkiS9UUELX0QVivILljtqEnazHyD0U5nuJdZmu5F7BKOlnU7i2tEowECWj_YHYyVZoxvj-qsF-MMSm8c/s1600/Vincent+van+Gogh+-+The+Starry+Night.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhQDoCyMony9zLzMgr64TcJ6CXlaPwwM1Q2-r_l7txVSuHsnkZ_940C8JN2iNxCkiS9UUELX0QVivILljtqEnazHyD0U5nuJdZmu5F7BKOlnU7i2tEowECWj_YHYyVZoxvj-qsF-MMSm8c/s1600/Vincent+van+Gogh+-+The+Starry+Night.jpg" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">"<a href="http://en.wikipedia.org/wiki/The_Starry_Night">The Starry Night</a>" by <a href="http://en.wikipedia.org/wiki/Vincent_van_Gogh">Vincent van Gogh</a></td></tr>
</tbody></table>2011年2月2日、Googleさんは<a href="http://www.googleartproject.com/">Google Art Project</a>というサイトをリリースしました。どういうサイトかっていうと<a href="http://www.moma.org/">MoMAニューヨーク近代美術館</a>や<a href="http://www.tate.org.uk/britain/">Tate Britain</a>など、世界中の美術館にある有名な絵画を渡り歩くことができるという趣旨のサイトで、誰もが一度は見たことがある絵画がー家に居ながらー非常にリアリスティックな解像度で見ることができる。恐ろしい時代になったもんです。<br />
<br />
ところがこのサイト、高解像度で見られるのは非常にいいことなのですが、<b>肝心のダウンロードができない</b>。せっかく高解像度なんですから保存できてもいいはずなのに、だけどできない。<a href="http://maps.google.co.jp/">Google Map</a>なんかは常に最新の情報に更新しなければならないので、ダウンロードに制限をかけているのは合理性があるんですけど、絵画とか全く更新する必要がないでしょう。著作権も切れてるのに。よく分からないです。<br />
<br />
ということで、なければ自分でなんとかしようという精神で、1時間くらいでちゃっちゃとPythonでスクリプト組んでダウンロードするようにしてみました。即興で組んだのでいろいろ粗い面はありますけど、そこら辺は勘弁してください。<br />
<br />
<pre class="prettyprint">#!/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())
</pre><br />
使い方は簡単で、<br />
<pre class="prettyprint">$ python artproject.py http://www.googleartproject.com/museums/tate/mariana-248
</pre>のようにURLを指定してあげると勝手に取得して結合してくれます。こんな感じで:<br />
<div class="separator" style="clear: both; text-align: center;"></div><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEip8_wJt8HuUgwR1fw6pFOkbh1Gcbu2XJeOkTaWxvLmToDkkznfSHIu3YD1lJeCyjHsWKStKkiTZNCp4IBPkbg1c07uLZhc8Wxu_6VsPoCHreBdBeDiwR_UiraQ8KUNCplUoGZHqF_ZTYs/s1600/Sir+John+Everett+Millais+-+Mariana.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEip8_wJt8HuUgwR1fw6pFOkbh1Gcbu2XJeOkTaWxvLmToDkkznfSHIu3YD1lJeCyjHsWKStKkiTZNCp4IBPkbg1c07uLZhc8Wxu_6VsPoCHreBdBeDiwR_UiraQ8KUNCplUoGZHqF_ZTYs/s1600/Sir+John+Everett+Millais+-+Mariana.jpg" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">"<a href="http://www.googleartproject.com/museums/tate/mariana-248">Mariana</a>" by <a href="http://ja.wikipedia.org/wiki/%E3%82%B8%E3%83%A7%E3%83%B3%E3%83%BB%E3%82%A8%E3%83%B4%E3%82%A1%E3%83%AC%E3%83%83%E3%83%88%E3%83%BB%E3%83%9F%E3%83%AC%E3%83%BC">John Everett Millais</a></td></tr>
</tbody></table>あぁそれにしてもお美しい……巷では<a href="http://www.salvastyle.com/menu_pre_raphael/millais_ophelia.html">オフィーリアちゃん</a>が有名だけど、やっぱり<a href="http://www.asahi.com/millais/detail/gallery.html">マリアナさん</a>のほうが綺麗だよ。絶対。<br />
<br />
<b>余談</b><br />
<ul><li>もともとの動機としては高解像度な絵画の壁紙が欲しくて、いろいろと探し回ったわけですが、どうも<a href="http://www.googleartproject.com/">Google Art Project</a>しか高解像度の絵画が載っていないので仕方なくといった次第です。</li>
<li>利用はどうやら<a href="http://www.google.com/accounts/TOS">Google's Terms of Service</a>に準じているようです。ざっと見した感じだとダウンロードを禁止する条項はないですし著作権も切れているので全く問題ないように思えるんですが、もし警告とかきたらチキンですから遠慮無く取り下げます。ご了承ください。</li>
<li>このネタを題材に時間があればHTMLのスクレイピング、パケット解析についての入門記事を書いてみたいです。</li>
<li>なんだかんだ言ってちゃんと実物見たほうが綺麗。ミレイの作品は以前東京に行ったときに<a href="http://www.bunkamura.co.jp/">Bunkamura</a>で見ました。Bunkamuraはそういう催し物よくやっていて、そこらへんが東京羨ましかったり</li>
</ul>rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com2tag:blogger.com,1999:blog-1062283300115822026.post-16122128434841646912011-04-19T01:28:00.001+09:002011-04-19T12:49:07.357+09:00並列環境におけるreduceとscanアルゴリズム<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiuw7Gy7MtArBak-6bAIdbx8aNrIvRyQtWLWJmnOWnYYz-1VJWN-C6ccct8J8v9ktin-vizjq49DuVjKpTqWfBfbI7xWAYigKqNOytSammS8Z-tBA9tLcQKKfzHX9TGRMGJt_hx2wUxHcg/s1600/fireking.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiuw7Gy7MtArBak-6bAIdbx8aNrIvRyQtWLWJmnOWnYYz-1VJWN-C6ccct8J8v9ktin-vizjq49DuVjKpTqWfBfbI7xWAYigKqNOytSammS8Z-tBA9tLcQKKfzHX9TGRMGJt_hx2wUxHcg/s1600/fireking.jpg" style="background-attachment: initial; background-clip: initial; background-color: initial; background-image: none; background-origin: initial; background-position: initial initial; background-repeat: initial initial; border-bottom-width: 0px; border-color: initial; border-left-width: 0px; border-right-width: 0px; border-style: initial; border-top-width: 0px;" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">"<a href="http://www.flickr.com/photos/cmyk1219/3064510251/">fireking</a>" by <a href="http://www.flickr.com/photos/cmyk1219/">cmyk1219</a></td></tr>
</tbody></table><span class="Apple-style-span" style="font-size: large;">1. 概要</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>と<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>はGPGPUなどの並列環境下において重要な役割を果たす関数だけど、あまりこれらの関数、特にscanについて言及している記事はあまり見かけない。なので今回は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>と<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>について調べた結果をまとめていきたいと思う。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">2. reduce, scan関数の概要</span><br />
<b>2.1. reduce</b><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>については<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>が分かればより理解しやすくなるので、ひとまずは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>から始めることにする。<br />
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>とは、端的に言うと、対象のリストを与えられた二項演算子で纏め上げる関数だ。とは言っても、いきなりそんなことを言われてもよく分からない。まずは例を見ていこう:<br />
<pre class="prettyprint">>>> 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</pre><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>において重要なのは二項演算子であり、上記の場合、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>は以下のような演算を行う。<br />
<pre class="prettyprint">(((1 + 2) + 3) + 4) + 5 = 15
(((1 - 2) - 3) - 4) - 5 = -13
(((1 * 2) * 3) * 4) * 5 = 120</pre>これにより、最初の演算は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">sum</span>関数、3番目の関数は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">prod</span>関数と等価であることが分かる。<br />
リストの型は何も数値に限られている訳ではない:<br />
<pre class="prettyprint">>>> reduce(operator.and_, [True, False, True])
False
>>> reduce(operator.or_, [True, False, True])
True</pre>これより、上記2つの演算は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">all</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">any</span>関数と機能的に等価であることが分かるのだけれど、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">all</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">any</span>関数のほうが効率的なので、実際はちゃんと<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">all</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">any</span>関数を使うべきだ。<br />
<br />
他にも<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>を用いて様々な関数、アルゴリズムを記述できるが、あまりにも多いので今回は後回しにして、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>を重点的に解説していく。<br />
<br />
<b>2.2. scan</b><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>関数は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>と比べて文献が少ない。言及している文章もあまり見かけないが、並列処理において<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>は重要な関数だ。残念ながらpythonで実装されていないけれど、仮に実装されていたとしたら以下のような働きをする:<br />
<pre class="prettyprint">>>> 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]</pre>きちんと書き下せばどのような働きをしているのか分かりやすい:<br />
<pre class="prettyprint">[1, 1+2, (1+2)+3, ((1+2)+3)+4, (((1+2)+3)+4)+5]</pre>結局のところ、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>は先頭から<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>していった結果をリストにまとめ、それを返す関数となる。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>は別名<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">prefix sum</span>とも呼ばれている。<br />
ちなみに、haskellでは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scanl1(scanl)</span>という名前で実装されている。<br />
<pre class="prettyprint">Prelude> :type scanl1
scanl1 :: (a -> a -> a) -> [a] -> [a]
Prelude> scanl1 (+) [1,2,3,4,5]
[1,3,6,10,15]</pre>また、<a href="http://www.nvidia.co.jp/object/cuda_home_new_jp.html">CUDA</a>ライブラリである<a href="http://code.google.com/p/thrust/">thrust</a>には<a href="http://thrust.googlecode.com/svn/tags/1.1.0/doc/html/group__prefixsums.html"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">inclusive_scan</span>(<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">exclusive_scan</span>)</a>という名前で実装されている。<br />
<pre class="prettyprint">#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}</pre><br />
<span class="Apple-style-span" style="font-size: large;">3. 並列による制限</span><br />
何故並列下において<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>が重要になってくるのか、それは、両方ともメニーコアな環境において効率的に動作するアルゴリズムが複数提唱されているからだ。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>は並列と相性がいい。つまり、既存のアルゴリズムをどうにかして<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>に落とし込むことができれば、並列化の恩恵を受けれられると。<br />
<br />
ただし、並列環境においては、演算にいくつかの制限が生じてくる。具体的には、集合Gと演算子"・"の組において、<br />
<ul><li>演算子"・"は交換則"a・b = b・a"を満たしていなければならない</li>
<li>演算子"・"は結合則"(a・b)・c = a・(b・c)"を満たしていなければならない</li>
</ul>後者に関してはあたり前の話で、これが成り立っていなければ並列に計算を行うことができない。前者については多少曖昧で、実際には満たさなくてもよいアルゴリズムも作ることはできるけど、コアレッシングなどの問題を考えると満たしていることが強く求められる。実際、thrustの<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>は交換則を満たしていることを前提としている。<br />
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>の場合はこの2つの制限で十分だけれど、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>の場合は上記2つの制限に加えてさらに2つの制限を満たしておくと都合がいい。具体的には、集合Gと演算子"・"の組において、<br />
<ul><li>ただ一つの単位元"e"が存在する</li>
<li>演算子"・"は任意の元"a"に対して"a・a<sup>-1 = </sup>a<sup>-1</sup>・a = e"を満たす、ただ一つの逆元"a<sup>-1</sup>"が存在する</li>
</ul>これら4つの制限を満たしている集合Gと演算子"・"の組を『<b><a href="http://ja.wikipedia.org/wiki/%E3%82%A2%E3%83%BC%E3%83%99%E3%83%AB%E7%BE%A4">アーベル群</a></b>』という。要するに、scan関数は、対象の型から成る集合と二項演算子の組がアーベル群を成していれば都合がいいということになる。<br />
<br />
具体例を出そう。整数と和の組<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">(Integer, (+))</span>は明らかに上記4つの制限を満たす。よって、この組は並列に<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>オペレーションを適用でき、かつ都合が良い。<br />
<pre class="prettyprint">1 + 2 == 2 + 1
(1 + 2) + 3 == 1 + (2 + 3)
1 + (-1) == (-1) + 1 == 0</pre><br />
真偽値と論理積の組<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">(Bool, (&&))</span>は上記2つの制限を満たすため並列に<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>オペレーションを実行可能だが、後半2つの制限は満たしていないため都合は良くない。<br />
<br />
3次元の回転行列と積の組<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">(Matrix, (*))</span>は交換則を満たしていないため―少なくともthrustでは―並列に<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce, scan</span>オペレーションを実行することはできない。<br />
<br />
<b>3.1. 都合が良いとは?</b><br />
とは言っても、一体アーベル群であればどこがどう都合がいいのか?それは、リスト中において局所的に<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>を行いたい場合に役に立つ。<br />
<br />
リスト<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">[a<sub>0</sub>, ... , a<sub>n</sub>]</span>にscanオペレーションを実行して<br />
<a href="http://formula.s21g.com/?%26%26%20%7B%5Crm%20scan%7D%5C;%5C;(%5Coplus)%5C;%5C;[a_0,%20%5Ccdots%20,%20a_n]%20%5C%5C%0A%26%26%20%3D%20[a_0,%20a_0%20%5Coplus%20a_1,%20%5Ccdots%20,%20a_0%20%5Coplus%20%5Ccdots%20%5Coplus%20a_n]%20%5C%5C%0A%26%26%20%5Cequiv%20[b_0,%20%5Ccdots,%20b_n]"><img src="http://formula.s21g.com/?%26%26%20%7B%5Crm%20scan%7D%5C;%5C;(%5Coplus)%5C;%5C;[a_0,%20%5Ccdots%20,%20a_n]%20%5C%5C%0A%26%26%20%3D%20[a_0,%20a_0%20%5Coplus%20a_1,%20%5Ccdots%20,%20a_0%20%5Coplus%20%5Ccdots%20%5Coplus%20a_n]%20%5C%5C%0A%26%26%20%5Cequiv%20[b_0,%20%5Ccdots,%20b_n].png" /></a><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"></span>を手に入れたとする。この場合、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">i < j</span>において<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">b<sub>j</sub>+b<sub>i</sub><sup>-1</sup></span>を計算すると、<br />
<a href="http://formula.s21g.com/?%26%26%20b_j%20%5Coplus%20b_i%5E%7B-1%7D%20%5C%5C%0A%26%26%20%3D%20(a_0%20%5Coplus%20%5Ccdots%20%5Coplus%20a_j)%20%5Coplus%20(a_0%20%5Coplus%20%5Ccdots%20%5Coplus%20a_i)%5E%7B-1%7D%20%5C%5C%0A%26%26%20%3D%20(a_0%20%5Coplus%20%5Ccdots%20%5Coplus%20a_j)%20%5Coplus%20(a_i%5E%7B-1%7D%20%5Coplus%20%5Ccdots%20%5Coplus%20a_0%5E%7B-1%7D)%20%5C%5C%0A%26%26%20%3D%20(a_0%20%5Coplus%20a_0%5E%7B-1%7D)%20%5Coplus%20%5Ccdots%20%5Coplus%20(a_i%20%5Coplus%20a_i%5E%7B-1%7D)%20%5Coplus%20(a_%7Bi+1%7D%20%5Coplus%20%5Ccdots%20%5Coplus%20a_j)%20%5C%5C%0A%26%26%20%3D%20a_%7Bi+1%7D%20%5Coplus%20%5Ccdots%20%5Coplus%20a_j"><img src="http://formula.s21g.com/?%26%26%20b_j%20%5Coplus%20b_i%5E%7B-1%7D%20%5C%5C%0A%26%26%20%3D%20(a_0%20%5Coplus%20%5Ccdots%20%5Coplus%20a_j)%20%5Coplus%20(a_0%20%5Coplus%20%5Ccdots%20%5Coplus%20a_i)%5E%7B-1%7D%20%5C%5C%0A%26%26%20%3D%20(a_0%20%5Coplus%20%5Ccdots%20%5Coplus%20a_j)%20%5Coplus%20(a_i%5E%7B-1%7D%20%5Coplus%20%5Ccdots%20%5Coplus%20a_0%5E%7B-1%7D)%20%5C%5C%0A%26%26%20%3D%20(a_0%20%5Coplus%20a_0%5E%7B-1%7D)%20%5Coplus%20%5Ccdots%20%5Coplus%20(a_i%20%5Coplus%20a_i%5E%7B-1%7D)%20%5Coplus%20(a_%7Bi+1%7D%20%5Coplus%20%5Ccdots%20%5Coplus%20a_j)%20%5C%5C%0A%26%26%20%3D%20a_%7Bi+1%7D%20%5Coplus%20%5Ccdots%20%5Coplus%20a_j.png" /></a><br />
となる。つまり、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">i+1</span>から<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">j</span>までの要素についての<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reduce</span>を僅か<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">O(1)</span>で求めることができる。<br />
<br />
これには大きな利点がある。例えば、大規模な時系列データを、その周りのデータから平均化(平滑化)したいとする。各要素ごとに計算を行うのは明らかに非効率であるので、まず<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>オペレーションを<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">(+)</span>に対して適用し、その後で差分を計算し、正規化を行ってやればよい。結局、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan (+)</span>のやっていることはf(x)の不定積分F(x)を求めることと本質的に等価となる。この考え方はCVでは、<a href="http://www.kec.jp/committee/johoshi/pdf/jyohoshi_210-3.pdf">積分画像</a>として様々な箇所で用いられる。<br />
<br />
そしてさらに、この操作はアーベル群であれば何でも良い……とは言っても、アーベル群であることは結構厳しい制限となる。<br />
<br />
例を出すと、リスト中の任意の範囲ーiからjの成分がすべてある条件式に適合しているかどうか確かめたいとする。<br />
関数型の考え方だと、まずは<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">map</span>関数を用いて(C++でなら<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><a href="http://www.boost.org/doc/libs/1_35_0/libs/iterator/doc/transform_iterator.html">transform_iterator</a></span>を用いて)リストの成分を<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Bool</span>に変え、その後で<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>を利用することを思いつくかもしれない。しかし、残念ながら<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">(Bool, (&&))</span>はアーベル群でないので、そのまま愚直に適用することができない。<br />
<pre class="prettyprint">Prelude> scanl1 (&&) $ map (<3) [1,2,3,2,1]
[True,True,False,False,False]
-- 最後の2要素については条件を満たしているのに、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">False</span>で埋もれてしまって分からない!
</pre>この場合は、単純にアーベル群である<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">(Integer, (+))</span>を<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">scan</span>に適用してやればよい。<br />
<pre class="prettyprint">Prelude> scanl1 (+) $ map (\x->if x<3 then 1 else 0) [1,2,3,2,1]
[1,2,2,3,4]
-- きちんと第3要素のみ条件から外れていることが分かる
</pre>Pythonっぽく書くとたぶんこんな感じ:<br />
<pre class="prettyprint">>>> import operator
>>> f = lambda x: 1 if x < 3 else 0
>>> scan(operator.add, [f(x) for x in [1,2,3,2,1]])
</pre><br />
<span class="Apple-style-span" style="font-size: large;">4. より詳しく知りたい方は</span><br />
scanアルゴリズムのGPGPU実装に関する詳細については、現在のところ『<a href="http://developer.download.nvidia.com/compute/cuda/1_1/Website/projects/scan/doc/scan.pdf">Parallel Prefix Sum (Scan) with CUDA(Mark Harris, NVIDIA)</a>』が一番分かりやすい。より詳しく知りたい方はそちらをどうぞ。<br />
<br />
というより、基本NVIDIAの出しているテキストは分かりやすい。すごいことです。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-35588746440107507602011-03-31T19:17:00.004+09:002011-03-31T22:22:06.264+09:00Moving<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgQ1jjTApN_2lWD0uNlCGsL9JGTXhsDPOUSST5UwpmQ5x2zXBVF835BUeeuxZ1y_QV5G0Z749jDlRJZXHNmcKN9h2wyJY9mjgW_BgOEIH2XVN1UbprAgRBxOZNGJeoh0d4IKwhDnR33u9g/s1600/moving.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgQ1jjTApN_2lWD0uNlCGsL9JGTXhsDPOUSST5UwpmQ5x2zXBVF835BUeeuxZ1y_QV5G0Z749jDlRJZXHNmcKN9h2wyJY9mjgW_BgOEIH2XVN1UbprAgRBxOZNGJeoh0d4IKwhDnR33u9g/s1600/moving.jpg" style="background:none; border:0;" /></a></div>本当は4/1に書く予定でしたが、そういえば4/1はエイプリルフールでしたので今日書きます。<br />
<br />
えーと、本日、<a href="http://www.phys.tohoku.ac.jp/">東北大学理学部物理学科</a>を卒業し、明日から<a href="http://www.is.tohoku.ac.jp/">東北大学大学院情報科学研究科</a>に入学します。所属する研究室はコンピュータビジョンなど、主に視覚的な情報について研究している研究室ですので、たぶんこれから書く記事はそれ系統の内容が多くなると思います。このブログのタイトルに少し近づきました。<br />
<br />
元々このブログは映像を作り、その技術的な情報についてログを残せるよう立ち上げたブログでした。確かに現在は映像の「え」の字もないブログに仕上がっていますが、もちろん今でも映像は好きですし、好きなプロダクションさんの映像はいつも観に行くようにしています。そして全く役に立たなかったかと言われればそういうわけではなくて、個人的に仕事を頼まれることもありましたし、Webサイトや印刷物のデザインをする際にも役立ちましたし、ここで得た知見のお陰で様々な方とお話することもできました。<br />
<br />
物理学科に入ったのも元々は相対論や量子力学の深遠な世界について知りたいために入りました。正直に言って完全に理解できたかというと怪しいところでありますし、理学から工学という他分野に移るわけで、一見すると今までの知識が役立たなくなるように見えます。が、実はそうでもなく、統計力学(*1)や相対論で使っている手法をコンピュータビジョンに輸入している論文があったり、案外様々な場面で役立ったりしているわけで(*2)。<br />
<br />
全然違う領域に見えても絡んでくる領域というものは少なからずあるものです。ということで、映像と物理がお互いに絡み合っている、コンピュータビジョンの道に進むことを決めました。<br />
<br />
これからもどうぞよろしくお願いします。<br />
<br />
p.s. 家族は無事でした、よかったよかった<br />
<br />
*1 ここで宣伝ですが、田崎先生の『<a href="http://www.gakushuin.ac.jp/~881791/statbook/">統計力学I, II</a>』はお薦めです。この本はぜひ理学系だけでなく工学系の方々にも読んでもらいたい本です。値段も2冊併せて7,000円くらいと安い。新品のゲームソフト1本買うくらいならこれを買いましょう<br />
*2 もちろん線形代数や基本的な偏微分も数式や概念を学んだりする上で役立ちました。一番学んで良かったのはブラケット記法とヒルベルト空間、あれは素晴らしいrezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-84040456103069887112011-03-11T20:51:00.001+09:002011-03-11T20:52:49.984+09:00自分は無事です自分は出張で京都にいましたので、無事です。もし身内の方が見ていましたらご安心ください。<br />
ただ暫くは京都に滞在することになると思います。早く交通機関が復旧してほしいです。<br />
<br />
ただ家族がどうなっているのか全く分からない。怖いです、不安です。<br />
<br />
[文句]<br />
<span class="Apple-style-span" style="font-size: large;">iPhoneではSoftbankの災害用伝言ダイヤルに『登録できない』!!!</span><br />
<br />
<span class="Apple-style-span" style="font-size: large;"></span>明らかにsoftbank側の怠慢、糞だ。糞すぎるrezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-49784359189520510732011-02-16T02:22:00.001+09:002011-02-16T02:27:23.503+09:00OMakeのContributionsに載ったった<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj69arVdw_t9FKKqjldD_FTAxuu4uVXkLMSOwUSrwzxL7h8YcNxqg0bcGHe8Kf8t3nuRwl1HNoig8osg7gkqx2gBfqtR5UUIW3npZwFpPOa1tlKnHRBSKZzV0MQYBrsfguF6mOaPSx7i3w/s1600/screenshot.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj69arVdw_t9FKKqjldD_FTAxuu4uVXkLMSOwUSrwzxL7h8YcNxqg0bcGHe8Kf8t3nuRwl1HNoig8osg7gkqx2gBfqtR5UUIW3npZwFpPOa1tlKnHRBSKZzV0MQYBrsfguF6mOaPSx7i3w/s1600/screenshot.jpg" style="background: none; border: 0;" /></a></div><div style="text-align: center;"><span class="Apple-style-span" style="font-family: monospace; white-space: pre-wrap;"><a href="http://omake.metaprl.org/contribs.html">The OMake build system: User Contributions</a></span></div><br />
<div style="text-align: center;">やったね!</div><br />
余談<br />
翻訳をすることは勉強にもなって非常にいいですし皆さんにもぜひやって欲しいのですが、<br />
その分、なんというか、日本の中へと篭ってしまうことにも繋がっている気がします。<br />
<br />
英語に臆することなく積極的にアウトプットできるよう、自分自身の心構えをもう少しきちんとしなければならんですね。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-84344027196649943072011-02-15T03:54:00.012+09:002011-03-31T19:31:40.851+09:00boost.GILで組む、ジェネリックな画像処理のアルゴリズム<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAqiDd-wKGUuOBR8cN-KRH9Y6T5trpttOEb4pxXoAR9KOENEfSF5AScP64L9wYuvXdL235MbGrF-7xgi8v3shoGbK-pR4hAjC54qUVAi-q-tg82z5xGOBB8_LuTOoK6aptStwILUrXsnI/s1600/human.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAqiDd-wKGUuOBR8cN-KRH9Y6T5trpttOEb4pxXoAR9KOENEfSF5AScP64L9wYuvXdL235MbGrF-7xgi8v3shoGbK-pR4hAjC54qUVAi-q-tg82z5xGOBB8_LuTOoK6aptStwILUrXsnI/s1600/human.jpg" style="background: none; border: 0;" /></a></div>"<a href="http://www.flickr.com/photos/gelinh/3279386782/">Image dans le néant</a>" by <a href="http://www.flickr.com/photos/gelinh/">gelinh</a><br />
<br />
<a href="http://opensource.adobe.com/wiki/display/gil/Generic+Image+Library">boost.GIL</a>は凄い!開発者の頭の良さがビシビシと伝わってくる!<br />
ということで、今回はGILに関しての紹介記事を書こうと思います。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">概要</span><br />
あなたは画像処理のエキスパート。顧客の依頼で、8bitのRGB画像を処理するアルゴリズムを記述していたとします。ところが対象となるデバイスの仕様を調べていた際に、実はRGBA画像にも対応させなければいけないことが分かりました。面倒だと思いながら書いていた矢先、さらにBGRやABGR,さらには16や24bitにも対応したアルゴリズムを記述しなければならないことが判明しました。なんということでしょう…これらの画像すべてに対してアルゴリズムを書くなんて、とてもじゃないですがやってられない。<i>やめてくれ!</i>って感じです。<br />
<br />
<a href="http://www.boost.org/">boost</a>に付いてくる<a href="http://opensource.adobe.com/wiki/display/gil/Generic+Image+Library">GIL</a>を用いることで、画像に対する操作をよりジェネリックに行うことが出来ます。具体的に言うと、あなたはGILの作法に従ってアルゴリズムを書き下すだけで、メモリに展開されている画像データの色空間、色深度、レイアウト、配置順に『<b>関わらず</b>』、任意のデータ配置に対して適用できるアルゴリズムを作ることができるのです。最初に提示した問題が一気に解決しちゃいました。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">長所と短所</span><br />
長所<br />
<br />
<ul><li>画像処理に対してジェネリックなアルゴリズムを記述でき、少ない記述量で全てをカバーできる。処理速度も普通に書き下した場合と同等(静的なポリモーフィズムしか使っていないので、コンパイル時に必要な計算はすべて終わっている)。</li>
<li>様々なアダプタが存在している(e.g. <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">subimage_view, rotated180_view</span>)ので、ばんばん遅延評価して一気に処理を行うことが可能。</li>
<li>普通の人が普通に書くよりも、GILのアルゴリズムに従って書いたほうが大抵の場合早い</li>
<li>ソースコードが比較的読みやすい。コード量も短くなる。</li>
<li>誰が作ったかわかんないもの利用したくねーという方も安心のAdobe製。画像処理の専門家が作ってるのでとても合理的。</li>
</ul><br />
短所<br />
<br />
<ul><li>テンプレートを多用しており『<b>非常に難解</b>』。</li>
<li>全てを抽象化しているために『<b>非常に難解</b>』。</li>
<li>C++と画像処理の両方に秀でていなければ記述できないので、使っている人が殆どいない。</li>
<li>英語のサイトですら文献が殆どない</li>
</ul>正直なところ、短所の部分がかなり厳しい。かなりC++をやりこんでいないと使えないということは、それだけでユーザ層を大分狭めます。仕方がないことですけど。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">Viewの概要</span><br />
GILでは、画像に対する操作はすべて『<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">View</span>』を起点にして行われます。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">View</span>は以下のような構造を持ち、ジェネリックな操作を行えるよう工夫してあります。<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><br />
</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">View</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">|-- Locator</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">| `-- Pixel</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">| |-- Channel</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">| `-- Layout</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">| |-- Color Space</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">| `-- Channel Mapping</span><br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">`-- dimensions</span><br />
<br />
<div style="text-align: center;">(<span class="Apple-style-span" style="font-size: x-small;">上図はMemory basedで、画素のメモリ配置がInterleavedな場合のView構造)</span></div><br />
以下のリストで、それぞれの名称がどのような役割を持っているのか説明します。<br />
<ul><li><b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">View</span></b>: (x, y)座標における画素値の他、特定のX、Y方向におけるイテレータや全ピクセルにわたって捜査できるxyイテレータなど、「様々な手法によって画素にアクセスできる手段」を提供する、一種の抽象的なクラスを表しています。画像の幅、高さを表す<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">dimensions</span>と、後述する<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Locator</span>を持ちます。(注: 画像のデータを実際に保持している『わけではありません』。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">View</span>はあくまで、対象の画素値へとアクセスするための『手段』を提供しているだけに過ぎないのです)</li>
<li><b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Locator</span></b>: 画像の2次元的な『位置』を表します。ランダムアクセスイテレータの2次元版のようなものと考えれば理解が早いと思います。具体的に言いますと、 <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">*loc </span>あるいは <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">*(loc + point2(dx, dy)) </span>などで現在の画素値や(dx, dy)だけ離れた位置の画素値を取得することは出来ますが、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Locator</span>の現在地(xy座標)や画像の大きさ(width, height)を取得することはできません。</li>
<li><b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Pixel</span></b>: <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Locator</span>によって返される画素値を表します。後述する<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Channel</span>と<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Layout</span>を持ちます。</li>
<li><b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Channel</span></b>: 画素の色深度を表します(8bit, 16bit, etc…)。各型の取りうる最大値、最小値も取得できます。</li>
<li><b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Layout</span></b>: 色空間のレイアウトを表します。ちょっと分かりにくいので、下の2つを読んだほうが早いでしょう。</li>
<ul><li><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><b>Color Space</b></span><b>:</b> 色空間を表します(e.g. RGB, CMYK, HSV, Gray, etc…)。</li>
<li><b><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Channel Mapping</span></b>: 実メモリ上での配置順を表します(e.g. RGB, BGR, etc…)。ただし、この情報はInterleave画像のみ用います(Planar画像だと配置順は関係ないため)。</li>
</ul></ul><br />
<span class="Apple-style-span" style="font-size: large;">アルゴリズムを用いて書いてみる</span><br />
それでは実際にGILに付属しているアルゴリズムを用いて実際に書き下してみます。以下の2つのアルゴリズムを用いてみましょう。<br />
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">boost::gil::transform_pixel(const View& src, const View& dst, F functor)</span><br />
対象のピクセルを受け取って別のピクセルを返す関数を、全ピクセルに対して適用します。この関数は、画素値と出力値が1対1対応である場合に用います(e.g. 明度、コントラスト, レベル補正, トーンカーブ, etc...)。<br />
<br />
<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">boost::gil::transform_pixel_positions(const View& src, const View& dst, F functor)</span><br />
周囲のピクセルの値を元に結果となるピクセルを返す関数を、全ピクセルに対して適用します。transform_pixelとは、受け取る引数がロケータであるという点が異なっており、周囲の画素値を元に実際の出力値を求める場合に有効です(e.g. Blur, Sovel, Laplacian, etc...)。<br />
<br />
実際にやってみれば、どんな感じか分かるでしょう:<br />
<br />
<pre class="prettyprint">#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;
}
</pre><br />
また、画像の出力例を下記に示します:<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAc067c6JHMX3xRF3Z_KdCdN8MAH6atjBQUQL1M99oKfueiP2POXiyUFTBPjVrrOiih-CZ3-BGbbjF6sngkx7OJGiBbn6_W2cHdk1vI_CKSt_06Db1F3EL-nSbQP4I6mDcvo1MnTX5OZY/s1600/out.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="100" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAc067c6JHMX3xRF3Z_KdCdN8MAH6atjBQUQL1M99oKfueiP2POXiyUFTBPjVrrOiih-CZ3-BGbbjF6sngkx7OJGiBbn6_W2cHdk1vI_CKSt_06Db1F3EL-nSbQP4I6mDcvo1MnTX5OZY/s400/out.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">左から、元画像, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">change_brightness, change_brightness_locator, x_gradient</span>の結果</td></tr>
</tbody></table><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">change_brightness_locator()</span>については、<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">for_each_pixel</span>用のアルゴリズムをlocatorを用いて適用させた場合にどうなるか確認するためだけに実装しましたので、実務的ではありません。また、値を型の持つ最大、最小値に丸めるために<a href="http://opencv.jp/"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">OpenCV</span></a>の<a href="http://opencv.jp/opencv-2svn/cpp/operations_on_arrays.html?#saturate_cast"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">cv::saturate_cast</span></a>を使用しています。<br />
<br />
main()についてはどんなことをやってるのか大体感覚的に理解できるはずですので、とりあえず3つのファンクタに関して説明を加えていきます。<br />
<br />
<pre class="prettyprint">typedef typename channel_type<Pixel>::type channel_type;</pre><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">channel_type<T>::type</span>で実際に使われているチャネルの型を取得できます。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">T</span>にはピクセル、ロケータ、ビューのどれを指定しても構いません。今回の場合ですと8bitのチャネルを指定しているので、型には<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">unsigned char</span>が入ります。<br />
<br />
<pre class="prettyprint">num_channels<Pixel>::value</pre>ピクセルを構成しているチャネルの数を取得できます。この場合ですと色空間がRGBなので値には3が入ります。この値はコンパイル時に決定されるので、最適化を用いることでループ展開が行われ、速度的なオーバーヘッドはありません。<br />
<br />
<pre class="prettyprint">typedef typename Locator::value_type pixel_type;</pre>ロケータが返すピクセルの型を取得できます。この場合ですと<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">rgb8_pixel_t</span>が入ります。<br />
<br />
<pre class="prettyprint">pixel_type src = loc(0, 0);</pre><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">()</span>演算子によって現在のピクセル値を取得しています。他にも<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">*</span>演算子でイテレータっぽくアクセスすることもできます。<br />
<br />
<pre class="prettyprint">explicit x_gradient_functor(Locator loc)
: left_(loc.cache_location(-1, 0)),
right_(loc.cache_location(1, 0)) {};</pre>ロケータの変動値をキャッシュしておくことで、高速化を図ります。<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Locator::cached_location_t</span>によってロケータのキャッシュ型を取得します。キャッシュを使ったアクセスには[]演算子を用います。<br />
仰々しく言っていますが、要は<span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">std::ptrdiff_t</span>みたいなものです。移動する値は同じなんだから記憶しておきましょうということです。<br />
<br />
<pre class="prettyprint">subimage_view(src, 1, 0, src.width()-2, src.height())</pre><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">subimage_view()</span>によって画像の領域を切り出しています。これはアダプタで、実際の評価は遅延的に行われます。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">まとめ</span><br />
大体こんな感じでアルゴリズムを組んでいくのがGIL流です。遅延評価とアルゴリズムを武器に、ジェネリックなアルゴリズムを作っていきましょう。<br />
<br />
<span class="Apple-style-span" style="font-size: large;">余談</span><br />
C++に関しての記事を投稿するのは正直に言って『少々怖い』です。大筋では合っていると思いますが、完全に合っているという自信があまりないからです。ですが、GILに関しては本当に、本当に参考となる文献が少ない。サンプルコードも多くない。こういった中でこの記事が少しでも参考程度になっていただければ幸いです。<br />
<br />
それと、ざっくり組んでしまったのでこうしたほうがいいよとか、間違い等あれば遠慮なく言ってもらえると助かります。<br />
<br />
(2/16) 追記: ソースコード若干の修正、説明や画像例の追加。<br />
(3/31) typo: e.q. -> e.g.rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-8624774289000338422010-12-18T18:25:00.010+09:002011-03-23T16:37:08.792+09:00画像をぼかす: いろんなボケ、リアルなボケ<div class="separator" style="clear: both; text-align: left;">『<i>画像をぼかす</i>』と一口に言っても、その方法にはいろんなものがあります。一番有名なのはPhotoshopの『ガウスぼかし』に見られるような<a href="http://en.wikipedia.org/wiki/Gaussian_blur">Gaussian Blur</a>ですが、残念ながらこれはカメラに見られるようなリアルなぼかしと比べるとどうしても見劣りします。原因はいくつか考えられて、一つはGaussian Blurに使われているガウス関数の勾配が穏やかなため、現実のレンズのシチュエーションを考慮していないという点があります。</div><div class="separator" style="clear: both; text-align: left;"><br />
</div><div class="separator" style="clear: both; text-align: left;">それじゃあ勾配を急にした関数ー<i>例えば<a href="http://en.wikipedia.org/wiki/Fermi%E2%80%93Dirac_statistics">フェルミ分布関数</a>のような</i>ーをカーネルに採用すれば良いのかというとそういうわけではありません。有名な<a href="http://ja.wikipedia.org/wiki/%E3%83%AC%E3%83%8A_%28%E7%94%BB%E5%83%8F%E3%83%87%E3%83%BC%E3%82%BF%29">Lena画像</a>で試してみましょう。</div><div class="separator" style="clear: both; text-align: left;"><br />
</div><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgBpokEL5zj4_XiyYaISfoENUeMHlnB5FoZi4QMqcD0abGTLUVzIV9i33ypxB1Bq-kDfcgD03Z44HRdwBn6AbTH-gmNUvrxKwBCMc_ktC9eLMzzyMo_DaXQEMh6lEMQU4-ZgzmJheVLcKQ/s1600/lena_normal_thumb.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="132" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgBpokEL5zj4_XiyYaISfoENUeMHlnB5FoZi4QMqcD0abGTLUVzIV9i33ypxB1Bq-kDfcgD03Z44HRdwBn6AbTH-gmNUvrxKwBCMc_ktC9eLMzzyMo_DaXQEMh6lEMQU4-ZgzmJheVLcKQ/s400/lena_normal_thumb.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">左が元画像、中央がカーネルにガウス関数を用いた画像、右がフェルミ分布関数を用いた画像。<br />
確かに中央よりは綺麗だけど、それでもまだリアルとは言えない</td></tr>
</tbody></table>それじゃあ他に何が違うのでしょう?端的に言うと、<b>黒も白も等価に扱って平均してしまうのがマズい</b>。<br />
<br />
画像をぼかすという行為は要は畳込み積分で、周囲のピクセル値を特定のウェイトで平均してやってることになります。で、普通にブラーしてしまうとピクセルの輝度に関わらず常に一定の割合で平均してしまいますので、なんか濁ったような画像が生成されてしまうと。<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj9ziZmSrlrJ-Y9YLWqu4lJ44jhMLkP3JjQP4nDtXKNlAY4AxOxnkEawTBdnasQOa95nr6atCD4ZLW2Yz700ynEnMJkvPir2tiWlIAWlPHc582osKQ5oJGnM_OTuagO2nY9fXK7IhdMzFE/s1600/lena_normalblur.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj9ziZmSrlrJ-Y9YLWqu4lJ44jhMLkP3JjQP4nDtXKNlAY4AxOxnkEawTBdnasQOa95nr6atCD4ZLW2Yz700ynEnMJkvPir2tiWlIAWlPHc582osKQ5oJGnM_OTuagO2nY9fXK7IhdMzFE/s400/lena_normalblur.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">通常のブラー。カーネルにはよりリアルに見せるため正六角形を用いている</td></tr>
</tbody></table>それではどのようなアプローチを用いればいいのでしょうか?レンズブラーを観察してみると、明るい箇所が強調されてぼかされていることが分かります。言い換えると、輝度の高いピクセルが優先的に混合されている、つまり予めRGBの全ピクセル値を適当な単調増加関数で写像してから畳込み積分を行って、その後対応する逆関数で戻してやればいい。要は<a href="http://www.hirax.net/dekirukana5/bokeboke1/">hirax.net</a>とやっていることは同じで、expしてからボカしてlogすると、こういうわけです。<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj0pd-ZCl9eHDsSgMakrHbvlS-0spMtOUqtMeHA4T8TAqQbtpoYb7vNs6Lfz2Gme6BfZ7JCaUO3wTuTlCaUDR1t5rC2NR9HMVeMKFFZV2aglnXv4Yh6MUfpAMQR8KXdUfmyDBNXg_hiYno/s1600/lena_realblur.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj0pd-ZCl9eHDsSgMakrHbvlS-0spMtOUqtMeHA4T8TAqQbtpoYb7vNs6Lfz2Gme6BfZ7JCaUO3wTuTlCaUDR1t5rC2NR9HMVeMKFFZV2aglnXv4Yh6MUfpAMQR8KXdUfmyDBNXg_hiYno/s400/lena_realblur.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">結果の画像。比較的リアルにボカされていることが分かる</td></tr>
</tbody></table><div class="separator" style="clear: both; text-align: center;"><br />
</div><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgx5Nk1wHWuo1_udqdtbWweq29wwnJjG87kEtHV0Ho7n1kOxt7Jq9HM9goglBjKL2YzgwZ3Sgx7l5MzoH63bMSpAYJ6X_5yKkTO-vuWfqyzZgWYw7BVHE6KDRQ0nBDJvbsZcb5yOpkoGn8/s1600/lena_thumb.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="131" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgx5Nk1wHWuo1_udqdtbWweq29wwnJjG87kEtHV0Ho7n1kOxt7Jq9HM9goglBjKL2YzgwZ3Sgx7l5MzoH63bMSpAYJ6X_5yKkTO-vuWfqyzZgWYw7BVHE6KDRQ0nBDJvbsZcb5yOpkoGn8/s400/lena_thumb.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">並べてみれば、違いがはっきり(クリックで拡大)</td></tr>
</tbody></table>Lenaさんだけでは少し分かりづらいので、適当な背景写真に適用してみましょう。<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjNF6HvHSmdq1anSgI_bmWEHvgWYSmraUZMCKttD_jNucnd1vai8dAARJ6Xq44GyaIh7gkRymsEjrPL7Ty1AEiujQIn2vRMh64-L1s22Lv0C1cP5AQyZZ_8-l6aU2X2iq_bmf37p0ByCBs/s1600/forest.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjNF6HvHSmdq1anSgI_bmWEHvgWYSmraUZMCKttD_jNucnd1vai8dAARJ6Xq44GyaIh7gkRymsEjrPL7Ty1AEiujQIn2vRMh64-L1s22Lv0C1cP5AQyZZ_8-l6aU2X2iq_bmf37p0ByCBs/s400/forest.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">森林画像</td></tr>
</tbody></table>まず、通常の手法でぼかしてしまうとどうなるでしょう?<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiUlCewYYKmzXKYAHIEuWVB5wzQjZT8qiDSTTYbryWQqPyYtVwBtDoflDG_24-U0w4C3GLTYWLr8dozEROa9gCd3T52u_tNgdB5sGh8O0C0D6_qBI78oJZBA-VJJ3F7XsfL2j7-PJPvaG8/s1600/forest_normalblur.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiUlCewYYKmzXKYAHIEuWVB5wzQjZT8qiDSTTYbryWQqPyYtVwBtDoflDG_24-U0w4C3GLTYWLr8dozEROa9gCd3T52u_tNgdB5sGh8O0C0D6_qBI78oJZBA-VJJ3F7XsfL2j7-PJPvaG8/s400/forest_normalblur.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">森林画像(通常のブラー)</td></tr>
</tbody></table>お世辞にも綺麗なボケとは言えません。それじゃあ今回の手法でぼかしてみます。<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEggqpwpthash1sq6FcvoFjoGVhwNKBTWoIe5slFUncsYNtgDRsrk4wcWO2onWRe9-wLTit14dOUiOPwuXyajscq75C-y5d9BLOy0Q18jgJG-pqHxn__TyZ-PgIES3BrcFLyZklfslHkSdQ/s1600/forest_realblur.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEggqpwpthash1sq6FcvoFjoGVhwNKBTWoIe5slFUncsYNtgDRsrk4wcWO2onWRe9-wLTit14dOUiOPwuXyajscq75C-y5d9BLOy0Q18jgJG-pqHxn__TyZ-PgIES3BrcFLyZklfslHkSdQ/s400/forest_realblur.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">森林画像(リアルブラー)</td></tr>
</tbody></table><b>なんということでしょう…</b>リアルなボケ画像ができました。<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjy88RfCUGHQ9ybZz5XB5u3teqAfZ56I4geD9zCjMc2djq5w6veDQAJ0xx2xy44BkPcmmqkS46QXZ1_8dURrAeAo_v5EKUaPWVmWCNfuDCkGzsVwrNkOYmfeowTcXLrVATLUV-nEtvk_Pw/s1600/forest_thumb.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="88" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjy88RfCUGHQ9ybZz5XB5u3teqAfZ56I4geD9zCjMc2djq5w6veDQAJ0xx2xy44BkPcmmqkS46QXZ1_8dURrAeAo_v5EKUaPWVmWCNfuDCkGzsVwrNkOYmfeowTcXLrVATLUV-nEtvk_Pw/s400/forest_thumb.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">比較画像(クリックで拡大)</td></tr>
</tbody></table>今回組んだOpenCVのソースコードを下記に示します。2.2でコンパイルしましたが、多分2.0以上だったら普通に通ると思います。<br />
<br />
<pre class="prettyprint">#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;
}
</pre><br />
<b>余談</b><br />
<ul><li>実質80行くらいで完成したので嬉しいです。可読性も多分結構読みやすい部類だと思う</li>
<li>C++のコードが何やってるか分からない方は、まず関数オブジェクトとSTLについて調べてみましょう。あれはいいもんです</li>
</ul><b>追記(11/03/23)</b><br />
<blockquote>どんなデバイスであるか、たとえば、それが入力デバイスであるのか、出力でバイスであるのか、はたまた、計算処理デバイスであるのか次第で「輝度値の単位系」をどのように保持(処理)するべきかは違うことでしょう。<br />
大切なことは、取り扱う「値」がどんな単位であるのかを見失わないことなのか、と思います。 </blockquote><blockquote><span class="Apple-style-span" style="font-family: monospace; white-space: pre-wrap;"><a href="http://hirax.net/techlog/2011/03/06.html">hirax.net::「コンピュータ画像処理」と「輝度値の単位系」</a></span> </blockquote>自分がこのような『任意の単調増加関数で良い』という具合で筆を進めたのには2つ理由があります。<br />
<br />
まず1つ目が、『<b>全てのデジタルカメラにおいて、ピクセルは輝度の対数の次元を取っているのだろうか?</b>』という単純な疑問です。<br />
<br />
きちんと物理的/工学的な意味合いを持つ画像処理にするためには、格納されているピクセルデータの次元がきちんと揃えられている必要があります。自分自身も普通に考えて輝度の対数が格納されているんだろうということで、任意ではなく指数関数に制限された関数系を使いましょうと運びたかったのですが、そうするための証拠が不足していた。<br />
<br />
もしかしたら自分の予想もつかない方法で処理されている可能性があるわけで、そういった証拠が十分にないことから、単調増加関数と少し曖昧とした意味合いを含ませたのです。<br />
<br />
もう一つの理由は『<b>今回の画像処理は物理的に正確である必要がない</b>』という理由です。これが大きかった。<br />
<br />
そもそも今回の画像処理はどういう目的があるのでしょうか?ぼけて普通よりもリアルに見えて面白い、それで終わりです。物理的に計測をするわけではなく、ざっくり『それっぽく』見えればそれでいいという類のものです。そのような目的の下では、下手に物理的正確さを求めて処理時間を長くするよりは、計算速度の早いフェイクのほうがより適しています。<br />
<br />
だから何も指数関数に制限する理由はあまりないわけで、より早くてよりそれっぽく見えれば、そちらのほうが優れていると、そう言えるわけです(ただ、今回の計算量からすると最も大きなオーダーとなっているのは畳み込み積分の箇所で、単調増加関数自体のコストはそれに比べれば微々たるものでしょう)。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-30128409231026839922010-10-27T02:03:00.002+09:002010-10-27T02:06:01.543+09:00OMakeマニュアル 日本語訳 PDFバージョンをリリースしました<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEijqs3shyphenhyphen73Z2rtYvdpQm8Y4wLuKzDuJerwt6JNegecElkrYFbqnKAHoli1zJq5vbJS6js25gxsO4IM3ZCHNhWlJV0s6WtyqR0NN1sFHfARYUW9uzR0Lf1Z__6V0U6L8FzZNowb8tgQMB8/s1600/cute.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEijqs3shyphenhyphen73Z2rtYvdpQm8Y4wLuKzDuJerwt6JNegecElkrYFbqnKAHoli1zJq5vbJS6js25gxsO4IM3ZCHNhWlJV0s6WtyqR0NN1sFHfARYUW9uzR0Lf1Z__6V0U6L8FzZNowb8tgQMB8/s1600/cute.jpg" style="background: none; border: 0;" /></a></div>"<a href="http://www.flickr.com/photos/clickflashphotos/3440590766/">Aint they cute?</a>" by <a href="http://www.flickr.com/photos/clickflashphotos/">Nicki Varkevisser</a><br />
<br />
題名の通りで、以前から取り組んでおりました『<a href="http://omake-japanese.sourceforge.jp/">OMakeマニュアル 日本語訳</a>』の全ページを一つに纏めたPDFバージョンをリリースしました。<br />
<br />
<a href="http://omake-japanese.sourceforge.jp/OMake.pdf">http://omake-japanese.sourceforge.jp/OMake.pdf</a><br />
<br />
よりダウンロードできます。電子書籍でまとめてご覧になりたい方や、職場などの環境で使いたいけれどネットが使えないのでごりごり印刷したい…といった場合に便利です。<br />
<br />
元々は<a href="http://www.shibu.jp/">渋川さん</a>の『<a href="http://sphinx-users.jp/cookbook/pdf/latex.html">LaTeX経由でのPDF作成</a>』という記事を拝見しまして、試験的にPDF版を作成してみたところ、思いのほか綺麗に製本されていることからきちんとリリースしてみようと思い立った次第です。殆ど素の状態から手を加えていないので暫定的ですが、後はきちんとスタイルを組んでいって、もうちょっと読みやすいレイアウトにしていけばいいかなって感じです。<br />
<br />
思えば個人的な暇潰しから始まったこのプロジェクトも現在はPDFにして約200ページにもなり、きちんと一つの事柄に集中して取り組めたことは、学部時代のいい経験になったと考えています。まあ、これからも暇があればちょくちょくレイアウト、文章含めて修正していきますので、どうぞ今後ともよろしくお願いします。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-74060543870891978942010-10-04T01:58:00.004+09:002010-10-04T16:28:17.124+09:00確率微分方程式を解くー幾何ブラウン運動、バシチェックモデルの場合<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEidekp8c2pod6sJqb2g8qeTb_KJMzcWntD0d8hlOyYKdX4OlsWBlyEqpal9A768Xx2jVr6i0_fW79vEhE2CJutpYdbaqgNmqE0_byOkcHdSaDTmh3yKrvI-ZFRPK1tznrpE-sd6Td2Zpec/s1600/New+York+Fashion+Week+Fall+2007.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEidekp8c2pod6sJqb2g8qeTb_KJMzcWntD0d8hlOyYKdX4OlsWBlyEqpal9A768Xx2jVr6i0_fW79vEhE2CJutpYdbaqgNmqE0_byOkcHdSaDTmh3yKrvI-ZFRPK1tznrpE-sd6Td2Zpec/s1600/New+York+Fashion+Week+Fall+2007.jpg" style="background: none; border: 0;" /></a></div>"<a href="http://www.flickr.com/photos/artcomments/382733093/">New York Fashion Week Fall 2007: Doo Ri</a>" by <a href="http://www.flickr.com/photos/artcomments/">Art Comments</a><br />
<br />
確率微分方程式という分野がある。どういうものかっていうと、要は確率過程を微分方程式によって表そうという試みであり、微分方程式の中に微小の確率変数が含まれている。具体的にはこんな感じ。<br />
<a href="http://formula.s21g.com/?%7B%5Crm%20d%7Dx(t)%20%3D%20a(x,%20t)%7B%5Crm%20d%7Dt%20+%20b(x,%20t)%7B%5Crm%20d%7DW_t"><img src="http://formula.s21g.com/?%7B%5Crm%20d%7Dx(t)%20%3D%20a(x,%20t)%7B%5Crm%20d%7Dt%20+%20b(x,%20t)%7B%5Crm%20d%7DW_t.png" /></a><br />
このdWはブラウン運動を表していて、時間がたつにつれてうねうねと確率的に動く。期待値と標準偏差はそれぞれ0と√tで、時間がたつにつれて動く場所もどんどん広がっていく。この運動はランダムウオークとか酔歩とも呼ばれていて、酔っ払ったおっさんがうねうねとあちこちに歩いているように見えることからこういう名前がついたそうな。上のような式で表すことができる確率過程を一般には『伊藤過程』と言って、他にもストラトノビッチ確率解析とかいうのもあるらしいけど、詳しくやっていないので分かりません。<br />
<br />
これの何が便利かっていうと確率過程を通常の微分方程式のように扱えるところで、確率変数をまるで微分方程式を解くかのように求めることができる。とは言っても通常の微分方程式にはない規則がいくつかあって、初めはすごく奇妙に思えるかもしれない。<br />
<br />
具体的には、以下のように確率変数x(t)とy(t)を定めた場合、<br />
<br />
<a href="http://formula.s21g.com/?%7B%5Crm%20d%7Dx(t)%20%26%3D%26%20a(x,%20t)%7B%5Crm%20d%7Dt%20+%20b(x,%20t)%7B%5Crm%20d%7DW_t%5C%5C%0Ay(t)%20%26%3D%26%20g(x,%20t)"><img src="http://formula.s21g.com/?%7B%5Crm%20d%7Dx(t)%20%26%3D%26%20a(x,%20t)%7B%5Crm%20d%7Dt%20+%20b(x,%20t)%7B%5Crm%20d%7DW_t%5C%5C%0Ay(t)%20%26%3D%26%20g(x,%20t).png" /></a><br />
このy(t)もまた伊藤過程に従い、確率微分方程式dyは以下のルールに従う。<br />
<br />
<a href="http://formula.s21g.com/?%7B%5Crm%20d%7Dy(t)%20%26%3D%26%20%5Cfrac%7B%5Cpartial%20g%7D%7B%5Cpartial%20t%7D%7B%5Crm%20d%7Dt%20+%20%5Cfrac%7B%5Cpartial%20g%7D%7B%5Cpartial%20x%7D%7B%5Crm%20d%7Dx%20+%20%5Cfrac%7B1%7D%7B2%7D%5Cfrac%7B%5Cpartial%5E2%20g%7D%7B%5Cpartial%20x%5E2%7D(%7B%5Crm%20d%7Dx)%5E2"><img src="http://formula.s21g.com/?%7B%5Crm%20d%7Dy(t)%20%26%3D%26%20%5Cfrac%7B%5Cpartial%20g%7D%7B%5Cpartial%20t%7D%7B%5Crm%20d%7Dt%20+%20%5Cfrac%7B%5Cpartial%20g%7D%7B%5Cpartial%20x%7D%7B%5Crm%20d%7Dx%20+%20%5Cfrac%7B1%7D%7B2%7D%5Cfrac%7B%5Cpartial%5E2%20g%7D%7B%5Cpartial%20x%5E2%7D(%7B%5Crm%20d%7Dx)%5E2.png" /></a><br />
普通だったら消えるはずの2乗の項が残ってしまう。それじゃあ(dx)^2はどのようにして計算するのかというと、以下の公式に従う。<br />
<br />
<a href="http://formula.s21g.com/?(%7B%5Crm%20d%7Dt)%5E2%20%26%3D%26%20%7B%5Crm%20d%7Dt%7B%5Crm%20d%7DW_t%20%3D%20%7B%5Crm%20d%7DW_t%20%7B%5Crm%20d%7Dt%20%3D%200%20%5C%5C%0A(%7B%5Crm%20d%7DW_t)%5E2%20%26%3D%26%20%7B%5Crm%20d%7Dt"><img src="http://formula.s21g.com/?(%7B%5Crm%20d%7Dt)%5E2%20%26%3D%26%20%7B%5Crm%20d%7Dt%7B%5Crm%20d%7DW_t%20%3D%20%7B%5Crm%20d%7DW_t%20%7B%5Crm%20d%7Dt%20%3D%200%20%5C%5C%0A(%7B%5Crm%20d%7DW_t)%5E2%20%26%3D%26%20%7B%5Crm%20d%7Dt.png" /></a><br />
なんと、(dW)^2がdtに変化してしまう!これを『伊藤の公式』という。なんでそうなるのか?それを知りたい方は専門書を漁ればいくらでも載っているので、そこらへんを参照すればいいと思う。ここが通常の微分方程式と違うところで、逆を言えばそこさえ守ればきちんと積分によって求めることができる。専門書だとこの箇所がすごく分かり辛い表現で説明されているけど、詳細をはしょればこんな感じなはず。<br />
<br />
それじゃあ、まず手始めに<a href="http://ja.wikipedia.org/wiki/%E7%A2%BA%E7%8E%87%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%E5%BC%8F">幾何ブラウン運動</a>と呼ばれるモデル<br />
<br />
<a href="http://formula.s21g.com/?%7B%5Crm%20d%7Dx(t)%20%3D%20%5Cmu%20x(t)%20%7B%5Crm%20d%7Dt%20+%20%5Csigma%20x(t)%20%7B%5Crm%20d%7DW_t"><img src="http://formula.s21g.com/?%7B%5Crm%20d%7Dx(t)%20%3D%20%5Cmu%20x(t)%20%7B%5Crm%20d%7Dt%20+%20%5Csigma%20x(t)%20%7B%5Crm%20d%7DW_t.png" /></a><br />
を求めてみる。これは株価の変動を表す最も簡単なモデルなんだけど、何故これが株価の変動を表しているのかっていうのは、とりあえず解いてみれば判明するはず。<br />
<br />
さて、普通の考えだったら『変数分離法を用いてxを分離させて、そのまま積分すればよくね?』って発想になるんだけど、何度も言っているようにこれは確率微分方程式なので、まずは伊藤の補題を用いて変数変換を行い、xを分離させてみることにする。y = ln(x)とおくと、伊藤の補題より<br />
<br />
<a href="http://formula.s21g.com/?%7B%5Crm%20d%7Dy%20%26%3D%26%20%5Cfrac%7B1%7D%7Bx%7D%7B%5Crm%20d%7Dx%20-%20%5Cfrac%7B1%7D%7B2x%5E2%7D%20(%7B%5Crm%20d%7Dx)%5E2%5C%5C%0A%20%26%3D%26%20%5Cfrac%7B1%7D%7Bx%7D(%5Cmu%20x%20%7B%5Crm%20d%7Dt%20+%20%5Csigma%20x%20%7B%5Crm%20d%7DW_t)%20-%20%5Cfrac%7B1%7D%7B2x%5E2%7D(%5Csigma%5E2%20x%5E2%20%7B%5Crm%20d%7Dt)%5C%5C%0A%20%26%3D%26%20(%5Cmu%20-%20%5Cfrac%7B1%7D%7B2%7D%5Csigma%5E2)%7B%5Crm%20d%7Dt%20+%20%5Csigma%20%7B%5Crm%20d%7DW_t"><img src="http://formula.s21g.com/?%7B%5Crm%20d%7Dy%20%26%3D%26%20%5Cfrac%7B1%7D%7Bx%7D%7B%5Crm%20d%7Dx%20-%20%5Cfrac%7B1%7D%7B2x%5E2%7D%20(%7B%5Crm%20d%7Dx)%5E2%5C%5C%0A%20%26%3D%26%20%5Cfrac%7B1%7D%7Bx%7D(%5Cmu%20x%20%7B%5Crm%20d%7Dt%20+%20%5Csigma%20x%20%7B%5Crm%20d%7DW_t)%20-%20%5Cfrac%7B1%7D%7B2x%5E2%7D(%5Csigma%5E2%20x%5E2%20%7B%5Crm%20d%7Dt)%5C%5C%0A%20%26%3D%26%20(%5Cmu%20-%20%5Cfrac%7B1%7D%7B2%7D%5Csigma%5E2)%7B%5Crm%20d%7Dt%20+%20%5Csigma%20%7B%5Crm%20d%7DW_t.png" /></a><br />
となる。(-1/2 σ^2)の項が出てくるところがミソで、普通に変数変換してしまうと出現しない。これで定数部分だけに分離できたので、そのまま愚直に積分を行うと、<br />
<br />
<a href="http://formula.s21g.com/?y(t)%20%26%3D%26%20y(0)%20+%20(%5Cmu%20-%20%5Cfrac%7B1%7D%7B2%7D%5Csigma%5E2)t%20+%20%5Csigma%20W(t)%20%5C%5C%0Ax(t)%20%26%3D%26%20x(0)%5Cexp%5Cleft(%20(%5Cmu%20-%20%5Cfrac%7B1%7D%7B2%7D%5Csigma%5E2)t%20+%20%5Csigma%20W(t)%20%5Cright)"><img src="http://formula.s21g.com/?y(t)%20%26%3D%26%20y(0)%20+%20(%5Cmu%20-%20%5Cfrac%7B1%7D%7B2%7D%5Csigma%5E2)t%20+%20%5Csigma%20W(t)%20%5C%5C%0Ax(t)%20%26%3D%26%20x(0)%5Cexp%5Cleft(%20(%5Cmu%20-%20%5Cfrac%7B1%7D%7B2%7D%5Csigma%5E2)t%20+%20%5Csigma%20W(t)%20%5Cright).png" /></a><br />
となる。これは対数正規分布の形となっているため、株価の変動を表す簡単なモデルとなっていることは容易に理解できると思う。<br />
<br />
それでは、次に<a href="http://ja.wikipedia.org/wiki/%E3%83%90%E3%82%B7%E3%83%81%E3%82%A7%E3%83%83%E3%82%AF%E3%83%BB%E3%83%A2%E3%83%87%E3%83%AB">Vasicekモデル</a>と呼ばれるモデル<br />
<br />
<a href="http://formula.s21g.com/?%7B%5Crm%20d%7Dr(t)%20%3D%20a(b-r(t))%7B%5Crm%20d%7Dt%20+%20%5Csigma%20%7B%5Crm%20d%7DW_t"><img src="http://formula.s21g.com/?%7B%5Crm%20d%7Dr(t)%20%3D%20a(b-r(t))%7B%5Crm%20d%7Dt%20+%20%5Csigma%20%7B%5Crm%20d%7DW_t.png" /></a><br />
について求めてみる。これは主に国債金利の利子率の時間的変動を表すモデルで、いったいこいつの何が重要なんじゃいっていうと、金融商品のリスクを推定する際に使われる。金融商品の価格は利子率によって変動するので、利子率の変動を推定することは金融商品のリスクを推定することに繋がる(*1)。つまり、今の自分のポートフォリオがどれくらいのリスクを持っているのか判断できる。そもそもこの記事を書くきっかけになったのはこいつが原因で、wikipediaさんにはただ『積分するとこれになりますよ』って書いてあるだけで、どうやって積分するのかについては書いていなかった。しょうがないので自分で手を動かして、ようやっと理解できたので公開してみようと。<br />
<br />
さて、まずはu=b-rとおいて変数変換すると、<br />
<br />
<a href="http://formula.s21g.com/?%7B%5Crm%20d%7Du%20%26%3D%26%20-%7B%5Crm%20d%7Dr%20%5C%5C%0A%26%3D%26%20-au%20%7B%5Crm%20d%7Dt%20-%20%5Csigma%7B%5Crm%20d%7DW_t"><img src="http://formula.s21g.com/?%7B%5Crm%20d%7Du%20%26%3D%26%20-%7B%5Crm%20d%7Dr%20%5C%5C%0A%26%3D%26%20-au%20%7B%5Crm%20d%7Dt%20-%20%5Csigma%7B%5Crm%20d%7DW_t.png" /></a><br />
となる。ただ普通にChain Ruleやってるように見えるけど、たまたま結果が同じだっただけで、どんな変数変換でもきちんと伊藤の補題を使わなきゃいけない。そんで、物理やってる人はピンとくるかもしれないけど、これは<a href="http://ja.wikipedia.org/wiki/%E3%83%A9%E3%83%B3%E3%82%B8%E3%83%A5%E3%83%90%E3%83%B3%E6%96%B9%E7%A8%8B%E5%BC%8F">ランジュバン方程式</a>の形となっているので、積分には同様の手法を利用できる。まぁ自分はピンとこなかったんですけど。それはさておき、v=e^(at)uとおいてさらに変数変換してやると、<br />
<br />
<a href="http://formula.s21g.com/?%7B%5Crm%20d%7Dv%20%26%3D%26%20ae%5E%7Bat%7Du%20%7B%5Crm%20d%7Dt%20+%20e%5E%7Bat%7D%7B%5Crm%20d%7Du%5C%5C%0A%26%3D%26%20ae%5E%7Bat%7Du%20%7B%5Crm%20d%7Dt%20+%20e%5E%7Bat%7D(-au%7B%5Crm%20d%7Dt%20-%20%5Csigma%7B%5Crm%20d%7DW_t)%5C%5C%0A%26%3D%26%20-%5Csigma%20e%5E%7Bat%7D%7B%5Crm%20d%7DW_t"><img src="http://formula.s21g.com/?%7B%5Crm%20d%7Dv%20%26%3D%26%20ae%5E%7Bat%7Du%20%7B%5Crm%20d%7Dt%20+%20e%5E%7Bat%7D%7B%5Crm%20d%7Du%5C%5C%0A%26%3D%26%20ae%5E%7Bat%7Du%20%7B%5Crm%20d%7Dt%20+%20e%5E%7Bat%7D(-au%7B%5Crm%20d%7Dt%20-%20%5Csigma%7B%5Crm%20d%7DW_t)%5C%5C%0A%26%3D%26%20-%5Csigma%20e%5E%7Bat%7D%7B%5Crm%20d%7DW_t.png" /></a><br />
となり、いい感じに余計な項が消えてブラウン運動だけになった。よってこれを積分してやると<br />
<br />
<a href="http://formula.s21g.com/?v(t)%20%26%3D%26%20v(0)%20-%20%5Csigma%20%5Cint_%7B0%7D%5E%7Bt%7De%5E%7Bas%7D%7B%5Crm%20d%7DW_s%5C%5C%0A%26%3D%26%20b%20-%20r(0)%20-%20%5Csigma%20%5Cint_%7B0%7D%5E%7Bt%7De%5E%7Bas%7D%7B%5Crm%20d%7DW_s%0A"><img src="http://formula.s21g.com/?v(t)%20%26%3D%26%20v(0)%20-%20%5Csigma%20%5Cint_%7B0%7D%5E%7Bt%7De%5E%7Bas%7D%7B%5Crm%20d%7DW_s%5C%5C%0A%26%3D%26%20b%20-%20r(0)%20-%20%5Csigma%20%5Cint_%7B0%7D%5E%7Bt%7De%5E%7Bas%7D%7B%5Crm%20d%7DW_s%0A.png" /></a><br />
あとはごちゃごちゃ整理してやると<br />
<br />
<a href="http://formula.s21g.com/?r(t)%20%3D%20r(0)e%5E%7B-at%7D%20+%20b(1-e%5E%7B-at%7D)%20+%20%5Csigma%20e%5E%7B-at%7D%5Cint_%7B0%7D%5E%7Bt%7De%5E%7Bas%7D%7B%5Crm%20d%7DW_s"><img src="http://formula.s21g.com/?r(t)%20%3D%20r(0)e%5E%7B-at%7D%20+%20b(1-e%5E%7B-at%7D)%20+%20%5Csigma%20e%5E%7B-at%7D%5Cint_%7B0%7D%5E%7Bt%7De%5E%7Bas%7D%7B%5Crm%20d%7DW_s.png" /></a><br />
これで<a href="http://ja.wikipedia.org/wiki/%E3%83%90%E3%82%B7%E3%83%81%E3%82%A7%E3%83%83%E3%82%AF%E3%83%BB%E3%83%A2%E3%83%87%E3%83%AB">wikipediaさんにあるのとおんなじ形</a>になった。めでたしめでたし。<br />
<br />
なんかブログに書くと導出するのは楽ちんなように思えるけど、結構悩んだし大変だった。どんなもんでも答えを見るとすぐ分かるし思いつきそうなんだけど、実際になにもない地点から思いつくかどうかっていうのはまた別な話で、それを考えると新しい学問を生み出すっていうのは凄まじいほどの労力が必要だ。実際、未だにどういう発想で量子力学を思いついたのか本気で理解できない。あんなのヒントがいくらあっても無理だ。そういうことを考えた一日でありましたとさ。<br />
<br />
*1: もちろんこの説明は適当でないが、詳しく説明してもあんまり意味がないのでこう書くことにする<br />
<br />
追記(2010/10/04)<br />
1/2の項が抜けていたので修正。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0tag:blogger.com,1999:blog-1062283300115822026.post-38315432638623925382010-09-28T18:31:00.002+09:002010-09-28T18:33:01.354+09:00マルコフ連鎖モンテカルロ法を実装してみよう<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjtCR3D_yv5pZd3AlbZpPI7cyhLUF4Pl8q4bf5-uAeNkPtWTqQ33mTEZqsG7wBlfy5o4uRmmjAbghtwAbUs1wXVApiCUPxZJ8SX9b4O-fmLd73mTsA9Qthz9JvZciFA0UaQLXdf7DLmQHY/s1600/chain.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjtCR3D_yv5pZd3AlbZpPI7cyhLUF4Pl8q4bf5-uAeNkPtWTqQ33mTEZqsG7wBlfy5o4uRmmjAbghtwAbUs1wXVApiCUPxZJ8SX9b4O-fmLd73mTsA9Qthz9JvZciFA0UaQLXdf7DLmQHY/s1600/chain.jpg" style="background:none; border:0" /></a></div>"<a href="http://www.flickr.com/photos/pedrosz/3642296149/">Chain Bridge, Budapest</a>" by <a href="http://www.flickr.com/photos/pedrosz/">szeke</a><br />
<br />
約2ヶ月ぶりということで、ご無沙汰しております。書きたいネタというのは結構あって書こう書こうとは思っているんですが、なにより書くとなるとあれもこれも伝えなきゃいけないみたいな思いになって、結局分量の多さから諦めてしまうというのが結構続いています。もう少し気を張らずに更新していきたいものです。<br />
<br />
さて、最近の自分はマルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo, MCMC)についておりました。何にも知らない方ににとってはよく分からないかもしれませんが、この手法はマルコフ連鎖が持つ簡明さとモンテカルロ法が持つ実用性が合わさった凄まじい手法なんです。そしてなによりとてもエレガント。<br />
<br />
とりあえず知らない人のためにてきとーに解説しますと、与えられた適当な関数から確率変数をサンプリングするための公式です。べつにこれじゃなくても確率変数のサンプリングにはいろんなアルゴリズムがあるわけですが、これを使うと関数の計算コストがやけに高い場合だとか、ピーク部分を集中的にサンプリングしてくれたりとかで結構便利で、金融や工学、物理の分野でも非常に役立っております。計算物理をやる人にとっては必須に近い手法といっても過言ではありません。<br />
<br />
そんなMCMCですが、実際に提案されている手法には様々なものがあります。不変分布から遷移確率を逆算するため、自由度が増えるから遷移確率は一意じゃなくなってしまうんです。今回は最も有名なMetropolis-Hastingsアルゴリズム(*1)を用いて、任意の関数に従って乱数のサンプリングを行うコードを練習がてらpythonで書いてみました。<br />
<br />
<pre class="prettyprint">#!/usr/bin/env python
#-*- coding:utf-8 -*-
import sys
import math
import random
def main():
f = lambda x: math.exp(-x*x / 10.0)
for x, fx in metropolis(f):
print x, fx
def 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 value
def _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_fx
if __name__ == '__main__':
sys.exit(main())
</pre>うーん・・・きれいだ(*2)・・・<br />
<br />
こいつのすごいところは、アルゴリズム自体は非常にシンプルで実装しやすいというところにあります。行数も上を見ると分かるのですがとっても少ないし、コア部分はたった25行くらいで構成されています。<br />
<br />
なのに使っている学問は結構複雑で、定常分布から遷移行列を導出するために詳細釣り合いの条件からMetropolis的制約を用いて導出、サンプリングを行っております。このアルゴリズム自体は結構知ってる人が多いかもしれませんが、実際に理論方面で理解している人はどうなんでしょう。<br />
<br />
んで、正規分布のガウス関数からサンプリングを行った結果がこれ。 <br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEirIA7cqJuOYVSYcBJsSO13kgN6xzwyOwoEaxyWD4HiCq9OIetxCTXjDp-pDzZVdStNPhbhcQUbZaWO5yU3wz3niHJIJd_VJZsfFpPoe5VRDGMaQCoQBB_JmSsP3rWDWVOw8x2FhhS3Lxg/s1600/test.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEirIA7cqJuOYVSYcBJsSO13kgN6xzwyOwoEaxyWD4HiCq9OIetxCTXjDp-pDzZVdStNPhbhcQUbZaWO5yU3wz3niHJIJd_VJZsfFpPoe5VRDGMaQCoQBB_JmSsP3rWDWVOw8x2FhhS3Lxg/s1600/test.png" /></a></div><br />
赤線が元の確率分布。きちんとサンプリングが行われているようです。今回は正規分布で計算を行いましたが、非負の関数であればなんでも構いません。<br />
<br />
あとはこれを研究テーマへと落とし込むことができれば、とりあえず一区切りつけるかなって感じです。<br />
<br />
*1: 提案確率関数が対称であるという制約を持つのがMetropolis。それを非対称の場合にも拡張したのがHastings。個人的にはMetropolis-HastingよりMetropolisのほうが綺麗で好きです。というか提案確率密度関数が非対称じゃないといけないようなシチュエーションがよく分からない<br />
*2: yieldとジェネレータって便利ですよね。計算速度も効率的だし、慣れてくるとすべて遅延評価で行わせたいくらいrezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com1tag:blogger.com,1999:blog-1062283300115822026.post-41258153325671387422010-07-23T01:49:00.009+09:002010-08-23T01:23:18.251+09:00ドルコスト平均法は賢い投資方法と言えるのか? 補足<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiEkPxHBRIEKIjtBl2cAUyR-KRPp0_qlIgen9yh8qUR80nLOZTBHUK_2Jd6O_08GXen3luVhshxWqqUZN3gvWsbR8G1VkIppIQWDrt_rU-dfoNWniDEVQfg4r13nBFhQDyrBHzo-QPLsKk/s1600/capitalism.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiEkPxHBRIEKIjtBl2cAUyR-KRPp0_qlIgen9yh8qUR80nLOZTBHUK_2Jd6O_08GXen3luVhshxWqqUZN3gvWsbR8G1VkIppIQWDrt_rU-dfoNWniDEVQfg4r13nBFhQDyrBHzo-QPLsKk/s1600/capitalism.jpg" style="background: none repeat scroll 0% 0% transparent; border: 0pt none;" /></a></div>"<a href="http://www.flickr.com/photos/wnorrix/262600004/">Pyramid of Capitalism</a>" by <a href="http://www.flickr.com/photos/wnorrix/">Warren Noronha</a><br />
<br />
<ul><li><a href="http://mglab.blogspot.com/2010/05/part1.html">ドルコスト平均法は賢い投資方法と言えるのか? part1</a></li>
<li><a href="http://mglab.blogspot.com/2010/05/part2.html">ドルコスト平均法ははたして賢い投資方法と言えるのか? part2</a></li>
</ul><br />
以前ブログの記事を読んでいただいた友人から感想を伺うと、「<i>なんかよく分からない</i>」といった感想が帰ってきました。基本的に説明が不足しているなぁと思った箇所は実際いくつもあって、金融の確率論を少しでもやると大体感覚的に理解できるのですが、普通の人はそんなのやらないのが当たり前ですので、基本的な概念について少し説明したいと思います。この分野に興味ない人でもちょこちょこと面白そうな話を挟んでいくので、暇つぶしに見てやってください。<br />
<br />
<h4>1. 確率変数、期待値、分散</h4><b>確率変数というのは試行によって出てくる値が確率的に異なってくる変数のこと</b>です。これは抽象的なので、もうすこし具体的にいきましょう。例えば、以下のコイントスゲームを考えます。最初あなたは1万円を所持していたとして、表であれば所持金が2倍となり、一方で裏であれば所持金が0.5倍、つまり半分になります。このような試行を3回行ったとして、最終的な資金額をXと表すと、このゲームは以下の式で表せます。<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?X%20%3D%2010000%20%5Ctimes%20X_1X_2X_3"><img src="http://formula.s21g.com/?X%20%3D%2010000%20%5Ctimes%20X_1X_2X_3.png" /></a></div><br />
ここでのX<sub>1</sub>はそれぞれの試行の確率変数に相当し、50%の確率で0.5, 50%の確率で2.0の値をとります。<br />
で、これだけではただ定式化しただけで、何にも面白いことは分かりません。ですが<b>こいつの期待値と分散を求めてやることで、このゲームは得なのか損なのか、またどれくらいのリスクがあるゲームなのか調べてやることができる</b>のです。<br />
<br />
例えば、Xの期待値をE[X]としてやりましょう。ここで、期待値の性質として以下を確認しておきましょう。まず、XとYが確率変数とすると、その和の期待値は単純にそれぞれの期待値の和で表せます。<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?%7B%5Crm%20E%7D[X%20+%20Y]%20%3D%20%7B%5Crm%20E%7D[X]%20+%20%7B%5Crm%20E%7D[Y]"><img src="http://formula.s21g.com/?%7B%5Crm%20E%7D[X%20+%20Y]%20%3D%20%7B%5Crm%20E%7D[X]%20+%20%7B%5Crm%20E%7D[Y].png" /></a></div><br />
また、XとYが独立(それぞれの試行がもう一方の試行に影響を及ぼさない)であるとすると、その積も同じような性質を持ちます。<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?%7B%5Crm%20E%7D[XY]%20%3D%20%7B%5Crm%20E%7D[X]%7B%5Crm%20E%7D[Y]"><img src="http://formula.s21g.com/?%7B%5Crm%20E%7D[XY]%20%3D%20%7B%5Crm%20E%7D[X]%7B%5Crm%20E%7D[Y].png" /></a></div><br />
この性質を用いると、Xの期待値は以下のように表せます。<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?%7B%5Crm%20E%7D[X]%20%26%3D%26%2010000%20%5Ctimes%20%7B%5Crm%20E%7D[X_1]%7B%5Crm%20E%7D[X_2]%7B%5Crm%20E%7D[X_3]%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%26%3D%26%2010000%20%5Ctimes%201.25%5E3%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%26%3D%26%2019531.25"><img src="http://formula.s21g.com/?%7B%5Crm%20E%7D[X]%20%26%3D%26%2010000%20%5Ctimes%20%7B%5Crm%20E%7D[X_1]%7B%5Crm%20E%7D[X_2]%7B%5Crm%20E%7D[X_3]%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%26%3D%26%2010000%20%5Ctimes%201.25%5E3%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%26%3D%26%2019531.25.png" /></a></div><br />
つまり、このゲームを行うことによってあなたは9500円くらいの利益を得ることが期待できます。もうこれはやるっきゃないと。<br />
<br />
このゲームが得であることは理解できましたが、ただこのゲームにどれくらいリスクがあるのかは不透明です。期待値が10000円プラスになるからといって、40%の確率で全財産を失ってしまうみたいなゲームですと、普通は行いたくないですよね。このリスクについて求めるためには分散を計算してやればいいのですが、これはそれなりに面倒です(*1)。前回ごっちゃごっちゃやってたのはこの面倒さによるものが大半ですので、今回は説明しません。とりあえず分散がリスクに対応しているということが大事です。<br />
<br />
世の中の投資行動には常にリスクとリターンが付きまといます。一般的にリスクとリターンは紙一重と言われており、リスクの高い行動はよりリターンも大きいというのが感覚的な理解です。ただ、このリスクとリターンを確率変数を用いて計算することによって、どの行動が一番リスクとリターンのバランスが取れている行動であるのか、判断することができます。<br />
<br />
で、どの戦略を用いれば合理的に資産を増やすことができるのかという問題ですが、実はこの戦略は既に知られており、<b>ある戦略を用いればリスクをリターンに変化させることで指数関数的に資産を増やすことができます。</b>その戦略について興味がある方は適当に調べると出てくると思います。<br />
<br />
<h4>2. なぜ乗法モデルを仮定しているのか?</h4>まず、株価の変動は確率的であるから確率変数として表せるという話は納得できると思います。では、株価などの資産モデルには様々なものがあるのに、なぜそれぞれの確率変数が独立と仮定した乗法モデル<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?S_N%20%3D%20R_1%20R_2%20%5Ccdots%20R_N%20S_0"><img src="http://formula.s21g.com/?S_N%20%3D%20R_1%20R_2%20%5Ccdots%20R_N%20S_0.png" /></a></div><br />
をわざわざ採用するのでしょうか?別に、例えば<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?S_N%20%3D%20S_0%20+%20R_1%20+%20R_2%20+%20%5Ccdots%20+%20R_N"><img src="http://formula.s21g.com/?S_N%20%3D%20S_0%20+%20R_1%20+%20R_2%20+%20%5Ccdots%20+%20R_N.png" /></a></div><br />
のような加法モデルでも構わないわけです。何故なんでしょう?<br />
<br />
まず独立と仮定した理由については、単に計算が楽だからです。正直独立と仮定しないと、ここまで綺麗な形がでるのか、そもそも解析的に解けるのかどうかすら分かりません。ただこれを仮定したことによって、ごちゃごちゃ計算できて結果的にそれなりに綺麗な数式へと展開することができたという次第です。<br />
<br />
乗法モデルを仮定した理由については、<b>金融業界の指数関数的な性格</b>によるものです。例えば、ローンなんかは年利15%ーフリーローンですとこんくらいでしょうかーとありますが、これは1年間につきローンが15%増加していくので、ふるまいとしては指数的になります。同様に、国債についても年利で計算するので指数的な計算となります。金利が低いといってほいほい借金をすると危険だという理由もこれによるものです。逆に、金持ちがより金持ちになる理由の一つでもあります。<br />
<br />
これは個人的な意見ですが、<b>金融の世界で乗法モデルが成り立っている理由は、人間が金に対して貧欲であるからだ</b>と思っています。はっきりと言いますが、世の中の結構な割合の人間は指数的なふるまいを理解できません。ですからマイホームを数十年ローンでぽんと買ってしまったり、リボ払いで決済をしたりしてしまいます。<br />
<br />
ちなみに、これらの金額はすべて年金公式と呼ばれる公式で算出できます。1期間後に支払いを始め、n期間にわたって金額Aを支払う場合、金利をrとおくと借入額Pは<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?P%20%3D%20%5Cfrac%7BA%7D%7Br%7D%20%5Cleft[%201%20-%20%5Cfrac%7B1%7D%7B%281+r%29%5En%7D%20%5Cright]"><img src="http://formula.s21g.com/?P%20%3D%20%5Cfrac%7BA%7D%7Br%7D%20%5Cleft[%201%20-%20%5Cfrac%7B1%7D%7B%281+r%29%5En%7D%20%5Cright].png" /></a></div><br />
もしくは<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?A%20%3D%20%5Cfrac%7Br%281+r%29%5EnP%7D%7B%281+r%29%5En%20-%201%7D"><img src="http://formula.s21g.com/?A%20%3D%20%5Cfrac%7Br%281+r%29%5EnP%7D%7B%281+r%29%5En%20-%201%7D.png" /></a></div><br />
と表せます。これはすごく役立つ公式なので、店員の言葉に騙されないでおくこともできるかもしれません。ついでに、ローンを借りて全額返した場合に損をする差額はnA-Pで計算できて、<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?nA%20-%20P%20%3D%20%5Cfrac%7B%281+r%29%5En%28nr-1%29+1%7D%7B%281+r%29%5En-1%7DP"><img src="http://formula.s21g.com/?nA%20-%20P%20%3D%20%5Cfrac%7B%281+r%29%5En%28nr-1%29+1%7D%7B%281+r%29%5En-1%7DP.png" /></a></div><br />
となります。こういう情報ってあんまり表沙汰にならないんですけれど。<br />
<br />
<b>追記(10/08/23)</b><br />
上の式は正直汚いからどんな振る舞いをするのか感覚的に理解し辛いですけど、利率が少ない、つまりrが微小のときは (1+r)^n ~ 1 + nr で近似できますから損失額は<br />
<br />
<div style="text-align: center;"><a href="http://formula.s21g.com/?nA-P%20%5Csimeq%20nrP"><img src="http://formula.s21g.com/?nA-P%20%5Csimeq%20nrP.png" /></a></div><br />
となります。ですから利率が低いときは、分割回数をちょっと多くしても線形的な負担にかならない。確かに損失額が増えることは確かですけど、一括払いで生活に支障がでるよりは、3, 4回くらいローンで払うっていう選択肢もありなのかもしれません(*2)。<br />
<br />
*1: 1次元のising模型の解析解を求めたり、極座標のラプラシアン導出するのよりは面倒じゃない<br />
*2: ただそこまでして本当に欲しいものかそれっていう疑問はありますし、30年ローンとかで家を購入するような場合には全く通用しません。そもそもこのご時勢に長期のローンを組むなんて狂気の沙汰です。狂っている。rezoolabhttp://www.blogger.com/profile/09706721371333005060noreply@blogger.com0