Syoyo Fujita's Blog

raytracing monte carlo

Month: January, 2005

休載

ド.. ドドド..  ドドド…. ドドドドド….
  ドド… ドドドド… ド…    ドド   ドド
    ド… ドド… ド…  ド…  
 ド… ド… ド…    ド…   ドド…

2/16 あたりまで、休載します。

  ドドド… ドド… ド…    ド…

ドド…     ドド….

  ドド… ドド.. ド…    ド…

Green’s intersection method

Green の交差判定法というのは以下のようです。

Simple, Fast Triangle Intersection, by Chris Green http://www.acm.org/tog/resources/RTNews/html/rtnv6n1.html#art6

 

ちなみに Green の交差判定法を含め、各種交差判定法の reference が載っている論文が

J. Amanatides and K. Choi, “Ray Tracing
Triangular Meshes”,
Proceedings of the Eighth Western
Computer Graphics Symposium,
April 1997, pp 43-52.

にあります。本論文では、プリュッカー座標系(Pluecker coordinate)による交差判定が各種交差判定法より 25
% 程度高速であると述べています。プリュッカー座標系は、レイが三角形の内部にあるかのテストが内積のみでできるので、
非常に数値計算がロバストであるという利点があるものの、プリュッカー座標では 6×3 = 18 floats
のストレージコストがかかるのと、t までの計算は効率的ですが、u, v まで求めるとなるとより extra
な計算量が必要になりそうなのが欠点でしょうか。

to be continued…

 

SIMD ray intersection, SIMD BSP traversal

リアルタイムレイトレの鬼と呼ばれる Ingo Wald 博士の博士論文が公開されている。

Ingo Wald,
Realtime Ray Tracing and Interactive Global Illumination
Ph.D. Thesis.
http://www.mpi-sb.mpg.de/~wald/PhD/index.html

氏がかかわってきたリアルタイムレイトレプロジェクト RTRT(Real Time Ray Tracing) や
SaarCOR の集大成でもありますが、とくに第七章では、それらのコアとなる SIMD ray/triangle 交差判定や
SIMD BSP トラバーサルのコードつきでの詳細な解説があります。

SIMD ray intersection

レイと三角形の交差判定には、重心座標をベースとした射影法(Projection Method)を採用しています。射影法は、
三角形の法線ベクトルの xyz 要素のうち、もっとも大きい平面へと三角形を射影することで、問題を 2
次元に落とし込んで処理を効率的に行うというものです。また、射影する平面は法線ベクトルの要素がもっとも大きい方向となるので、
determinant の計算でゼロ割が生じることもないため、より計算はロバストになります。

射影法はまた、各種係数を前計算しておくことができるため、交差判定計算も単純になります。博士論文では、
各三角形毎に交差判定に必要な前計算の係数データは 48 byte(キャッシュラインにのせるためにパディングを含む)となっています。

struct TriAccel
{
// first 16 byte half cache line
// plane:
float n_u; // == normal.u / normal.k
float n_v; // == normal.v / normal.k
float n_d; // constant of plane equation
int k; // projection dimension

// second 16 byte half cache line
// line equation for line ac
float b_nu;
float b_nv;
float b_d;
int pad; // pad to next cache line

// third 16 byte half cache line
// line equation for line ab
float c_nu;
float c_nv;
float c_d;
int pad; // pad to 48 bytes for cache alignment purposes
};

これはもちろん三角形の頂点データとは別に必要となるので、非常にポリゴン数が多いシーンでは、逆にメモリネックとなり得ます。
博士論文によると、ポリゴン数が多いシーンでは、通常よく用いられている extra なデータを必要としない Mueller
式に切り替えているとのこと。

射影法は Graphics Programming Methods の mlt や kd-tree
の記事のコードでも用いられており、そこでは Green’s intersection method となっていました。Wald
の博士論文ではネタ元は書かれていませんでした。うーん、この射影法のネタ元は Green って人の発案なのかな?

to be continued…

Solid Harmonics

さて、今となっては球面調和関数を調べることは
ほとんどなくなってしまったのだが、有益な変換式をいまさらながら
見つけました。

通常球面調和関数は球面座標で与えられますが、
デカルト座標での (x, y, z) による解析的な式で表現することができる
ことはよく知られています。


