Syoyo Fujita's Blog

raytracing monte carlo

Month: September, 2009

Radeon HD 5870 has 544 DP GFlops

(via twitter)

http://pc.watch.impress.co.jp/img/pcw/docs/317/309/html/kaigai12.jpg.html

– double 精度で 544 Gflops(およそ  Core i7 の 10 倍)
– IEEE754-2008 互換

Spec を見る限りでは、「おぉ、ブラボー!」
GPGPU による、lucille のようなオフラインレイトレレンダラのアクセラレーションに本格的に使えそうです.

CPU での演算と比べて GPU で計算させても double 精度でも誤差が発生せず、
(double 精度でもきちんと IEEE754 互換であれば)
かつ CPU の 10 倍ほどの(理論ピークでの比較ですが)浮動小数点パフォーマンスがあるわけですから、
オフラインレンダラのような精度が求められつつも高速化したいというアプリに向いていると思います.
shared memory も GeForce より多いサイズを積んでいるようなので,
たとえば CUDA ベースのレイトレカーネルも Radeon HD 5870 で問題なく動きそう.

ただ、ATI はどちからというと既存グラフィックス性能重視、NV は汎用演算重視らしいので、
CUDA レイトレのように比較的分岐やランダムアクセスが多いカーネルは、
もしかしたら Radeon ではうまくパフォーマンスが出ない可能性もあります.

いずれにせよ、まずは OpenCL 対応を!… という状況ですね.
(待てないなら Brooks? でやる手もある… のかな)

World’s first SIMD ray-triangle intersection with ARM/NEON

I’ve wrote possibly world’s first SIMD ray-triangle intersection code with ARM/NEON instruction.

The code runs finely on iPod touch 3G(ARM Cortex) and iPhone 3GS.
The performance on iPod touch 3G is around 348 cycles per isect4().
# of assembly instruction of isect4() is around 140, thus the CPI = 348/140 = 2.48, which is not an efficient number for me(CPI should be near to 1.0).

Note that there’s still a room to optimize the code.
e.g. use fmad instruction instead of fadd/fmul combination.

Related articles

– SIMD ray-triangle intersection in Intel/AVX
http://lucille.atso-net.jp/blog/?p=649

– SIMD ray-triangle intersection in OpenCL
http://lucille.atso-net.jp/blog/?p=910

[Google groups] SIMD-cafe for SIMD lovers

なぜ SIMD が重要か

プロセッサのクロック数は頭打ち…
マルチコア化ではコアが増えるたびにコア間の通信がネックになる問題がある.
専用プロセッサを作るにも, 固定したアルゴリズム&量が出ないとビジネスにならないので,
やはり CPU 的にソフトなプロセッサが求められる.
一方で、エコ時代であるので省電力が求められる.
そんな HW の制約の中で、SW(アプリ)は富豪的プログラミングで HW の向上に頼るだけでは対応しきれなくなってきています.
特に認識処理、グラフィックス、データベース処理でパフォーマンスを出していくことが求められている.

このような、現在のプロセッサ(半導体)の技術的壁やプロセッサデザインの収益モデル(ビジネスモデル)を考えると、
特に SIMD 演算器とその活用が重要になってきます.

実際, x86/AVX(wider SIMD), iPhone 3GS(ARM/NEON), Atom/SSE が登場してきています.
アプリを SIMD 最適化する必要性は今後また重要になってくるでしょう.
(SIMD 最適化の最初のブームは、 MMX/SSE/SSE2 の出てきた 2000 年頃でしょうか)

しかし、どんなにコンパイラが賢くても、効率的な SIMD コードを自動で吐いてくれるのはまだまだ難しい.

そこで、プログラマがいかにアルゴリズムを SIMD フレンドリーにするか、
SIMD 演算器をいかにまく使ってプログラムを最適化していくかが、
今また大きく求められている、そんな気がしています.

また、ときどき SIMD 野郎と話すことがあるのですが、
結構みんな「ほかに SIMD のことを濃く話すことができるひとが居なくて…」
と言われることもあり, SIMD 野郎のための google groups を作ってみました.

