Tuesday, October 23, 2012

Interesting ideas in CACM (Vol.55, No.9)


CACM (Vol.55, No.9)に 2 つほど面白い話があったので紹介する.

1. DieHard

多くの C++ で書かれたソフトウェアのバグは多少の buffer overflow である.そこで memory allocator が random に allocation のサイズを拡大する.また,allocator は heap の allocation を randomize する.すると,この種のバグはdeterministic ではなく,Heisenberg bug になる.

これは devloper からすると悪夢である.バグが再現できたりできなかったりするので,それを修正するのが困難になるからだ.

しかし,ユーザからすれば,バグは起きないこともある.クラッシュしたら再起動すれば,運が良ければクラッシュしないかもしれない.

この論文は DieHard というタイトルで Microsoft research から入手できる.しかしバグを隠すことは Seatbelts や Airbags のアナロジではないと個人的には思う.面白いアイデアとは思うが,プログラマとしては deterministic なバグの方が好みであるし,できればバグは隠すよりも直したいと思う.

http://queue.acm.org/detail.cfm?id=2333133

2. CriptDB

暗号化したデータをサーチしたり,圧縮したデータを圧縮したまま処理する方法というのがあると聞いていったいどうするのだろうと思ったことがある.詳細は私にはわからないが,一方向関数でも対応がある程度明らかなものでは可能のようだ.たとえば,英語の Databse があるとする,これをドイツ語に翻訳するのが暗号化と思えば,ドイツ語の単語でサーチができる.たとえば,山(Mountain)という単語を探そうとすると,このキーワードを暗号化した Bergによって探索する.話はそんなに単純ではないが,原理はこれに類したものであるようだ.

http://dl.acm.org/citation.cfm?doid=2330667.2330691

Monday, October 22, 2012

鉄道ファンのためのグラフ理論: マルコフ行列の中の著者達番外編


この話は マルコフ行列の中の著者達の番外編である.鉄道ファンとグラフ理論の関係は何かという話である. (English version)

同僚の Dietger がグラフの全ての edge をたどる問題を知っているかと聞くので,Hamiltonian path のことかと答えると,Hamiltonian path は全ての nodeを一度通るもので,そうではなく全ての edge を通るものだということだった.

私はVertex と Edge の Dual を考えれば良いではないかと言ったのだが,Triangle mesh ではVertex と Face は Dual になるが,Edge はそうではない.(図1, 2, 3)考えてみると,点は辺を介して接続され,面も辺を介して接続されている.では辺は何を介して接続されているのだろうか.面か点である.面と点の接続関係は辺によるものなので,辺はある意味この面・点とは異なる.だから面と点を入れかえることはまだできるが,辺は何かと入れかえられないと考えた.

Figure 1: Duality of face and vertex: faces to a dual grah.
Figure 2: Duality of face and vertex: a dual graph to faces.

Figure 3: Edge's duality?

これはどんな問題を考えているのかと尋ねたところ,彼の友人は鉄道の線路の管理をしているということであった.つまり線路を定期的に全て検査しなくてはならない.効率良く線路を全てまわるのはどうすれば良いのかという問題をグラフの問題として考えたものである.結局,これは中国人の郵便配達問題と呼ばれていることがわかった.中国人の数学者(Mei Ko Kuan)がこの問題に関して研究していたことにちなむらしい.これは Euler path 問題,つまり一筆書き問題とかかわりが深い.

私は昔,PTT (Programming tools and techniques) という集まりでこれに類する話を聞いたことがある.この会はソフトウェア関係の会であるが,鉄道好きの人も多数いらした.葛西さんという方が,この会で発表したのは,最長の片道切符はどんなものかというものであった.これは JR (Japan Railway) で買うことができる最長の切符である.http://www.swa.gr.jp/lop/index.html

referred from http://www.swa.gr.jp/lop/lop_res.html
彼の計算では線形計画法と総当たり法の二つで求めた解が一致している.ただし,総当たり法ではグラフの簡略化や分割など様々なテクニックが利用されている.

葛西氏は日本語ではあるが,アルゴリズムに関しても詳しく説明されている.2000年当時,彼の使った線形計画法のソフトウェアは商用の高価なものであったが,2005 年には計算機の性能の向上,そして,GLPK (GNU Linear Programming Kit) の登場によって家庭でもこのような計算が可能になっているそうである.その方法も示されている.

