Ceres Solver で遊んでみる

仕事に空きができたので、面白そうだけどフォローできていなかったライブラリやツールの調査をしている。
で、Googleがオープンソースで公開している"Ceres Solver"が凄そうだ。非線形の最小二乗法ソルバーでgoogle社内の様々なプロジェクトでも利用されているようだ。
ただwebを探しても情報が結構少なく、チュートリアルもものすごく単純なケースか、Bundler Adjustmentのような複雑な例しかない。以下のような2D図形でゴニョゴニョするサンプルコードなんかがあるといいんだけど見つからない。


そこで、習作としてひとまず2Dポリゴンの位置合わせなんかに使えないかなと実験してみた。

アルゴリズムは単純なものにした。位置合わせ先と元の2Dポリゴンの境界線上で適当にサンプリングして、得られた点群同士を 単純なICPのちょっと変形版で位置合わせを行うというものだ。
パラメータの自由度は 平行移動(x,y)、スケール(x,y)、xy平面での回転の5つ。

誤差関数は以下のような感じ。

誤差関数

template<typename T>
void Transform(T x, T y, T tx, T ty, T r, T sx, T sy, T& xNew, T& yNew) {
    xNew = sx * (x * cos(r) - y * sin(r)) + tx;
    yNew = sy * (x * sin(r) + y * cos(r)) + ty;
}

struct PointResidual {
    PointResidual(const Point& model, const Point& data)
        : model_(model), data_(data) {}

    template <typename T>
    bool operator()(const T* const t, const T* r, const T* s, T* residual) const {
        T x, y;
        Transform(T(data_.x()), T(data_.y()), t[0], t[1], r[0], s[0], s[1], x, y);

        T dx = T(model_.x()) - x;
        T dy = T(model_.y()) - y;
        residual[0] = sqrt(dx*dx + dy*dy);
        return true;
    }

private:
    Point model_;
    Point data_;
};

自動微分を使うために、operator()内の演算をすべてジェネリック型で行ったが、数値微分とあまり結果が変わらなかった。何か間違っているのだろうか?
そして、自動微分なんて技術があるなんて初めて知った。用意された演算子や関数を使って関数書くだけで厳密な微分を行ってくれる。ceresだと普通にdoubleで計算してるみたいなコードで書ける。すごい。

データ群の設定は以下の様な感じ。
ICPに詳しくないのでここの論文を参考にしてみたが、論文では変位した変換元の点に対して、変換先の最近点を対にするだけで良いようだがこれではスケール成分が0になって変換先の1点に収束する現象が生じてしまった。
仕方ないので逆方向の変換先の点に対して、変位した変換元の点も対として用いるようにした。一般的にはどうやって対応しているんだろう?

データ入力

    while (true) {
        ...

        Problem problem;
        std::vector<Point> dataT;
        for (auto ptD : data.outer()) {
            // 変換元の点を現在の変換行列で変位
            double x, y;
            Transform(ptD.x(), ptD.y(), param[0], param[1], param[2], param[3], param[4], x, y);
            Point ptT(x, y);
            dataT.push_back(ptT);

            // 変換先の点群から最も近い点を検索する
            size_t iNearest;
            FindNeareset(ptT, model, iNearest);

            // コスト関数に登録
            CostFunction* cost_function =
                new AutoDiffCostFunction<PointResidual, 1, 2, 1, 2>(
                    new PointResidual(model.outer()[iNearest], ptD));
            problem.AddResidualBlock(cost_function, NULL, param, &param[2], &param[3]);
        }

        for (auto ptM : model.outer()) {
            // 変換先の点に対し、変位した変換元の点群の中で最も近い点を検索する
            size_t iNearest = 0;
            FindNeareset(ptM, dataT, iNearest);

            // コスト関数に登録
            CostFunction* cost_function =
                new AutoDiffCostFunction<PointResidual, 1, 2, 1, 2>(
                    new PointResidual(ptM, data.outer()[iNearest]));
            problem.AddResidualBlock(cost_function, NULL, param, &param[2], &param[3]);
        }

        // 解いてみる
        Solver::Options options;
        options.linear_solver_type = ceres::DENSE_QR;
        options.minimizer_progress_to_stdout = false;
        Solver::Summary summary;
        Solve(options, &problem, &summary);

        // 変位行列が前回と変わらなかったら終了, そうでなければ繰り返し
        ...
    }

結果は以下のようになった。太い灰色線が求めるべき変換先のポリゴン、緑色線が変換元のポリゴン、赤い線がceresで求めた変換行列を緑色ポリゴンに施したもの。サンプリング点数は数百点ほど。

例1

真値: 平行移動=(-20,12), 回転角=30, スケール=(3,0.8)
計算値: 平行移動=(-20.1121,11.9887), 回転角=29.8778, スケール=(2.99513,0.801453)

例2

真値: 平行移動=(-20,10), 回転角=30, スケール=(3,0.8)
計算値: 平行移動=(-20.3742,10.0149), 回転角=29.8033, スケール=(3.01447,0.79894)

結構単純な方法だが結構良い精度で求まっているようだ。
ただしサンプリング点数を少なくすると、やはり変なところに収束してしまう。例2のサンプリング点数を50点ほどに少なくすると以下のようになる。

例3