SIMD cafe
http://groups.google.co.jp/group/simd-cafe

SIMD 野郎の参加をお待ちしております.

The effect of importance sampling

In this post, I’d like to let you know the power of importance sampling.

First of all, see the images rendered with importance sampling and without importance sampling when choosing random ray direction.

pt_is.png
(Path tracer with importance sampling)

pt_nois.png
(Path tracer without importance sampling)

Both use same the number of samples per pixel, thus the rendering time are identical for both images.
But you can see much more visual noise in the non importance sampled image.

For the image with importance sample, I’ve drawn the samples according to the following probability function,


(i.e. generate samples proportional to cosine attenuation)

to generate random ray direction.
This importance sampling technique is widely used in many G.I. renderer.

For non importance sampled image, I’ve used following naive probability function


(i.e. uniformly sample point on the hemisphere)

to generate random ray direction, then multiply cosine attenuation to the contribution.

Why so different?

Let me briefly explain why so the result image is different with and without importance sampling.

Importance sampling for the light integral.

Mathematically, Things what I’ve did in the above importance sampled image is reduced to computing the following integral.


(w_i is sampled with the PDF function p(w) = cos(theta) sin(theta) / PI)

Imagine if L(w) was always 1(e.g. uniformly lit sky), then the result of the integral become

(1 + 1 + 1 + … + 1) * (PI/N) = PI,

