Syoyo Fujita's Blog

raytracing monte carlo

Month: September, 2005

amd64 dual core 対応

これからの CPU のマルチコア化、64bit 化の流れを見据え、
lucille の開発環境に AMD のデュアルコア CPU を追加しました。

特に外部 GPU を利用する必要もないので、
今回は内臓グラフィックスを持つ、キューブベアボーンである ST20G5 に、
最近のお気に入りである debian(64bit の amd64 版)をインストールしました。

http://www.debian.org/ports/amd64/index.ja.html

現在リリースされているインストーラでは、ネットワークをうまく認識してくれなかったので、
daily のインストーラを利用しました。


http://cdimage.debian.org/pub/cdimage-testing/daily/amd64/current/

それでも、ACPI 関連ではカーネルパニックに、フレームバッファ関連では
ウィンドウの表示に失敗したので、起動オプションに以下を指定して回避させました。

# install acpi=off debian-installer/framebuffer=false
 
現在の linux では、kernel 2.6.12 以上でないと、AMD のデュアルコア CPU を、

CPU が二つあると認識してくれないようです。
(cat /proc/cpuinfo で processor が 2 個出てくると認識されていることになる)

64bit バイナリとして luculle をコンパイルしてみましたが、今のところ 64bit-ize しても
問題なく動いているようです。
その昔に G5 の 64bit モードで動くように、コードをコントリビュートしていただいた
ことがありましたので、それがまだまだ効いているのでしょう。

基本的には、64bit 化するときには、sizeof(void *) != sizeof(int) であることに
注意してコーディングしていれば OK だと思います。

パフォーマンス測定

chrysler シーンを、
o AMD X2 3800+(2.0GHz) 64bit + 2 threads(デュアルコアをフルに使う)
o AMD X2 3800+(2.0GHz) 64bit + single
thread(デュアルコアの片側のコアだけつかう)
o 既存の P4 32bit 2.4 GHz(シングルスレッド、ハイパースレッド非搭載)

でレンダリングした結果です。

o AMD DualCore 2.0GHz 2 threads: 49.7 sec
o AMD DualCore 2.0GHz single thread: 68.7 sec
o P4 2.4 GHz single thread: 112.55 sec

2 threads では、 single thread の約 38 % の高速化となりました。
理論値では処理速度は 2 倍が期待できますが、実際は 80 % の向上が見込めれば
十分だと思います。今はまだそれの半分の効率ですね。並列処理の部分は
まだまだアルゴリズムの改良、チューニングの余地がありそうです。

trac : Integrated SCM and Project Management

Subversion 実践入門
http://ssl.ohmsha.co.jp/cgi-bin/menu.cgi?ISBN=4-274-06613-4

を見ていたら、trac という SCM ツールがあることを知りました。

http://www.edgewall.com/trac/

wiki と subversion とバグトラッキングが一緒になっているのが特徴のようです。

インストールして試してみましたが、デザインもすっきりしていて、
なかなかよさそうでした。

subversion のコード表示では、リビジョン間での差分なども視覚的に表示できるので、
どこが変わったとかも web 上でみることができます。

自分でも lucille のコードを書いていて、どこがどう変わったのかも
覚えなくなっているときがあるので、もう少し使ってみて、使えそうならこの trac で
lucille のサイトからコードまでを一元管理できればと考えています。

グリッドハッシュ

lucille では主にレイトレの空間データ構造には一様グリッドを用いており、
一度にガバっと三次元ボクセル(たとえば voxel[128][128][128])という形でメモリを確保しています。

これだと、もっとも単純なのですが、物体の無いカラのボクセルでもメモリを
確保してしまうため、メモリを無駄に消費してしまうという欠点があります。

ひとつの解として、ハッシュを用いてカラのボクセルはメモリ確保しないという
方法があります。

グリッドハッシュは少し調べて試してみたことがあったのですが、
Arnold も空間データ構造にグリッドハッシュを利用していると聞き、
lucille でも、本格的に lucille に組み込んでデフォルトの空間データ構造に
しようかと考えています。