ICPなんかをもっと勉強しないと実用するのは厳しそうか。
Ceres Solver自体はいろんな用途で使えそうだ。

UE4のビルドでハマった

UE4を githubから取得してVS2015でビルドしているのだが、ある日リポジトリから最新版を取得したところビルドが通らなくなった。
Oculus RiftがらみのDLLやIntelのライブラリがリンクできなくなって、いろいろ試行錯誤したが一向に治らない。
webで検索しても同じ症状の人がちらほら見つかるのだが、解決まで至っていない。

で、結局原因は単純なもので、UE4のルートフォルダのSetup.batとGenerateProjectFiles.batを再実行していないことだった。
ソースの更新内容によっては、これらの手順を再度踏まないといけないようだ。
最初にビルドしてから時間が空いてたから、そんな手順があったってことすっかり忘れてた。

Visual Studio 2015で gflagsを x64/Unicode文字セットのアプリケーションから利用する

google先生の作った gflagsというコマンドライン引数を扱うライブラリを タイトルの環境で使えるようにしたのだけど、とてもつかれた・・・ また同じ思いはしたくないのでメモ。ちなみにまだ簡単なサンプルでしか動作確認していないので注意。

[CMakeのインストール]
・cmakeがインストールされてなければ、CMake公式サイトから win32/x86のバイナリをダウンロードしてインストール

[gflagsのビルド]
1. githubから ソースを入手
2. VS2015の各種ツールのパスの通ったコンソールから、gflagsの CMakeLists.txt があるフォルダに移って、"cmake ."を実行
3. gflags.sln が生成されたはずなので、これをVS2015で開いてビルドすると、ひとまずx86でマルチバイト文字セットのライブラリができる

[gflagsをx64でビルド]
1. gflags.slnを開いた状態で、"ビルド->構成マネージャー"メニューを実行し、右上のアクティブソリューションプラットフォームの欄で新規作成を選び x64プラットフォームを作成する。
2. gflags_staticとgflags_nothred_staticプロジェクトのプロパティを開き、ライブラリアン->その他のオプションの項目で MachineX86オプションが明示的に指定されているのでこれらを除去。
3. これでx64のライブラリがビルドできるはず。

[gflagsをUnicode文字セットのプロジェクトからリンクする]
1. gflagsをUnicode文字セットでビルドするのはソースをガシガシ書き換える必要がありそうなので諦めて、gflags自体はマルチバイト文字セットのままにする。
2. Unicode文字セットのプロジェクトでは普通に gflagsをリンクし、google::ParseCommandLineFlags()で渡す引数をマルチバイトに変換して渡す。
3. これで利用できるはずだが、VS2015ではロケールの扱い周りでバグがあるらしく(?)、std::cout、std::wcoutにどちらかしか日本語出力ができなかったり、mbstowcs系の関数がうまく動かなかったりでハマった。 ひとまず Win32 APIの MultiByteToWideChar()を使ったり, coutへの出力は諦めて対応。コードはざっと以下のような感じ。

#include "../../../gflags/include/gflags/gflags.h"
#pragma comment(lib, "../../../gflags/x64/release/gflags_static.lib")
#pragma comment(lib, "shlwapi.lib") // __imp_PathMatchSpecAのリンクに必要

#include <iostream>
#include <vector>
#include <windows.h>

DEFINE_bool(big_menu, true, "Include 'advanced' options in the menu listing");
DEFINE_string(language, "japanese", "default language");

bool MBS2WCS(const std::string& s, std::wstring& sResult)
{
    int nBufSize = MultiByteToWideChar(CP_ACP, 0, s.c_str(), -1, NULL, 0);
    std::vector<wchar_t> vBuf(nBufSize);
    if (MultiByteToWideChar(CP_ACP, 0, s.c_str(), -1, &vBuf[0], nBufSize) == 0) {
        _ASSERTE(false);
        return false;
    }
    sResult = std::wstring(&vBuf[0]);
    return true;
}
bool WCS2MBS(const std::wstring& s, std::string& sResult)
{
    int nBufSize = WideCharToMultiByte(CP_ACP, 0, s.c_str(), -1, NULL, 0, NULL, NULL);
    std::vector<char> vBuf(nBufSize);
    if (WideCharToMultiByte(CP_ACP, 0, s.c_str(), -1, &vBuf[0], nBufSize, NULL, NULL) == 0) {
        _ASSERTE(false);
        return false;
    }
    sResult = std::string(&vBuf[0]);
    return true;
}

int _tmain(int argc, TCHAR* argv[], TCHAR* envp[])
{
    setlocale(LC_ALL, "");
    std::locale loc("");
    std::wcout.imbue(loc);

    std::vector<char*> argvA(argc);
    std::vector<std::vector<char>> buf(argc);
    for (int i = 0; i < argc; ++i) {
        std::string s;
        WCS2MBS(argv[i], s);
        buf[i].resize(s.size() + 1);
        strcpy_s(&buf[i][0], buf[i].size(), s.c_str());
        argvA[i] = &buf[i][0];
    }

    gflags::SetUsageMessage("gflagsライブラリのテスト");
    char** pA = &argvA[0];
    google::ParseCommandLineFlags(&argc, &pA, true);

    std::wcout << FLAGS_big_menu << std::endl;

    // coutに日本語が出力できないのでwcoutに出力
    std::wstring s;
    MBS2WCS(FLAGS_language, s);
    std::wcout << s << std::endl;

    return 0;
}

