2004年03月31日
light mapping
vray 1.5 で搭載されるという light-mapping(ライトマッピング) が面白そうなので、ちょこっと更新。
http://vray.info/newsread.asp?ID=81
いわゆるフォトンマップ以前の、フォトンマップの基礎となったライトマッピング(ライトからレイトレーシングを行う手法)とは異なる手法で、視点からのレイトレーシングになります。
視点からのレイの各跳ね返りでは、照明値をフォトンマップのように三次元空間データ構造に保存するようで、視点からフォトンを放出するフォトンマップと言えますね。
ライトマッピングの利点は、
o フォトンマップでは、ライトが複数ある場合、それぞれでフォトンの放出をセットアップしなければならないが、ライトマッピングでは視点からのみなので単純
o ライトマッピングでは、自己発光体、非物理的ライト、スカイライトなどのどの種類のライトでも扱える。それに対して、フォトンマッピングでは、再生成(reproduce)可能な照明効果のみしか扱えない(たとえば、自乗減衰しない点光源は扱えない)
o ライトマップはシーンのカドや小さな物体の周りの照明も正しく扱える。それに大して、フォトンマップは密度推定(density estimation)などの手法が必要になるが、これはときどき間違った結果になる。
o 多くの場合、ライトマッピングは直接に可視化できる(directly visualized)(つまりファイナルギャザーを必要としない!?)
サンプル画像を見る限り、ライトマッピングは非常によい結果が見受けられます。
またスピードもフォトンマッピングと同等で、かつアニメーションにも容易に適用できるとの事。
視点からのトレース時に、どのように照明値を計算するのかは分かりませんが、ある意味視点からのインポータインスフォトンマッピングのようなものになるのでしょうか?
コースティクスもきれいに表現できるのでしょうか?
もしコースティクスがうまく表現できないとしても、コースティクスはフォトンマップで行い、
そのほかの間接照明効果ではライトマッピングを用いるということができるでしょう。
多くの場合では、ライトマッピングは直接に可視化できるらしいので、間接照明の表現にグローバルフォトンマップで必要となるファイナルギャザーがいらないというのは大きなアドバンテージな気がします。
メトロポリス光輸送とも組み合わせることができたりするかな。
とまあフォトンマップに比べてこれでもかってくらいアドバンテージが述べられています。
個人的にはフォトンマップ本でメトロポリス光輸送がけっこうボロクソに書かれていて悲しかったので(近年のラージステップ MLT 法の提案などで、結構フォトンマップ本で書かれている欠点の多くはすでに解消されていると思います)、ポストフォトンマップとなるレンダリングアルゴリズムを支持したいところですね。
2004年03月22日
BSP ツリー実装
いつまでたってもレイトレの空間データ構造が一様グリッドとオクツリーだけでは、
たとえば非常に広大な地平面があるシーンなどでは効率が落ちるので、
BSP ツリーを実装しました。
Graphics Gems を参考にした、単純な軸平行(Axis-aligned) BSP ツリーです。
同じシーンを一様グリッドと BSP ツリーでレイトレにかかった時間を測定。
一様グリッドのトラバースでは、 altivec を無効にしています。
chrysler.rib:
BSP: 73.762986 秒
GRID: 96.403030 秒
chinahou.rib:
BSP: 6.590258 秒
GRID: 9.106296 秒
BSP の実装は非常にストレートで最適化も何も行っていませんが、
一様グリッドよりも効率がよい....
今まで一様グリッドを必死になって最適化して
きたのは無駄だったかなぁ...
BSP ツリーもこれから SIMD で高速化したり、
ツリーの構築のときの分割面の決定をよりよいもの
に改善するなどすれば、さらなる向上が期待できそうです。
たとえば、 Igno Ward 氏らによるコヒーレントレイトレーシング
SIMD高速化 参照
では、軸平行 BSP ツリーのトラバースを SIMD で行って一次レイの高速化を図っています。
また、一様グリッドまたはオクツリーと BSP の混合手法というのもありますね。
最適な BSP の構築というのは非常によく研究されている分野で、
突き詰めるのは難しそうですが、
もう少し構築の改善とトラバースの SIMD 化を行っていこうかと思います。
BSP ツリーは要素が多い時のメモリ量と構築時間の問題というのもあるらしいので、
非常にポリゴン数が膨大な時に、BSP ツリーでは
どれくらいメモリ量と構築時間がかかるかというのも測定しなければなりませんね。
2004年03月21日
rover と線形合同法
MLT では、画像面のピクセルをランダムにサンプルするプロセス(レンズサブパス変異)があります。
これを、ただ単に擬似乱数でサンプルしようとしても、
すべての画像面のピクセルを一様にまんべんなく取り上げていくことはなく、
サンプルされない領域とサンプルされすぎる領域に分かれてしまい、
非常に効率が悪いというのは良く知られています。
この問題には、層別化や LDS を使うというのも考えられていますが、
Eric Veach 氏の博士論文では、これに rover(ローバー)と呼ばれる
メモリ管理法からの応用を利用することを述べられています。
(11.4.4 節)
rover というのは、各ピクセルが一度だけ表れるようなランダムなピクセル順を
生成するインデックスアルゴリズムのことのようです。
(rover memory とかで検索すると火星探査機の rover がヒットしてしまう)
各ピクセルにはサンプルされた回数を記録するカウンタを用意しておき、
カウンタがある値を超えたら(サンプルされすぎたら)、
rover によるインデックスを参照して次のランダムのピクセルを
採用するようにします。
rover では、各ピクセルが一度だけ表れるようなランダムなピクセル順を生成するのに、
線形合同法(linear congruential) の低位ビットを利用しています。
各ピクセルが一度だけ表れるというのは linear pixel shuffling でも同じですね。
(linear pixel shuffling と線形合同法は式も似ていますしね)
線形合同法
線形合同法は、擬似乱数を生成するアルゴリズムの古典的な方法で、たとえば
gcc では以下のように実装されているらしいです。
(たぶん最近では Mersenne Twister に代わっているでしょうか?)
unsigned int seed = 1;
int
randLC()
{
seed = 1103515245 * seed + 12345;
return seed;
}
という非常に単純なアルゴリズムです。
線形合同法では、低位 n ビットを取り出すと、それが周期 2^n の
数列を成すことが知られています。
またこのとき、各周期 2^n では、おのおの重なる数字は現れません。
この性質を利用する事で、各要素の数字が必ず 1 度しか現れないランダムな
数列を生成することが可能になります。
ためしに 2^8 = 256 個の数列を作って重なりがないかテストするコードを書いてみました。
/* lc.c */
#include <stdio.h>
#define GENSIZE 256
unsigned int seed = 1;
int table[GENSIZE];
int
randLC()
{
seed = 1103515245 * seed + 12345;
return seed;
}
void
init()
{
int i;
for (i = 0; i < GENSIZE; i++) {
table[i] = 0;
}
}
void
test()
{
int i;
int r;
int t;
int count = 0;
for (i = 0; i < GENSIZE; i++) {
r = randLC();
t = r & (GENSIZE-1); /* extract lower bits */
if (table[t] != 0) {
printf("boo! t = %d\n", t);
}
table[t] = 1;
}
for (i = 0; i < GENSIZE; i++) {
if (table[i] == 1) count++;
}
printf("count = %d\n", count);
}
int
main()
{
init();
test();
}
実行結果。
$ gcc lc.c $ ./a.out count = 256
と、重なりが無い事が確認できます。
Eric Veach の博士論文では、さらに各ピクセルの内部(サブピクセル)では、
Poisson minimum-disk パターンというのを用いているそうです。
2004年03月16日
アンビエントオクルージョンシェーダ
ずっと手つかずで忘れていた RenderMan シェーダ言語の
G.I. 関数である occlusion() 関数を、 lucille に移行しました。
ノイズがあるのは LDS ではなく乱数を用いているからです。
サンプル生成関数
MC, QMC に関わらず、C シェーダや内部実装で、
サンプル点生成のインターフェイスとなる関数をそろそろ実装しようといけないなぁ
と感じています。
mental ray では、一次元・多次元にかかわらずサンプル点の取得には、
mi_sample()
という1 つの API ですべて行えます。
これは乱数ではなく、常に LDS がサンプル点として返ってくるようです。
サンプル数が少ない(数十以下)と QMC ではエイリアシングが出やすい傾向にあるので、
求めるサンプル点が少ない場合は乱数に切り替えるようにするということもできるでしょう。
次の実装は trace() かな。
ヒルベルト曲線順を組み込んでみた。
ヒルベルト曲線順で一様グリッドデータを構成するように変更してみました。
chrysler シーンでテスト。
結果は... 既存のコードと実行時間何も変わらず....
うーむ、残念。
変更した部分は丁度プロファイラ Shark が検出した最も時間のかかっている部分に直に関連する部分なんだけどなぁ...
というかヒルベルト曲線順のインデックスを求めるのに 20 秒程度(128*128*128)かかり
逆に全体では遅くなってしまいました。
コードの構成自体にも問題があるのかもしれませんが、
メモリキャッシュに絡む最適化は難しいですね...
Mac OS X 10.3.3 OpenGL 変更点
Mac OS X 10.3.3 アップデートがリリースされました。
Mobility Radeon 9600 上での変更点は
拡張に GL_ARB_point_parameters が追加されただけで、
OpenGL のバージョンが 1.3 から 1.4 になっただけでした。
GLSL はまだまだ Mac にくるのは遠そうですね...
当分は Cg を使った方がよさそうです。
ヘッダファイルにも GLSL 関連の拡張の記述がありませんでしたし。
2004年03月15日
Linear Pixel Shuffling
Linear Pixel Shuffling は、画像のピクセルをプログレッシブかつ一様に擬似乱数的にサンプリングしていく手法です。
この手法は Peter G. Anderson 博士により発案されたとのこと。
http://www.cs.rit.edu/~pga/abstracts.html
http://www.kcg.ac.jp/kcg/35/lps.html
同氏は j 言語
と呼ばれるなんとも不思議な言語を作ってもいるようです。
Linear Pixel Shuffing については、同氏によりいくつか論文が発表されていますが、最初の導入には、
Peter G. Anderson,
Linear Pixel Shuffling for Image Processing, an Introduction
Electron. Imag. 2: 147-154(1993).
が参考になります。
よりすすんで解説したものについては、
Peter G. Anderson, Jonathan S. Arney, and Kevin Ayer
Linear Pixel Shuffling (I): New Paradigms for New Printers
Presented at the IS&T Conference, on Non-Impact Printing, Vancouver, 2001
などが参考になります。
基本は フィボナッチ数的な考えでピクセルをサンプリングしていくというもののようです。
(フィボナッチ格子点 の生成と似ているようです)
ディゾルブ効果や誤差拡散、サブピクセルサンプリングなどの、
ピクセルをランダム(かつ一様に)に選びたいが、
そのときに同じピクセルを二度サンプルしないようにしたい、
というような手法に使えます。
プログレッシブに一様にサンプルしていく、
つまり 25% でサンプル生成を止めたら
画像全体のちょうど 1/4 の部分を一様に埋めるようになる、
という非常に有益な性質があります。
この手法はまた、論文でも述べられているように LDS とみなすことで QMC にも使えますし、
FPrime のようなプログレッシブレンダラにおいてサンプル点を生成するのにも使えます。
それ以外にも空間検索やニューラルネットワークなどの多岐にわたって応用があるという
スゴイ手法。しかも生成は簡単な式で求めることができます。
(なんか graphics gems で紹介されてそうなネタですね)
メトロポリス光輸送でのピクセルサンプリングの部分にも使えそうです。
2004年03月13日
virtualight
virtualight がフリーになるようです。
なんかもともとフリー的なソフトウェアだった気が...
virtualight のホームページの機能リストを見ると、最終最新リリースは 2002 年 10 月と 1 年半ぐらい前ですが、
すでに何気に結構先端的なレンダリング技術がいろいろと実装されています。
グローバルイルミネーション、プログラマブルシェーダ、準モンテカルロなどなど。
しかし、virtualight がフリーになるということは結構どうでも良くて、驚愕したのはフリー版 virtualight の説明文。
今までのレンダリング機能よりもさらにすすんだものを提供します。... Metropolis Light Transport ...
な、なんと...
メトロポリス光輸送は結構 lucille の売り込みネタになるかと思っていましたが、すでにやられていたとは...
やっぱり世界って広いわ。
きっとこれからの各種レンダラ戦争(?)で、フォトンマップの次にくるレンダリングアルゴリズムはメトロポリス光輸送に違いない。
2004年03月12日
情報理論とグラフィックス
Entropy-based adaptive sampling を読み直しています。
この著者のホームページを訪れて知ったのですが、著者である Jaume Rigau 氏は、情報理論(information-theory)をコンピュータグラフィックスへ応用する論文を数多く発表されています。
情報理論のグラフィックスへの応用については、同ホームページからもダウンロードできる、
M. Sbert, M. Feixas, J. Rigau, F. Castro, and P.P. Vázquez.
Applications of Information Theory to Computer Graphics.
Proceedings of 5th International Conference on Computer Graphics and Artificial Intelligence, 3IA'02 (Limoges, France), pp. 21-36, 2002.
が概観をつかむものとして参考になります。
いくつか応用例を挙げると、情報理論において情報量を測る尺度であるエントロピー(entropy, wikipedia 参照) を用いて、
o レイトレにおける適応サンプリングによるサブピクセルサンプリング
o フォームファクタ計算のための遮蔽計算
などの興味深い話題を取り扱う方法が述べられています。
(また 3Dモデルの最適な視点を見つけるのにエントロピーを用いるという面白い例も示されています。
これはゲームなどでカメラの視点をオートマチックに動かすのに使えたりできそうです)
リファインの基準
サブピクセルサンプリングや遮蔽の計算(予測?)を効率的に行いたい場合、
そのリファインを決定するための基準となる物差しに何を使うかが問題になります。
つまり、" 1 ピクセル 1 レイサンプルの結果ではこれくらい情報量(画像品質)が足りないから、
あと何レイサンプル必要だ " というような評価をどう求めるのか? ということです。
サブピクセルサンプリングの既存研究では、
例えばさらにサブピクセル領域を細かくサンプリングを行うのを決定するための基準として、
ピクセルサンプルのコントラストや Z 値を用いるというのがありますが、
あくまでこれらの尺度は経験上うまくいく定量であって、
これを頑健な尺度、つまりどんな入力の状況に対してもウソを付かないような尺度、
を求めたいということがあります。
情報理論で用いられている情報量を測る尺度であるエントロピーを用いる事で、
頑健な尺度かどうかは分りませんが、しかし
既存の手法よりもより良い結果が得られているようです。
エントロピーの計算自体は log を取るぐらいの簡単な式で求められます。
また Rigau 氏は、エントロピーと似た f-情報量(f-divergence) と呼ばれる、
これもまた情報理論で用いられる尺度を用いて
リファインの基準の計算に用いる論文を発表されています。
J. Rigau, M. Feixas, and M. Sbert.
Refinement Criteria Based on f-Divergences.
Rendering Techniques 2003 (14th Eurographics Workshop on Rendering, Leuven, Belgium, June 25-27), pp. 260-269.
Per H. Christensen and Daniel Cohen-Or (eds.), Association for Computing Machinery, New York (NY), USA, 2003, ISBN 1-58113-754-0.
エントロピーのような情報理論で用いられている情報量を測る尺度のグラフィックスへの応用というのは
いろいろとありそうでおもしろそうです。
特に、QMC やグローバルイルミネーションとの連携がこの情報理論とうまくくっつきやすい気がします。
情報理論はサンプリング、数論、LDS とも関わりがあるようですし。
あと、エントロピー使ってシェーダアンチエイリアシングとかできそうだな。
GLSL とかだと log() もハードウェアでサポートされているし。
2004年03月10日
SBR 2004 設営中
SBR 2004( SIGGRAPH 論文実装レース 2004) の参加者を受け付けるッ!!!
SBR 2004 のサイトは以下を予定しています。
http://lucille.sourceforge.net/cgi-bin/sbr2004/wiki.cgi
plone も便利ですが、wiki も便利ですね。
2004年03月09日
Efficient Illumination by High Dynamic Range Images テストレンダリング
Efficient Illumination by High Dynamic Range Images の実装のテストレンダリング画像。
シェーディング点あたり 128 サンプル, 512x512, 1x1 アンチエイリアシングで約 1 分。
オラオラ!! テメ〜、実装しろやオラ!!
Sir!! yes Sir!!!
てな感じで実装させたやつ。
サンプル点の生成は 128 個で約 1 分、 300 個で約 2 分 30 秒なので、
品質、サンプル点の生成速度のどちらも Structured Importance Sampling よりもいいような気がします。
実装の詳細や、ソースはもう少ししたら実装した人に説明してもらおうと思います。
ちなみにテストレンダリングに用いたレンダラの名前は Caesar(カエサル) だそうです。
lucille を参考にしてくださっているのですが、
lucille よりも高速っぽいのでこっちも結構がんばらないとなぁ。
Efficient Illumination by High Dynamic Range Images ではまた、
サンプリングが決定的であるので、
サンプル点をサブピクセル数分のレイヤに分けて、効率を挙げる手法も提案されています。
これは、まずサンプル点をレイヤに分けます。
このとき、各レイヤのサンプル点の分布をはなるべくもとのサンプル点の分布を保持したままにします。
個々のサブピクセル毎に各レイヤを割り当てることで、各シェーディング点で用いるサンプル点を減らすということです。
たとえば、サブピクセル数(アンチエイリアシング)を 2x2 = 4 とします。
もとのサンプル数 128 を、 32 個のサンプル数を持つ 4 つのレイヤを作成します。
各サブピクセルでは、 32 サンプルでシェーディングを行います。
結果としては ピクセルあたり 128 サンプルになりますが、
各サブピクセルでのシェーディングでは 32 サンプルしか用いないことになります。
単純に考えて毎サブサンプルあたり 128 サンプルするのに比べて 4 倍高速になります。
SIS でもこの手法が使えるかと思います。
ちなみに私は今は少佐です。全国ランキング 80 位ぐらいか...
2004年03月08日
リレンダリング、リライティング
FPrime というレンダラを知りました。
モデラーと連携して高速にプレビュー、リレンダリングしていくインターフェイスは素晴らしいですね。
(ところで G2 はどうなっちゃうのでしょうか...)
というわけで、リレンダリング(re-rendering, 再レンダリング)、リライティング(re-lighting, 再ライティング)ネタについて。
ここらへんの研究も多く行われているようです。
多くの場合、我々はレンダリングのプレビューを行う時に、サーフェスのプロパティ(反射率や色)やライトのプロパティ(位置や色)は頻繁に変更しますが、視点は固定のままになるかと思います。
このとき、画像面(カメラ)から見えるサーフェスの情報は 1 度だけ計算しておけばよいので、ピクセルにそれらの情報を保存しておけば、後々の計算では隠面消去を繰り返す必要がなく、効率的なシェーディングの再計算役立つというのが分ります。
このような手法の基本はやはり G-buffer でしょうか。
Takafumi Saito and Tokiichiro Takahashi.
Comprehensible rendering of 3-d shapes.
In Computer Graphics Annual Conference Series 1990, pages 197–206, August 1990.
いわゆるフレームバッファに色以外の情報(法線などのジオメトリ情報)を保存しておくという手法です。
同じ意味での G-buffer って Brazil とかでも使われているんですね。
次は、多分 RenderMan のリレンダラ irma で使われている感じがする手法。
Reid Gershbein, Pat Hanrahan
A Fast Relighting Engine for Interactive Cinematic Lighting Design
SIGGRAPH 2000, Pages: 353 - 358, 2000.
より映画産業方面でのリライティング、リレンダリングについての論文。今になって読み返して知ったのですが、OpenGL 使ってるんですね。
主要な点は、OpenGL でのレンダリングに適したものになるように、
o シャドウマップによるライティング
o BRDF のテクスチャへの分解
でしょうか。
以下論文 abstract 日本語訳。
プロシージャルシェーダを用いる複雑なシーンを、インタラクティブに映画的なライティングデザインを行うための新しい手法について提案します。
画像に表示されるサーフェスのジオメトリと光学的な情報を保存するために、デープ-フレームバッファ(deep-framebuffer)を用います。
ジオメトリ情報は向きのある点の集合として表現され、光学的な情報は双方向反射分布関数(BRDF)として表現されます。
BRDF は、サーフェスの見た目が空間的に変化するプロシージャルに定義されたサーフェステクスチャ関数により生成されます。
デープフレームバッファ情報は、OpenGL グラフィックスパイプラインを用いたマルチパスアルゴリズムでレンダリングされます。
物理的に正確な反射モデルと、映画産業で用いられる非リアルな反射モデルの両方を取り扱えるように、BRDF を独立な要素に分解しグラフィックスハードウェアのライティングユニットとテクスチャユニットの両方に割り当てます。
似たような分解はライティング分布のコントロールにも用いられます。こららの技法を用いることで、ライティングの計算は既存の手法に比べて 2,500 倍高速に評価することができます。これは、ライティングが変化する場合に静的な環境では 20 フレームでレンダリングできることになります。このときシーンは百万のオブジェクトが含まれ、数十種類のプロシージャルに定義されたサーフェスプロパティと数十のライトで構成されています。
ところで、FPrime がプログラマブルシェーダに対応しているのかは分りませんが
(というかそもそも lightwave でプログラマブルシェーダって扱えるのかな)、
プログラマブルシェーダに対しても、コンパイラにより変数間の依存関係を抜き出しておき、
変更部分だけを再評価して高速にリレンダリングを行う手法である
specializing shaders(日本語訳はシェーダ分離化が妥当?)が提案されています。
Brian Guenter, Todd B. Knoblock, Erik Ruf.
specializing shaders.
SIGGRAPH '95. Pages: 343 - 350, 1995.
これは、なんか奇特な人により 全文日本語訳があったり します。
基本的に、上記すべての手法はローカルイルミネーションのみへの対応です。
グローバルイルミネーションに対応していて、かつプログラマブルシェーダにも対応っていうのはあったりするのでしょうか。
(FPrime を開発した会社はスタンフォード大とも繋がっているようで、
同社のファーレンダラサスクアッチの技術が SIGGRAPH 2003 に採択されていましたね。
もしかしたら SIGGRAPH 2004 の論文に FPrime の技術が採択されていたりするかな...)
ところで specializing shaders を出したのはマイクロソフトです。
なんでコンパイラ技術ネタがマイクロソフトから?と思っていたのですが、
まさかもしかして HLSL のコンパイラの作成はこの頃から始まっていたのでしょうか...
三次元ヒルベルト曲線順パフォーマンス 2
セルデータのメモリ構成をブロックに分けて、キャッシュの改善があるかテストしてみました。
/* main_block.c */
#include "mapping.c"
#ifdef WIN32
#include <windows.h>
#else
#include <sys/time.h>
#endif
#define CELLSIZE 128
#define WIDTH 4
#define BLKSIZE (CELLSIZE / WIDTH)
#define SHIFTSIZE 2 /* log(2, WIDTH) */
#define MASK (WIDTH - 1)
#define NITER 100
typedef float data_t[3];
double
get_seconds()
{
double seconds;
#ifdef WIN32
LARGE_INTEGER freq, timer;
QueryPerformanceFrequency(&freq);
QueryPerformanceCounter(&timer);
return (double)timer.QuadPart / (double)freq.QuadPart;
#else /* linux, MacOS-X */
struct timeval t;
gettimeofday(&t, NULL);
seconds = (double)t.tv_sec + (double)t.tv_usec / (1000 * 1000);
#endif
return seconds;
}
float
hilb_test(int *table, data_t *dat)
{
int i, j, k, n;
int index;
int blkindex;
int pos;
float ret = 0.0f;
/* Z->Y->X order */
for (i = 0; i < CELLSIZE; i++) {
for (j = 0; j < CELLSIZE; j++) {
for (k = 0; k < CELLSIZE; k++) {
blkindex = ((k >> SHIFTSIZE) *
BLKSIZE * BLKSIZE) +
((j >> SHIFTSIZE) * BLKSIZE) +
(i >> SHIFTSIZE);
index = (k & MASK) * WIDTH * WIDTH +
(j & MASK) * WIDTH + (i & MASK);
pos = table[blkindex *
WIDTH * WIDTH * WIDTH +
index];
ret += dat[pos][0] + dat[pos][1] +
dat[pos][2];
}
}
}
return ret;
}
float
xyz_test(data_t *dat)
{
int i, j, k;
float ret = 0.0f;
for (k = 0; k < CELLSIZE; k++) {
for (j = 0; j < CELLSIZE; j++) {
for (i = 0; i < CELLSIZE; i++) {
ret += dat[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][0] +
dat[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][1] +
dat[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][2];
}
}
}
return ret;
}
float
zyx_test_blocked(data_t *dat)
{
int i, j, k;
int index;
int blkindex;
float ret = 0.0f;
for (i = 0; i < CELLSIZE; i++) {
for (j = 0; j < CELLSIZE; j++) {
for (k = 0; k < CELLSIZE; k++) {
blkindex = ((k >> SHIFTSIZE) *
BLKSIZE * BLKSIZE) +
((j >> SHIFTSIZE) * BLKSIZE) +
(i >> SHIFTSIZE);
index = (k & MASK) * WIDTH * WIDTH +
(j & MASK) * WIDTH + (i & MASK);
ret += dat[blkindex * WIDTH * WIDTH * WIDTH +
index][0] +
dat[blkindex * WIDTH * WIDTH * WIDTH +
index][1] +
dat[blkindex * WIDTH * WIDTH * WIDTH +
index][2];
}
}
}
return ret;
}
void
test(int num)
{
int i, j, k;
int blkindex;
int index;
float ret;
int *table = NULL;
data_t *dat1D = NULL;
data_t *dat3D = NULL;
data_t *dat3Dblk = NULL;
double start, end;
Point p, pdec;
Hcode h;
printf("Generating data... ");
fflush(stdout);
table = (int *)malloc(sizeof(int) * CELLSIZE * CELLSIZE * CELLSIZE);
dat1D = (data_t *)malloc(
sizeof(data_t) * CELLSIZE * CELLSIZE * CELLSIZE);
dat3D = (data_t *)malloc(
sizeof(data_t) * CELLSIZE * CELLSIZE * CELLSIZE);
dat3Dblk = (data_t *)malloc(
sizeof(data_t) * CELLSIZE * CELLSIZE * CELLSIZE);
/* build test data */
for (k = 0; k < CELLSIZE; k++) {
for (j = 0; j < CELLSIZE; j++) {
for (i = 0; i < CELLSIZE; i++) {
p.hcode[0] = i;
p.hcode[1] = j;
p.hcode[2] = k;
blkindex = ((k >> SHIFTSIZE) *
BLKSIZE * BLKSIZE) +
((j >> SHIFTSIZE) * BLKSIZE) +
(i >> SHIFTSIZE);
index = (k & MASK) * WIDTH * WIDTH +
(j & MASK) * WIDTH + (i & MASK);
h = H_encode(p);
table[blkindex * WIDTH * WIDTH * WIDTH +
index] = h.hcode[0];
dat1D[h.hcode[0]][0] = (float)k;
dat1D[h.hcode[0]][1] = (float)j;
dat1D[h.hcode[0]][2] = (float)i;
dat3D[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][0] = (float)i;
dat3D[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][1] = (float)j;
dat3D[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][2] = (float)k;
blkindex = ((k >> SHIFTSIZE) *
BLKSIZE * BLKSIZE) +
((j >> SHIFTSIZE) * BLKSIZE) +
(i >> SHIFTSIZE);
index = (k & MASK) * WIDTH * WIDTH +
(j & MASK) * WIDTH + (i & MASK);
dat3Dblk[blkindex * WIDTH * WIDTH * WIDTH +
index][0] = (float)i;
dat3Dblk[blkindex * WIDTH * WIDTH * WIDTH +
index][1] = (float)j;
dat3Dblk[blkindex * WIDTH * WIDTH * WIDTH +
index][2] = (float)k;
}
}
}
printf("OK.\n");
if (num == 1) {
ret = 0;
start = get_seconds();
for (i = 0; i < NITER; i++) {
ret += hilb_test(table, dat1D);
}
end = get_seconds();
printf("--- Hilbert curve encoded ---\n");
printf("traverse time = %f\n", (end - start) / (double)NITER);
printf("ret = %f\n", ret);
} else if (num == 2) {
ret = 0;
start = get_seconds();
for (i = 0; i < NITER; i++) {
ret += xyz_test(dat3D);
}
end = get_seconds();
printf("--- normal array(X->Y->Z order) ---\n");
printf("traverse time = %f\n", (end - start) / (double)NITER);
printf("ret = %f\n", ret);
} else {
ret = 0;
start = get_seconds();
for (i = 0; i < NITER; i++) {
ret += zyx_test_blocked(dat3Dblk);
}
end = get_seconds();
printf("--- blocked array(Z->Y->X order) ---\n");
printf("traverse time = %f\n", (end - start) / (double)NITER);
printf("ret = %f\n", ret);
}
}
int
main(int argc, char **argv)
{
int n;
if (argc < 2) {
printf("usage: %s N\n", argv[0]);
printf("N = 1: test Hilbert curve encoded traverse.\n");
printf("N = 2: test X -> Y -> Z order traverse.\n");
printf("N = 3: test Z -> Y -> X order traverse.\n");
exit(-1);
}
n = atoi(argv[1]);
test(n);
}
ヒルベルト曲線順(Z->Y->X オーダー相当)の table と、
Z->Y->X オーダーの通常配列に対して、
128x128x128 セルを4x4x4 のブロックごとに分けて
メモリ配置するように変更しています。
実行結果。
Macintosh% gcc -O2 main_block.c Macintosh% ./a.out 1 Generating data... OK. --- Hilbert curve encoded --- traverse time = 0.385858 ret = 39864975360.000000 Macintosh% ./a.out 2 Generating data... OK. --- normal array(X->Y->Z order) --- traverse time = 0.110044 ret = 39864975360.000000 Macintosh% ./a.out 3 Generating data... OK. --- blocked array(Z->Y->X order) --- traverse time = 0.447105 ret = 39864975360.000000
通常配列の X->Y->Z がぶっちぎりなのは置いておいて、
前回での Z->Y->X オーダーのヒルベルト曲線順が
traverse time = 0.851813
、通常の Z->Y->X オーダーが
traverse time = 0.784434
でしたので、ヒルベルト曲線順(Z->Y->X)、通常配列(Z->Y-X)のどちらでも 2 倍程度の向上が見られました。
ブロック化後の配列に対して Z->Y->X オーダーでアクセスするよりも、
ヒルベルト曲線順で Z->Y->X オーダーでアクセスする方が効率的なようです。
主要なものを早い順にまとめてみると、
1. 通常配列 X->Y->Z
2. ヒルベルト曲線順 Z->Y->X(ブロック化後)
3. 通常配列 Z->Y->X(ブロック化後)
になります。なので table 参照によるヒルベルト曲線順でもそれなりに効率的にふるまうようです。
しかしこれは、
CPU のタイプやキャッシュメモリの大きさ、コンパイラの最適化にも関係することなので、
今回の効率の結果が必ずしも絶対であるとは云えません。
あとは実際にレイトレのトラバースに組み込んでテストしてみる必要がありますね。
今の lucille では空間構造のポリゴンデータを断続してメモリに配置しているので、
一次元にまとめて配置しなおすなどデータ構造をキャッシュフレンドリーにするだけでも
効果が出そうです。
2004年03月07日
本到着
ゴゴゴゴゴゴゴ... ゴゴゴゴ... ゴゴ...
某所から本が届きました。なんだろう?
ゴゴゴゴゴゴゴ... ゴゴゴゴ... ゴゴ... ゴゴゴ...
一体、一体なんだと云うのだァ〜〜〜〜!!!!!!!!!!!!!!!!!!!!!!
落ち着け、落ち着くんだ。素数を数えて落ち着くんだ。2... 3.. 5... 7... 9.. いや違う、 11...
うおおおおおおおおおおおおおおおおおおおお!!!!
ドォ〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜ン!!!!!
バア〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜ン!!!
実践ッ!! CG へのッ!!! 誘いィィィ〜〜!!!
謹呈... うれしすぎて涙ッ!!!
絶版により、もう一生手に入れることができないだろうと思っていた本が、出版社の好意で入手することができました。
書籍って形はやっぱりいいですね。ホントいい本は一生ものです。
これを糧に、がんばって RenderMan レンダラの実装していって恩を返していきたいと思います。
おまけ。
今まで使っていた全コピしてファイリングしていたやつ。
ヒルベルト曲線 追記
hacker's delight にも、一章を割いてヒルベルト曲線についての記述があります。
hacker's delight では 2 次元でのヒルベルト曲線についての演算方法がメインに記述されています。
章の最後に書かれているヒルベルト曲線の応用の節では、やはりコンピュータグラフィックスへの適用について述べられています。
画像の圧縮と、レイトレ時のスキャンオーダーです。
今まではレイトレの一様グリッドやオクツリーのトラバースに三次元ヒルベルト曲線順でやるとパフォーマンスがあがるかなーと注目していたのですが、hacker's delight の方では、レイを投げるときのピクセルのスキャンオーダー(二次元ヒルベルト曲線)について述べられています。
通常視点からの一次レイを投げるときにはスキャンラインのように、左から右、上から下とピクセルを辿って投げることになりますが、このときに画像面のピクセルをヒルベルト順にスキャンしていって投げるように変更すればシーンデータベースへのヒットの局所性が保たれてページング(キャッシュ)が有効にできるよ、というもの。
ヒルベルト曲線順でのポリゴンラスタライズ( ヒルベルト曲線順とメモリキャッシュ 参照) と同じ考えですね。
でもこれは、間接反射がなく一次レイが優勢なローカルイルミネーション系のレイトレ(古典的なレイトレやリアルタイムレイトレ)でのみ効果がありそうです。
2004年03月06日
三次元ヒルベルト曲線順パフォーマンス
以下のテストプログラムを作成して、ヒルベルト曲線順のパフォーマンスを計測してみました。
#include "mapping.c"
#ifdef WIN32
#include <windows.h>
#else
#include <sys/time.h>
#endif
#define CELLSIZE 128
#define NITER 100
typedef float data_t[3];
double
get_seconds()
{
double seconds;
#ifdef WIN32
LARGE_INTEGER freq, timer;
QueryPerformanceFrequency(&freq);
QueryPerformanceCounter(&timer);
return (double)timer.QuadPart / (double)freq.QuadPart;
#else /* linux */
struct timeval t;
gettimeofday(&t, NULL);
seconds = (double)t.tv_sec + (double)t.tv_usec / (1000 * 1000);
#endif
/* todo : MacOS-X version code here */
return seconds;
}
float
hilb_test(int *table, data_t *dat)
{
int i, j, k, n;
int index;
float ret = 0.0f;
/* Z->Y->X order */
for (i = 0; i < CELLSIZE; i++) {
for (j = 0; j < CELLSIZE; j++) {
for (k = 0; k < CELLSIZE; k++) {
index = table[i * CELLSIZE * CELLSIZE +
j * CELLSIZE + k];
ret += dat[index][0] + dat[index][1] +
dat[index][2];
}
}
}
return ret;
}
float
xyz_test(data_t *dat)
{
int i, j, k;
float ret = 0.0f;
for (k = 0; k < CELLSIZE; k++) {
for (j = 0; j < CELLSIZE; j++) {
for (i = 0; i < CELLSIZE; i++) {
ret += dat[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][0] +
dat[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][1] +
dat[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][2];
}
}
}
return ret;
}
float
zyx_test(data_t *dat)
{
int i, j, k;
float ret = 0.0f;
for (i = 0; i < CELLSIZE; i++) {
for (j = 0; j < CELLSIZE; j++) {
for (k = 0; k < CELLSIZE; k++) {
ret += dat[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][0] +
dat[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][1] +
dat[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][2];
}
}
}
return ret;
}
void
test(int num)
{
int i, j, k;
float ret;
int *table = NULL;
data_t *dat1D = NULL;
data_t *dat3D = NULL;
double start, end;
Point p, pdec;
Hcode h;
table = (int *)malloc(sizeof(int) * CELLSIZE * CELLSIZE * CELLSIZE);
dat1D = (data_t *)malloc(
sizeof(data_t) * CELLSIZE * CELLSIZE * CELLSIZE);
dat3D = (data_t *)malloc(
sizeof(data_t) * CELLSIZE * CELLSIZE * CELLSIZE);
/* build test data */
for (k = 0; k < CELLSIZE; k++) {
for (j = 0; j < CELLSIZE; j++) {
for (i = 0; i < CELLSIZE; i++) {
p.hcode[0] = k;
p.hcode[1] = j;
p.hcode[2] = i;
h = H_encode(p);
table[p.hcode[2] * CELLSIZE * CELLSIZE +
p.hcode[1] * CELLSIZE + p.hcode[0]] =
h.hcode[0];
dat1D[h.hcode[0]][0] = (float)k;
dat1D[h.hcode[0]][1] = (float)j;
dat1D[h.hcode[0]][2] = (float)i;
dat3D[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][0] = (float)i;
dat3D[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][1] = (float)j;
dat3D[k * CELLSIZE * CELLSIZE +
j * CELLSIZE + i][2] = (float)k;
}
}
}
if (num == 1) {
ret = 0;
start = get_seconds();
for (i = 0; i < NITER; i++) {
ret += hilb_test(table, dat1D);
}
end = get_seconds();
printf("--- Hilbert curve encoded ---\n");
printf("traverse time = %f\n", (end - start) / (double)NITER);
printf("ret = %f\n", ret);
} else if (num == 2) {
ret = 0;
start = get_seconds();
for (i = 0; i < 100; i++) {
ret += xyz_test(dat3D);
}
end = get_seconds();
printf("--- normal array(X->Y->Z order) ---\n");
printf("traverse time = %f\n", (end - start) / (double)NITER);
printf("ret = %f\n", ret);
} else {
ret = 0;
start = get_seconds();
for (i = 0; i < 100; i++) {
ret += zyx_test(dat3D);
}
end = get_seconds();
printf("--- normal array(Z->Y->X order) ---\n");
printf("traverse time = %f\n", (end - start) / (double)NITER);
printf("ret = %f\n", ret);
}
}
int
main(int argc, char **argv)
{
int n;
if (argc < 2) {
printf("usage: %s N\n", argv[0]);
printf("N = 1: test Hilbert curve encoded traverse.\n");
printf("N = 2: test X -> Y -> Z order traverse.\n");
printf("N = 3: test Z -> Y -> X order traverse.\n");
exit(-1);
}
n = atoi(argv[1]);
test(n);
}
結果。
Macintosh% gcc -O2 main.c Macintosh% ./a.out 1 --- Hilbert curve encoded --- traverse time = 0.229896 ret = 39864975360.000000 Macintosh% ./a.out 2 --- normal array(X->Y->Z order) --- traverse time = 0.117733 ret = 39864975360.000000 Macintosh% ./a.out 3 --- normal array(Z->Y->X order) --- traverse time = 0.784434 ret = 39864975360.000000
通常の三次元配列での X->Y->Z 順のトラバースは最もキャッシュヒットが高くなり、
通常の三次元配列での Z->Y->X 順のトラバースは最もキャッシュミスが発生しやすくなっているはずです。
ヒルベルト曲線順では、Z->Y->X 順のセルのトラバース相当を行っています。
明らかに Z->Y->X 順でトラバースするよりはヒルベルト順でトラバースしたほうが効率は良さそうで、
テストでは 4 倍程度高速になっています。
しかし、xyz の位置からヒルベルト曲線順で一次元に再配置されたインデックスを求めるのにテーブル参照を
行っているので、結局計時にはそのテーブルへのメモリアクセスも含まれることになります。
試しに hilb_test() でループの順番を i -> j -> k から k-> j -> i に変更すると、
traverse time = 0.851813
と最悪になってしまいます(table のメモリアクセスがキャッシュミスを起こして時間がかかっている)。
これは、xyz の位置から直接高速にインデックスを求める関数とかあると解決できるのですが、
そのような関数はあるのでしょうか?
もしくは table を数値計算での行列計算の時のようにブロックに分けることでも
table へのメモリアクセスの効率化が図れるかもしれません。
2004年03月05日
Algorithms for Division Free Perspective Correct Rendering
パースペクティブ補正を、除算なしで行うというもの。
B. Barenbrug, F.J. Peters, C.W.A.M. van Overveld,
Algorithms for Division Free Perspective Correct Rendering
SIGGRAPH Eurographics Graphics Hardware Workshop, pp.7-13, 2000.
中点アルゴリズム(midpoint algorithm)を基本とするもののようで、除算するものに比べて、 2 倍程度高速になるようです。除算なしなので、除算器がないハードウェアでも動くというものです。
ネタ的にはインパクトがあり有益なのですが、しかし最近では携帯などの組み込み CPU でも除算器が搭載されつつあり、また携帯に ハードウェア 3D エンジンも載るくらいの御時世なのですよね...
自分で 3D エンジンを書くときに使えるかも。でもそれぐらいしか役えないかなぁ(今では ATI か NVIDIA の GPU ではラスタライズにこの方法が用いられていたりするのでしょうか)。
以下論文 abstract 日本語訳
平面ポリゴンのパースペクティブ補正レンダリングの良く知られた実装では、ピクセル毎に除算が必要になります。このような除算は、シリコンゲートやクロックサイクルの観点からみて高価な演算であるので、回避されるのが望ましいでしょう。本論文では、除算を必要としない効率的な中点アルゴリズム(midpoint algorithms)の一種を提案します。
これらのアルゴリズムはピクセル毎に、いくつかの加算演算以外には何も必要としません。
これらの手法がスキャンラインアルゴリズムやミップマップのアルゴリズムに組み込めることを示します。
ソフトウェア実装による実験では、ポリゴンが非常に小さくないときに除算なしアルゴリズムは 2 倍高速になりました。
これらのアルゴリズムはハードウェアで実装されたときにより有益になるでしょう。
三次元ヒルベルト曲線順インデックス生成
ヒルベルト曲線順とメモリキャッシュ からの続きで、レイトレのトラバースなど、グラフィックスで有効でありそうな三次元ヒルベルト曲線順について。
三次元ヒルベルト曲線の生成方法にはいろいろあるようです。
また、その生成される三次元ヒルベルト曲線の構造は実に 1536 パターンの種類 [1] があるとのこと。
実際にはこの構造パターンの違いというのは我々の場合には重要ではなく、三次元のセルの位置をヒルベルト曲線順で一次元の配列にインデックス化できれば十分です。
というわけで、ソースコードが公開されている多次元ヒルベルト曲線のコードを見つけました(他にも探すといろいろと見つかります)。
J.K.Lawder. Calculation of Mappings Between One and n-dimensional Values Using the Hilbert Space-filling Curve. Research Report BBKCS-00-01 (formerly JL1/00), August 2000
このコードで、三次元、4x4x4 セルのヒルベルト曲線順でのインデックスを求めるサンプルプログラムを組んでみました。
/* main.c */
#include "mapping.c"
#define CELLSIZE 4
Hcode table[CELLSIZE * CELLSIZE * CELLSIZE];
int
main(int argc, char **argv)
{
int i, j, k;
Point p, pdec;
Hcode h;
for (k = 0; k < CELLSIZE; k++) {
for (j = 0; j < CELLSIZE; j++) {
for (i = 0; i < CELLSIZE; i++) {
p.hcode[0] = i;
p.hcode[1] = j;
p.hcode[2] = k;
h = H_encode(p);
table[h.hcode[0]] = p;
}
}
}
for (i = 0; i < CELLSIZE * CELLSIZE * CELLSIZE; i++) {
p = table[i];
printf("%d : (%d, %d, %d)\n",
i, p.hcode[0], p.hcode[1], p.hcode[2]);
}
}
実行結果。
$ gcc main.c $ ./a.out 0 : (0, 0, 0) 1 : (0, 1, 0) 2 : (1, 1, 0) 3 : (1, 0, 0) 4 : (1, 0, 1) 5 : (1, 1, 1) 6 : (0, 1, 1) 7 : (0, 0, 1) 8 : (0, 0, 2) 9 : (1, 0, 2) 10 : (1, 0, 3) 11 : (0, 0, 3) 12 : (0, 1, 3) 13 : (1, 1, 3) 14 : (1, 1, 2) 15 : (0, 1, 2) 16 : (0, 2, 2) 17 : (1, 2, 2) 18 : (1, 2, 3) 19 : (0, 2, 3) 20 : (0, 3, 3) 21 : (1, 3, 3) 22 : (1, 3, 2) 23 : (0, 3, 2) 24 : (0, 3, 1) 25 : (0, 3, 0) 26 : (0, 2, 0) 27 : (0, 2, 1) 28 : (1, 2, 1) 29 : (1, 2, 0) 30 : (1, 3, 0) 31 : (1, 3, 1) 32 : (2, 3, 1) 33 : (2, 3, 0) 34 : (2, 2, 0) 35 : (2, 2, 1) 36 : (3, 2, 1) 37 : (3, 2, 0) 38 : (3, 3, 0) 39 : (3, 3, 1) 40 : (3, 3, 2) 41 : (2, 3, 2) 42 : (2, 3, 3) 43 : (3, 3, 3) 44 : (3, 2, 3) 45 : (2, 2, 3) 46 : (2, 2, 2) 47 : (3, 2, 2) 48 : (3, 1, 2) 49 : (2, 1, 2) 50 : (2, 1, 3) 51 : (3, 1, 3) 52 : (3, 0, 3) 53 : (2, 0, 3) 54 : (2, 0, 2) 55 : (3, 0, 2) 56 : (3, 0, 1) 57 : (3, 1, 1) 58 : (2, 1, 1) 59 : (2, 0, 1) 60 : (2, 0, 0) 61 : (2, 1, 0) 62 : (3, 1, 0) 63 : (3, 0, 0)
とうまく生成できているようです。
ここらへん OpenGL でサクッと可視化するプログラムが作れるといいですね。
Homage to Hilbert
三次元ヒルベルト曲線で思い出したのが、 SIGGRAPH '98 のエレクトロニックシアターのムービー
当時見た時はこれがヒルベルト曲線だったとは知りませんでしたが、なんかスゲームービーだなーと感じていました。
(Video math festival '98 でも採択されているようです。homage to Hilbert は含まれていませんが、Video math festival のトレーラーを こちら で見る事が出来ます)
ってこれって OpenGL でレンダリングだったんですね。しかも作ったのは Nelson Max かよ!!
こういうのを見ると思うのは、数学的な事柄はビジュアライズすると本当に理解しやすいということです。レンダリングというプロセスもある意味数学的なことのビジュアライズですね。
[1] J. Alber and R. Niedermeier. On multi-dimensional Hilbert indexings. In Proceedings of the 4th Annual International Computing and Combinatorics Conference (COCOON), number 1449 in Lecture Notes in Computer Science, pages 329-338. Springer-Verlag, 1998
2004年03月04日
浮動小数点計算誤差の推定 2
前回 のコードをレイと三角形の交差判定に組み込んでみました。
/* * Estimates catastrophic cancellation of floating point arithmetrics. * * This is a implementation of * * Hiroshi Suzuki, Hajime Ohiwa, * "A program for estimating catastrophic cancellation of floating-point * arithmetrics", * IPSJ SIGNotes High Performance Computing, No.053, 1994. * http://www.ipsj.or.jp/members/SIGNotes/Eng/12/1994/053/index.html * * * Note that program might be wrong... */ #includestatic int estimate(float a, int aDigit, float b, int bDigit, float c) { int validDigit; unsigned int *xp1; unsigned int *xp2; unsigned int *xp3; int xe1; // exponent part int xe2; // exponent part int xe3; // exponent part int lack; // number of lacked digits xp1 = (unsigned int *)&a; xp2 = (unsigned int *)&b; xp3 = (unsigned int *)&c; // extract expornent part xe1 = (*xp1 >> 23) & 0xff; xe2 = (*xp2 >> 23) & 0xff; xe3 = (*xp3 >> 23) & 0xff; if (xe1 == xe2) { if (aDigit < bDigit) { validDigit = aDigit; } else { validDigit = bDigit; } lack = xe1 ^ xe3; validDigit -= lack; if (validDigit < 0) { validDigit = 0; } } else { if (xe1 > xe2) { int biased; lack = xe1 ^ xe3; biased = xe1 > 127 ? xe1 - 127 : 0; biased += bDigit; if (aDigit < biased) { validDigit = aDigit; } else { validDigit = biased; } validDigit -= lack; if (validDigit < 0) { validDigit = 0; } } else { int biased; lack = xe2 ^ xe3; biased = xe2 > 127 ? xe2 - 127 : 0; biased += aDigit; if (bDigit < biased) { validDigit = bDigit; } else { validDigit = biased; } validDigit -= lack; if (validDigit < 0) { validDigit = 0; } } } return validDigit; } class FloatEst { public: FloatEst() { val = 0.0f; validDigit = 23; } FloatEst(float f) { val = f; validDigit = 23; } ~FloatEst() {}; FloatEst &operator =(float f) { val = f; validDigit = 23; return *this; } FloatEst operator +(const FloatEst &f) const { FloatEst ret; float z; z = this->val + f.val; ret.validDigit = estimate(this->val, this->validDigit, f.val, f.validDigit, z); ret.val = z; return ret; } FloatEst operator -(const FloatEst &f) const { FloatEst ret; float z; z = this->val - f.val; ret.validDigit = estimate(this->val, this->validDigit, f.val, f.validDigit, z); ret.val = z; return ret; } FloatEst operator *(const FloatEst &f) { FloatEst ret; /* no valid digits are lacked. choose lesser one */ if (this->validDigit < f.validDigit) { ret.validDigit = this->validDigit; } else { ret.validDigit = f.validDigit; } ret.val = this->val * f.val; return ret; } FloatEst operator /(const FloatEst &f) const { FloatEst ret; /* no valid digits are lacked. choose lesser one */ if (this->validDigit < f.validDigit) { ret.validDigit = this->validDigit; } else { ret.validDigit = f.validDigit; } ret.val = this->val / f.val; return ret; } float val; int validDigit; }; std::ostream &operator <<(std::ostream &s, const FloatEst &f); std::ostream &operator <<(std::ostream &s, const FloatEst &f) { return (s << "val = " << f.val << ", validDigits = " << f.validDigit); } /* ray - triangle intersection arithmetric */ static int isect(FloatEst org[3], FloatEst dir[3], FloatEst v0[3], FloatEst v1[3], FloatEst v2[3], FloatEst *t, FloatEst *u, FloatEst *v) { FloatEst one = 1.0f; FloatEst e1[3], e2[3], tvec[3], pvec[3], qvec[3]; FloatEst det, inv_det; e1[0] = v1[0] - v0[0]; e1[1] = v1[1] - v0[1]; e1[2] = v1[2] - v0[2]; e2[0] = v2[0] - v0[0]; e2[1] = v2[1] - v0[1]; e2[2] = v2[2] - v0[2]; pvec[0] = dir[1] * e2[2] - dir[2] * e2[1]; pvec[1] = dir[2] * e2[0] - dir[0] * e2[2]; pvec[2] = dir[0] * e2[1] - dir[1] * e2[0]; det = e1[0] * pvec[0] + e1[1] * pvec[1] + e1[2] * pvec[2]; std::cout << "det = " << det << std::endl; if (det.val > 1.0e-6f) { tvec[0] = org[0] - v0[0]; tvec[1] = org[1] - v0[1]; tvec[2] = org[2] - v0[2]; (*u) = tvec[0] * pvec[0] + tvec[1] * pvec[1] + tvec[2] * pvec[2] ; std::cout << "u = " << (*u) << std::endl; if ((*u).val < 0.0 || (*u).val > det.val) return 0; qvec[0] = tvec[1] * e1[2] - tvec[2] * e1[1]; qvec[1] = tvec[2] * e1[0] - tvec[0] * e1[2]; qvec[2] = tvec[0] * e1[1] - tvec[1] * e1[0]; (*v) = dir[0] * qvec[0] + dir[1] * qvec[1] + dir[2] * qvec[2]; std::cout << "v = " << (*v) << std::endl; if ((*v).val < 0.0 || (*u).val + (*v).val > det.val) return 0; } else { return 0; } inv_det = one / det; (*t) = e2[0] * qvec[0] + e2[1] * qvec[1] + e2[2] * qvec[2]; (*t) = (*t) * inv_det; (*u) = (*u) * inv_det; (*v) = (*v) * inv_det; std::cout << "t = " << (*t) << std::endl; std::cout << "u = " << (*t) << std::endl; std::cout << "v = " << (*t) << std::endl; return 1; } using namespace std; int main(int argc, char **argv) { FloatEst org[3], dir[3], v0[3], v1[3], v2[3], t, u, v; org[0] = 0.0f; org[1] = 0.0f; org[2] = 0.0f; dir[0] = 0.0f; dir[1] = 0.0f; dir[2] = 1.0f; v0[0] = 10.0f; v0[1] = -10.0f; v0[2] = 10.0f; v1[0] = -10.0f; v1[1] = -10.0f; v1[2] = 10.0f; v2[0] = 0.0f; v2[1] = 10.0f; v2[2] = 10.0f; isect(org, dir, v0, v1, v2, &t, &u, &v); }
結果は以下の通り。
det = val = 400, validDigits = 4 u = val = 100, validDigits = 1 v = val = 200, validDigits = 7 t = val = 10, validDigits = 4 u = val = 10, validDigits = 4 v = val = 10, validDigits = 4
多分重要なのは det, t, u, v の精度でしょう。
ということで、サンプルの演算
レイ: 原点 (0, 0, 0) -> 方向 (0, 0, 1)
三角形 (10, -10, 10), (-10, -10, 10), (0, 10, 10)
では 2 進数で 4 桁、つまり 10 進では 1 / 2^4 = 0.0625 の解像度が安全に表現できる下限ということになる...のでしょうか?
いくつかランダムにサンプルデータを生成してテストを回してみる必要もありますね。
2004年03月03日
浮動小数点計算誤差の推定
浮動小数点計算の誤差は、レイトレの三角形の交差判定コードを書いた時から気になっていましたが、 MLT を実装していくうちでやっぱりちゃんと対処していかないといけないと考えるようになりました。
まあ今まではそれほど数値計算の誤差が致命的になるような状況になかったので、適当に機械エプシロン
#define EPS 1.0e-6f
で乗り切れるようなことが多かったのですが、 MLT やパストレーシングでの距離の 2 乗の計算部分で浮動小数点値が膨大になったり微小になったりする(フォトンマップ本でもちょこっとここらへんが指摘しされており、power heuristic 法を使うとうまくいくとか)のが問題になりそうなのでいろいろと考えてみることに。
まあ double を使えばいいのでしょうが、せっかく lucille はレイトレの交差判定に SSE(float 精度 x 4 同時実行)を使っているので、
double を使うのは気が引けるし、SSE2(double 精度 x 2) は 2 命令同時実行なので SSE に比べればそれほどパフォーマンスが得られるわけでもないのでこれも気が引けるし...
交差判定の演算はもう決まっているので、float でも精度をうまく考慮すればうまくいけるのではないかと思うのです。
とりあえず float(浮動小数点) の計算で最も考慮すべきことは、
o 近い値の減算は大きな桁落ちを引き起こす
o float の場合、値間の開きは高々 1.2e-7(= 2^(-23)) の解像度
という感じでしょうか。たとえば (0, 0) と (100, 0) の距離の逆 2 乗が 1.0/100^2 = 1.0e-6 なので、 1.2e-7 の解像度というのは結構少ないのかもしれません(もっと float は精度がある印象だったのですが)。
誤差の推定
で、演算を行っていく時に、結局どれくらいの精度があればいいのか推定できる手法というのがないかと探したところ、それらしきものが。
鈴木 弘, 大岩 元
浮動小数点演算における精度見積りアルゴリズムとその評価
情報処理学会 研究報告 「ハイパフォーマンスコンピューティング」No.053, 1994.
この手法では、丸め誤差、情報落ち、桁落ちのうち、桁落ちのみを推定するものです。
論文を読んで初めて知ったのですが、乗算と除算では桁落ちって生じないのですね。
論文の手法では、C++ の演算子オーバーロードを用いて演算後の有効桁数を測ると云うものです。
といわけで(忌まわしき)C++ を久しぶりに使って(C++ は書かなくなるとホント忘れてしまいますね)とりあえずこの手法を実装してみました。
プログラムは、 2 進での演算後の浮動小数点の有効桁数を求めます。
main() 内のサンプル演算は、推定ではいくつか桁落ちが示されていますが、
実際には桁落ちが起こらない演算です。
あとはこのコードを用いて交差判定の演算を行ってみれば誤差の推定値が得られます。
論文での比較によると、推定値は、実際の精度よりも悪い値を返すようです(最悪の条件が判断できるのでそれはそれでいいかと)。
桁落ちが浮動小数点演算の誤差のすべてではありませんが、この推定値の算出は結構役に立つかもしれません。
/* * Estimates catastrophic cancellation of floating point arithmetrics. * * This is a implementation of * * Hiroshi Suzuki, Hajime Ohiwa, * "A program for estimating catastrophic cancellation of floating-point * arithmetrics", * IPSJ SIGNotes High Performance Computing, No.053, 1994. * http://www.ipsj.or.jp/members/SIGNotes/Eng/12/1994/053/index.html * * * Note that program might be wrong... */ #includestatic int estimate(float a, int aDigit, float b, int bDigit, float c) { int validDigit; unsigned int *xp1; unsigned int *xp2; unsigned int *xp3; int xe1; // exponent part int xe2; // exponent part int xe3; // exponent part int lack; // number of lacked digits xp1 = (unsigned int *)&a; xp2 = (unsigned int *)&b; xp3 = (unsigned int *)&c; // extract expornent part xe1 = (*xp1 >> 23) & 0xff; xe2 = (*xp2 >> 23) & 0xff; xe3 = (*xp3 >> 23) & 0xff; if (xe1 == xe2) { if (aDigit < bDigit) { validDigit = aDigit; } else { validDigit = bDigit; } lack = xe1 ^ xe3; validDigit -= lack; if (validDigit < 0) { validDigit = 0; } } else { if (xe1 > xe2) { int biased; lack = xe1 ^ xe3; biased = xe1 > 127 ? xe1 - 127 : 0; biased += bDigit; if (aDigit < biased) { validDigit = aDigit; } else { validDigit = biased; } validDigit -= lack; if (validDigit < 0) { validDigit = 0; } } else { int biased; lack = xe2 ^ xe3; biased = xe2 > 127 ? xe2 - 127 : 0; biased += aDigit; if (bDigit < biased) { validDigit = bDigit; } else { validDigit = biased; } validDigit -= lack; if (validDigit < 0) { validDigit = 0; } } } return validDigit; } class FloatEst { public: FloatEst() { val = 0.0f; validDigit = 23; } FloatEst(float f) { val = f; validDigit = 23; } ~FloatEst() {}; FloatEst &operator =(float f) { val = f; validDigit = 23; return *this; } FloatEst operator +(FloatEst &f) { FloatEst ret; float z; z = this->val + f.val; ret.validDigit = estimate(this->val, this->validDigit, f.val, f.validDigit, z); ret.val = z; return ret; } FloatEst operator -(FloatEst &f) { FloatEst ret; float z; z = this->val - f.val; ret.validDigit = estimate(this->val, this->validDigit, f.val, f.validDigit, z); ret.val = z; return ret; } FloatEst operator *(FloatEst &f) { FloatEst ret; /* no valid digits are lacked. choose lesser one */ if (this->validDigit < f.validDigit) { ret.validDigit = this->validDigit; } else { ret.validDigit = f.validDigit; } ret.val = this->val * f.val; return ret; } FloatEst operator /(FloatEst &f) { FloatEst ret; /* no valid digits are lacked. choose lesser one */ if (this->validDigit < f.validDigit) { ret.validDigit = this->validDigit; } else { ret.validDigit = f.validDigit; } ret.val = this->val / f.val; return ret; } float val; int validDigit; }; std::ostream &operator <<(std::ostream &s, const FloatEst &f); std::ostream &operator <<(std::ostream &s, const FloatEst &f) { return (s << "val = " << f.val << ", validDigits = " << f.validDigit); } using namespace std; int main(int argc, char **argv) { FloatEst a, b, c, d, e, f; a = 1.5655f; b = 1.5f; f = 16.0f; // val valid digits c = a - b; // 0.0625 19 cout << "c = " << c << endl; d = a - c; // 1.5 19 cout << "d = " << d << endl; e = a - d; // 0.0625 15 cout << "e = " << e << endl; e = f + e; // 16.0625 19 cout << "e = " << e << endl; }
「実践ＣＧへの誘い」（絶版本）の寄贈のご案内
急げ!!!!
レンダラ、RenderMan やるなら必携本。
これほどまでに絶版が残念でならない本はありません。
ちなみに私は図書館で全ページをコピーしたクチです。