http://www.cs.princeton.edu/~ah/publications/multipole/html/node30.html

リアルタイムグラフィックスでの利用を考えると、
評価が容易なデカルト座標での解析的な式が好まれるでしょう。

んで、実は 5 次以上とかになるとあまりネットや文献ではこの
デカルト座標での解析式は載っていないので、
(高次になると数値誤差も問題になりますしね)
当時はもんもんとしながらどこかに無いものかと探し回っていたのですが、
なんと Mathematica でお手軽に計算する方法がッ!!


http://mathforum.org/epigone/comp.soft-sys.math.mathematica/kigruzhon

虚数部は単に捨てて実数部だけ使えばいいはず。
ヤベー、マジすげーよ。Mathematica。ひれ伏します。

まあ、今となっては生真面目に p+q+s=l and p-q=m となるような
(p, q, s) の組を総当りで探すプログラムを書くのでもよかったんだなぁと
思いを馳せる日々であった。

SIMD ray/box intesection 2

前回 SIMD でレイとボックスの交差判定を行うコードを紹介しましたが、

http://lucille.sourceforge.net/blog/archives/000372.html

flipcode のほうで、 Inf
となるコーナーケースの問題にも対応したよりロバストなバージョンが公開されていました。

http://www.flipcode.com/cgi-bin/fcarticles.cgi?show=65014

しかし、このコードは SSE では動くが( min(NaN, Inf) = Inf ),
altivec に移植しようとなると altivec は java モードであるため、
min(NaN, Inf) = NaN となってしまうので注意が必要である。  


http://www.simdtech.org/apps/group_public/email/altivec/200412/msg00006.html

というわけで、AltiVec 版に書き直してみた。

#include <stdio.h>
#include <math.h>

// AtiVec version of ray/box intersection from flipcode
static vector float vplus_inf;
static vector float vminus_inf;
static vector float vzero;
static vector float vnan;

#define elem(v, i) *(((float *)&(v)) + (i))

typedef struct _aabb_t {
  vector float bmin;
  vector float bmax;
} aabb_t;

typedef struct _ray_t {
  vector float pos;
  vector float inv_dir;
} ray_t;

typedef struct _ray_segment_t {
  float t_near, t_far;
} ray_segment_t;