実行するとこんな感じ。

>gflag_test.exe -language 日本語
1
日本語

>gflag_test.exe -language English -nobig_menu
0
English

江添亮さんがC++11の解説本を公開してくれました!

C++標準化委員会のメンバーで、最近はクッキーを増やすゲームのおそらく日本で一番面白いレビューを書いてくれた江添さんが、C++11の解説本をなんとコピーレフトで公開してくれました。

http://ezoeryou.github.io/cpp-book/C++11-Syntax-and-Feature.xhtml で見ることができます。

まだ途中までしか見ていませんが、簡潔で、かつ、プログラマの視点で重要な部分の説明に重点を置いた本ですね。
さらには仕様の背景にある考えも説明してくれていてありがたいです。

重要で知っておくべきなのに、複雑だしあまり面白くないのがADLのあたりだと思いますが、この本ではかなり分りやすく書いてあるのでC++プログラマは必見だと思います。
「オーバーロードの解決があいまいです」ってエラーがでてウガーッ!!!ってなると思いますが、ADLの仕様を知っておくと対処が楽になりますね。

とりあえず、auto、ラムダ式はあたりは当然ですが、地味なところでscoped enumと生文字列リテラルを一刻もはやく使って幸せになりたいです(;_;)

elevatedで使われている技術

2009年にスペインのiq氏がelavatedという衝撃的なデモを発表した。
わずか4KBのプログラムで、リアルな山岳地帯の風景を3Dでリアルタイム描画し、音楽も鳴らすというものだ。


「たった4KBでこんなことができるわけがない! チートだ!」と叫びたくなるプログラムだが、製作者(iq氏)がいろんな資料を公開してくれていて、 それを読むと「天才ならこのぐらいはできるかもな」ぐらいには思えるようになる。

iq氏の資料はココにまとまっている。

elavatedの概要のpdfを読むと、

• 地形は2D格子で、各格子点の高さはパーリンノイズで求める。(パーリンノイズについての詳細)
• 通常のポリゴンベースの描画ではなく、ピクセルシェーダを使ったレイマーチン法でレンダリング。(地形とのレイマーチン法についての詳細)
• カメラの位置、方向、パスもピクセルシェーダ内で算出
• 地形の法線ベクトルは隣のピクセルとの差分から算出するのではなく、パーリンノイズの式を微分して解析的に求める。
• 影はフェイクで、ローパスフィルタをかけた地形の拡散光を影の成分とする。

といった手法を使っているようだ。

さらにはelevatedのShaderToy版もあり、簡易版ながらコードを見たり、いじることもできる。
魔法がとけてよかった。

MSXML6のgettext()が異常に遅い件

MSXML6をXMLの読み込みに使っているが、国土地理院のページで配布されている標高メッシュのXMLを読み込ませたところ、これがとにかく遅い！

問題のXMLは以下のように values要素が約84万個大量に並んでいるもので、容量はおよそ100MB程度である。

<jps:CV_GridValuesMatrix>
  <jps:values>
    <jps:memberValue>
      <type>その他</type>
      <alti>-9999.00</alti>
    </jps:memberValue>
  </jps:values>

  <jps:values>
    <jps:memberValue>
      <type>その他</type>
      <alti>-9999.00</alti>
    </jps:memberValue>
  </jps:values>
    ...
</jps:CV_GridValuesMatrix>

これを以下のようなコードで読み込ませてみると、数10分かかっても処理が終わらない。

MSXML2::IXMLDOMNodePtr p = pCV_GridValuesMatrix->firstChild;
while (p) {
  if (p->GetnodeName() == _bstr_t(L"jps:values")) {
    MSXML2::IXMLDOMNodePtr pAlti = p->firstChild->firstChild;
    double d = _wtof(pAlti->Gettext());
  }

  p = p->nextSibling;
}

プロファイルをとってみると、Gettext()関数が処理の大半を占めている。
同じフォーマットで10MB程度のXMLは数秒で終わることから、どうやらGettext()がファイルサイズに対してO(N^2)のような処理オーダーになっている印象を受ける。
これは奇妙だ。Gettext()の処理内容を想像すると定数時間で済みそうなものなのに!

試行錯誤して、Gettext()の変わりにfirstChild->nodeValue を呼ぶようにしたら処理は数秒で終わるようになった。たったこれだけで数千倍も処理時間が違う。

MSXML2::IXMLDOMNodePtr p = pCV_GridValuesMatrix->firstChild;
while (p) {
  if (p->GetnodeName() == _bstr_t(L"jps:values")) {
    MSXML2::IXMLDOMNodePtr pAlti = p->firstChild->firstChild;
    double d = static_cast(pAlti->firstChild->nodeValue);
  }

  p = p->nextSibling;
}

ってことで、いまいち理由がわからないが巨大なXMLを扱うときには Gettext()は呼び出したらアウトのようだ。

ちなみに今回の現象で悩んだ際にMSXML関連の情報をあさってみたが、以下のようなコードで遅いと言っている人もいた。

MSXML2::IXMLDOMNodeListPtr l = pCV_GridValuesMatrix->GetchildNodes();
for (int i = 0; i < l->Getlength(); ++i) {
  MSXML2::IXMLDOMNodePtr p = l->Getitem(i);
  ...
}