ボクセルハッシュ

視覚的な解説は、


http://www-evasion.imag.fr/Publications/2005/ZTKHRF05/spatial_subdivision.pdf

が参考になります。

通常のハッシュアルゴリズムとほとんど同じです。

まず、構築時では、物体をボクセルに登録します。
物体位置 (x, y, z) からハッシュ値を計算します。

hashval = hash(x, y, z) % HASH_TABLESIZE

この値でハッシュテーブルをルックアップします。

voxel = hashtable[hashval]

この voxel の場所にポリゴンを追加します。
しかし、ハッシュテーブルの数は有限で通常ボクセル分割数よりも少ないので、
異なる (x, y, z) でも同じハッシュ値と衝突してしまう可能性があります。
この場合はリンクリスト(チェイン)を作成します。
ボクセルには自身が対応する (x, y, z) の情報も保持しておきます。

found = false;
while (voxel) {
    if (voxel->x == x && voxel->y == y
&& voxel->z == z) {
        物体をこのボクセルに登録
        found = true;
        break;
    }

    voxel = voxel->next;
}

if (found!) {
    newvoxel = new voxel;
    物体と位置情報をこのボクセルに登録

    // リンクリストを作成
    newvoxel->next = hashtable[hashval];
    hashtable[hashval] = newvoxel;
}

トラバーサルのときも同じようにします。
まず 3DDDA でトラバース位置を計算、(x, y, z) からハッシュ値を計算し、
ハッシュテーブルからボクセルをルックアップします。
同じ (x, y, z) の値を持つボクセルでなかったら、リンクをたどって同じ (x, y, z)
の値を持つボクセルを見つけるまで繰り返します。

ボクセルハッシュの問題は、このリンクをたどる操作は、
そのままでは単純な線形探索になってしまうということです。
同じハッシュ値のテーブルに複数のボクセルが登録されてしまうと、
トラバーサル時のパフォーマンスが低下してしまいます。
これは、ボクセルの分割数が多くなると顕著に

とはいえ、たとえばリンクリストのボクセルを、
さらに kd-tree 化して空間データ構造
を作成するようにすれば、効率がよくなりそうです。  

ハッシュ値の計算

どのハッシュスキームにも言えることですが、
(x, y, z) のペア(今回は一様グリッドなので、x, y, z は整数でよい)から、
ハッシュ値を計算するハッシュ関数は、
なるべく計算が高速で、かつ生成されるハッシュ値にかたよりがないのが望ましいです。

私が知っているハッシュ関数は、以下の二つです。

Wyvill, Brian, 3D Grid Hashing Function,
Graphics Gems, p. 343-345, code: p. 733-734,
http://www.acm.org/pubs/tog/GraphicsGems/gems/Hash3D.c

M. Teschner, B. Heidelberger, M. Mueller, D. Pomeranets, M.
Gross,
Optimized Spatial Hashing for Collision Detection of
Deformable Objects
,”
Proc. Vision, Modeling, Visualization VMV’03,
Munich, Germany,
pp. 47-54, Nov. 19-21, 2003
http://cg.informatik.uni-freiburg.de/publications.htm

少しテストした感じでは、前者の分布は、後者にくらべるとあまりよくないようです。
また、後者は単純な線形合同法のような感じなので、非常に効率よくハッシュ関数を
計算できるという利点があります。

空間データ構造にハッシュを用いるのは、後者の論文のように、
軟体の衝突判定ではよく用いられているようです。

M. Teschner, S. Kimmerle, B. Heidelberger, G. Zachmann, L.
Raghupathi, A. Fuhrmann, M.-P. Cani, F. Faure, N.
Magnenat-Thalmann, W. Strasser, and P. Volino,
Collision Detection for Deformable
Objects
http://eg04.inrialpes.fr/Programme/STAR/STAR5.html
State of the Art Reports, EUROGRAPHICS 2004.

特に軟体のように毎フレーム物体の位置が変化する場合、
bsp のように物体を基準に空間を切るのではなく、
一様グリッドのように空間を基準に空間を切る(物体の形状に非依存)
ほうが効率がよいからのようです。

