Syoyo Fujita's Blog

raytracing monte carlo

Month: May, 2005

ERPT 実装 : 点を画像面へ射影する。

経路の変異を実装しているのであるが、はたと気づいたことがある。

コースティクス摂動では、光源からの方向を摂動させる。
これがたとえば LDE の経路だったばあい、
摂動した方向に光源からレイを飛ばして交点を求め、その交点から
視点を繋ぎ、またその交点 – 視点間のレイが画像面のどのピクセル位置
に対応するのか求めなければならない。

まあふつーに幾何を解いていけばよいのでしょうが、
なんだか計算するの面倒ですなー、と思っていたら、
David Cline の A practical 〜 には親切御丁寧にも
このための変換行列が付録に書かれていました。

こいつはグレートだぜ、、、
Cline 先生、、、すばらしすぎです、、、

という訳で、ワールド座標での点(視点から最初にヒットする点)を
画像面に射影する行列の実装は以下のようになるでしょうか。

// Create the matrix which transforms a point in world space onto
// its pixel coordinate.
void
matrxi_project_point_onto_the_image_plane(double M[16])
{
        double C[16], V[16], T[16], TC[16];
        vector_t eye;

        vcpy(&eye, &render->eyepos);
        memcpy(T, render->T, sizeof(double) * 16);

        double w = render->width;
        double h = render->height;

        double tw = 45.0 * M_PI / 180.0;        // horizontal field of view
        double th = 45.0 * M_PI / 180.0;        // vertical field of view

        C[ 0] =    1.0; C[ 1] =   0.0; C[ 2] =   0.0; C[ 3] = 0.0;
        C[ 4] =    0.0; C[ 5] =   1.0; C[ 6] =   0.0; C[ 7] = 0.0;
        C[ 8] =    0.0; C[ 9] =   0.0; C[10] =   1.0; C[11] = 0.0;
        C[12] = -eye.x; C[13] = eye.y; C[14] = eye.z; C[15] = 1.0;

        V[ 0] = w / 2.0 * tan(tw / 2.0);
        V[ 1] = 0.0;
        V[ 2] = 0.0;
        V[ 3] = 0.0;

        V[ 4] = 0.0;
        V[ 5] = h / 2.0 * tan(th / 2.0);
        V[ 6] = 0.0;
        V[ 7] = 0.0;

        V[ 8] = w / 2.0;
        V[ 9] = h / 2.0;
        V[10] = 1.0;
        V[11] = 1.0;

        V[12] = 0.0;
        V[13] = 0.0;
        V[14] = 0.0;
        V[15] = 0.0;

        // M = V . T . C
        mmulmat(TC, T, C);
        mmulmat(M, V, TC);
}

render->T はカメラ行列ですが、これは gluLookAt() で計算させれば簡単に求まります。

行列 M を点 p にかけてあげればピクセル座標が求まる(はず。まだ未検証)。

to be continued…

Advertisements

ERPT 実装: エネルギー再分布ステップの枠組みと摂動ルーチン

さて、そろそろ ERPT のコアとなる部分を実装していくことにしよう。

エネルギー再分布

パストレにより、まずはデポジットする e_d が求まったはずであるので、
あとはこれをおのおのの変異した経路の対応するピクセルに加算していけばよいはずだ。

e_d はパストレの結果画像バッファから引っこ抜いてきた値を e_ave として、
変異させる回数で割れば求まる。

e_d の計算のコードはこんな感じになるだろうか。

for (j = 0; j height; j++) {
        for (i = 0; i width; i++) {

                vector_t e_ave;
                vector_t e_d;

                // get average energy of this pixel from
                // pass1's path tracing result.
                e_ave.x = render->e_ave[3 * (i + j * render->width) + 0];
                e_ave.y = render->e_ave[3 * (i + j * render->width) + 1];
                e_ave.z = render->e_ave[3 * (i + j * render->width) + 2];

                // e_d is the deposition energy for each mutated
                // path Y.
                e_d.x = e_ave.x / (double)render->nmutations;
                e_d.y = e_ave.y / (double)render->nmutations;
                e_d.z = e_ave.z / (double)render->nmutations;

                energy_redistribution_sampling(i, j, e_d);
        }
}

render->e_ave はパストレの結果の画像バッファをコピーしたものである。
render->nmutation が変異の回数(式 7 の k)である。

energy_redistribution_sampling() は、図 7 相当の処理を行うルーチンであり、
以下のようになる。