int
ray_box_intersect(const aabb_t *box, const ray_t *ray, ray_segment_t *rs)
{
  // you may already have those values hanging around somewhere
  const vector float plus_inf = vplus_inf;
  const vector float minus_inf = vminus_inf;

  // use whatever's apropritate to load.
  const vector float box_min = box->bmin;
  const vector float box_max = box->bmax;

  const vector float pos = ray->pos;
  const vector float inv_dir = ray->inv_dir;

  // use a div if inverted directions aren't available.
  vector float l1 = vec_madd(vec_sub(box_min, pos), inv_dir, vzero);
  vector float l2 = vec_madd(vec_sub(box_max, pos), inv_dir, vzero);

  // the otder we use for those min/max is vital to filter out
  // NaNs that happens when an inv_dir is +/- inf and
  // (box_min - pos) is 0. inf * 0 = NaN.

  /* In AltiVec min(x, NaN) = NaN,
   * so we modify it so that being min(x, NaN) = x.
   */
  vector bool int not_nan1 = vec_cmpeq(l1, l1);
  vector float y = vec_sel(plus_inf, l1, not_nan1);
  const vector float filtererd_l1a = vec_min(plus_inf, y);

  vector bool int not_nan2 = vec_cmpeq(l2, l2);
  y = vec_sel(plus_inf, l2, not_nan2);
  const vector float filtererd_l2a = vec_min(plus_inf, y);

  /* In AltiVec max(x, NaN) = NaN,
   * so we modify it so that being max(x, NaN) = x.
   */
  y = vec_sel(minus_inf, l1, not_nan1);
  const vector float filtererd_l1b = vec_max(minus_inf, y);

  y = vec_sel(minus_inf, l2, not_nan2);
  const vector float filtererd_l2b = vec_max(minus_inf, y);

  // now that we'ew back on our feet, test those slabs.
  vector float lmax = vec_max(filtererd_l1a, filtererd_l2a);
  vector float lmin = vec_min(filtererd_l1b, filtererd_l2b);

  // unfold back. try to hide the latency of the shufps & co.
  // (a, b, c, d) -> (b, c, d, a)
  const vector float lmax0 = vec_sld(lmax, lmax, 4);
  const vector float lmin0 = vec_sld(lmin, lmin, 4);

  // lmax(a, b) = (min(a0, b0), a0, a1, a2)
  const vector unsigned char permtable0 = (vector unsigned char)(
          0x00, 0x01, 0x02, 0x03, 0x14, 0x15, 0x16, 0x17,
          0x18, 0x19, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f
  );
  const vector float s1 = vec_min(lmax, lmax0);
  lmax = vec_perm(s1, lmax, permtable0);

  // lmin(a, b) = (max(a0, b0), a0, a1, a2)
  const vector float s2 = vec_max(lmin, lmin0);
  lmin = vec_perm(s2, lmin, permtable0);

  // low(a, b, c, d) | high(e, f, g, h) = (c, d, g, h)
  const vector unsigned char permtable1 = (vector unsigned char)(
          0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f,
          0x18, 0x19, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f
  );

  const vector float lmax1 = vec_perm(lmax, lmax, permtable1);
  const vector float lmin1 = vec_perm(lmin, lmin, permtable1);

  // lmax(a, b) = (min(a0, b0), a0, a1, a2)
  const vector float t1 = vec_min(lmax, lmax1);
  lmax = vec_perm(t1, lmax, permtable0);

  // lmin(a, b) = (max(a0, b0), a0, a1, a2)
  const vector float t2 = vec_max(lmin, lmin1);
  lmin = vec_perm(t2, lmin, permtable0);

  // store lmin[0] to rs.t_near
  // store lmax[0] to rs.t_far

  rs->t_near = elem(lmin, 0);
  rs->t_far  = elem(lmax, 0);

  // bool ret = mm_comige_ss(lmax, mm_setzero_ps()) &&
  // mm_comige_ss(lmax, lmin);

  bool ret = (rs->t_far >= 0.0) && (rs->t_far >= rs->t_near);


  return ret;

}

int
main(int argc, char **argv)
{
  aabb_t box;
  ray_t  ray;
  ray_segment_t rs;

  const float flt_plus_inf = -logf(0);
  elem(vplus_inf, 0) = flt_plus_inf;
  elem(vplus_inf, 1) = flt_plus_inf;
  elem(vplus_inf, 2) = flt_plus_inf;
  elem(vplus_inf, 3) = flt_plus_inf;

  elem(vminus_inf, 0) = -flt_plus_inf;
  elem(vminus_inf, 1) = -flt_plus_inf;
  elem(vminus_inf, 2) = -flt_plus_inf;
  elem(vminus_inf, 3) = -flt_plus_inf;

  box.bmin = (vector float)(-1, -1, -1, 0);
  box.bmax = (vector float)(1, 1, 1, 0);

  ray.pos = (vector float)(-1, -1, -1, 0);
  elem(ray.inv_dir, 0) = flt_plus_inf;
  elem(ray.inv_dir, 1) = flt_plus_inf;
  elem(ray.inv_dir, 2) = -1.0;
  elem(ray.inv_dir, 3) = 0.0;

  vzero = (vector float)(0.0, 0.0, 0.0, 0.0);

  vnan = vec_madd(vplus_inf, vzero, vzero);

  ray_box_intersect(&box, &ray, &rs);

  printf("t near = %f, t far = %f\n", rs.t_near, rs.t_far);

}

コンパイルは

$ gcc -faltivec raybox_altivec.c -lmx

たぶん SSE のときと同様に動くはず。

入力としてレイの方向はあらかじめ逆数であたえるため、コードのように逆数がゼロ割により Inf になる可能性がある( (0,
0, -1) という単純なレイの方向ベクトルの逆数を計算するだけで容易に Inf が現れる) 。

とはいえ、altivec で逆数(除算)を計算する場合、vec_re() + ニュートンラフソン法を用いるのが通例である。
vec_re(0.0) の結果は Inf であるので問題ないが、結局ニュートンラフソンの部分で 0 * Inf となり NaN
となってしまい、交差判定部分での NaN の処理というのは無駄になってしまう。