これを上記のXMLでやると、このトラバースだけで数10分かかっても終わらない。
これはおそらく DOMNodeListは線形リストの実装になっていて、ランダムアクセスにO(N)の時間がかかるためだろう。
MS公式のFAQにも「こういうトラバースはしないで!」という情報があったが、だったら初めからGetitem()なんて関数用意しなきゃいいのに。

猿でも分かった～型消去技法とは

ひさびさにC++の話。

C++の技法に型消去(Type Erasure)というものがありますね。なんかかっこいい名前だし、いろんなすごいライブラリで使われてるっぽい話はよく聞くので、「Type Erasure? なにそれ? やだ、こわーい」となる人もいるでしょう。(ぼくもそう)
でも型消去はトモダチ、怖くない!

唐突ですが、継承関係はないけどとにかく"f"という関数を持ったクラス群があったとしましょう。

struct A {
  void f() { cout << "A!"; }
};
struct B {
  void f() { cout << "B!"; }
};

役にたつかどうかは分からんけど、こういったものを保持しておいて後からfを呼び出すことができるクラスholderが作れるか考えてみましょう。
これは簡単で、テンプレートを使って以下のように書けますね。

template<typename T>
struct holder {
  holder(const T& obj) : m_obj(obj) {}

  void f() { m_obj.f(); }

  T m_obj;
};

holder<A> a(A());
holder<B> b(B());
a.f();
b.f();

さてここでholderはテンプレートクラスですが、テンプレートでないクラスにすることはできるでしょうか? もしこれができれば、以下のようにとにかくfという関数をもったあらゆるオブジェクトをobjsという配列に溜め込んで、その後適切にfを呼び出すようなことができます。

vector<holder> objs;

...

for (size_t i = 0; i < objs.size(); ++i) {
  objs[i].f();
}

驚くべきことに、これが実に単純な方法で実現できます。仮想関数とテンプレートを組み合わせて使って型に関する情報を派生クラス側に追い出すだけです。

struct holder {
  virtual void f() = 0;
};

template<typename T>
struct holder_sub : public holder {
  holder_sub(const T& obj) : m_obj(obj) {}

  virtual void f() { m_obj.f(); }

  T m_obj;
};

使ってみましょう。