驚くことに1960年頃から手計算で片道切符の旅に取り組んでいる方がいらして,その方の経路はほとんどこの計算結果と同じというのも驚いた.人間の脳がすごいのか,この人が特殊なのか,とにかくすごいことである.

ところで,世界には鉄道ファンは多数いるはずである.ドイツにも鉄道ファンのテレビがあると聞く.そういう番組では,きっとこのような試みがあるであろう.多分.しかし,日本は地理的条件からグラフの分割というものが可能である.(グラフ的には島を結ぶ辺が極端に少ないので,葛西氏はこれを利用して計算量を減らしている.)ドイツの鉄道の最長片道切符,あるいはヨーロッパの最長片道切符はどういうものか,ご存知の方は御一報下さい.

Wednesday, October 10, 2012

マルコフ行列の中の著者達: どの著者がもっとも人々に影響を与えたのか? (6)


数学では番号だけ考えていれば良いが,実際にはこの番号は何かに対応している.何かとは区別さえできればなんでもよいのであるが,具体的な例を示した方がイメージがわくと思うので,2つほど例を示そう.

例1:点が著者を示す場合

以下の4人の著者を考える.

  • William Shakespeare
  • Lewis Carroll
  • Raymond Smullyan
  • Martin Gardner

これらの著者がそれぞれどのような影響を受けたかは様々な見解があるだろうが,とりあえず,Shakespeare は Carroll に影響を与え,Carroll は Smullyan とGardner に影響を与えたとしよう.そのグラフは図 6のようになるだろう.

Figure 6. Graph example 1. Each node is an English author.

例2:点が駅名を示す場合

点が駅名を示す場合を考える.

  • Weinmeisterstr
  • Alexanderplatz
  • Hackescher Markt
  • Jannowitzbruecke

これらの駅間が隣あっている場合,それらの駅には関係があるものとして辺で接続しよう.するとそのグラフは図 7 のようになる.ところでBerlin の市内電車はよく工事をしていて,ある区間が不通であったり,困ったことに一方通行しかない場合があったりする.一方通行の場合には,グラフは有向グラフとなるであろう.
Figure 7. Graph example 2. Each node is a train station.

ここで一つ注意して欲しい.点の例が何であっても,これらのグラフは同じ形をしている.同じ形をしているグラフは全て数学では同じであって区別しない.つまり,図 4, 5, 6, 7 (図4,5 は以前の blog を参照)は全て同じグラフである.数学ではこの形のグラフから何が言えるかということに関して考える.

著者と駅名が同じ形で示されたということは奇妙なことかもしれない.しかし数学は世界からパターンをみつけだし,それらに共通のことを考える学問である.そして私はいつも驚くのであるが,この世界には同じようなパターンが出現することがいくつもある.新しく出会ったパターンが何かに似ている場合,同じ問題が隠されているかもしれない.その時,初めて出会った新しい問題にもかかわらず,その問題を数学で解くことができることもあるのだ.もちろん,いつも解けるとは限らないし,似ているのに実際は違っていたということもある.最も基本的なパターンは数えられるものに共通するパターンである.人も犬も街も星も教会も数えられるということでは同じものである.だから足し算を一度習えば,どんな数えられるものも足すことができるのである.グラフは関係だけを示す数学のパターンである.点と辺だけをみていると味気ないかもしれないが,これが,Web page 間の関係,人間関係,街と街を継ぐ道路,電話を継ぐネットワーク,石油のパイプライン,などに共通してでてくるパターンである.

では,この関係をどうやって記述すれば良いのだろうか.記述することは重要である.特に同じものを同じものとして記述できれば「比較」ということが可能になる.2つのグラフが同じかというのはまたちょっと難しい問題だが,ここでは点は一意に並べられるものとする.その場合,同じ関係を記述したら同じものになれば比較ができる.次回はそのような方法を述べよう.

Tuesday, October 9, 2012

How to get the camera parameters from OpenGL perspective/modelview matrix? (2)

モデルビュー行列からカメラのパラメータを計算する.