void
energy_redistribution_sampling(int x, int y, vector_t e_d)
{
        bool     hit;
        int      k;
        double   u, v;
        vector_t dir;

        path_t   path;
        double   X_f;
        char     str[1024];

        if (fabs(e_d.x) nmutations; k++) {

                u = (double)x + randomMT();
                v = (double)(render->height - 1.0)
                  - (double)y + randomMT();

                get_eyedir(&dir, &render->eyepos, &render->eyedirpos,
                           &render->du, &render->dv, render->T, u, v);

                // create a path, x, int the current pixel.
                hit = path_trace(&path, &render->eyepos, &dir);

                // RGB -> Luminance
                X_f = CIEcol(path.radiance);

                // X_f(x) > 0.0
                if (hit && X_f > 0.0) {
                        equal_deposition_flow(&path,
                                              path.radiance,
                                              render->nchains,
                                              e_d);
                }
        }
}

まず、e_d の輝度値がゼロのとき、つまりピクセルが完全に遮蔽されているような領域だったり、
背景にヒットするような位置である場合は、処理を打ち切る。
(たまたまサンプル数が不十分で、光源へのシャドウレイが全部遮蔽する場合でも
e_d がゼロになる可能性があるが、、、まあそのときはしょうがないということにしよう)

e_d に輝度があれば、まず適当にパストレを行って初期の経路を生成する。
あとはこの経路が同じくどこかに遮蔽されたり背景にヒットするなどされずに
非ゼロな輝度を持つのであれば、 equal_deposition_flow() へと進む。

equal_deposition_flow() はこんな感じ。

void
equal_deposition_flow(path_t *x, vector_t e, int m, vector_t e_d)
{
        int    i, j;
        int    num_chains;
        double e_luminance, e_d_luminance;

        e_luminance    = CIEcol(e);
        e_d_luminance  = CIEcol(e_d);

        num_chains = (int)floor(randomMT() + e_luminance / (m * e_d_luminance));

        render->total_chains += num_chains * render->nchains;
        render->chain_events++;

        for (i = 0; i nchains; j++) {
                        // mutate_path(y, x);
                        // relative_path_density_(y, x);
                        // ごにょごにょ...
                }
        }
}

というわけで、ERPT(MLT でもそうだが) の核心のひとつである、
経路の変異処理、mutate_path()、を実装する部分までやってきた。

変異処理

ERPT では、変異はコースティクス摂動とレンズ摂動の二つしかない
(少なくとも論文ではこの二つしか使っていない)。
つまりはピクセルの位置や光源の方向をちょっとだけずらしてあげる。

また、各変異において、経路の長さは一定で変わらない。

MLT は bidirectional mutation により経路の長さが変わったりするので、
ERPT のほうが変異処理はかんたんである。

まずは、コースティクス摂動とレンズ摂動のために、レイの方向を摂動
(つまりちょっとだけずらす)させる必要があるので、このための基礎コード
を実装することにします。

ERPT でも、摂動処理の基本は Eric Veach の MLT の論文をもとにしている。

これを行うためのコードが、すでに ERPT の第一著者が大学のテクニカルレポートと
して出版している “A practical introduction to Metropolis Light Transport”
に付録として載っている。こちらはパストレを用いて MLT の入門的な実装を
解説している。ERPT と似ている部分も多く、
摂動部分についての解説は ERPT にそのまま流用できるはずだ。

David Cline and Parris K. Egbert,
“A Practical Introduction to Metropolis Light Transport”
BYU Technical Report, 2005
http://rivit.cs.byu.edu/a3dg/publications.php

というわけでこれをそのまま利用すればよい。簡単だ。

void
pixel_offset(double *x, double *y, double r1, double r2)
{
        double phi = randiomMT() * 2.0 * M_PI;
        double r   = r2 * exp(-log(r2 / r1) * randomMT());

        (*x) += r * cos(phi);
        (*y) += r * sin(phi);
}