また、ハッシュをよく使う分野、たとえば通信や暗号などの研究分野にも、
もしかしたら計算が早く、グラフィックスにも使えるような
効率的なハッシュ関数があるのかもしれませんね。

オクツリーハッシュ

応用として、オクツリーをハッシュ化するオクツリーハッシュがあります。

Michael S. Warren and John K. Salmon,
“A parallel hashed oct-tree N-body
algorithm”,

In Supercomputing ’93, pages 12-21, Los Alamitos, IEEE
Comp. Soc, 1993
http://www.thesalmons.org/john/pubs/references.html

この論文では、ハッシュ値の計算には x, y, z のビットインターリーブを用いています。
子のオクツリーボクセルのトラバースを効率よく行うためのようです。
既存の CPU ではビットインターリーブを計算する専用の命令はないため、
以下のテーブルルックアップによる方法を使うのがベストでしょうか。

http://lucille.sourceforge.net/blog/archives/000365.html
http://lucille.sourceforge.net/blog/archives/000363.html

オクツリーハッシュも時間があれば試してみたいですね。

 

wavelet noise and modified noise 2


http://lucille.atso-net.jp/blog/archives/2005/08/wavelet_noise_a.html

からの続きで、ここで実際に Perlin noise と wavelet noise を FFT 変換してみようと思う。

FFT 変換は、octave や mathematica などの数値計算ソフトでサポートされているが、
ここでは自分で作成したノイズプログラムを FFT への入力としたいから、
C 言語のライブラリであるほうが都合がよい。

今回は FFT 変換を行うライブラリとして有名な fftw を利用した。
(ちなみに octave では内部で fftw が使われている)

http://www.fftw.org/

クロスプラットフォームであるのも魅力で、cygwin のパッケージなどにも含まれている。
最新版はバージョン 3 で、バージョン 2 からは API 周りなども違っているので注意が必要である。
ここではバージョン 3 を用いた。

フーリエ像の生成は、

フーリエ変換したい元の 2D 画像を作成 -> 離散フーリエ変換(FFT)

という手順で FFT を行うことで得られる。

もちろん入力を離散画像とすることにより、
ナイキスト限界である 2 ピクセル周期以下の高周波は正確に表現できなくなるので、
入力の関数は 2 ピクセル周期以上としなければならない。

つまり、画像の解像度を 256 としたばあい、

for (i = 0; i < 256; i++) {
        image[i] = noise(0.5 * (i + 0.5));
}

と、noise() への引数入力は 1/2 倍以下にする必要がある。
(ほんとはちょうど 1/2 がナイキスト限界になるので、
これよりすこし周期を大きくする必要があると思う)

もちろん、 Perlin ノイズは周期 1 でかつ整数点ではゼロの値となるので、
いずれにせよピクセルサイズに対して 1/2 倍にする必要がある。

Wavelet noise では、今回は、ちょうど 2D 画像と同じ大きさのノイズテーブルを作成している。
論文の中では、ノイズのテーブルを埋める初期の乱数列には、ガウスノイズ gaussianNoise() を
使うことが推奨されており、これは、[-1, 1] の正規分布に従う乱数列となっている。

この正規分布乱数列の生成には、たとえば論文で述べられているように、
Box-Muller 法を用いて求めることができる。

http://en.wikipedia.org/wiki/Box-Muller_transform

しかし、正規分布は無限のサポート [-∞, ∞] を持つため、生成される値は
必ずしも [-1, 1] の範囲にはない。そのため結果を [-1, 1] へ
クランプしなければならない。

これはどうしたらよいかと迷ったが、平均 0, 分散 1 の
正規分布 N(0, 1) の性質から、4 で割ることにした。
( N(0, 1) の場合、x 軸の値は [-4, 4] の間に 100% の確率で分布するため。
正確には +- 3.9 付近らしい)
本当はここらへんはもう少し考慮する必要があるかもしれない。
( [-1, 1] に分布するような分散 a を持つ N(0, a) にしなければ正しくないなど)