なのであらかじめ交差判定ルーチンに渡す前にレイ方向の逆数が NaN であれば Inf
にするなど対処を施しておく必要がありそうです。

ところで、 -logf(0) で +Inf を得るというのはよくつかうテクニックなのかな?

bit interleaving 3

以前、Hacker’s Delight(ハッカーのたのしみ)
において、
ビットインターリーブ(bit interleave)がグラフィックスに応用があるがその元ネタはどこ?
と述べましたが、

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

実は Graphics
Gems I
でちゃんと取り上げられていました。

Hacker’s Delight において、ヒルベルト曲線のグラフィックスへの応用
の部分の解説部分で引用されていたのが
Graphics Gems II であったので、
その時点でもしや? と気づくべきでした。

Graphics Gems I でのビットインターリーブは、
クアッドツリー(四分木)とオクツリー(八分木)のアドレッシングへの応用であり、
私の Z 曲線順のなんかよりもはるかに詳細に解説されています。

それにしても、やはり Graphics Gems はすごいなぁ、なんか大体思いつくと
だいたいはすでにそこに載っているようなものだったりする。
今でも内容は色あせていないし。

Graphics Gems は「こりゃスゲーよ。Graphics Gems にはッ!!、
全巻買わなければ!! と思わせる『スゴ味』があるッ!」
という感じでとりあえず全巻所有しているのだが、
その正統なる後継である? Game Programming Gems や GPU Geoms
には日本語版があるにもかかわらず、
まったく持って食指が動かないのはなぜだろうか?

さらにその Game Programming Gems に続く Graphics Programming Methods

『スゴ味』を感じて大変よいのだが、続編がでてくれるのかどうかどうか怪しくてちょっと残念。

valgrind, x86 linux プログラムデバッガ & プロファイラ

==1278== Use of uninitialised value of size 4
==1278==    at 0x806D014: ri_vector_mul (in /home/syoyo/work/lucille/src/lsh/lsh
)
==1278==    by 0x4242C99E: ???

Matt Pharr 大先生も御用達? の valgrind というツールがあります。

http://valgrind.kde.org/

このツールは、メモリリークの検出、スレッドエラーの検出、
キャッシュのプロファイル、ヒープのプロファイルなどと、
非常に高性能かつ多機能なメモリ関連ツールです。

本ツールの大きな特徴としては、
たとえば dmalloc のように特別になんらかのライブラリとリンクしたり、
ソースコードを書き換える必要はなく、すでにコンパイルされたプログラムを
そのまま処理できる点でしょう。

もちろん、デバッグ情報付きでコンパイルされていれば、
たとえばメモリエラーが起きたときには、
そのエラーの起きた関数名と行番号が出力されます。

x86 linux 限定であるという制約はあるものの、
Mac OS X の malloc_history, leaks と並び非常に使えるツールです。
(PowerPC/linux への移植版も行われているようです)

Windows でも、colinux で動かせば問題なく使えます。

ためしに lucille に valgrind をかけて見ました。

...
==1278== More than 30000 total errors detected.  I'm not reporting any more.
==1278== Final error counts will be inaccurate.  Go fix your program!
==1278== Rerun with --error-limit=no to disable this cutoff.  Note
==1278== that errors may occur in your program without prior warning from
==1278== Valgrind, because errors are no longer being displayed.
...
==1278== ERROR SUMMARY: 30000 errors from 6 contexts (suppressed: 2 from 1)
==1278== malloc/free: in use at exit: 6136545 bytes in 37165 blocks.
==1278== malloc/free: 98898 allocs, 61733 frees, 13423847 bytes allocated.
==1278== For a detailed leak analysis,  rerun with: --leak-check=yes
==1278== For counts of detected errors, rerun with: -v

うーん、メモリリークや未使用変数の利用がいっぱい…

 

valgrind の GUI フロントエンドとして、
キャッシュプロファイルやコールグラフの可視化を行う
KCachegrind というのもあります。


http://kcachegrind.sourceforge.net/cgi-bin/show.cgi/KcacheGrindIndex

Thai Statue