for any random number and for any N(# of MC samples).

The light integral with ordinal sampling.

On the other hand, The case of non-importance sampled image is formulated as follows.


(w_i is sampled with the PDF function p(w) = sin(theta) / (2 PI) )

Again imagine if L(w) always returned 1, then the integral will give

(cos(u0) + cos(u1) + cos(u2) + … + cos(uN)) * (2 PI / N) =
(0.34 + 0.56 + 0.90 + … + 0.12) * (2 PI / N) ~= PI

which is nearly equal to PI when N goes large, but not exactly goes to PI as in the importance sampled case.
There’s variance(error) in this case.

Conclusion

As a conclusion, the reason why 2 images are so different is due to the variance of cos(theta).
The variance of cos(theta) turned into more visual noise in non-importance sampled image.

Caustic Graphics ray traced global illumination demonstration

(via Twitter)

Caustic Graphics ray traced global illumination demonstration
http://www.pcper.com/comments.php?nid=7801

おー、Caustic でパストレっぽいこともできるのね.
ただ、パフォーマンスはちょっと微妙な感じ.
2M ~ 10M rays/sec くらいという感じかな.
ただ、実際の製品のパフォーマンスの 1/10~1/20 と言われている FPGA で動かしてこのレベルということであれば、まあ妥当な数字かもしれません.

しかし相変わらずいつもビデオでデモしているときの場所が大学の研究室っぽいのはなにか狙っているのだろうか.
アジェンダが、ホワイトボードに手書きだし 🙂

でもこれで $10M くらいレイズしているんだよなー、うらやましいぜ!

Review of Cloudy: From the perspective of a renderer writer.

cloudy_with_a_chance_of_meatballs.jpg

I’ve watch the Cloudy(Cloudy with a chance of meatballs) and I’d like to give visual impression from the perspective of a render writer.

Overall, Cloudy was beautifully rendererd compared to previous motion picture Monster-house(Both used same renderer).

There’s no visually recognizable Monte-carlo noises in Cloudy
(noise in global illumination, noise in motion blur, etc) even seen though on the IMAX quality(4k?) screen.
It’s perfect. I’ve convinced raytracing-based renderer is now practical!

Motion blur

Motion blur would be the one of the biggest challenge in Cloudy, because rendering motion blur effect in raytracing-based render requires a lot of time to render(sometimes 10x than still images).

I cannot see any visual Monte-carlo noise(which is one of the biggest artifact in raytracing-based renderer).
I think they did well on rendering beautiful motion-blur effect.

G.I./Raytracing

There’s no G.I. noise in Cloudy, awesome!
But I could see perfect reflection and refraction effect(glass material), which sometimes gives distinct appearance compared to characters and buildings(toon-like appearance).

And… I have to say there wasn’t no caustics effect onto the table from the Cognac glass in the restaurant scene.
(I know its OK, because the table are self-emitted).

Hairs

Rendering hairs are also one of the most challenging factor in Cloudy.
Hairs was rendered beautifully.
I also couldn’t see any visual artifacts for the hairs(monte-carlo G.I. noise, motion blur noise).

Clothes

I could see a kind of BRDF-based cloth apparance, but overall cloth seems are rendererd with simple appearance model.

Effects(Particles and smokes)

As far as I know, effects are rendered SPIW’s in-house render.
Effects give TOTALLY distinct appearance from characters and environments: it has very realistic appearance.
I think, personally, it’d be better if effects have somewhat toon like appearance like characters.

[Ja]

Cloudy を観てきました.

ストーリーは、私はアメリカナイズされたアニメものは Pixar 作品ですら好まないので「へー」という感じでしたが、
子供が観たらいろいろ楽しいのかも.
まさかあそこまでのフードディザスタームービー(食べ物災害映画)になるとは!

さて、期待のレンダリングの出来ですが、、、
モーションブラーも G.I. もヘアーもよくできています。
IMAX 3D で観てきましたが、あれだけの解像度でも視覚的なモンテカルロノイズはまったくありませんでした.
(逆に奇麗すぎてもうちょっと映画特有なフィルムノイズ的なものとかあったほうがいいのでは、と思うくらい)

「あー、やっぱりもうレイトレレンダラでいいんじゃね?ってかこれからはレイトレレンダラっしょ」

そんな感じですね. もうここまでレイトレベースでクオリティ出るなら、 REYES の時代は終わりに近いんじゃないかなー.
(レンダリング時間が解決されれば)

ただ、いかにもレイトレなエフェクト(窓やグラスの反射や屈折)はいかにもレイトレ!って感じで、
ちょっとキャラクターや背景と比べて浮いている感は否めませんでした.
glossy 反射やフロストグラス的な表現があるとよかったかも.

あとはパーティクルと煙も、ちょっと浮いている感じでした.
パーティクル系は別途のインハウスレンダラを使っているとのことなので、そのせいなのかもしれませんが、
もうちょっと各エフェクトとキャラクターの質感との間で「なじむぞッ!!」感を出すとよかったのでは。

あとはレストランシーンで、ワイングラスからのコースティクスが無かったのが気になりますねー.
テーブルが発光するようなもののように見えたので(テーブルにライトが埋め込まれている)、
だからコースティクスが無いだけだったのかもしれませんが.

エンドクレジットにて、レンダラ野郎を確認しました。
顔ぶれはほとんど変わっていませんでしたね。

Watching Cloudy gathering 9/19(Sat)

せっかくなんで、某スタジオの某レンダラの出来を確認するのを裏の目的として,
Cloudy(くもりときどきミートボ-ル) 鑑賞オフをしたいと思います.

N 年の歳月を費やした? モーションブラーの出来はどうか!ヘアーの出来はどうか!
ぜひともご確認ください.
観賞後、ご意見あれば直接某レンダラ開発チームにレポートしますんで 😉

時間と場所ですが、
9/19(土) 109シネマズ川崎 IMAX シアターになります.

残念ながらまだスケジュールは決まっていないようです。
また、インターネットからの予約も3日前からしかでません。

上映時間は夕方あたりを選択して、観賞後は(勝手に)某レンダラの出来を語る懇親会にしたいと思います.

参加希望の方は以下の URL にて、参加表明をお願いします.
http://atnd.org/events/1542

おまけ

109 シネマズ川崎で、ちょうど今、ダークナイトも IMAX 上映してます.

OpenCL ray-triangle intersection kernel.

openclraycastsc.png

I’ve wrote possibly World’s first ray-triangle intersection kernel in OpenCL.
I am now working with porting full raytracing kernel(traversal, shader, isector, etc) into OpenCL.
This article is just exposing the small part of my professional(?) OpenCL work to you.

Here’s the sample code of ray-triangle intersection kernel in OpenCL.


//
// Copyright 2009 Syoyo Fujita.
//
typedef struct _ray4_t
{
    float4 rox, roy, roz;
    float4 rdx, rdy, rdz;
    float4 t;
} ray4_t;

typedef struct _triangle4_t
{
    float4 v0x, v0y, v0z;
    float4 e1x, e1y, e1z;
    float4 e2x, e2y, e2z;
} triangle4_t;

typedef struct _triangle_t
{
    float v[3][4];  // (x,y,z,w) * 3
} triangle_t;

static inline float4
mycross4(float4 a, float4 b, float4 c, float4 d)
{
    return ((a * c) - (b * d));
}

static inline float4
mydot4(float4 ax, float4 ay, float4 az, float4 bx, float4 by, float4 bz)
{
    return (ax * bx + ay * by + az * bz);
}

static int4
isect4(
    float4 *t_out,
    float4 *u_out,
    float4 *v_out,
    ray4_t ray,
    triangle4_t tri)
{
    const float4 px = mycross4(tri.e2z, tri.e2y, ray.rdy, ray.rdz);
    const float4 py = mycross4(tri.e2x, tri.e2z, ray.rdz, ray.rdx);
    const float4 pz = mycross4(tri.e2y, tri.e2x, ray.rdx, ray.rdy);

    const float4 sx = ray.rox - tri.v0x;
    const float4 sy = ray.roy - tri.v0y;
    const float4 sz = ray.roz - tri.v0z;

    const float4 vone  = (float4)(1.0f);
    const float4 vzero = (float4)(0.0f);
    const float4 veps  = (float4)(1.0e-6f);

    const float4 det = mydot4(px, py, pz, tri.e1x, tri.e1y, tri.e1z);
    const float4 invdet = (float4)(1.0f) / det;


    const float4 qx = mycross4(tri.e1z, tri.e1y, sy, sz);
    const float4 qy = mycross4(tri.e1x, tri.e1z, sz, sx);
    const float4 qz = mycross4(tri.e1y, tri.e1y, sx, sy);

    const float4 u = mydot4(sx, sy, sz, px, py, pz) * invdet;
    const float4 v = mydot4(ray.rdx, ray.rdy, ray.rdz, qx, qy, qz) * invdet;
    const float4 t = mydot4(tri.e2x, tri.e2y, tri.e2z, qx, qy, qz) * invdet;

    int4 mask = (fabs(det) > veps) & ((u+v)  vzero) & (v > vzero)
              & (t > vzero) & (t < ray.t);


    (*t_out) = t;
    (*u_out) = u;
    (*v_out) = v;

    return mask;

}

__kernel
void
main(
    __global uint *out)
{
    int x = get_global_id(0);
    int y = get_global_id(1);
    int w = get_global_size(0);
    int h = get_global_size(1);

    ray4_t ray4;

    ray4.rox = 2.0f * ((x - w) / (float)w) + 1.0f;
    ray4.roy = 2.0f * ((y - h) / (float)h) + 1.0f;
    ray4.roz = (float4)(0.0f);

    ray4.rdx = (float4)(0.0f);
    ray4.rdy = (float4)(0.0f);
    ray4.rdz = (float4)(-1.0f);

    ray4.t = (float4)(1.0e+30f);



    triangle_t tri;
    triangle4_t tri4;

    tri.v[0][0] = -0.5f;
    tri.v[0][1] = -0.5f;
    tri.v[0][2] = -1.0f;
    tri.v[1][0] =  0.1f;
    tri.v[1][1] =  0.75f;
    tri.v[1][2] = -1.0f;
    tri.v[2][0] =  0.5f;
    tri.v[2][1] = -0.75f;
    tri.v[2][2] = -1.0f;

    tri4.v0x = (float4)(tri.v[0][0]);
    tri4.v0y = (float4)(tri.v[0][1]);
    tri4.v0z = (float4)(tri.v[0][2]);
    tri4.e1x = (float4)(tri.v[1][0] - tri.v[0][0]);
    tri4.e1y = (float4)(tri.v[1][1] - tri.v[0][1]);
    tri4.e1z = (float4)(tri.v[1][2] - tri.v[0][2]);
    tri4.e2x = (float4)(tri.v[2][0] - tri.v[0][0]);
    tri4.e2y = (float4)(tri.v[2][1] - tri.v[0][1]);
    tri4.e2z = (float4)(tri.v[2][2] - tri.v[0][2]);

    float4 t, u, v;
    int4 mask = isect4(&t, &u, &v, ray4, tri4);

    //int r = mask.x & 0xff;
    //int g = mask.y & 0xff;

    int r = (int)(t.x * 255);
    int g = (int)(u.x * 255);
    int b = (int)(v.x * 255);

    r = mask.x & r;
    g = mask.x & g;
    b = mask.x & b;

    out[x+y*get_global_size(0)] = (uint)(r | (g<< 8) | (b << 16) | (255<<24));
}

Similar articles

World’s first SIMD ray-triangle intersection in AVX?
http://lucille.atso-net.jp/blog/?p=649

AMD OpenCL on windows too slow…?

cltestsc.JPG

I’ve wrote a simple OpenCL test program which simply fills pixels and display it through OpenGL.
When I run it on ATI Stream SDK 2.0 beta on Windows,
I got terribly slow performance: around 8 secs per frame!
(On SnowLeopard OpenCL(CPU version) it runs around 0.0001 secs per frame)

What’s wrong? Because AMD version is a beta version?
I’d like to discuss this problem but I don’t know where is the good place to ask…

FYI, Here’s OpenCL kernel used in the test program.


__kernel
void
main(
   __global uint *out,
   uint col) // not used.
{
   int x = get_global_id(0);
   int y = get_global_id(1);

   out[x+y*get_global_size(0)] = (uint)(x | (y << 8) | (255 << 16) |
(255<<24));
}

[Ja]

単純に 256×256 画像のピクセルをフィルするだけの OpenCL プログラムを書いてみたのですが、
ATI StreamSDK 2.0 beta(Windows, CPU 利用)の OpenCL 実装だと、8 秒近くもかかってとてつもなく遅い!
(SnowLeopard だと 0.0001 秒なのに…)

追記

Visual Studio 2008 を使っているのですが、デバッグ -> デバッグ開始だととても遅くなる模様(Release ビルドでも).
デバッグ -> デバッグなしで実行だとそれなりの速度でうごきました.
んー、なにか変な実行時デバッグチェックが、OpenCL ドライバかどうか分かりませんが、入っているのですかね.
ちなみに「デバッグなしで実行」でも、0.05 secs くらいかかりました.
(SnowLeopard の 0.0001 秒に比べるとまだ遅い)

Kai-fu Lee quits Google and launch his venture fund.

http://www.reuters.com/article/technologyNews/idUSTRE5850BQ20090906

http://www.bloomberg.co.jp/apps/news?pid=90003011&sid=aiAVEbSlU7ek&refer=jp_asia

お、李開復 Google 辞めてファンド作るのか。
んー、まさに私がやりたいことだなー.
私の場合: 稼いだマネーでレンダラファンド+レンダラ財団立ち上げ.

李開復の100 億円のファンドは、良い目標になるけど、
とりあえずは自分のほうは 5 億でファンド or 財団を立ち上げられないかなーと考えています.
…まあまずは 5 億稼ぐネタを探すところからですかね 😉

参考資料

ビル・ゲイツ、北京に立つ
http://www.amazon.co.jp/dp/4532313031
MSRA とか李開復とかあたりがいくらか書かれています.
SIGGRAPH ネタもほんの少しだけどあるよ.