fftw によるフーリエ変換

FFTWを使ったサンプルスクリプト
http://kuma2.isl.titech.ac.jp/~tomoya/blog/archives/000008.html

を参考にし、fftwを利用して Perlin noise と wavelet noise のフーリエ変換を
行うコードを angelina に記述した。

http://lucille.atso-net.jp/svn/angelina/noise/

FFTW が出力するフーリエ像では象限の定義が異なるのか、
中心に低周波(DC component)がこない配置となっているようである。
そのため、中心に低周波が来るように、
フーリエ像を画像ファイルに出力する段階で象限の並べ替えを行っている。

以下が結果画像である。


2D Perlin ノイズの画像とそのフーリエ像


2D Wavelet ノイズの画像とそのフーリエ像

どうも Perlin ノイズと Wavelet ノイズでは基本となる周期が異なるようで、
フーリエ像も同スケールの周波数帯となってはいないが、
(wavelet noise のほうがなめらかっぽいのでフーリエ像における周波数帯も低周波部分に寄っている)
とりあえず論文のようなフーリエ像が得られていることがわかるのでよしとする。

wavelet noise のフーリエ像にはかすかに縦線と横線が見える(輝度スケールを大きくすると分かる)が、
これは画像生成時点か画像の境界に、非連続部分の周波数が現れているのかもしれない。

また、Perlin noise のフーリエ像にも、その形からして周期 1 以下(画面外側)の周波数がすこしカットされているので、

ノイズ画像を作るときに、もうすこししっかりとスケールを設定し、またノイズを数値積分計算して
ナイキスト周波数以下の周波数となるようにより正確なピクセル値を決定する必要がありそうだが、
いまのところこれ以上詳細に踏み込む必要もないのでここでやめにする。

次回は multi-band 化の実装と考察の予定である(Olano’s modified noise
も取り上げねばならない)

To be continued…

64bit 環境とコンパイラ, Intel or AMD?

64bit環境で試すコンパイラ –
アセンブルコードに見るその実力

最近は 64bit CPU も出てきましたし、デュアルコア CPU も出てきたので、
そろそろ lucille も 64bit + デュアルコアという環境で開発を行っていきたいと
思っています。

そんなわけで、上記のページはとても参考になります。
かいつまんでみると、

– 64bit 化で 32bit に比べて 10% 程度の向上(コンパイラの最適化有効で)
– デュアルコアは理論どおり 2 倍くらいのパフォーマンスアップ
– SSE2(とあるけど、実際には float x 4 の SSE 演算なのかな?) を使う以外は、
だいたい Intel よりも AMD のほうが性能がよい

AMD は、SSE 演算は遅いがそれ以外は Intel よりも演算性能がよい、
というのを聞いたことがありますが、だいたいその通りのようですね。

ILM などの多くのプロダクションでは、レンダーファームには AMD のプロセッサが
使われているようです。レンダリングの計算性能がよいからなのでしょうか。

私はいままで AMD の CPU を使ったことがないのですが、
デュアルコア化は AMD のほうが設計がよいらしいので、
これを機会に実際に AMD の CPU も使ってみて、いろいろと試してみたいですね。

MQOtoRIB

前回 MQOtoRIB はもう公開されていないようだ、、、と書きましたが、
作成していただいた作者である Soichi.H. 先生から
パワーアップされて再公開されているとのご連絡をいただきました。

http://homepage.mac.com/soichi_h/renderer/

の PolyConv です。

これらエクスポータを最大限活用していただくためにも、
lucille の RIB 読み込み部分やパーサなどもしっかり作りこんでいかなければなりません。
(いまのところパースエラー処理もまともにできないダメダメなので)

SBR 2005 観覧会

前回告知の通り、9/23 日に SBR 2005 観覧会を開かせていただきます。


http://lucille.atso-net.jp/sbr2005/index.php?SBR%202005%20%B4%D1%CD%F7%B2%F1

SIGGRAPH 2005 の論文、未来の SIGGRAPH 論文ネタにご興味のあるかたは、ぜひともご参加ください。