void
angular_offset(vector_t *N, double theta1, double theta2)
{
        vector_t U, V;
        vector_t axisZ, axisY;

        axisZ.x = 0.0; axisZ.y = 0.0; axisZ.z = 1.0;
        axisY.x = 0.0; axisY.y = 1.0; axisY.z = 0.0;

        // Make a UVN coordinate system from N
        if (fabs(N->x) x = N->x + r * cos(phi) * U.x + r * sin(phi) * V.x;
        N->y = N->y + r * cos(phi) * U.y + r * sin(phi) * V.y;
        N->z = N->z + r * cos(phi) * U.z + r * sin(phi) * V.z;

        vnormalize(N);
}

pixel_offset() がレンズ摂動(画像面のサンプル位置 x, y をずらす)、
angular_offset() がコースティクス摂動(角度をちょっとずらす)となる。

randomMT() は [0.0, 1.0) の範囲で実数ランダム値を返す。

pixel_offset() の r1, r2 は摂動量を調整するパラメータである。
MLT の元論文では、r1 は 0.1、r2 は画像面の面積(横ピクセル数 x 縦ピクセル数)となっている。
A practical 〜 では r1 は 0.1、r2 は横ピクセル数の 10 % に設定している

ちなみに、MLT のレンズ部分経路変異(lens subpath mutation) に相当する処理を行う場合は、
A practical 〜 では pixel_offset() により大きな値を設定、r1 は 1.0、r2 は横ピクセル数の 25 %、
することで代用しているとのこと。

angular_offset() の theta1, theta2 は、MLT も A practical 〜 でも、
theta1 は 0.0001、theta2 は 0.1 に設定しているとのこと。

では、実際に経路を摂動させてみよう…

to be continued…

ERPT 実装その1

ERPTを実装中。

とりあえずはいろいろコードを流用して saku saku っとパストレ部分を完成。
まずはディフューズ面のみ。

絵は 36 samples/pixel(=6×6 stratified sampling)。
まだまだノイズは多いですね。
ピクセルフィルタとか最適化とかまったく行っていないからこんなものでしょうか。

あとは mutation 処理と path density の計算をやれば ERPT の完成。
このまま問題なく行ければ結構簡単やも。

経路のデータ構造はこんな感じで作っておくとよいだろうか。
(ほとんど Graphics Programming Methods の MLT コードと同じだけど)

/* path vertex information */
typedef struct _path_vertex_t
{
        vector_t p;             // position;
        vector_t indir;         // Incident direction
        vector_t outidr;        // Outgoing direction
        int      type;          // type of vertex. 'E', 'D', 'S' or 'L'
        vector_t bsdf;          // BSDF(RGB)
        double   G;             // Geometric coefficient
} path_vertex_t;

/* record of light path */
typedef struct _path_t
{
        path_vertex_t verts[MAX_PATH_VERTICES]; // Path vertices
        int           length;                   // The number of the vertices

        int           px, py;                   // Pixel position

        vector_t      radiance;                 // Eval of this path(RGB).

} path_t;

ERPT のコードはこちらにおいてあります。

http://lucille.atso-net.jp/svn/angelina/lighttransport/erpt/

コンパイルには fltk と boost が必要です。
Mac OS X(メイン環境), linux, Windows(cygwin) で動くはず。

Energy Redistribution Path Tracing

David Cline, Justin Talbot and Parris K. Egbert
Energy Redistribution Path Tracing
ACM Transactions on Graphics (SIGGRAPH 2005 proceedings)
http://rivit.cs.byu.edu/a3dg/publications.php

メトロポリス光輸送(MLT)のような考え方を、
パストレーシングのコンテキストに導入した手法。
略して ERPT。

メトロポリス光輸送も経路のサンプリングのために
パストレーシングなどを用いますが、
こちらはメトロポリス法の考えをパストレーシングの
中に取り込んでいます。
(論文では両者を組み合わせたハイブリッドな手法と言っています)

それにより、パストレーシングの持つ利点を生かせ、
逆にメトロポリス光輸送の持つ欠点を克服することが
できるようになっています。

本手法の利点は、

  • スタートアップバイアスがない
  • ピクセル単位で層別化しやすい
    (メトロポリス光輸送は画像面で大域的に動くので、
    ピクセル単位でコントロールしずらい)
  • パストレ自体がエルゴード性を満たすので(ちょっとこれは実際のシーンでは疑問だけど…)、
    MLT のように複雑な経路の変異(mutation)を行わなくてもよい。

です。

基本がパストレなので、MLT よりも実装は容易となっています。

今回じっくり論文を読み、大体しくみがわかりました。
これなら実装できそうですし、実際に正しく動くのか検証もしてみたいので、
SBR 2005 に ERPT を引っさげて本格出走します。
(まだ公開はされていない wavelet importance sampling も実装してみたいですが) 

ERPT のより詳細な論文解説は SBR 2005 で。


http://lucille.atso-net.jp/sbr2005/index.php?Energy%20Redistribution%20Path%20Tracing

ATI XENOS

XBOX 360 の GPU ですが、ブロック図がちょっと公開されています。

また、48 シェーダ演算器というのは、16 ピクセルパイプ相当ではなく、
真に 4 要素ベクトル演算器 が 48 本(48 ピクセルパイプ相当)が正しいとのこと。

http://fcj.s18.xrea.com:8080/modules/news/article.php?storyid=1451
http://www.firingsquad.com/features/xbox_360_interview/

ひょえー、す、すごいなぁ、、、
現在の GPU が 16 ピクセルパイプ + 6 頂点パイプぐらいですから、
優に 2 倍ですか….

ということで、48 G shader op/sec の shader op は 4 fp 演算
が正しそうです。となるとやはりシェーダ内部は 1 GHz クロックの
倍速で動いているのかな?

48 * 4 * 2(fmadd) = 384 Gflops ???

さらに、うまくプログラムすれば 1 クロックで
2 シェーダを同時実行できるとのこと。うーん、どういう構成なんだろ。
NV40 みたいに fp 縱積み 2段ということなのでしょうか?
(この場合は fmadd がなくて、2 段(fadd + fmul)を通って fmadd 相当を行うのかな)

PS3

PS3 もでましたね。

こ、こんどは 2Tflops ですか、、、
GPU である RSX に関しては、シェーダパイプ数は明らかに
まだされていないようですが、1.8T flops の性能があるとのこと。
ひょえ~、、、ってありえねーよ。

XBOX 360 の 1 T flops もそうだけど、
やっぱり ROP での比較・ブレンドやテクスチャのフィルタリングに加えて、
ラスタライザでの頂点パラメータの補間とかも flops に入れての数字の
ようです。

RSX は Geforce 6800 Ultra の二倍程度の性能だそうです。
GPGPU 関連だったかで nvidia が fmadd のシェーダで 40 – 50 Gflops
の実測値を得たという
スライドを見た気がします(実際 gpubench の結果もそれくらいでした)。
なので、RSX の実効値は 80 – 100 Gflops 程度と見るのが妥当でしょうか。

FF7 の技術デモ(たぶんリアルタイムだよね?)も発表で流れていたようです。

http://media.ps3.ign.com/articles/614/614940/vids_1.html

いやー、時代の流れを感じました。すごいなー、今じゃリアルタイムで動いちゃう
んですね。

ちなみに、RSX は 22.4 GB/s(read+write) しかフレームバッファとの帯域がありません。
(XBOX 360 は 256 GB/s (read+write) のスペックになっていますが、これは
4x MSAA の 4 倍展開時の実効値なので、実際は 64 GB/s(read+write)になります)。
その代わり eDRAM の少ないメモリではなく GDDR3 256 MB もの大きなメモリ
を扱うことができる利点があります。

PS3 のゲームは、PS2 みたいにマルチパスで何度も書き込みんで
ディテールを出すのではなく、1 回のシェーダでなるべく全部演算し
てディテールを出していくスタイルになりそうです。

とはいえ、これからの時代、ゲームでリアリティを出すには、水とかの透明感がカギと
なってくると思うのですが、これは基本的に複数描画+アルファブレンドと
なり帯域を多く消費する傾向になるので問題になりそうな感じがしますが、、、

まあでもそこらへんはやはり杞憂で、プログラマブルシェーダのチカラで
なんでもきっと解決のでしょう。

Digital baby

先月くらいの CG world でも紹介されていましたが、
映画「レモニースケットの世にも不幸せな物語」では、
いくつかのショットで実写とみまちがうほどの
赤ちゃん(サニー)の CG レンダリングが行われています。


http://www.panoscan.com/PanoPress/2005Press/Lemony/LemonySnicket.html

サニーのレンダリングには、CG world に書いてあるように、
大域照明+サブサーフェススキャタリング+HDRI が使われています。
(放射照度キャッシュも使われたみたい)

特に、サニーの顔は、サブサーフェスの計算を行うために、
ポイントでレンダリングしたそうです。
ちなみにレンダラは PRMan のようです。

映画館で見てきましたが、マジでぜんぶ本物だと思いました。
…ってかぜんぜん CG だってわかんねーよ!

これも大域照明とサブサーフェススキャタリングのなせる業ということでしょうか。

ついでに、パペットモデルも使われたようです。

http://dvdtimes.co.uk/content.php?contentid=56902

XBOX 360

XBOX 360 が公表されました。

http://www.itmedia.co.jp/games/articles/0505/13/news012.html

おどろいたのは、システム全体で 1T flops という性能。
うーん、すごいなぁ、けど、そんな数字ってありえるのかな?

というわけで少し考えてみました。

まず、CPU は、9 G 内積演算で 115.2 Gflops です。

3 component(3 fmadd) の内積専用命令があると考えれば、

3.2(GHz) * 3(プロセッサ) * 4(SIMD) * 3(fmadd)  = 115.2 G
flops

と計算どおりです。

んで、 GPU は、48 way のシェーダパイプとありますが、
これは今の GPU で用いられているピクセルパイプで言えば 12 本相当
という意味でしょう(12 * 4 = 48)。

統一シェーダ化により、4 要素固定のピクセルパイプではなく、
より柔軟な 3 要素のベクトル演算とか、2 要素の
ベクトル演算を動的に構成していくアプローチを取っていそうです。
そのため、48 way という表記を
用いているのではないでしょうか。

シェーダパイプは秒間 48 G シェーダ命令(スカラ命令?)を
実行できるようです。

GPU はコアクロック 500 MHz ですが、シェーダパイプは
倍速の 1 GHz で動いるのでしょう。

NV みたいに、1 シェーダパイプ内に縦積みで 2 fp 演算器があるとすると、

48(inst) * 1(GHz) * 2(fp ALUs) * 2(fmadd) = 192
Gflops

CPU とのトータルで 1 T にはまだまだですね。

1 シェーダパイプ内に 3 fp, 4 fp 演算器があるとすると、

48(inst) * 1(GHz) * 3(fp ALUs) * 2(fmadd) = 288
Gflops
48(inst) * 1(GHz) * 4(fp ALUs) * 2(fmadd) = 384
Gflops

うーん、やっぱりまだ足りない。

というわけで、どうもテクスチャユニットでの fp テクスチャを引くときの
フィルタリングとか、ROP(Raster OPeration unit)での fp ブレンドとかの
flops も考慮に入れて、トータルで 1 T って言っているのかなぁと思います。

EUROGRAPHICS 2005 Papers

EUROGRAPHICS 2005 の論文リストが Tim Rowley 大先生のところで公開されています。

http://myweb.hinet.net/home7/hks/Papers2005/eg2005Papers.htm

Kinuwaki 先生の放射照度キャッシュの論文が採択されています。

ここ最近は、放射照度キャッシュやらファイナルギャザーやらアンビエントオクルージョンなどの、
半球上からの寄与率の計算を効率的に行う手法の提案が活発に行われている感じが
します。

○○○○はファーストしか認めん!?

、、、いやあ、壮絶なるラストの座の争奪戦になってもおかしくないところでした、、、

Graphics Hardware 2005 に出した論文の採択通知が来ました。
三人の共著ですが、もちろん私は第一著者ではなく、ラストオーサーとなる予定です。
(それに、実際元になったアイデアは私の発案でないので)

仕事としての論文ですので、pdf をすぐに公開できそうにないのが
残念ですが、いずれなんらかの形でここのページなどで解説できればと思います。

どのようなものかというと、

  • 2 のべき乗、これは美学
  • 世界でもっとも美しい並びのパターンを作った
  • それが本論文の寄与(contribution)だ

という感じです。

すでに完成されたハードウェアネタは内部ではいっぱいあるのですが、
利権などもあり、なかなか外に出す機会がないんですよね。

もったいないので、これからもコンスタントに
出していけたらいいなぁと思います。

グランドスラム達成に一歩近づいたかな?

グランドスラムとは?

SIGGRAPH, EUROGRAPHICS, EGSR, GH, PACIFIC GRAPHICS, とあとなんか 2 学会くらい、

と世界最高峰の 7 大グラフィックス学会すべてに論文を通すことを言う。。。

過去、日本人でグランドスラムを達成したのは、、、、いるのかな?
Nishita 先生と Dobashi 大教授くらい?

○○○○はファーストしか認めない!?

ちなみに、○○○○はファーストしか認めない、という風潮がありますね。

○ン○○とか、○ン○○とか、、、、

でも、わたしはラストが大好きです。