Thai Statue も lucille でレンダリングしてみました。512×512 ピクセル、2×2 アンチエイリアシング、Structured Importance Sampling を用いてシェーディング点あたり 300 サンプル。メッシュデータは簡略化せずに元の 1000 万三角形です。

うーん、やはり 512×512 じゃ 1000 万三角形の細かさは見えませんね。やはり最小でも 1024×1024 でレンダリングしないといけませんね。

Jim Blinn’s Corner 日本語版

jim1.jpgjim2.jpg

 

Jim Blinn のコラム本が日本語となって出版されている。

Jim Blinn’s Corner
日本語版 (1) A Trip Down the Graphics Pipeline

http://www.amazon.co.jp/exec/obidos/ASIN/4274065367

Jim Blinn’s Corner日本語版 (2) Durty Pixels
http://www.amazon.co.jp/exec/obidos/ASIN/4274065650

Jim Blinn’ Corner といえば、今回日本語版を読むまでは、
SIGGRAPH 2002 のときに無償で配られていた
IEEE Computer Graphics and Applications 誌でちらりと見たきりである。
(その号での記事は “Visualize Whired 2×2 Matrix”)

今読み返して見ると、この号には、 Andrew Glassner’s Note
(これは Jim Blinn’s Corner と同じく、Andrew Glassner により
長らく連載されているコラムである)
 と Geometric Algebra: A
Computational Framework For Geomtrical
Applications part2
という興味深い記事も掲載されていた。

当時、Jim Blinn ってブリンシェーディングのひと?
程度の知識しかなかった私としては、2×2 行列の可視化も GA も
なんじゃらほいという内容でまったく理解できなかったが。。。

Jim Blinn’s Corner の内容は、サブピクセルサンプリング、
ポリゴンクリッピング、同次座標による透視変換など、
普段 OpenGL の API をたたくようなレベルでは気にしないが、
実際にオフラインレンダラやソフトウェアリアルタイムレンダリングエンジン、
GPU 設計のときに皆が詰まるような問題を取り上げている。

CG の奥深くを知るにあたっては非常に有益かと思います。

書籍の紹介によれば、Jim Blinn 独特の言い回しで解説とのこと
であるが、やはり日本語からはそのうまみ?をあまり味わうことが
できない気がします。やはり原書も買うべきでしょうか。

Xorshift RNGs

G. Marsaglia. Xorshift RNGs. Journal of Statistical Software,
8(14) :1 6, 2003
http://www.jstatsoft.org/v08/i14/xorshift.pdf

George Marsaglia 氏により 2003 年に考案された、xor
とシフトを使うだけの超高速な擬似乱数生成器(Random Number Generator, RNG)です。周期は 2^k-1(k
= 32, 64, 96, 128, 160, 192)。ランダム性のテストにも十分合格するとのこと。たとえば、周期が
2^128-1 の場合のルーチンは以下のようになり、乱数値の計算部分はわずか 1 行である。

unsigned long xor128(){
  static unsigned long x=123456789,y=362436069,z=521288629,w=88675123;
  unsigned long t;
  t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) );
}

Xorshift 法については、擬似乱数研究で有名な Pierre L’Ecuyer 氏によりその理論の解析が行われている。

F. Panneton and P. L’Ecuyer, On the Xorshift Random
Number Generators
, 2004
http://www.iro.umontreal.ca/~lecuyer/papers.html

また、Richard Brent 氏による Xorshift 法による擬似乱数ルーチンのノートと、少なくとも
(2^4096-1)までの周期まで拡張が施された C 言語の実装 xorgen を以下のページから手に入れることができる。

Note on Marsaglia’s xorshift RNGs

http://web.comlab.ox.ac.uk/oucl/work/richard.brent/pub/pub218.html

Some Uniform and Normal Random Number Generators
http://web.comlab.ox.ac.uk/oucl/work/richard.brent/random.html

暗号処理に xorshift 法を使う、 xse というのもあります。

http://sandbox.emboss.co.nz/algorithm.php?CipherID=5

ただし、ニュースグループでのやりとりを見るとあまり評判はよくないみたいです。


http://www.derkeiler.com/Newsgroups/sci.crypt/2003-11/1355.html