Security-Enhanced Versions of CRT Functions

Visual C++ 2005 Express(beta)  で lucille のコンパイルテストをしていると、

warning C4996: ‘fopen’ が古い形式として宣言されました。
c:\program files\microsoft visual studio 8\vc\include\stdio.h(235)
: ‘fopen’ の宣言を確認してください。

というような警告が出ていました。
ついに VC では標準 C ライブラリ関数すらサポートしなくなるのか?…
と思いながら華麗にスルーしていたのですが、そろそろいい加減原因を調べてみるかという気分となり、
よくよく C4996 のエラーコードで MSDN を検索すると、以下のようなお返事が。


http://whidbey.msdn.microsoft.com/library/default.asp?url=/library/en-us/dv_vccomp/html/926c7cc2-921d-43ed-ae75-634f560dd317.asp

つまりは、いくつかの関数については、セキュリティを高めた標準 C ライブラリ関数の
拡張版があるのでそっちを使ってくれ、という意味での警告ということでした。

セキュリティを高めた標準 C ライブラリ関数の代替関数のリストは、以下の MSDN にあります。

http://whidbey.msdn.microsoft.com/library/default.asp?url=/library/en-us/dv_vccrt/html/f87e5a01-4cb2-4379-9e8f-d4693828c55a.asp

たとえば、fopen() の場合ですと、fopen_s() を替わりに使うことが推奨されています。

http://whidbey.msdn.microsoft.com/library/default.asp?url=/library/en-us/dv_vccrt/html/c534857e-39ee-4a3f-bd26-dfe551ac96c3.asp

セキュリティを高めたバージョンの標準 C ライブラリ関数の拡張版は、VC 固有のものでしょうか?

というわけで、これからは、VC 2005 向けには、
ワーニングを抑制するために、コンパイラオプションで _CRT_SECURE_NO_DEPRECATE を定義するか、
以下のように ifdef で切るなりしてコードを書いていくのがよさそうです。

#if defined(_MSC_VER) && (_MSC_VER >= 1400)  /* VC 2005 */

#else /* others */

#endif

Spiral バケットレンダリング

lucille では、Z 字順、ヒルベルト曲線順でバケットをレンダリングするコードを書いてみたのですが、

http://lucille.sourceforge.net/blog/archives/000356.html (ヒルベルト曲線順)
http://lucille.sourceforge.net/blog/archives/000116.html (なぜバケットレンダリング順序が大切か?)

vray などの商用レンダラのように Spiral 順(真ん中から渦を巻くように外側に)
にレンダリングするようにしたいなーと思っていたら、

Production Rendering : Design and Implementation
http://www.amazon.co.jp/exec/obidos/ASIN/1852338210

http://bookweb.kinokuniya.co.jp/guest/cgi-bin/booksea.cgi?ISBN=1852338210

に、ビンゴなコードがあったので実装してみました。

http://lucille.atso-net.jp/svn/lucille/src/spiral.c

うーん、結構良いですね。実装も簡単ですし。

非正方形な画面サイズでも対応できるので、ヒルベルトや z 曲線順のように
(一時的に)バケットが画面からはみ出ることなく、
バケットを画像面上で連続してスイープさせることができます。

SBR 2005 観覧会のお知らせ

長らくお待たせしていましたが、
SBR 2005(SIGGRAPH 2005 論文実装レース) の観覧会の日程
が決まりましたので告知させていただきます。

期日: 9/23(金曜祝日)
場所: 東京都目黒区駒場あたりを予定
内容: SIGGRAPH 2005 論文解説および実装者による実装解説

今年(セカンド・ステージ)は、なんと SIGGRAPH 2005 論文採択者
自らの論文実装解説などがご聴講いただける予定です。

発表者の内容および観覧会の詳細は、SBR
2005
にて後日お知らせいたします。
(数日中に専用ページを作成します)

SIGGRAPH の論文や最新グラフィックス技術にご興味のあるかたは、
ぜひともご観覧にいらしてください。

とりいそぎ告知まで。