modelview matrix は model に作用する matrix も含んでいるので,シーンがスケール,移動,回転されていたりするとその効果が含まれている.これをmatrixの情報のみを使って分離する方法はない.そこで,ここではシーンに対するtransformation matrix は Identity matrix, つまり全てはカメラの効果としてカメラのパラメータを抽出する.もし,シーンに対する効果が分離されているようなレンダラと OpenGL レンダラを併用する場合には,この transformationmatrix を undo する必要がある.

modelview matrix は以下のような成分になっている.この成分の詳しいことに関しては私の以前の blog(http://shitohichiumaya.blogspot.de/2011/01/what-matrix-glulookat-generates-1.html)を参照のこと.
\begin{eqnarray*}
 \left[
  \begin{array}{cccc}
   x_x & y_x & z_x & 0 \\
   x_y & y_y & z_y & 0 \\
   x_z & y_z & z_z & 0 \\
   -(\vec{x} \cdot \vec{e}) &
    -(\vec{y} \cdot \vec{e}) &
    -(\vec{z} \cdot \vec{e}) & 1 \\
  \end{array}
 \right]
\end{eqnarray*}
4行を \(a,b,c\)で置きかえるとこの行列は以下のように書ける.
\begin{eqnarray*}
 \left[
  \begin{array}{cccc}
   x_x & y_x & z_x & 0 \\
   x_y & y_y & z_y & 0 \\
   x_z & y_z & z_z & 0 \\
   a   & b   & c   & 1 \\
  \end{array}
 \right]
\end{eqnarray*}
この行列には既にカメラの basis があるので,問題は視点の位置である.これは4行の1-3列だけを取り出して行列表示する.(ここでは行列とベクトルをカラムベクトル表示 \(\vec{x}, \vec{y}, \vec{z}, \vec{e}\)している.)
\begin{eqnarray*}
 \left[
  \begin{array}{ccc}
   & & \\
   \vec{x} & \vec{y} & \vec{z}\\
   & & \\
  \end{array}
 \right]
 \left[
  \begin{array}{c}
   \\
   \vec{e} \\
   \\
  \end{array}
 \right]
 & = &
 \left[
  \begin{array}{c}
   a \\
   b \\
   c \\
  \end{array}
 \right]
\end{eqnarray*}
であるから,
\begin{eqnarray*}
 \left[
  \begin{array}{c}
   \\
   \vec{e} \\
   \\
  \end{array}
 \right]
 & = &
 \left[
  \begin{array}{ccc}
   & & \\
   \vec{x} & \vec{y} & \vec{z}\\
   & & \\
  \end{array}
 \right]^{-1}
 \left[
  \begin{array}{c}
   a \\
   b \\
   c \\
  \end{array}
 \right]
\end{eqnarray*}
が視点の位置となる.従ってプログラムの実装例は以下のようになる(C++).

/// get camera parameters from
/// the OpenGL modelview matrix
///
/// \param[out] eyepos  eye position
/// \param[out] viewdir viewing directio
/// \param[out] updir   up direction vector
void gl_get_camera_parameters_from_modelview_matrix(
    Vector3d & eyepos,
    Vector3d & viewdir,
    Vector3d & updir)
{
    GLdouble mat[16];
    glGetDoublev(GL_MODELVIEW_MATRIX, mat);

    Vector3d xdir(0.0, 0.0, 0.0);
    Vector3d ydir(0.0, 0.0, 0.0);
    Vector3d zdir(0.0, 0.0, 0.0);

    xdir[0] = mat[0];  ydir[0] = mat[1];  zdir[0] = mat[2];
    xdir[1] = mat[4];  ydir[1] = mat[5];  zdir[1] = mat[6];
    xdir[2] = mat[8];  ydir[2] = mat[9];  zdir[2] = mat[10];

    // This is a, b, c components.
    Vector3d bvec(-mat[12], -mat[13], -mat[14]);

    Matrix33d basis_mat(xdir[0], xdir[1], xdir[2],
                        ydir[0], ydir[1], ydir[2],
                        zdir[0], zdir[1], zdir[2]);
    // This matrix should not be singular.
    // invert() gives matrix inverse.
    basis_mat.invert();

    eyepos  = basis_mat * bvec;
    viewdir = -zdir;
    updir   =  ydir;
}


ここでも簡単のためにエラーチェックをしていないが,個人的にはいくつかのassertion はしておきたい.

まとめ

ここでは,perspective 行列からカメラの perspective parameter を取り出す方法と,modelview 行列からカメラの位置や向きを取り出す方法について述べた.このようなことを計算することはあまりないかもしれないが,どうやって計算するのか興味のある人には面白いかもしれない.

How to get the camera parameters from OpenGL perspective/modelview matrix? (1)


Abstract

この話はコンピュータグラフィクスとプログラミングの話である.以前 gluLookat 関数の作る Matrix に関して書いた(http://shitohichiumaya.blogspot.de/2011/01/what-matrix-glulookat-generates-1.html).この記事を書いた動機はOpenGLとは独立したレンダラ(ray tracer など)を書いており,しかし OpenGL と画像をoverlay したいということがあったからである.そういう意味では特殊な話であるがこの記事は意外に参照されている.今回,我々の顧客が特殊なシステムを利用していて,この逆演算をする必要がでてきた.つまり OpenGL のProjection/Modelview matrix からカメラのパラメータを知りたいという要求である.これについてここに書いておく.

はじめに

世の中にはいろいろな可視化システムがあるもので,可視化システムのカメラのパラメータが不明,つまりカメラがどこにあるのかわからない,というものがあるようだ.話を聞くと実は OpenGL に依存した大規模なシステムで,あるソフトウェアのレイヤーではカメラのパラメータがわからないという状況のようだ.しかし,OpenGL を使っているので,projection matrix と model view matrix にはアクセスできるということである.そこでこれらの matrix からカメラのパラメータを抽出する方法というものを尋ねられた.

これはちょっと特殊な状況であると思うが,OpenGL の Projection/Modelviewmatrix が与えられた場合に,カメラのパラメータを抽出するにはどうするかというパズルとして考えれば面白いかもしれないと思い,ここに blog として記す.

パースペクティブ行列からカメラのパラメータを計算する.

まずは Perspective matrix に関連したパラメータを取りだそう.

gluPerspective 関数は以下のような関数である.詳しくは OpenGL の manual を参照のこと.

void gluPerspective(GLdouble fovyd,  GLdouble aspect,
                    GLdouble zNear, GLdouble zFar);

  • fovyd: field of view y (degree)
  • aspect: aspect 比
  • zNear: 視点に近い方のクリッピング平面までの距離
  • zFar:  視点から遠い方のクリッピング平面までの距離

OpenGL のマニュアルによると,このパラメータによって与えられる
perspective matrix は以下のようになる.
\begin{eqnarray*}
 \left[
  \begin{array}{cccc}
   \frac{fr}{asp}  & 0  & 0 & 0  \\
   0               & fr & 0 & 0  \\
   0               & 0  & \frac{z_f + z_n}{z_n - z_f} & -1 \\
   0               & 0  & \frac{2 z_f z_n}{z_n - z_f} & 0  \\
  \end{array}
  \right]
\end{eqnarray*}
ここで,\(fr = \frac{1}{\tan(\frac{fovy}{2})}\) であり,\(fovy\) はradian で示した \(fovyd\) である.\(asp\) は aspect ratio であり,\(z_n\) は視点に近い方のクリッピング平面までの距離,\(z_f\) は視点から遠い方のクリッピング平面までの距離である.perspective matrix の要素を次のように示すと,
\begin{eqnarray*}
 \left[
  \begin{array}{cccc}
   aa & 0  & 0  & 0  \\
   0  & bb & 0  & 0  \\
   0  & 0  & cc & -1 \\
   0  & 0  & dd & 0  \\
  \end{array}
  \right]
\end{eqnarray*}
\begin{eqnarray*}
 asp  &=& \frac{bb}{aa}\\
 fovy &=& 2 \arctan(\frac{1}{bb})\\
 z_f  &=&\frac{c-1}{c+1} z_n
\end{eqnarray*}
ここで,\(kk = \frac{c-1}{c+1}\)と置くと,
\begin{eqnarray*}
 z_n &=& \frac{dd(1-kk)}{2k}\\
 z_f &=& \frac{dd(1-kk)}{2}
\end{eqnarray*}
が得られる.ただし,\(z_n \neq 0\) とした.

これを元にした実装 (C++) は以下のようになる.

/// get perspective camera parameters
/// from the OpenGL projection matrix
///
/// \param[out] fovy_rad     field of view in radian
/// \param[out] aspect_ratio aspect ratio
/// \param[out] clip_min     clipping min distance
/// \param[out] clip_max     clipping max distance
void gl_get_camera_parameter_from_perspective_matrix(
    double & fovy_rad,
    double & aspect_ratio,
    double & clip_min,
    double & clip_max)
{
    GLdouble mat[16];
    glGetDoublev(GL_PROJECTION_MATRIX, mat);

    GLdouble const aa = mat[0];
    GLdouble const bb = mat[5];
    GLdouble const cc = mat[10];
    GLdouble const dd = mat[14];

    aspect_ratio = bb / aa;
    fovy_rad     = 2.0f * atan(1.0f / bb);

    GLdouble const kk = (cc - 1.0f) / (cc + 1.0f);
    clip_min = (dd * (1.0f - kk)) / (2.0f * kk);
    clip_max = kk * clip_min;
}

簡単のために division by zero のチェックは省いている.もちろん正しくOpenGL が設定されていればこれでも動くが,個人的には assert は入れておきたい.

Friday, October 5, 2012

0 とプラスの関係


8 歳の T 君が 7 について学んでいる.7 はどんな加算によって作られるかというものである.ここで利用している教材の一つは図 1 にあるような Zahlenhaus というものである.2つの部屋があり,それぞれに何人かがいる.全体で何人が一つの家にいますか? ということで加算というものを考えるものである.

Figure 1. Zahlenhaus: 7 = 4 + 3
Figure 2. Zahlenhaus: Questions
そこには,

  •  7 = 3 + ?
  •  7 = 4 + ?

のような問題(図2)があり,? を埋めるのである.T 君は上記の質問にはまったく問題なく答える.

しかし,次の質問がわからないという.

  • 7 = 0 + ?

ところで,学ぶ時,何がわからないのかをわかっているというのはとてもやりやすい.私自身,数学の本を読んでいてどこでわからなくなったのかをみつけるのに苦労することがある.そしてわからない部分がわかれば道が見えることが多い.どこがわからないかわからないと,どこまで本を戻ればいいのか見当がつかないからである.

さて,Zahlenhaus に戻ろう.左の部屋には誰もいない.右の部屋には何人いる?と尋ねると,2 + 5 人とか3 + 4 人という.私はこれには困った.つまり,T 君は

  • 7 = 0 + 2 + 5

と答えたのである,これは数学的にはまったく正しい.一つの部屋をまとめて数えなくてはいけない理由は特にないし,そういう仮定を明確に言ったわけではない.だいたい最初に 7 はいくつといくつ? というように聞いているのは練習の意味が強い.7 は 7 である.と答えて間違いはない.これはゲームだと思ってもらった方がいいかもしれない.

しかし,それぞれの部屋に一つだけ数字を割り当てるという暗黙のルールによって 7 = 0 + 2 + 5 は間違いとされる.Zahlenhaus には部屋が 3つないからである.もし,Zahlenhaus に部屋が 3つあれば,これは正しくなる.しかし,ある計算が部屋の数で間違いだったり正しかったりするのは逆に混乱するのではないだろうか?

私はよくチャールズ・ドジソン(ペンネーム: ルイス・キャロル)の本「鏡の国のアリス」を思い出す.

王とアリスが出あった時の会話である.(残念ながら日本語は no- の形を翻訳するのが難しいので英文をつけておく.)
「...その道を見よ,そしてわしに告げるがよい.もしそなたがどちらか片方でも見ることができたのであれば.」
「道には誰もいないのが見えるわ.(I see nobody on the road.)」アリスは答えました.
「そなたのような目が持てたならな.」王はいらだったような声で言いました.「『誰もいない』のが見えるとな!(To be able to see Nobody!)」 [Quote 1]
Zaehlenhause の片方の部屋には誰もいない.誰もいないということを記すにはどうするのか,人類は長い年月が必要だった.誰もいないということを 0 人がいると記すのは簡単なことなのだろうか.たとえば,0 人がいないということは,0人ではないということであるから誰かいることにはならないか.0 人はいる.ないということ(0人)がその部屋にある.もし,これが彼の問題であったらこれはなかなか難しい.しかし,話をしていくうちに 0 がそこにあるというのは問題ないらしいことがわかった.誰もいないということを記す方法として 0 を書くということがわかったのかどうかはわからないが,誰もいない部屋を機械的に 0 と呼ぶというのはよいらしい.

では一つの部屋の人数を一つの数字で示すということが難しいのだろうか.「普通」一つの数にまとめられるだけまとめるのは,後に比較をするのが簡単だからである.1 + 1 + 1 + 1 + 1 + 1 + 1 と 1 + 1 + 1 + 1 + 1 + 1 とを比較する時,7 と 6 を比較するのでは 7 と 6 を比較する方が人間には簡単である.だから0+ x を x と書くのが「普通」である.「普通」なのはある機能として簡単になるからである.「普通」というものと「正しい」というのは時に異なり,「正しい」かどうかの方が数学では重要である.

私は思った.またアリスだ.私はアリスほど数学的な本はないと友人に説明しようとしたことがある.
アリスが白の女王に尋ねられる.「足し算はできるかえ?一たす一たす一たす一たす一たす一たす一たす一たす一たす一たすは何だい?」
「わからないわ,何回一と言ったのか聞き逃したの」とアリスは答えました.
「この子は足し算ができないね.」赤の女王が割り込んで言いました.「引き算はどうだい.8から9を引いたらいくつかえ?」
「8から9は引けないわ.あたりまえでしょう」アリスはすばやく答えました.「でも...」
「この子は引き算もできないね.」白の女王が言いました.「割り算はどうかね.パンをナイフで割ったらどうなるかね?」 [Quote 2]
どんどん分解していくと,白の女王の問題になってしまう.しかし,ある意味白の女王は正しい.どんな自然数も1の足し算でできているのだ.それをなんとか説明できないものだろうかと思った.

私は T 君に,2 ユーロのアイスを買いたい時,1 ユーロ持っていたらいくら母親にねだるんかを尋ねた.すると,1 ユーロと答えた.では,1 ユーロも持っていなかったら? つまり 0 ユーロ持っている時は,彼は 2 ユーロと答えた.それでは 7 ユーロ欲しい時,なにも持っていない,つまり 0 ユーロ持っている時は?と尋ねるとやはりわからないという.では 1 ユーロ持っている時は,6 ユーロと答えた.つまり以下が T 君の理解である.

  • 2 = 1 + 1
  • 2 = 0 + 2
  • 7 = 1 + 6
  • 7 = 0 + ? わからない.

私は困ってしまった.T 君はどういうモデルでこのような答えを出すのだろうか?結局,

  • 1 ユーロ欲しい時,1 ユーロ欲しいと言う.わかった.
  • 2 ユーロ欲しい時,2 ユーロ欲しいと言う.わかった.
  • 3 ユーロ欲しい時,3 ユーロ欲しいと言う.わかった.
  • では,7 ユーロ欲しい時は? 7 ユーロ欲しいと言う.ついにわかった.

T 君はこの後いくつかの 0 + x の質問に正しく答えたのであるが,今回私は彼が何をわかっていなかったのかわからなかった.だからこれで良いのかちょっと心配である.ところで,もしいつか数学を仕事に使う時があなたに来たら,その時あなたは数字を使うことはないだろう.現在,数字はコンピュータのものである.

Quote 1

Quote 1 は Lewis Carroll による Through the Looking Glass and What Alice Found There. Chapter 7, The Lion and the Unicorn. から.また,テキストは The Annotated Alice, the definitive edition, edited by Martin Gardner, Penguin books, page 234, ISBN-13: 978-0140289299 による.翻訳は私のものであり,間違いがあればそれは私による.

Quote 2

Quote 2 は Lewis Carroll による Through the Looking Glass and What Alice Found There. Chapter 9, Queen Alice, Through the Looking Glass and What Alice Found There. から.また,テキストは The Annotated Alice the definitive edition, pp. 265-266, Penguin books, ISBN-13: 978-0140289299. による.翻訳は私のものであり,間違いがあればそれは私による.