vector<shared_ptr<holder> > objs;
objs.push_back(new holder_sub<A>(A());
objs.push_back(new holder_sub<B>(B());

for (size_t i = 0; i < objs.size(); ++i) {
  objs[i]->f();
}

ついでに、使うときに派生クラスを見せるのはかっこ悪いので以下のようなラッパーをかませましょう。

struct anyf
{
  template<typename T> anyf(const T& a) { m_p.reset(new holder_sub<T>(a)); }

  void f() { m_p->f(); }

  shared_ptr<holder> m_p;
};

するとこんな風に使えます。

vector<anyf> objs;
objs.push_back(A());
objs.push_back(B());

for (size_t i = 0; i < objs.size(); ++i) {
  objs[i].f();
}

anyfにはfをもった要素ならなんでも突っ込めるし、そうでなければちゃんとコンパイルエラーになってくれます。
これが型消去技法です。きわめて単純なコードなのに、魔法みたいな挙動が実現できてますね!

これまでのサンプルの関数fは非常につまらない関数でしたが、元の型にキャストするような関数に適用するとどうでしょうか?
そう、boostの超便利クラスanyみたいなものが出来そうです。anyも基本的にはこの単純な型消去技法を使って実装されているようです。

Boost.Geometryの設計がすごい

最近のboost C++ライブラリは数値計算や幾何計算のほうも充実してきているようで、先日Boost.Geometryの2Dブーリアンを使ってみたらそのパフォーマンスと頑強性に驚かされた。
さらにこのBoost.Geometryは設計が素晴らしい。そのDesign Ratonaleを読んでみたところとても面白かったので、より多くの人に読んでもらいたいなと思って和訳してみた。
わかりやすい題材を元に、traits, タグディスパッチ、コンセプト、メタ関数などのテンプレート周りの独特の技法が少しずつ登場してくるので、これらの概要を把握するにもとても良い資料だと思う。

原文はココ。

なお、「ジェネリック」、「特殊化」といったジェネリックプログラミング独特の用語が頻出するので、なじみのない方は あらかじめここで概要をつかむといいかもしれない。

---

設計の理論的基礎(Design Rationale)

仮に2つの点の間の距離を計算するC++プログラムを、あなたが必要になったとする。 あなたは構造体を定義するかもしれない。

struct mypoint
{
    double x, y;
};

そしてアルゴリズムを記述する関数も。

double distance(mypoint const& a, mypoint const& b)
{
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}

とてもシンプルで、使うにも十分だ。しかし、ジェネリックではない。ライブラリとして、もっと先に進んだデザインにするべきだ。 上記のデザインで利用可能なのは、2次元の点の、しかもmypointという構造体の、しかも直交座標系に限ってである。 ジェネリックなライブラリは以下のような距離の計算ができるべきである:

• mypoint構造体だけでなく、あらゆる点のクラス、構造体に対して
• 2次元以上の空間でも
• 地球や球面などの他の座標系においても
• 点や線分、その他の幾何形状との間で
• doubleよりも高い精度でも
• 平方根の計算を避けつつ(平方根の計算は相対的にコストの高い関数で、距離の比較では必要がないので、使いたくない)

こことそれに続くセクションでは、ステップ・バイ・ステップでよりジェネリックなこのような設計をしていく。

テンプレートの利用(Using Templates)

distance関数はテンプレート関数へと変更できる。これは簡単で、mypoint型以外の他のポイント型についても距離を計算することができるようになる。 異なる2つのポイント型も許容させるために、2つのテンプレート引数を追加する。

template <typename P1, typename P2>
double distance(P1 const& a, P2 const& b)
{
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return std::sqrt(dx * dx + dy * dy);
}

このテンプレートバージョンは多少は良い。しかし十分ではない。

メンバ変数x やy に直接アクセスできない C++クラスを考えてみよう。したがって、この節はさっさと終えて、次に進もう。

Traitsの利用(Using Traits)

我々はジェネリックなアプローチを取り、distance関数の入力にあらゆるポイント型を許可する必要がある。 メンバ変数x やy に アクセスする代わりに、traitsシステムを用いて、数段階の間接参照を追加する。関数は次のようになる:

template <typename P1, typename P2>
double distance(P1 const& a, P2 const& b)
{
    double dx = get<0>(a) - get<0>(b);
    double dy = get<1>(a) - get<1>(b);
    return std::sqrt(dx * dx + dy * dy);
}

この改修されたdistance関数は、点の座標にアクセスするため、次元をテンプレート引数として受け取るジェネリックなget関数を呼び出している。 このget関数は次のように定義されたtraitsシステムへと向かう:

namespace traits
{
    template <typename P, int D>
    struct access {};
}

以下は我々のmypoint型に特殊化されたもので、getというstatic関数で実装されている:

namespace traits
{
    template <>
    struct access<mypoint, 0>
    {
        static double get(mypoint const& p)
        {
            return p.x;
        }
    };
    // same for 1: p.y
    ...
}

traits::access<mypoint, 0>::get(a) を呼び出すことで我々のx座標が返ってくる。ナイスと思う? この関数の呼び出し方法はよくライブラリに用いられるものとしてはあまりに冗長である。我々は次の短い関数(free extra function)を追加できる:

template <int D, typename P>
inline double get(P const& p)
{
    return traits::access<P, D>::get(p);
}

これによって、この節の最初のdistanceアルゴリズム内で記述されているように、traits::accessが特殊化されたあらゆるポイント型に対して、 get<0>(a)を 呼び出すことができる。 したがってx()のようなメソッドを持つクラスに対して使いたいなら、 次元0に対してx()を、次元1に対してy()を、 返すようなスタティックなget関数を通して構造体にアクセスするような特殊化があれば良い。

次元への非依存性(Dimension Agnosticism)

今、我々はあらゆる構造体やクラスの2次元の点の間の距離を計算可能となった。 しかし、我々は3次元も扱いたい。なので、次元について非依存にする必要がある。 これは我々のdistance関数をややこしくする。 我々はforループを使って次元を歩き通すことができる。 だがforループは余分なパフォーマンスが必要となる。 我々はさらなるテンプレートの利用によって、より複雑だがテンプレートマニアを魅了させるようなdistanceアルゴリズムを作ることができる:

template <typename P1, typename P2, int D>
struct pythagoras
{
    static double apply(P1 const& a, P2 const& b)
    {
        double d = get<D-1>(a) - get<D-1>(b);
        return d * d + pythagoras<P1, P2, D-1>::apply(a, b);
    }
};

template <typename P1, typename P2 >
struct pythagoras<P1, P2, 0>
{
    static double apply(P1 const&, P2 const&)
    {
        return 0;
    }
};

distance関数は次元数を指定して、pythagoras構造体を呼び出す:

template <typename P1, typename P2>
double distance(P1 const& a, P2 const& b)
{
    BOOST_STATIC_ASSERT(( dimension<P1>::value == dimension<P2>::value ));

    return sqrt(pythagoras<P1, P2, dimension<P1>::value>::apply(a, b));
}

次元は次のように定義された別のtraitsクラスを使って参照される:

namespace traits
{
    template <typename P>
    struct dimension {};
}

これはmypoint構造体で特殊化される必要がある。

値を発行するだけなので、Boost.MPLの class boost::mpl::int_を派生させると便利である:

namespace traits
{
    template <>
    struct dimension<mypoint> : boost::mpl::int_<2>
    {};
}

前述のget関数のように、ライブラリはdimensionメタ関数も持つ。

template <typename P>
struct dimension : traits::dimension<P>
{};

なぜこの宣言が役に立つかを説明しよう。今、我々は次元に対して非依存となった。 distance関数は3次元やそれ以上の次元の点を扱えるようになったが、 2次元の点aと3次元の点bに対する演算をコンパイル時のアサーションによって防ぐことができる。

座標値型(Coordinate Type)

これまで我々は座標値はdouble型と決めてかかってきた。もし我々の点がinteger型だったらどうする?

traitsクラスを簡単に加えることができるし、実際にそうする。しかしながら、2つの整数値座標の距離はdoubleのままにできる。 また、それとは別にして、デザインゴールは平方根を避けることだった。 続く別の節で、これらのケースを扱う。しかし今は2つの整数値座標を扱うが結果はdoubleとして返すこととしよう。 座標値型を定義するには、もう1つ別のtraitsクラスcoordinate_typeを追加する。 これはライブラリのユーザーによって以下のように特殊化されるべきである:

namespace traits
{
    template <typename P>
    struct coordinate_type{};

    // specialization for our mypoint
    template <>
    struct coordinate_type<mypoint>
    {
        typedef double type;
    };
}

access関数でgetフリー関数を定義したように、プロキシを同様に加える。より長いバージョンは後ほど登場するが、短いバージョンは次のようになる:

template <typename P>
struct coordinate_type :    traits::coordinate_type<P> {};

今再びdistance関数が修正できる。依然としてdouble型で結果を返すので、我々は pythagorasクラスだけ修正すればよい。 結果はインプットと同じ座標値型を返すべきである。しかしインプットは2つの異なる型を受け取ることができる。 それらは異なる座標値型かもしれない。これはあまりないことであるが我々はジェネリックなライブラリをデザインしており、 これらの特異なケースも扱えるべきである。 我々は2つのうち精度の高い方を座標値型として選ぶ必要がある。 このことについてはここでは説明しない。あまりに長くなることと、幾何とは関係がないからだ。 今は、適切な型を選ぶメタ関数としてselect_most_preciseがあることと仮定しよう。

我々の計算クラスはこのようになる:

template <typename P1, typename P2, int D>
struct pythagoras
{
    typedef typename select_most_precise
        <
            typename coordinate_type<P1>::type,
            typename coordinate_type<P2>::type
        >::type computation_type;

    static computation_type apply(P1 const& a, P2 const& b)
    {
        computation_type d = get<D-1>(a) - get<D-1>(b);
        return d * d + pythagoras <P1, P2, D-1> ::apply(a, b);
    }
};

異なる幾何タイプ(Different Geometries)

あらゆる型と座標値の点をサポートする次元に非依存なシステムを設計できた。 まだいくつかの調整があるが後ほど説明しよう。 今からはどのようにして、点とポリゴンや点と線分との距離を計算するかを見ていこう。 これらの公式はより複雑でデザインへの影響もより大きいものとなる。 我々は次のような別の名前を持つ関数を加えたくはない:

template <typename P, typename S>
double distance_point_segment(P const& p, S const& s)

我々はジェネリックにしたい。distance関数は扱う幾何タイプを知ること無しに呼び出される必要がある。 なので名前はdistanceである必要がある。 我々はまたオーバーロードすることもできない。なぜなら同じテンプレートシグネチャを持つ場合に曖昧になってしまうからだ。 これには2つの解決策がある:

• tag dispatching(タグディスパッチ)
• SFINAE

我々はタグディスパッチを選んだ。traitsシステムによくフィットしたからだ。 初期のバージョンのBoost.GeometryではSFINAEを使っていたが、このような大きなデザインに対していくつかの 欠点が見つかったため、タグディスパッチへと変更された。

タグディスパッチによってdistanceアルゴリズムは入力パラメータの幾何タイプを調べることができる。 distance関数はこのように変更される:

template <typename G1, typename G2>
double distance(G1 const& g1, G2 const& g2)
{
    return dispatch::distance
        <
            typename tag<G1>::type,
            typename tag<G2>::type,
            G1, G2
        >::apply(g1, g2);
}

tagメタ関数を参照し、dispatch::distanceのapplyメソッドを呼び出している。 tagメタ関数はまた別のtraitsクラスであり、ポイント型ごとに特殊化される。 それらはこのようになる:

namespace traits
{
    template <typename G>
    struct tag {};

    // specialization
    template <>
    struct tag<mypoint>
    {
        typedef point_tag type;
    };
}

coordinate_systemやgetのようなフリーメタ関数は:

template <typename G>
struct tag : traits::tag<G> {};

Tags (point_tag, segment_tag, etc)はdispatch構造体を特殊化する目的の空の構造体である。 distanceに対するdispatch構造体とその特殊化は次のような隔離された名前空間に全て定義される:

namespace dispatch {
    template < typename Tag1, typename Tag2, typename G1, typename G2 >
    struct distance
    {};

    template <typename P1, typename P2>
    struct distance < point_tag, point_tag, P1, P2 >
    {
        static double apply(P1 const& a, P2 const& b)
        {
            // ここにすでにやってきたような
            // pythagorasを呼び出す
            ...
        }
    };

    template <typename P, typename S>
    struct distance
    <
        point_tag, segment_tag, P, S
    >
    {
        static double apply(P const& p, S const& s)
        {
            // ここで 2,3次元での
            // 点と線分の距離を計算する
            // 別の関数を参照する
            ...
        }
    };

    // ここには以下のようなさらなる特殊化が
    // あるかもしれない
    // 点-ポリゴン, ボックス-円, etc.

} // namespace

これでdistanceアルゴリズムはジェネリックになり、異なる幾何タイプもサポートできるようになった。 1つの欠点は、点-線分、線分-点の2つの特殊化をする必要があることだ。 これは後の節で解決される。以下は異なるポイント型、幾何タイプ、次元に対するサンプルである:

point a(1,1);
point b(2,2);
std::cout << distance(a,b) << std::endl;
segment s1(0,0,5,3);
std::cout << distance(a, s1) << std::endl;
rgb red(255, 0, 0);
rbc orange(255, 128, 0);
std::cout << "color distance: " << distance(red, orange) << std::endl;

核心部分に再び戻る(Kernel Revisited)

traits名前空間でそれぞれ定義されたtraitsクラスcoordinate_typeを使ってきた。 これは実際は必要ではない。なぜなら唯一の違いは名前空間だからだ。しかし今や我々は他の幾何タイプを持つようになった。これは根本的な問題だ。 我々は、あらゆる幾何タイプ、点、線分、ポリゴンなどに対応した、タグディスパッチを用いて実装されたcoordinate_typeメタ関数を 呼び出すことができる:

template <typename G>
struct coordinate_type
{
    typedef typename dispatch::coordinate_type
        <
            typename tag<G>::type, G
        >::type type;
};

dispatch名前空間の内部ではこのメタ関数は2回実装される。ジェネリックバージョンと点に特殊化されたものである。 点に特殊化されたものはtraitsクラスを呼び出す。ジェネリックバージョンは点の特殊化されたものを呼び出す。 再帰的なメタ関数定義のようなものである:

namespace dispatch
{

    // Version for any geometry:
    template <typename GeometryTag, typename G>
    struct coordinate_type
    {
        typedef typename point_type
            <
                GeometryTag, G
            >::type point_type;

        // Call specialization on point-tag
        typedef typename coordinate_type < point_tag, point_type >::type type;
    };

    // Specialization for point-type:
    template <typename P>
    struct coordinate_type<point_tag, P>
    {
        typedef typename
            traits::coordinate_type<P>::type
            type;
    };
}

さて、別のメタ関数point_typeを呼び出している。ここでは詳しく説明しないが、これが全ての幾何タイプ、幾何を形作る点のtypedefに対して呼び出すことを可能にしている。

これと同様のものを メタ関数dimensionや続いて登場するメタ関数coordinate sysytemに適用する。

座標系(Coordinate System)

ここまでは直交座標系と仮定してきた。しかし我々は地球は平面でないことを知っている。 直交座標系で2つのGPS座標の距離を計算することは無意味である。 したがって再び設計を拡張する。traitsシステムを再び使うことで、 それぞれのポイント型に対して座標系を定義する。それから座標系に合った計算を作成する。

座標系は座標値型と似ている。メタ関数で、dispath関数をもつあらゆる幾何タイプを呼び出し、 その点の特殊化に向かい、最終的にはtraitsクラスを呼び出し、座標系を持った型を定義する。 ここではその全ては示さない。定義の一部は以下のとおりである:

struct cartesian {};

template<typename DegreeOrRadian>
struct geographic
{
    typedef DegreeOrRadian units;
};

直交座標系(cartesian)は単純だ。地理座標系(geographic)では、度(degree)またはラジアン(radian)のいずれであるかが選択できる。

distance関数は、座標系に対応した計算方法を選択し、それからdistanceに対するdispath構造体を呼び出すように変更される。 我々は座標系に特殊化された計算方法をストラテジ(strategy)と呼ぶ。新しいバージョンのdistance関数はこのようになる:

template <typename G1, typename G2>
double distance(G1 const& g1, G2 const& g2)
{
    typedef typename strategy_distance
        <
            typename coordinate_system<G1>::type,
            typename coordinate_system<G2>::type,
            typename point_type<G1>::type,
            typename point_type<G2>::type,
            dimension<G1>::value
        >::type strategy;

    return dispatch::distance
        <
            typename tag<G1>::type,
            typename tag<G2>::type,
            G1, G2, strategy
        >::apply(g1, g2, strategy());
}

strategy_distanceは 異なる座標系に対する特殊化をもった構造体である。

template <typename T1, typename T2, typename P1, typename P2, int D>
struct strategy_distance
{
    typedef void type;
};

template <typename P1, typename P2, int D>
struct strategy_distance<cartesian, cartesian, P1, P2, D>
{
    typedef pythagoras<P1, P2, D> type;
};

再びpythagorasが登場していて、ここではストラテジとして定義されている。 distanceはapplyメソッドを呼び出すことで関数をdispatchする。

そしてこれが重要なステップとなる。球面および地理座標系に対して、別のストラテジを実装することが可能になる。 球面座標系はhaversine公式がある。したがってディスパッチするtraits構造体はこのように特殊化される

template <typename P1, typename P2, int D = 2>
struct strategy_distance<spherical, spherical, P1, P2, D>
{
    typedef haversine<P1, P2> type;
};

// struct haversine with apply function
// is omitted here

地理座標系では我々はいくつかの別の距離計算方法がある。 Andoyer法は高速で正確だ。そしてVincenty法は遅いがより正確だ。 また、他にも精度の劣るアプローチがいくつか存在する。

座標系ごとにデフォルトのストラテジが1つ定義されている。 他のストラテジも同様に使えるようにするため、我々は設計を再び修正し、 distanceアルゴリズムが第3引数としてストラテジを受け取るようにオーバーロードを 追加する。

この新しいオーバーロードのdistance関数は、ストラテジがdistance関数の外側で構築できるという利点も持つ。 なぜならもし内側で構築されると、構築時のパラメータを持つことができないからである。Andoyer法やVincenty法、haversine公式では 地球の半径をパラメータとして受け取るコンストラクタを持つことで意味をなす。

オーバーロードされたdistance関数はこのようになる:

template <typename G1, typename G2, typename S>
double distance(G1 const& g1, G2 const& g2, S const& strategy)
{
    return dispatch::distance
        <
            typename tag<G1>::type,
            typename tag<G2>::type,
            G1, G2, S
        >::apply(g1, g2, strategy);
}

ストラテジは2つの点を引数として受け取るメソッドを持つ必要がある。 これはスタティックなメソッドである必要はない。ストラテジは変数を受け取り、メンバ変数に格納するための コンストラクタを定義するかもしれない。このようなケースではスタティックなメソッドでは不便だろう。 これは通常のメソッド(const性を持った)として実装できる。

ここでは全ての実装を列挙しない。Vincenty法は大半のページは数学で埋め尽くされている。しかしアイディアは理解できるだろう。 我々はdistanceを以下のように呼び出すことができる:

distance(c1, c2)

ここで c1 と c2 は直交座標系での点である。または:

distance(g1, g2)

ここで g1 と g2 は地理座標系の点であり、 デフォルトのストラテジ(例えばAndoyer)を呼び出している。または:

distance(g1, g2, vincenty<G1, G2>(6275))

ここではストラテジは明示的に指定され、半径値で構築されている。

ポイントコンセプト(Point Concept)

ここまでの節で述べられた5つのtraitsクラスはポイントコンセプト(Point Concept)を形作る。 traits名前空間で特殊化が実装されたあらゆるポイント型は、有効な型として受け付けられるだろう。 Point Conceptは以下から成り立つ:

• traits::tagの特殊化
• traits::coordinate_systemの特殊化
• traits::coordinate_typeの特殊化
• traits::dimensionの特殊化
• traits::accessの特殊化

最後のものはgetメソッドと、(オプションで)setメソッドを持ったクラスである。 前の4つは型定義と値の宣言をするメタ関数である(MPLに従う)。

さて、我々は今、次元の数への非依存性、座標系への非依存性を得た。 このデザインはあらゆる座標値型、異なる幾何タイプを扱うことができる。 さらには独自のストラテジを指定でき、2つの点の次元が異なる場合は(アサーションによって) コンパイルされないようになる。さらには2つの点の座標系が異なる場合もコンパイルされない(特殊化が存在しないため)。 ライブラリはポイント型に課せられたコンセプトを満たしているかチェックできる。 これはやがてConcept Checkingの節で扱われる。

返り値型(Return Type)

我々はstd::sqrt の呼び出しはいつでも必須というわけではないと約束した。 なのでdoubleに変換可能な2乗の値を持つdistance result構造体を定義する。 しかしながら、pythagorasにのみ定義される。 球面での距離関数は平方根の計算をしない。そのため比較的高価な平方根の計算を避けることが必要でないためだ。 それらは単純に距離を返す。

distance result構造体はストラテジに依存する。そこでストラテジのメンバ型を作成する。 result構造体は次のようになる:

template<typename T = double>
struct cartesian_distance
{
    T sq;
    explicit cartesian_distance(T const& v) : sq (v) {}

    inline operator T() const
    {
        return std::sqrt(sq);
    }
};

平方根を計算することなしに自分自身と他の結果を比較する演算子を持つかもしれない。

それぞれのストラテジはその内部に返り値型を持つべきである。例えば:

typedef cartesian_distance<T> return_type;

あるいは:

typedef double return_type;

といったように、直交座標系(pythagoras)や球面座標系それぞれに。

再びdistance関数は新しい返り値型に対応するために修正される。 ストラテジによってオーバーロードし、これはそれほど複雑ではない。

template < typename G1, typename G2, typename Strategy >
typename Strategy::return_type distance( G1 const& G1 , G2 const& G2 , S const& strategy)

ストラテジを持たないものを除いて、我々はストラテジや座標値タイプなどを選択する必要がある。 それをやるには一行では長すぎるので我々は別のメタ関数を追加する:

template <typename G1, typename G2 = G1>
struct distance_result
{
    typedef typename point_type<G1>::type P1;
    typedef typename point_type<G2>::type P2;
    typedef typename strategy_distance
        <
            typename cs_tag<P1>::type,
            typename cs_tag<P2>::type,
            P1, P2
        >::type S;

    typedef typename S::return_type type;
};

そしてdistance関数を修正する:

template <typename G1, typename G2>
inline typename distance_result<G1, G2>::type distance(G1 const& g1, G2 const& g2)
{
    // ...
}

もちろんdispatchの特殊化においてもこのような結果を返す。 それらはストラテジをテンプレートパラメータとしてどこでも持ち、できるかぎり冗長性のないバージョンとなっている。

まとめ(Summary)

この設計の理論的基礎では、Boost.Geometryがタグディスパッチ、コンセプト、traits、そしてメタプログラミングによって 一歩ずつ設計されている。 我々はデザインを示すため、よく知られたdistance関数を使った。

Boost.Geometryはここで述べられたようなものと、もういくつかのテクニック(automatically reversing template arguments, tag casting, and reusing implementation classes or dispatch classes as policies in other dispatch classes.)で設計されている。

Copyright © 2011 Barend Gehrels, Bruno Lalande, Mateusz Loskot Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

---

英語は苦手なので、間違いがあったらツッコんでください。