Comments
Description
Transcript
グローバルイルミネーション入門
グローバルイルミネーション入門
藤田将洋 Syoyo Fujita
更新 平成 16 年 9 月 11 日
注意
1
本書は未完であり、また書き始めたばかりである。そのためまだ途中までしか書いていなかった
り、構成がよくなかったり、ところどころ間違った記述が修正されずに残っていることに注意して
ほしい。
導入
2
本文では、グローバルイルミネーションの入門から、より詳しいアルゴリズムの解説までを行な
い、グローバルイルミネーションの包括的かつ技術的な知識を与えることを目的としている。
グローバルイルミネーション (global illumination)1 は、三次元コンピュータグラフィックスにお
ける、現在主流かつ研究がさかんなレンダリングアルゴリズムの総称もしくは研究領域である。グ
ローバルイルミネーションの特徴は、写実的な CG 画像を作成したり、物理学的に正確な光のふ
るまいによるシミュレーションによる CG 画像を生成することが可能なことである。
2.1
グローバルイルミネーションとは?
図 2 に、グローバルイルミネーションを用いないで直接光のみでレンダリングした画像と、グ
ローバルイルミネーションを用いて間接的な照明も考慮してレンダリングした画像の比較を示す。
グローバルイルミネーションの威力を理解するのは明らかだと思う。
基本的なグローバルイルミネーションのアルゴリズムの研究は二十年以上前から行なわれてきて
いるが、近年とくにこの分野に注目が集まってきている。その理由は、もともとグローバルイルミ
ネーションは膨大な計算量が必要であったが、コンピュータの処理性能の向上や、新しい高速なア
ルゴリズムの提案により、ゲームなどのインタラクティブグラフィックスにおいてグローバルイル
ミネーションを用いることができるようになってきていることが挙げられる。多くのゲームは映像
のリアリティを求めており、グローバルイルミネーションはリアリティを与える手法であるため、
両者が結びつく。
1 大域照明、広域照明とも訳される
1
(a)
(b)
図 1: (a) グローバルイルミネーションなし。 (b) グローバルイルミネーションあり。シーンモデ
ル (Sponza Atrium) の著作権は Marco Davoic に帰属する
映画や TV などの CG 映像の制作に時間をかけることができる分野では、すでにその多くにグ
ローバルイルミネーションが用いられており、リアリティの向上に貢献している2 。
グローバルイルミネーション
ローカルイルミネーション
物理シミュレーション
知覚的、経験則
遅い
非常に高速
やや困難
容易
アルゴリズムの根幹
処理速度
アートディレクション
表 1: グローバルイルミネーションとローカルイルミネーションの比較
ただし、グローバルイルミネーションには、計算量以外にも大きな問題がある。それは物理的な
正確性に基づいてアルゴリズムが記述されているので、アルゴリズムは数学的、物理学的に複雑
で、深い知識がないと理解するのが困難であるということである。
本文は、より技術的なグローバルイルミネーションの解説、とくにグローバルイルミネーション
のアルゴリズムの解説を行なう。また同時に、それらを理解するために必要な数学的、物理学的知
識についても解説したいと思う。
本文の構成は以下のようになっている。まず最初は、グローバルイルミネーションアルゴリズム
の歴史を見ていく。これはグローバルイルミネーションのアルゴリズムの論文が、 SIGGRAPH な
どのコンピュータグラフィックスの学会で発表された年度をもとにしている。その後、年代順に提
案されたグローバルイルミネーションのアルゴリズムの詳細を解説する。次に、グローバルイルミ
ネーションの計算を高速化することが可能な、より少ないサンプル数で収束が早い準モンテカルロ
2 ただしフル CG 映画などでは、それでもまだすべてをグローバルイルミネーションでレンダリングするには時間がか
かりすぎるので、ローカルイルミネーションが用いられているのがほとんどである
グローバルイルミネーション入門 2
法について解説する。
2.2
グローバルイルミネーションのススメ
さて、本書ではグローバルイルミネーションの技術に興味があるのだけれども、そのまえにグ
ローバルイルミネーション、ひいてはコンピュータグラフィックスのレンダリングアルゴリズムと
いう分野のすばらしさを 3 点にまとめて紹介し、広く本分野へ興味を引かせたいと思う。
まず、レンダリングは取り組みやすい。結果が画像として現れるので、自分のやった作業の結果
がわかりやすい。そしてとくにグローバルイルミネーションは写実的できらびやかな画像を生成で
きるので、取り組むモチベーションが湧いてくる。そしてこれらはプログラムを書くだけで生成で
きるのである! 大規模な実験機器や特殊なグラフィックスハードウェアや、ソフトウェアなどは必
要ない。そのため最初の理論を理解するだけであれば、レンダリング時間にはそれほど気をかける
ことはないだろうから、ノート PC でも趣味感覚でゆっくりお気軽に行なうことができる。
また、レンダリングはシミュレーションであるので客観的な仕事である。キャラクタモデルの
制作やデザインの仕事は個人の美的センスが問われるが、レンダリングはそのようなことは無い。
100 人がレンダラを書けば、だいたいが同じ結果の画像が生成されるのが期待される。レンダラに
主観的な要素を取り入れるのはナンセンスである。レンダラはデータに客観的であるべきである
し、それを超えて画像に主観性を与えるのは、アーティストに委ねられるべきである。よって、レ
ンダラの制作において、レンダラ書き自身の美的センスという才能的な要素を問われることはない
ので、知識さえあればだれでもレンダラの制作を行なう事ができる3 。
そして最後に、レンダリングは学問的に興味深くそして奥深い分野である。高校時代に、いまま
での学校での教育は将来役に立つのか? という疑問は誰しも抱くものであろう。レンダラが生成す
るきらびやかな画像の裏側は、数学のオンパレードである。とくに中学・高校で習う (たぶん習う
であろう) 代数・ベクトル・行列がそのままでてくる。ベクトルの内積という訳の分からなかったも
のがシェーディングの中核で用いられ、きれいな絵と結びつく。またグローバルイルミネーション
ではさらに物理の知識も含まれてくる。レンダリングは、高校時代に感じた疑問に応えてくれる貴
重な分野の一つであろう。コンピュータグラフィックスに深く関わるようになると、コンピュータ
グラフィックスは物理学、統計学、金融工学などいろいろな分野からの応用で成り立っているとい
うことが判るかと思う。そのためコンピュータグラフィックスの研究では幅広い分野の研究を知っ
ておく必要があり、そのたびに知的探究心が刺激される。
2.3
レンダリングとは?
レンダリングは、現実世界におけるカメラとシーン (世界、環境) を、コンピュータの中でシミュ
レートすることである。こんな説明ができるだろう。
真っ暗で何もない世界がある。世界を照らすには、まずライトが必要だ。ライトを配置
3 近年はアニメ調や鉛筆画風などのアーティステックなレンダリング画像を生成するノンフォトリアリスティックレンダ
リングというアルゴリズムの研究がさかんであるが、これはそのアート的な処理を一連の手続きで与えるという意味でや
はり客観的なシミュレーションとして考えることができよう。
グローバルイルミネーション入門 3
しよう。そして、ライトで照らされる物体が必要だ。机やいすや壁などを配置しよう。
そして、カメラを配置し、シャッターを押して世界をフィルムに保存する。このシャッ
ターが押された時にフィルムに飛び込んでくる光を考えてみよう。
まず、光はライトから発せられる。その発生のしかたを考えねばならない。光が物体
に当たったら、その反射のふるまいを考えねばならない。金属のようにするどく反射
する場合もあるだろうし、ざらざらした紙のように均一に反射する場合もあろう。ビ
ロードの服だったら、ある方向の向きに反射の特徴が表れるだろう。鏡があったら映
り込みが生じるし、ガラスがあったら光は透過・屈折していく。光は物体に遮られるか
もしれない。物体に遮られるかどうかの判定も必要だ。最終的に光は何回か反射して、
そのわずかないくつかはカメラに飛び込んでくるだろう。これをフィルムに記録する。
うぉぉぉぉぉ、これだけでも膨大な計算量になるのは想像に難くないッ!!!。しかし、だッ!!
実は、光源は隣の部屋で点いていて、カメラは隣のライトのない真っ暗闇の部屋にあ
るかもしれない。このときは、真っ黒な画像ができあがる。明らかに今までの計算は
無駄であった! もしカメラがライトで照らされていない方向を向いているのであるのを
知っていれば、何も全部計算しなくてもよかったのに...
こんなことが無いように、特にグローバルイルミネーションでは、いかにカメラに入射してくる
重要な光の経路のみを見つけ出すか、いかに物理的に正確でありつつも無駄な計算をはしょって高
速化を図るか、が研究の中心である。
アルゴリズムの歴史
3
主要なグローバルイルミネーションアルゴリズムの年表と遷移を図 ?? に示す。アルゴリズムが
発表されたときのそのアルゴリズムの形態を、大きくオフラインとリアルタイムで分けた。オフラ
インはオフライン向け (映画などのようにレンダリングに時間をかける分野) のアルゴリズム、リ
アルタイムはリアルタイム向け (ゲームや GPU で処理できる分野) のアルゴリズムを示している。
近年はリアルタイム向けのグローバルイルミネーションアルゴリズムの提案が行なわれていること
が見て取れると思う。図中の矢印は、アルゴリズムの発展関係を示している。たとえば、分散レイ
トレーシングは、レイトレーシングを基としておりその発展版であるため、レイトレーシングから
分散レイトレーシングへと矢印を引いてある。
3.0.1
レイトレーシング法
最初のグローバルイルミネーションアルゴリズムは、レイトレーシング法 (ray tracing)4 である。
最初のレイトレーシング法は Appel ら [2] により提案された5 。これは再帰の無い一次レイ (視点
から発せられるレイ) のみ、かつ直接光の影のみを考慮していた。その後 Whitted [24] はこのレ
イキャスティング法を再帰的に適用し、反射・屈折透過を表現した6 。Whitted は [24] で global
4 光線追跡法とも訳される
5 Appel
らのレイトレーシング法は、現在ではレイキャスティング法 (ray casting) として知られている。
により提案されたレイトレーシング法は、その後の各種のレイトレーシング法の拡張手法と区別するため、
現在ではクラシックレイトレーシング (classic ray tracing) 法とも呼ばれている。
6 Whitted
グローバルイルミネーション入門 4
オフライン
リアルタイム
レイトレーシング
1980
Whitted1979
分散レイトレーシング
パストレーシング
Cook1984
Kajiya1984
ラジオシティ
Cohen1984,Nishita1985
1990
双方向パストレーシング
Lafortune1993,Veach1994
フォトンマッピング
Jensen1993
メトロポリス光輸送
Veach1997
2000
PRT
Sloan2002
ウェーブレット
ライティング
Ng2003
図 2: 主要なグローバルイルミネーションアルゴリズムの年表と遷移
illumination という言葉を用いてレイトレーシング法を説明している7 。
Whitted のレイトレーシング法の提唱により、ローカルイルミネーションからグローバルイル
ミネーションへと方法論が進化した。Whitted のレイトレーシング法は、今までは分けて考えて
いた遮蔽 (occlusion)・反射屈折効果・影 (shadow) を統一的に表現することを示したという点でも
重要である。
3.0.2
分散レイトレーシング法
Whitted のレイトレーシング法では、各レイがサーフェスで跳ね返るときの分岐は最大でも、反
射と、屈折もしくは透過の 2 本であるために、完全な鏡面反射 (perfect specular reflection)・完全
7 私が知る限りでは、global
illumination という単語が出てきたのはこの Whitted の論文が最初である。
グローバルイルミネーション入門 5
Computer Graphics
Volume18, Number 3
July 1984
Figure 8. 1984.
図 3: 分散レイトレーシングによる画像。[7] より引用。
な屈折 (perfect refraction)・ハードシャドウ (hard shadow) などのくっきりしたレンダリング画像
しか表現できなかった。そこで Cook ら [7] は、1 本のレイがサーフェスで跳ね返るときに複数の
レイを分散して新しく飛ばすことで、ぼやけた反射や屈折 (blurred reflection and refraction)・ソフ
トシャドウ (soft shadow) を表現することができる分散レイトレーシング (distribution ray tracing)
法が提案された8 (図 3 参照)。
145
またこの手法では、被写界深度 (depth of field)・モーションブラー (motion blur) などのカメラ
効果も統一した枠組みで表現することが可能になっている。
通常、分散レイトレーシング法と以下で述べるパストレーシング法はモンテカルロレイトレーシ
ング法 (Monte Carlo ray tracing) とも呼ばれることが多い。
3.0.3
ラジオシティ法
ラジオシティ(Radiosity) 法は、Goral ら [10] および Nishita ら [15] により提案された。ラジオシ
ティ法では、伝熱工学の理論を元に、光をエネルギーとして捉える。シーン内の物体をパッチ (微
小面積) 要素へと分割し、それらパッチ間でエネルギー伝達を行うことにより光輸送の計算を行う。
8 オリジナルの Cook らの論文の題名は Distributed Ray Tracing であるが、後に複数のコンピュータで並列にレイト
レーシング処理を行うという意味での distributed ray tracing と区別するために、現在では Cook らの手法は distribution
ray tracing と呼んで区別されている。
グローバルイルミネーション入門 6
ラジオシティ法では、レイトレーシング法では表現が難しい柔らかな間接光の表現が可能になっ
ているため、建築物のレンダリングや照明シミュレーションなどの分野で広く使用されている。
3.0.4
パストレーシング法
パストレーシング法 (path tracing)9 は、Kajiya[12] によりレンダリング方程式の提案と同時にそ
のレンダリング方程式の解決として提案された手法である。パストレーシング法では、レンダリン
グ方程式をモンテカルロランダムウォークで解くことにより、各レイがサーフェスで跳ね返るとき
には 1 本しか反射レイをしないレイトレーシングである。そのため分散レイトレーシングのよう
に各レイの跳ね返りで指数関数的にレイの本数が増えていくことを抑えることができる。
3.0.5
メトロポリス光輸送
メトロポリス光輸送 (Metropolis light transport, MLT) は、都市の名前ではありません。... 現在
では、パストレーシング法は多くの商用レンダラのグローバルイルミネーションアルゴリズムとし
て用いられている。
3.0.6
前計算放射輝度伝達 (PRT) 法
前計算放射輝度伝達法 (Precomputed Radiance Transfer, PRT 法) は、Sloan [20] らにより提案
された手法である。PRT 法では、放射輝度の伝達をあらかじめシミュレーションしておき、また
同時に球面調和関数 (spherical harmonics) でそれら伝達の結果を展開しておくことで、前述の手法
とは異なりインタラクティブに変化するシーンをレンダリングすることができる。
PRT 法では、低周波の照明分布において遮蔽も考慮された正確な照明積分を表現する。
3.1
照明アルゴリズムのまとめ
’70 年代、 ’80 年代初頭で提案されてきた、ローカルイルミネーションモデル、レイトレーシ
ング法、ラジオシティ法などの各種レンダリングアルゴリズムは、Kajiya のレンダリング方程式
[12](付録 A 参照) により一旦まとめられ、それらが統一された概念で記述できることが示される
ことになる。
グローバルイルミネーションの研究は、Kajiya のレンダリング方程式とその後に提案された光
輸送方程式 [3, 23] を基礎とした物理学的な観点から ’90 年代を通じて発展してきた。
グローバルイルミネーションアルゴリズムが、既存のローカルイルミネーションアルゴリズムと
区別される点をまとめると、以下の 3 点に集約される。
• 間接照明 (indirect illumination)
間接照明は、直接サーフェスに入射してくるのではなく、壁などに何度か反射してからサー
フェスに入射してくる光のことである。図 2b を見ても分るように、間接照明はリアルな画
像生成において非常に重要な要素である。
9 経路追跡法とも訳される
グローバルイルミネーション入門 7
• 遮蔽 (occlusion)
遮蔽 (および影) は、光が壁や物体などに遮られて生じる現象である。遮蔽は、遮蔽の無いレ
ンダリング結果と比べて、視覚的に深みを与え、物体の存在感を際立たせる要素として重要
である。遮蔽の計算は解析的に行うのは非常に困難であるので、通常はモンテカルロレイト
レーシングなどのモンテカルロ法を用いて確率論的に計算することが多い。また物体が動く
動的なシーンでは前もって遮蔽の変化を正確に計算することは困難である。正確な遮蔽の計
算は現在のチャレンジングな研究の 1 つとなっている。
• 物理学に基づく照明計算
正確なシミュレーションによりレンダリングを行うには、物理学に基づいて単位を取り決め
て照明計算を行うことが必要となる。グローバルイルミネーションの研究はこの物理学に基
づく照明計算と共に発展してきた [?]。
プログラミングツールと準備
4
本書でのレンダリングアルゴリズムの実装では、フリーソフトウェアのプログラミングツールを
使用します。レンダラは本質的にはデータを受け取って画像を生成するプログラムなので、コマン
ドライン環境を前提として、OS やウィンドウシステムには極力依存しないようにします。UNIX
系 OS(Linux, FreeBSD, Solaris) や Mac OS X にはコマンドライン環境が提供されていますので問
題はありませんが、Windows にはありません。Windows を開発環境に使用されたい人は cygwin
をインストールされることをお勧めします。cygwin は Windows 上で UNIX 環境を提供するソフ
トウェア群で、gcc などの開発環境も簡単にインストールできます。
本書のプログラムのコンパイルには gcc(GNU Compiler Collection) を推奨します。gcc はほぼす
べての OS に移植されているコンパイラです。
本書のプログラムはこの gcc でコンパイルの確認をしています。しかし基本的にコンパイラ依存
のコードは書かないようにしているので、他のコンパイラでのコンパイルには大きな問題は無いと
思います。
4.1
画像の保存
レンダラの仕事は、シーンデータをファイルやモデリングソフトから受け取って画像を生成する
ことです。まずは画像だけを保存するプログラムを作成します。ここで PNG や Jpeg などの形式
のフォーマットで画像を保存するプログラムを自分で 1 から組むのはとても大変です10 。そこでこ
のチュートリアルでは、プログラムが簡易なフォーマットである PPM(Portable PixMap) を採用
することにする。PPM は簡単なアスキー形式のフォーマット (バイナリフォーマットのバージョ
ンもある) である。PPM は読み書きのプログラムを簡単に作成することができるが、欠点として
ファイルサイズが巨大 (数 M バイト程度) になってしまうので注意して欲しい。PPM はだいたい
どのような画像ビューワでも見ることが可能である11 。以下は PPM で画像を保存するサンプルプ
10 libpng(URL をここに) や libjpeg(URL をここに) などのフリーの画像フォーマットライブラリを利用する手もある
がここでは取り扱わない
11 たとえば、Windows だと IrfanView(URL をここに) などが挙げられます
グローバルイルミネーション入門 8
ログラムです。
/*
* ppmsave.c
*/
#include <stdio.h>
int
main(int argc, char **argv)
{
FILE
int
int
*fp;
i, j;
width = 256;
int
const char filename[] =
height = 256;
"image.ppm";
fp = fopen(filename, "w");
if (fp == NULL) {
printf("Cannot create file : %s\n", filename);
exit(-1);
}
fprintf(fp, "P3\n");
/* magic number
fprintf(fp, "%d %d\n", width, height);
fprintf(fp, "255\n");
/* max pixel value
*/
*/
for (j = 0; j < height; j++) {
for (i = 0; i < width; i++) {
fprintf(fp, "%d %d 255 ", i, j);
}
}
fclose(fp);
}
このプログラムを実行すると、image.ppm というファイルが作成され、そのファイルを画像ビュー
ワで開くと図 4 のような画像が確認できる。
グローバルイルミネーション入門 9
図 4: PPM 画像出力例
4.2
OpenGL による PPM ビューワ
OpenGL になじみのある読者は、PPM 画像の表示は画像ビューワを利用しなくとも以下の
OpenGL(+GLUT) プログラムで表示を行なうことも可能である。
#if defined(__APPLE__) && defined(__MACH__)
#include <GLUT/glut.h>
#else
#include <GL/glut.h>
#endif
#include <stdio.h>
#include <stdlib.h>
static int
static GLubyte
width, height;
*image;
void
display()
{
glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
glClear(GL_COLOR_BUFFER_BIT);
glRasterPos2i(0, 0);
glDrawPixels(width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
グローバルイルミネーション入門 10
glFlush();
glutSwapBuffers();
}
void
key(unsigned char key, int x, int y)
{
if (key == 27) exit(0); /* ESC key */
}
void
reshape(int w, int h)
{
glViewport(0, 0, (GLint)w, (GLint)h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluOrtho2D(0.0, 1.0, 0.0, 1.0);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
}
int
main(int argc, char **argv)
{
FILE
*fp;
char
int
int
buf;
i, j;
maxval;
int
r, g, b;
if (argc < 2) {
printf("Usage: %s filename.ppm\n", argv[0]);
exit(-1);
}
fp = fopen(argv[1], "r");
if (fp == NULL) {
printf("Can’t open file: %s\n", argv[1]);
exit(-1);
グローバルイルミネーション入門 11
}
/* file has a correct magic number? */
if (fgetc(fp) != ’P’ || fgetc(fp) != ’3’ || fgetc(fp) != ’\n’) {
printf("File is not a ppm format.\n");
exit(-1);
}
/* skip comments. */
while((buf = fgetc(fp)) == ’#’) {
while(fgetc(fp) != ’\n’);
}
/* skip until newline. */
ungetc(buf, fp);
/* get image width and height. */
fscanf(fp, "%d %d\n", &width, &height);
/* get max value. */
fscanf(fp, "%d\n", &maxval);
image = (GLubyte *)malloc(width * height * 3 * sizeof(GLubyte));
for (j = height - 1; j >= 0; j--) {
for (i = 0; i < width; i++) {
fscanf(fp, "%d %d %d ", &r, &g, &b);
image[(i + j * width) * 3
] = (GLubyte)r;
image[(i + j * width) * 3 + 1] = (GLubyte)g;
image[(i + j * width) * 3 + 2] = (GLubyte)b;
}
}
glutInit(&argc, argv);
glutInitWindowPosition(0, 0);
glutInitWindowSize(width, height);
glutCreateWindow(argv[0]);
glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
glutKeyboardFunc(key);
glutDisplayFunc(display);
glutReshapeFunc(reshape);
グローバルイルミネーション入門 12
glutMainLoop();
}
レイトレーシング実装
5
まずはレイトレーシングの実装から始めよう。
5.1
直線と球との交点の計算
レイの進む方向の直線が球と交差するかを求める必要がある。これは中学で習う直線と円 (球)
の交わりを計算することと同じである。直線と球の交点は、代数的に求めることもできるが、グラ
フィックスではベクトル計算が主になることから、ここではベクトル幾何の考えから直線と球の交
点を求めることにしよう。
−
→
−
→
レイは、原点 p と方向ベクトル d (¯ここでは計算を簡単にするために、この方向ベクトルは長
→¯¯
¯−
さが 1 の単位ベクトルとする。つまり ¯ p ¯ = 0) の 2 つの変数で定義することができる (図 5)。
→
−
d
O
図 5: レイの表現
−
→
このレイがたどる軌跡によってできる直線上の任意の点 x は、パラメータ t を用いて、
−
→ −
→
→
−
x = p +td
(1)
−
→
と表現することができる。ここで、先ほど述べたように d は単位ベクトルであるので、パラ
−
→
−
→
メータ t はレイの原点 p から点 x までの距離と考えることができる。
球は、
¯−
→¯¯2
¯→ −
(2)
¯ x − c ¯ = r2
−
→
と表現することができる。ここで c は球の中心、r は球の半径である。直線の式 1 を、球の式 2
に代入すると、
r2
¯−
→
−
−
→¯¯2
¯→
= ¯( p + t d ) − c ¯
³−
→
→
−
−
→´ ³ −
→
−
→
−
→´
=
(p +td)− c
(p +td)− c
³−
→
→ −
−
→ ´³ −
→
−
→ −
→´
=
td +(p − c ) td +(p − c )
¯−
¯
¯−
¯
−
→ −
→
¯→¯
¯ →¯
= ¯ d ¯ t2 + 2 d · cpt + ¯cp¯
グローバルイルミネーション入門 13
よって t に関する二次方程式
¯−
¯
¯−
¯
→ −
−
→
¯→¯ 2
¯ →¯
¯ d ¯ t + 2 d · cpt + ¯cp¯ − r2 = 0
が得られる。ここで二次方程式の解の公式を当てはめると、
t=
−
→ −
→
− d · cp ±
r³
´
− −
→
→´2 ¯¯−
→¯¯ ³ −
→
d · cp − ¯ d ¯ |cp|2 − r2
¯
¯−
¯→¯2
¯d¯
交差判定を行ない、レイがヒットしたところは青色で、ヒットしないところは黒色でレンダリン
グするプログラムである。
1
/*
* raytrace1.c
*/
2
3
4
#include <stdio.h>
5
6
static double sphere_position[3] = {0.0, 0.0, -2.0};
7
static double sphere_radius
= 0.5;
8
11
int
intersect(double ray_origin[3],
double ray_direction[3]
12
{
9
10
13
14
15
double oc[3];
/* unit vector
*/ )
/* the vector from sphere center to ray origin */
double b, c, d;
double sr2;
16
18
oc[0] = ray_origin[0] - sphere_position[0];
oc[1] = ray_origin[1] - sphere_position[1];
/* x */
/* y */
19
oc[2] = ray_origin[2] - sphere_position[2];
/* z */
17
20
21
22
23
24
25
b = 2.0 * (ray_direction[0] * oc[0] +
ray_direction[1] * oc[1] +
ray_direction[2] * oc[2]);
sr2 = sphere_radius * sphere_radius;
c = oc[0] * oc[0] + oc[1] * oc[1] + oc[2] * oc[2] - sr2;
26
27
28
d = (b * b - 4.0 * c);
if (d > 0.0) return 1;
/* hit */
29
グローバルイルミネーション入門 14
return 0;
30
31
/* no hit */
}
32
33
34
35
36
37
38
39
40
41
42
43
int
main(int argc, char **argv)
{
FILE
*fp;
int
int
i, j;
width = 256;
int
double
double
height = 256;
sx, sy;
hw, hh;
double
const char filename[] =
org[3], dir[3];
"image.ppm";
44
45
46
fp = fopen(filename, "w");
if (fp == NULL) {
printf("Cannot create file : %s\n", filename);
exit(-1);
47
48
49
}
50
51
52
53
54
55
hw = (double)width
/ 2.0;
hh = (double)height / 2.0;
dir[0] = 0.0;
dir[1] = 0.0;
dir[2] = -1.0;
56
59
fprintf(fp, "P3\n");
fprintf(fp, "# raytrace1.c \n");
fprintf(fp, "%d %d\n", width, height);
60
fprintf(fp, "255\n");
57
58
/* magic number
/* comment
/* max pixel value
*/
*/
*/
61
62
63
64
65
66
67
68
for (j = 0; j < height; j++) {
for (i = 0; i < width; i++) {
sx = ((double)i - hw) / hw;
sy = ((double)j - hh) / hh;
org[0] = sx;
org[1] = sy;
org[2] = 0.0;
69
70
if (intersect(org, dir)) {
グローバルイルミネーション入門 15
fprintf(fp, "0 0 255 ", i, j);
71
} else {
72
fprintf(fp, "0 0 0 ", i, j);
73
}
74
}
75
}
76
77
fclose(fp);
78
79
}
図 6: raytrace1.c の実行結果
実行すると、図 6 のようなレンダリング結果となる。
6
パストレーシング
パストレーシングでは、視点からのレイトレーシングを行なう。レイがサーフェスに当たって反
射するときには、1 本のレイのみしか発生させない。そしてレイが光源もしくは背景に当たったと
きに反射を終了する。つまり全体的に見て、レイのパス (path, 経路) は 1 本のみに限定される。パ
ストレーシングは、Kajiya [12] により、レンダリング方程式の提案と同時にそれを解く手法とし
て一つの論文内で提案された。
グローバルイルミネーション入門 16
6.1
パストレーシング実装
では、まずはインポータンス (importance) や確率密度関数 (Probability Density Function, PDF)
などのコトは考えずに、またほとんどの係数やパラメータを決めうちとし、単純な設定でパスト
レーシングの実装を行なっていこう。サーフェスはすべてディフューズ面とし、カメラの位置やジ
オメトリの頂点などはワールド座標で指定し座標変換も行なわない。ライトは無限遠の球体とし
よう。つまりレイがどのオブジェクトにも当たらずに背景に飛んでいく場合に、それをライトに当
たったと見なすことにする。
まず、レイはカメラから発せられる。カメラの原点からあるピクセルを通ってレイは放たれる
(図 7)。
図 7: さて、どの方向へレイを反射させようか...
うんぬん....
レイがディフューズ面に当たったら、ロシアンルーレットでレイをさらに再帰的に反射させて追
跡を続けるか、もう追跡を諦めるかをディフューズ反射係数 Kd を用いて決定する。今回 Kd は
0.5 とした。
int
russian_roulette()
{
const double Kd = 0.5;
double r;
/* 50% reflectivity */
r = randomMT(); /* Mersenne Twister random number generator */
if (r < Kd) {
/* accepted */
return 1;
}
/* rejected */
return 0;
}
グローバルイルミネーション入門 17
randomMT() は、メルセンヌ・ツイスター法 (Mersenne Twister)12 による乱数生成関数であり、
[0.0, 1.0) 13 の範囲の浮動小数点数を返す。randomMT() の実装は以下のようになる。
#include "mt19937ar-cok.c"
/* Mersenne Twister random number generator */
#define randomMT genrand_real2
グローバルイルミネーションでは、乱数の質はレンダリング画像の品質に直接結びつくので優良
な乱数生成器を選ぶのは重要である。ただし、今は乱数の質はそれほど問題ではないので、この部
分は通常のライブラリ関数である rand(), drand48() などを用いてもよい。
?
?
?
図 8: さて、どの方向へレイを反射させようか...
続いて、レイがサーフェスに当たった後、どの方向へ反射したレイを飛ばすかについての戦略を
考えよう (図 8)。ディフューズ面は、半球上のどの方向へも光を均一に跳ね返す。つまり半球上の
どの方向も等確率でサンプルすればよい。しかし、式... でも判るように、コサインで減衰する。こ
こはコサインの減衰を考慮して反射方向をサンプルした方が効率がよい。コサインによる減衰も考
慮して、ディフューズ面での反射方向をサンプルするコードは以下のようになる。
static void
sample_cosweight(double outdir[3], const double normal[3])
{
int
i;
double r[2];
double cost, sint, phi;
double v[3];
double basis[3][3];
ortho_basis(basis, normal);
r[0] = randomMT();
r[1] = randomMT();
12 現在最良の乱数生成アルゴリズムとして知られている手法である。http://これこれ/。ただしその道の研究者に云わ
せるとほんのちょっと質に問題があるときがあるらしい。
13 半開区間であることに注意。つまり 1.0 は含まれない。
グローバルイルミネーション入門 18
cost = sqrt(r[0]);
sint = sqrt(1.0 - r[0]);
phi = 2.0 * M_PI * r[1];
v[0] = cos(phi) * sint;
v[1] = sin(phi) * sint;
v[2] = cost;
/* outdir = Basis . v */
for (i = 0; i < 3; i++) {
outdir[i] = v[0] * basis[0][i]
+ v[1] * basis[1][i]
+ v[2] * basis[2][i];
}
}
ortho basis() は、法線方向から正規直交基底行列を計算する。つまりはローカル座標での z
方向を法線方向へと合わせるための回転行列である。サンプル点はローカル座標で生成されるの
で、これを法線方向へと傾かせるために必要となる (図 9)。
図 9: ローカル座標で生成されたサンプル点を、サーフェスの法線方向を z 軸とした空間へ回転さ
せて合わせる。
グローバルイルミネーション入門 19
...
双方向パストレーシング
7
双方向パストレーシング (Bi-directional Path Tracing) は、視点からのレイトレーシングと光源
からのライトトレーシングを組み合わせた手法である。双方向パストレーシングは、Lafortune ら
[13] と Veach ら [22] によりそれぞれ独立して発表された。双方向パストレーシングでは、視点か
らのレイトレーシングの利点 (鏡面反射など) とライトトレーシングの利点 (コースティクスなど)
の両方を表現することが可能である。...
7.1
ノート
双方向パストレーシングは、Lafortune ら [13] と Veach ら [22] によりほぼ同じ時期に異なる学
会でそれぞれ独立して発表された。
双方向パストレーシングの文献や、オープンに公開されているコードというのは書籍でも web
でもあまり見受けられないように思われる。双方向パストレーシングの理論は、手法の提案者でも
ある Eric Veach の博士論文の第 10 章で詳細に解説されている。RenderPark は、双方向パスト
レーシングが実装されている数少ないレンダラであり、そのソースコードが参考になるであろう。
Advanced Global Illumination [9] には、双方向パストレーシングのコードフラグメント (コードの
一部) が掲載されている。ただしそのコードに対する解説は書中では行なわれていない。
...
トーンマッピング
8
現在の CRT モニタや液晶モニタなどのフレームバッファは、色の各要素を [0, 255] の 8 ビット
で表現している (低ダイナミックレンジ)。しかし高ダイナミックレンジ画像は非常に色の表現が
広いため、このような 8 ビットのデバイスでは単純に高ダイナミックレンジ画像を表示すること
はできない。そこで、トーンマッピング (tone mapping)14 とよばれる、高ダイナミックレンジ画像
を、うまく低ダイナミックレンジデバイスで表示するための手法が必要になる。
8.1
ノート
Tumblin-Rushmeier トーンマッピングオペレータ [21] は、コンピュータグラフィックスで初め
て用いられたトーンマッパーである。
14 階調変換とも訳される
グローバルイルミネーション入門 20
8.2
ローカルイルミネーション
ローカルイルミネーション (local illumination)15 アルゴリズムは、サーフェスに直接入射してく
る光のみを考慮してシェーディング (色付け) 処理を行う。
Z バッファレンダリング [5]、REYES アルゴリズム [6] などの隠面消去アルゴリズムと [1, 4, 11, 19]
などの経験的シェーディングモデルとの組み合わせがローカルイルミネーションアルゴリズムの代
表である。つまりは現在の GPU ハードウェアで行われている処理の多くがローカルイルミネー
ションの部類に入る。ローカルイルミネーションでは、光源から直接来る光しか考慮しない。その
ためローカルイルミネーションによるレンダリングは非常に高速に行うことができるが、間接光を
表現できないので、視覚的に豊かな表現を行うにはデザイナによるアートディレクションが必要に
なる。また影や遮蔽を表現するには [25, 8] などのアルゴリズムが別途必要になる。
また、ローカルイルミネーションの多くで用いられるシェーディングモデルは、コンピュータでの
計算に適するように現実世界の観測・経験結果に近似式を当てはめたものである。加えて 0 − 255 な
どのディスプレイハードウェアに合わせた直感的な単位を用いて照明計算を行うことが多い。その
ため物理学に基づいた照明計算を行うように設計されていないことから、物理的に正確なシミュー
レションが必要なレンダリングにローカルイルミネーションアルゴリズムを用いることは難しい。
しかしローカルイルミネーションアルゴリズムは、ライトや物体を動かした時の照明効果の変化
を捉えやすく (間接的な照明効果の依存関係がない)、照明値には直感的な単位を用いているので、
ライティングや色調をコントロールしやすいという利点がある。そのためスキルのあるアーティス
トは局所照明アルゴリズムを好む場合もある。
現在でも、多くのプロダクションでは局所照明アルゴリズムに基づいたレンダラを用いてレンダ
リングを行っている。
サンプリングとアンチエイリアシング
9
コンピュータグラフィックスは、本質的にサンプリングの問題であると言ってよいだろう。いか
にして少ないサンプルで誤差なく目的のレンダリング画像を計算するかがカギとなる。レンダリン
グにおけるサンプリングの問題は、80 年代にかけて多くの研究が行なわれている。
シーンは連続的で無限の周波数を持つ関数であるが、最終的にラスタ型ディスプレイに画像を表
示するために、ピクセルの色を決めねばならない。たとえばレイトレーシングは本質的にポイント
サンプリングである。そして
9.1
ノート
元の信号を復元するには、少なくとも元の信号の周波数の 1/2 以上の周波数でサンプリングす
る必要がある、というナイキストの定理 (Nyquist theorem) は、Harry Nyquist により [16, 17] で
初めて発表され定式化された。... アンチエイリアシングの研究は Mitchell ...
15 局所照明とも訳される
グローバルイルミネーション入門 21
10
準モンテカルロ法
モンテカルロ (Monte Carlo, MC) 法が乱数を用いて積分計算などの問題を解くのに対して、準
モンテカルロ (quasi-Monte Carlo, QMC) 法は、低食い違い量列 (low discrepancy sequence) と呼
ばれる決定論的 (deterministic) な数列を用いて積分計算を解く手法です。
準モンテカルロ法は、低次元の問題においてモンテカルロ法に比べて求めたい値に収束するのが
早く (最近では高次元での効率的な準モンテカルロ法も開発されているようです)、また低食い違
い量列を求めるコードは乱数を生成するコードよりも単純でハードウェアでの実装に向いていると
いう特徴があります。
aa aaa bbb ccc
ここでは、準モンテカルロ法の数学的な基礎や、各種低食い違い量列の説明と実装を説明し、準
モンテカルロ法を実際にコンピュータグラフィックスのレンダリングへの応用について解説します。
準モンテカルロ法は、モンテカルロ法に比べて高速であり少ないサンプル数で同じ品質を実現す
ることができるため、モンテカルロ法を利用している部分を準モンテカルロ法で置き換えることを
強く勧めます。本文がその手助けになれば幸いです。
以降、問題を簡単にするため、積分の区間は [0, 1]n (n は次元数) の超立方体の中で行うことと仮
定します。
10.1
食い違い度
食い違い度 (discrepancy) とは、ある数の集合 PN = x0 , . . . , xN −1 (N は数の個数、xi は i 番目
の数の値) がどれくらい一様に (等しく) 分布しているかというのを数学的に測るものさしです。
¯
¯
¯
A(J ∗, P) ¯¯
∗
¯
DN (P) = sup ¯λ(J ∗) −
¯
s
N
J ∗∈I
ここで、s は次元数で、J ∗ = [0, u1 ) × [0, u2 )×, . . . , ×[0, un ) =
Qs
i=1
ui ∈ [0, 1)s = I s になります。
この J ∗ はつまり、s 次元の超立方体 [0, 1)s 内で i 次元目の辺の長さが ui である部分空間 (subspace)
の取り得るすべての組み合わせを表します。たとえば 2 次元の場合、[0, 1)2 = I 2 の正方形内に含
まれ、[0, u1 ) × [0, u2 ) で表されるいろいろな組み合わせの長方形 (ヨコ u1 、タテ u2 ) になります。
ここで注意しておきたいのは、I も J ∗ も半開区間であることです (右側が閉じていない)。この
ことはこの後に J ∗ の面積を考えるときに重要になります。
sup は、上限 (supremum) を表し、sup f (x) は f (x) が取りうる値のうちで最大 (上限) のもの
を定義します。sup はほとんど max と似ていますが、max が定義できないものにも、sup が定義
できるのがあります。たとえば、f (x) = −1/x2 の関数では、f (x) < 0 となりますが、f (x) はゼ
ロにはならないので、max f (x) = 0 は定義できませんが、取り得る値の上限が 0 ということで、
グローバルイルミネーション入門 22
sup f (x) = 0 は定義することができます。
図 10: 大きな正方形は I 2 を表す。長い破線の領域に含まれる点の個数は 4 個なので、A = 4。短
い破線の領域に含まれている点の個数は 3 個なので A = 3
A(J ∗, P は、空間 I s に散らばっている点群 P のうち、部分空間 J ∗ に含まれる点の個数を表し
ます (図??を参照)。
λ(J ∗) は J ∗ のルベーグ測度 (ルベーグの意味での面積) になります。これはつまり J ∗ の面積
Qs
Q
Qs
|J ∗A| = i=1 ui (ここで、 は総積記号であり、 i=1 ui = ui u2 · · · un と定義されます) というこ
とになります。
なぜここで普通に面積 |J ∗| としないでルベーグ測度で λ(J ∗) とするのでしょうか。それは、J ∗
が半開区間であるからです。
結論としては、この場合の半開区間の範囲のルベーグの意味での面積は今までの考えでいう面積
と同じであるのですが、より厳密にこの半開区間の面積 (ルベーグ測度) を考えておくことにしま
す。
我々が小学校で習った意味での面積は、どちらも閉区間の閉じた範囲において与えられます (図
11の右)。
たとえば 2 次元で考えると、[0, u1 ] × [0, u2 ] で囲まれる面積は、ヨコの長さ u1 、タテの長さ u2 の
積、u1 u2 になります。しかし、今回の場合 J ∗ の張る区間の片方が開区間で閉じていないため (図
11の左)、今までの考えでの面積を厳密には与えることができません ([0, u1 ) × [0, u2 ) の面積を定
義することができない)。
しかし、ルベーグ測度という概念で面積というものを考えると、このような閉じていない空間において
も面積を与えることができるようになり、[0, u1 )×[0, u2 ) のルベーグ測度は λ([0, u1 )×[0, u2 )) = u1 u2
グローバルイルミネーション入門 23
(a)
(b)
図 11: 左: 図形は閉じている。右: 図形の上と右は閉じていない。
となり、閉区間の場合の面積と等しくなります。よって、J ∗ のルベーグ測度というのは、J ∗ が
閉区間の場合の面積 |J ∗| と同等になります。
よって、点群 P の star discrepancy D∗ (P) とは、あらゆる y i , i = 1, . . . , n の組み合わせでできる
部分空間 A に含まれる点群 P の個数を点群の総数 N で割った値と、部分空間の面積 |A| の差が、
最も大きくなるときの値、ということになります。
定理 1 もし 0 ≤ x1 ≤ x2 ≤ · · · ≤ xN ≤ 1 であれば、
∗
DN
(x1 , x2 , . . . , xN )
¯
¯
¯
1
2n − 1 ¯¯
¯
=
+ max xn −
2N 1≤n≤N ¯
2N ¯
証明 1 x0 = 0, xn+1 = 1 とします。P を点群 x1 , . . . , xn とすると、
¯
¯
¯ A([0, u); P )
¯
∗
¯
DN (P ) = max
sup
− u¯¯
¯
0≤n≤N xn ≤u≤xn+1
N
¯n
¯
¯
¯
= max
sup
¯ − u¯
0≤n≤N xn ≤u≤xn+1 N
¯ ¯n
¯´
³¯ n
¯
¯ ¯
¯
= max max ¯ − xn ¯ , ¯ − xn+1 ¯
0≤n≤N
N
N
¯¶
µ¯
¯ ¯¯ n − 1
¯
¯n
¯
= max max ¯ − xn ¯ , ¯¯
− xn ¯¯
1≤n≤N
N
N
¯
¯
¯
¯
1
2n
−
1
¯.
=
+ max ¯xn −
2N 1≤n≤N ¯
2N ¯
ここで、n = A([0, u); P ) は部分空間 A に含まれる点群 P の個数です。
定理 1 Koksma の不等式 f が [0, 1] で有界変動 Vf をもてば、
が成り立ちます。
¯Z
¯
N
¯ 1
¯
1 X
¯
¯
∗
f (u) du −
f (xn )¯ ≤ V (f )DN
(S)
¯
¯ 0
¯
N n=1
(3)
グローバルイルミネーション入門 24
証明 1 Koksma の不等式の証明 x1 ≤ x2 ≤ · · · ≤ xN とし、また x0 = 0, xN +1 = 1 とします。
式 (3) の右辺の絶対値の中の総和は、summation by parts の公式、
n
X
j=m
aj bj = [Sn bn − Sm−1 bm ] +
n−1
X
j=m
Sj (bj − bj+1 ), SN =
N
X
an .
n=1
を用いて、
N
X
f (xn ) =
n=1
N
X
n=1
=
1 · f (xn )
N
−1
X
[N f (xN ) − 0 · f (1)] +
= N f (xN ) −
N
−1
X
n=1
n=1
n (f (xn ) − f (xn+1 ))
n (f (xn+1 ) − f (xn ))
= N f (xN +1 ) − N (f (xN +1 ) − f (xN )) −
= N f (xN +1 ) −
= N f (xN +1 ) −
= N f (1) −
N
X
n=0
N
X
n=1
N
X
n=0
N
−1
X
n=1
n (f (xn+1 ) − f (xn ))
n (f (xn+1 ) − f (xn ))
n (f (xn+1 ) − f (xn ))
n (f (xn+1 ) − f (xn ))
よって、
N
N
X
n=1
f (xn ) = f (1) −
N
X
n
(f (xn+1 ) − f (xn )) .
N
n=0
右辺の絶対値の中の積分は、部分積分を用いて、
Z 1
Z 1
f (u)u0 du
f (u) du =
0
0
= [f (u)u]10 −
= f (1) −
= f (1) −
Z
1
Z
u
0
Z
1
f 0 (u)u du
0
df (u)
du
du
1
u df (u)
0
になります。よって、式 (3) の右辺の絶対値の中は、
Z
1
0
N
1 X
f (u) du −
f (xn )
N n=1
µ
Z
=
f (1) −
=
Z
1
0
1
0
!
¶ Ã
N
X
n
(f (xn+1 ) − f (xn ))
u df (u) − f (1) −
N
n=0
u df (u) −
N
X
n
(f (xn+1 ) − f (xn ))
N
n=0
グローバルイルミネーション入門 25
となり、ここで、
Z
1
u df (u) =
0
N Z
X
n=0
および、
xn+1
u df (u),
xn
Z
f (xn+1 ) − f (xn ) = [f (u)]xxnn+1 =
を用いて整理することで、
Z 1
N
1 X
f (u) du −
f (xn )
N
0
n=1
1 df (u),
xn
N Z
X
=
xn+1
xn+1
xn
n=0
N Z xn+1
X
=
xn
n=0
となります。
³
u−
u−
n´
df (u)
N
n
df (u).
N
定理 (1) から (2 行目に注目)、0 ≤ n ≤ N となる固定された n および各 xn ≤ uxn+1 に対して、
¯
n ¯¯
¯
∗
(x1 , . . . , xn )
¯u − ¯ ≤ DN
N
が成り立つので、
N Z xn+1 ¯
N Z xn+1
X
X
n ¯¯
¯
∗
df (u)
DN
¯u − ¯ df (u) ≤
N
x
x
n
n
n=0
n=0
∗
DN
=
N Z
X
n=0
∗
= DN
N
X
xn+1
1 df (u)
xn
(f (xn+1 ) − f (xn ))
n=0
また、有界変動 V (f ) は、
V (f ) = sup
I
N
X
∆J=∈I,i=0
(f (xi + ∆) − f (xi ))
で定義されていることにより (J は I をあらゆる分割方法の可能性により分割した領域とし、その
分割数を N 、その幅を ∆ とする)、
N
X
(f (xn+1 ) − f (xn )) ≤ V (f )
n=0
となるので (この不等式は、右辺は N 個の荒い分割で測った変動量であるのに対し、左辺はより精
細な間隔で測った上にその最大 (上限) を取った変動量であることをイメージすれば、理解し易い
と思います)、よって最終的に、
¯Z
¯
N
¯ 1
¯
1 X
¯
¯
f (u) du −
f (xn )¯
¯
¯ 0
¯
N
n=1
≤
∗
DN
≤ V
N
X
(f (xn+1 ) − f (xn )
n=0
∗
(f )DN
となります。 ■
グローバルイルミネーション入門 26
10.1.1
k 次元の場合
定理 1 Koksma-Hlawka の不等式 f が [0, 1]k で、Hardy-Krause の意味での有界変動 (the vari-
ation in the sense of Hardy-Krause)16 Vf をもてば、
¯Z
¯
N
¯
¯
1 X
¯
¯
∗
f (x) dx −
f (xn )¯ ≤ Vf DN
(S)
¯
¯ [0,1]k
¯
N
n=1
が成り立ちます。
...
10.2
Hardy-Klause の意味での有界変動
Hardy-Klause の意味での有界変動は、Vitali の意味での有界変動 (the variation in the sense of
Vitali) から導きだされます。
10.3
Vitali の意味での有界変動
関数がどれくらい変化しているか (どれくらい起伏の程度があるか) を測った量を変動量 (varia-
tion) と呼びます。変動量が有限の値である (無限にならない) とき、これを有界変動と呼びます。
たとえば、y = 3 のような関数は全く起伏がないので、変動量はゼロになります。
10.4
一様格子分布 (regular grid, uniform grid)
空間を一様に分布する配置として、通常我々は空間を一様の大きさに分割するような点の配置が
思い浮かぶと思います。たとえば一次元の場合、サンプル点が N 個、一様格子分布にもとづいた
数列を Punif orm とすると、
Punif orm =
0 1
N −1
, ,...,
N N
N
となる数列です。
直感的に考えて、一様格子分布は、一様に分布しているのだから最も食い違い度が低い数列ではな
いのかと思ってしまいますが、実際に食い違い度を測るとそうでないことが分かります。一様格子
分布のスター食い違い度は、
1
∗
s
Dunif
orm = O(f rac1N )
になります17 。ここで s は次元数です。つまり一次元に限れば、一様格子分布のスター食い違い度
は O(1/N ) になるので O(logs N/N ) となるような低食い違い列よりも勝ることがありますが、一
様格子分布のスター食い違い度は次元数が増えるほど次元爆発 (dimensional explosion) が起こ
り大きくなってしまうのが分かります。
16 この説明
17 これの証明
グローバルイルミネーション入門 27
10.5
一様分布乱数列 (uniformly distributed random series)
モンテカルロ法では、次元爆発を避けるために乱数を用いました。これは、スター食い違い度を
∗
測ることで確認するこができます。一様分布乱数列 Drandom
のスター食い違い度は、
µ
¶
log log N
∗
Drandom
=O
2N
になります (確率 1 で)。
証明
一様に分布している乱数列 xi から、以下により新しい乱数 ξi を定義します。
ξi =
(
a 2
c
4
もし xi が [0, 1) で一様に分布しているのであれば、 ξi の期待値と分散は、
E(ξi ) = λ(A), V 2 (ξi ) = λ(A) − λ(A)2 ≤
1
4
ここで λ(A) は A のルベーグ測度 (ルベーグの意味での面積) です。
もし xi が独立乱数であれば、生成される ξ1 , ξ2 , . . . , ξN 乱数もまた独立になり、期待値 E(ξ) と分
散 V (ξ) をもつ等しい分布になります。
10.6
van der Corput 列
低食い違い量列は、根基逆関数 (radical inverse function) を元にして生成します。18 10 進法の
整数 n を、b 進法 ( b ≥ 2 ) で表し、このときの j + 1 桁目の数字を aj (n) とします (式の定義が
j = 0 から始まるのですが、日本語で 0 桁目というのは変なので j + 1 桁目という表現を使いまし
た)。すると、10 進法での n は、
∞
X
n=
aj (n)bj
j=0
で表すことができます。たとえば 10 進数での 6 は 2 進数で表すと 110 になるので、 a(0) =
0, a(1) = 1, a(2) = 1 になります。
このとき、基数 b の根基逆関数は、以下で定義されます。
Φb (n) =
∞
X
aj (n)b−j−1
j=0
つまり、b 進法での Φb (n) は、 n を小数点で対称に折り返したものになります。また、n ≥ 0 の
すべての整数 n で Φb (n) ∈ [0, 1) になるという性質があります。図 2 に基数 2 のときの van der
Corput 列を示します。
18 vandercorput
列のオリジナルの論文
グローバルイルミネーション入門 28
n
2 進法での Φ2 (n)
10 進法での Φ2 (n)
0
1
2
0
0.1
0.01
0.000
0.500
0.250
3
4
0.11
0.001
0.750
0.125
5
6
7
0.101
0.011
0.111
0.625
0.375
0.875
表 2: 基数 2 の van der Corput 列
10.7
Halton 列
Halton 列では、van der Corput 列を次元数分並べたベクトルになります19 。s 次元の Halton
列は、以下のように定義されます。
xsn = (Φb1 (n), Φb2 (n), . . . , Φbs−1 (n))
このとき、根基逆関数の bi には重ならない素数を使います。例えば 2 次元で、 b1 = 2, b2 = 3 と
した時の、 Halton 列を表 3 に示します。
1 次元の Halton 列というのは、 van der Corput 列と同じです。
n
2 進法での Φ2 (n)
3 進法での Φ3 (n)
(Φ2 (n), Φ3 (n))
0
1
0
0.1
0
0.1
(0.000, 0)
(0.500, 1/3)
2
3
0.01
0.11
0.2
0.01
(0.250, 2/3)
(0.750, 1/9)
4
5
6
0.001
0.101
0.011
0.11
0.21
0.001
(0.125, 4/9)
(0.625, 7/9)
(0.375, 1/27)
7
0.111
0.101
(0.875, 10/27)
表 3: 2 次元の Halton 列 (b1 = 2, b2 = 3)
19 halton
列の論文
グローバルイルミネーション入門 29
10.8
Hammersley 列
Hammersley 列は、 第一次元の要素を n/N (N は生成するサンプル数。 n ≤ N − 1) にする以
外は、Halton 列と似ています20 。
xsn = (n/N, Φb1 (n), Φb2 (n), . . . , Φbs−2 (n))
Hammersley 列は、N に依存するので、前もって生成するサンプル数が決定されていなくては
ならないという欠点があります。その点、 van der Corput 列や Halton 列は、N に依存しないの
で、インクリメンタルに計算を行うことができます。その代わり、Hammersley 列は、 Halton 列
よりも少しよいスター食い違い度を持ちます。
Hammersley 列で最も特徴的な点の並び方は、2 次元の時で基数が b1 = 1 、かつ生成するサン
プル数が N 2 のときです。このときのみ、 Hammersley 列の点群は、空間 N で等分してできる格
子上に置かれることになります。
n
n/N
2 進法での Φ2 (n)
(n/N, Φ2 (n))
0
1
0/8
1/8
0/8
0.1
(0 , 0.000)
(1/8, 0.500)
2
3
4
2/8
3/8
4/8
0.01
0.11
0.001
(2/8, 0.250)
(3/8, 0.750)
(4/8, 0.125)
5
6
5/8
6/8
0.101
0.011
(5/8, 0.625)
(6/8, 0.375)
7
7/8
0.111
(8/8, 0.875)
表 4: N = 8 のときの 2 次元の Hammersley 列 (b1 = 2)
準モンテカルロ法はとくに金融工学の分野で広く用いられています.
10.9
ノート
準モンテカルロの研究は、Harald Niederreiter 博士が多大なる貢献をしている。[14] は、この
分野におけるバイブルとも云える書籍である。
有限体理論からの低食い違い量列
多くの低食い違い量列は、 Faure 列 ( Faure sequence) を元に
生成されてきましたが、 Niederreiter 氏と Xing 氏により、有限体 (finite field) に基づいた全く新
しい低食い違い列である Niederreiter-Xing 列 (Niederreiter-Xing sequence, NX 列)[26] が考案さ
れました。
NX 列は、既存の LDS の中でもかなり優秀なもののようで、それなりの高次元でも使えるよう
です。
NX 列の実装は [1] で入手することができます。
20 hammersley
の論文
グローバルイルミネーション入門 30
この NX 列の理論の元になっている有限体とは、足したり掛けたりの演算を行っても有限の位
数 (要素数) しか取りえない世界 (体) のことです。有限で閉じているので、コンピュータでの演算
に適しており、実際有限体は楕円曲線、暗号、符号理論、低食い違い量列などと幅広く応用されて
いるようです。
10.10
一様分布以外の生成
参考文献
[1] 湯前祥二, 鈴木輝好著, モンテカルロ法の金融工学への応用
[2] Alexander Keller, Strictly Deterministic Sampling Methods in Computer Graphics, (mental
images technical report, 2001) in ”Monte Carlo Ray Tracing”, SIGGRAPH ’2003 Course
44, San Diego, July, 2003.
[3] T. Wong, W. Luk, and P. Heng, Sampling with Hammersley and Halton points, Journal of
Graphics Tools, vol. 2, no. 2, pp. 9-24, 1997
11
数学的予備知識
11.1
放射測定
11.1.1
放射輝度
コンピュータグラフィックスにおいて、放射輝度はもっとも重要な単位である。なぜなら、放射
輝度はレンダリングされる画像のピクセル値と直接に対応するものであるからである。
11.1.2
測度論的放射輝度
測度論からのアプローチによる放射輝度の扱いは
11.2
ノート
放射測定については、 Nicodemus が先駆者であり多大なる仕事をしている。Nicodemus による
放射測定の仕事は、NIST((米国) 標準技術局) から [?] にまとめられている。放射測定の考えは、
Kajiya がレンダリング方程式 [12] を導出するためにコンピュータグラフィックスの分野に初めて
導入された。
11.3
随伴法
通常、光のシミュレーションというと、フォトン (光子) を光源から放出し、サーフェスで反射
し、最終的にカメラやフィルム面などの検出器 (detector) へと飛び込んでくる光の経路をシミュ
グローバルイルミネーション入門 31
レーションする。しかし、レイトレーシングなどで見てきたように、カメラから光源へと逆方向の
流れで光の経路を計算しても同じ結果が得られる。この対称性を数学的に記述し証明するのが、随
伴法 (Adjoint method) である。
11.4
インポータンスサンプリング
インポータンスサンプリング (importance sampling)21 は、結果への寄与が大きい重要な (important)
部分を、集中的にサンプリングする手法である。グローバルイルミネーションなど、サンプリング
の問題へと帰着することの多いコンピュータグラフィックスの分野では、インポータンスサンプリ
ングは重要な概念である。なぜならより少ないサンプル数で、より正確なレンダリング結果を得る
ことができるからである。
便宜を図るため、d 次元超立方体 [0, 1)d での積分
I=
Z
f (x) dx
[0,1)d
を考えよう22 。これは、適当な関数 p(x) を用いて、
I=
Z
f (x)
p(x)dx
p(x)
と変形して表現することもできる。ここで、p(x) は [0, 1)d 上で定義される確率密度関数 (prob-
ability density function, PDF) とすると、我々が示したいことが容易になる。p(x) から独立して
サンプルされた Xi を用いた I のモンテカルロ推定 Iˆ は
n
となる。
11.4.1
1 X f (Xi )
Iˆ =
n i=1 p(Xi )
ノート
インポータンスサンプリングに関する統括的な解説かつ奥深い洞察については、 Owen らの [18]
が参考になるだろう。
11.5
ロシアンルーレット
ロシアンルーレットのコンピュータグラフィックスへの導入は、Arvo ら [3] により初めて行なわ
れた。
一次元の積分
21 加重サンプリング、重点的サンプリング、要点抽出法などとも訳される
22 変数変換や区間の写像を行なえば、多くの積分はこの形式で表現することができる
グローバルイルミネーション入門 32
I=
Z
1
f (x)dx
0
を考えよう。ここで f を水平方向に P でスケーリングし、垂直方向に 1/P でスケーリングす
ることを考える (P ≤ 1)。このときの積分値 Irr は、
Irr =
となる。I = Irr ...
Z
P
0
1 x
f ( )dx
P P
ロシアンルーレットによる、分散の評価は、... が行なっている。
11.6
レンダリング方程式
3D CG におけるレンダリングの照明計算式は、レンダリング方程式 [12] で定式化を行うことが
できる。サーフェス上のある点 x から、方向 ω へと出射していく放射輝度 L(x, ω) は以下のよう
に定義される [12]。
L(x, ω) = Le (x, ω) +
Z
L(x, ω 0 )f (ω, ω 0 ) cos θ dω 0
(4)
Ω
ここで積分区間 Ω は、サーフェス上の点 x の法線方向の半球になる。dω 0 は微小立体角であり、
これは球座標 (θ, φ) で示すと sin θdθdφ となる。f は点 x での BRDF (Bidirectional Reflectance
Distribution Function、双方向反射分布関数)、cos θ の θ は点 x の法線と ω 0 との角度である。
Le (x, ω) はサーフェスから発光する放射輝度で、これは通常光源となるサーフェスでのみ非ゼロ
になる。
この式は図 12 で図示することができる。
=
+
図 12: Illustration of the rendering equation.
つまりこの積分の示すところは、ある点からある方向へと出射する放射輝度は、法線方向の半球
にあらゆる方向から入射してくる放射輝度に BRDF 積を掛けたものの積分を取るということを示
している。
式 (4) は、積分演算子23 T を用いて、
L = Le + T L
(5)
と簡潔に表記する事ができる。ここで T は、
23 operator
の訳語は物理学分野では演算子、数学分野では作用素と訳される。ここでは物理学の方を用いることにした
グローバルイルミネーション入門 33
0
(T ◦ L)(x, ω ) :=
と定義される。
Z
L(x, ω 0 )f (ω, ω 0 ) cos θ dω 0
Ω
式 (5) は両辺に L があるためそのままでは評価を行うことができない。レンダリング方程式
は第二種フレッドホルム積分方程式として捉えることができる。エネルギー保存則が成り立つと
すれば、積分演算子のノルムは |T | < 1 と仮定できるので、式 (5) にノイマン級数展開
P∞
1 + z + z 2 + · · · = i=0 z i (z < 1)24 を適用し,
L = Le + T Le + T (T Le ) + · · · =
∞
X
i=0
T i Le
1
1−z
=
(6)
を得る。これにより右辺は未知関数 L の代わりに既知関数 Le で置き換えられるので、 i を適
当な次数で打ち切ることで L の評価を行うことができる。
12
参考書籍
ここでは、グローバルイルミネーションを知る・学ぶ上で有益な書籍を紹介したいと思う。Advanced
Global Illumination12.
24 つまりは初項
1, 公比 z の等比級数の公式である
グローバルイルミネーション入門 34
解析概論。別名高木本。古典的な初等数学の入門書。世の大学の数学科の学生は必携だとか。
グローバルイルミネーション入門 35
参考文献
[1] Arthur Appel. The notion of quantitative invisibility and the machine rendering of solids.
In Proceedings of ACM National Conference, pages 387–393, 1967.
[2] Arthur Appel. Some techniques for shading machine renderings of solids. In Proceedings of
the Spring Joint Computer Conference, pages 37–45, 1968.
[3] James Arvo and David Kirk. Particle transport and image synthesis. In Computer Graphics
(SIGGRAPH ’90 Proceedings), pages 63–66, 1990.
[4] W. Jack Bouknight. A procedure for generation of three-dimensional half-toned computer
graphics presentations. Communications of the ACM, 13(9):527–536, 1970.
[5] Edwin E. Catmull. A Subdivision Algorithm for Computer Display of Curved Surfaces. PhD
thesis, Dept. of CS, U. of Utah, December 1974.
[6] R. Cook, L. Carpenter, and E. Catmull. The reyes image rendering architecture. In Proceedings of SIGGRAPH 1987, pages 95–102, 1987.
[7] Robert L. Cook. Shade trees. In Proceedings of SIGGRAPH 1984, pages 223–231, 1984.
[8] Franklin C. Crow. Shadow algorithms for computer graphics. Computer Graphics (SIGGRAPH ’77 Proceedings), 11(2):242–248, July 1977.
[9] Philip Dutré, Philippe Bekaert, and Kabita Bala. Advanced Global Illumination. A.K. Peters
Ltd., 2003.
[10] Cindy M. Goral, Kenneth E. Torrance, Donald P. Greenberg, and Bennett Battaile. Modeling the interaction of light between diffuse surfaces. In Proceedings of SIGGRAPH 1984,
pages 213–222, 1984.
[11] Henri Gouraud. Continuous shading of curved surfaces. IEEE Transactions on Computers,
20(6):623–628, 1971.
[12] James T. Kajiya. The rendering equation. In Proceedings of SIGGRAPH 1986, volume
20(4), pages 143–150, 1986.
[13] Eric P. Lafortune and Yves. D. Willims. Bi-directional path tracing. In Third International
Conference on Computational Graphics and Visualization Techniques(Compugraphics ’93),
pages 145–153, December 1993.
[14] Harald Niederreiter.
Random Number Generation and Quasi-Monte Carlo Methods.
Philadelphia, Pa. : Society for Industrial and Applied Mathematics, 1992.
[15] Tomoyuki Nishita and Eihachiro Nakamae. Continuous tone representation of 3-d objects
taking account of shadows and interreflection. Computer Graphics (SIGGRAPH ’85 Proceedings), 19(3):23–30, July 1985.
グローバルイルミネーション入門 36
[16] Harry Nyquist. Certain factors affecting telegraph speed. Bell System Technical Journal,
page 324, 1924.
[17] Harry Nyquist. Certain topics in telegraph transmission theory. A.I.E.E. Trans., 47:617,
1928.
[18] Art B. Owen and Yi Zhou. Safe and effective importance sampling. volume 000, pages
00–00, 1998.
[19] Bui-T. Phong. Illumination for computer generated pictures. Communications of the ACM,
18(6):311–317, June 1975.
[20] Peter-Pike Sloan, Jan Kautz, and John Snyder. Precomputed radiance transfer for real-time
rendering in dynamic, low-frequency lighting environments. In Proceedings of SIGGRAPH
2002, pages 527–536, 2002.
[21] J Tumblin and H. E. Rushmeier. Tone reproduction for realistic images. volume 13, pages
42–48, 1996.
[22] Eric Veach. Bidirectional estimators for light transport. In Fifth Eurographics Workshop
on Rendering, pages 147–162, June 1994.
[23] Eric Veach. Non-symmetric scattering in light transport algorithms. In Eurographics Rendering Workshop 1996 Proceedings, pages 82–91, June 1996.
[24] Turner Whitted. An improved illumination model for shaded display. Communications of
the ACM, 23(6):343–349, 1980.
[25] Lance Williams. Casting curved shadows on curved surfaces. In Computer Graphics (SIGGRAPH ’78 Proceedings), volume 12, pages 270–274, August 1978.
[26] Chaoping Xing and Harald Niederreiter. A construction of low-discrepancy sequences using
global function fields. volume 73, pages 87–102, 1995.
グローバルイルミネーション入門 37
索引
確率密度関数, 17
クラシックレイトレーシング, 4
グローバルイルミネーション, 1
双方向パストレーシング, 20
めとろぽりすひかりゆそう, 7
メルセンヌ・ツイスター, 18
レイキャスティング, 4
レイトレーシング, 4
ローカルイルミネーション, 21
38
13
更新履歴
• 2004 年 7 月 30 日
– レイトレーシング実装の節を追加
• 2004 年 7 月 6 日
– 放射測定の節を追加
– 準モンテカルロ法:ノートを加筆
• 2004 年 7 月 5 日
– アルゴリズムの歴史を加筆
• 2004 年 7 月 4 日
– プログラミングツールと準備の節を追加
• 2004 年 7 月 1 日
– 準モンテカルロ法の節を追加 (今まで書きかけだったのをマージ)。
– フッターを追加。
• 2004 年 6 月 29 日
– 参考書籍の節を追加。
– 更新履歴を取リ始める。
グローバルイルミネーション入門 39