今までの画像は、ライティングがいかにもうそ臭いものでした。
この点も修正して、実写らしい見た目になるようにがんばって見ましょう。
今回のプログラムは、次のものです。
まぁ、いつものように適当にファイルが入っています。
APP WIZARD から出力されるフレームワークのファイルは紹介を省かせていただきます。
今回から、オブジェクトのためのクラスを別ファイル prim.h、prim.cpp に分けました。
| render.cpp | レイトレ用の描画関数群 |
| render.h | レイトレ用の描画関数群 |
| prim.cpp | オブジェクトのデータのための関数 |
| prim.h | オブジェクトのデータのための関数 |
| mainDlg.cpp | ダイアログを管理するクラスのメソッドが書かれたファイル |
| mainDlg.h | ダイアログを管理するクラスのヘッダ |
あと、実行ファイル、リソースファイル、プロジェクトファイルが入っています。
モンテカルロ法とは、乱数を使って計算を何度も行い、その結果の合計を試行回数で割った平均（アンサンブル平均）を正式な計算の結果として使うというものです。
レイトレースに関していえば、本来は、レイが物体にあたったときは、BRDFに従ってさまざまな方向からの光を寄せ集めてレイに入射する光を計算しなくてはなりません。
それはめんどくさいので、大まかに反射、屈折および拡散の現象に分けて、それらの反射率に従って１つのレイを飛ばして、色を求めます。
このとき、反射、屈折、拡散のどの種類の反射をするのかは乱数を振るまでわからないので、ロシアンルーレット式の方法とも呼ばれます。
さらに、拡散反射の場合には、四方八方の光を調べなくてはならないのですが、ここにも確率的な考え方を取り入れます。一度の試行では、拡散反射についてでもレイを１方向だけしか飛ばしません。ただし方向をランダムにすることによって全方位からの光を取り入れます。
たとえば、ちょっとつるつるな物質を考えます。たとえば、2:8の割合で光を薄く反射する物質だとします。
鏡面反射率 : 0.2 拡散反射率 : 0.8 屈折反射率 : 0.0
この物質にレイがあたったとき、まず０～１の範囲で乱数を計算します。 乱数の値が 0.2 よりも小さな時は、次のレイは反射ベクトルに従って反射します。 乱数の値が 0.2～1.0の時は、拡散反射の計算を行います。 つまり、てんでばらばらな方向にレイを飛ばして多重散乱の計算を行います。
あと、今まで出てこなかった材質について考えなくてはなりません、ライトのような光を発する物質です。 このような物質へレイがあたった場合には、そこから先はレイは反射せずに、ライトの色をレイに入射する色に与えて、計算を終了します。
というか、今回は、Blinn-Phongシェーディングは捨て去ってしまいたいので、レイの計算が終わる方法は（レイが物体にあたらないで空間を突き抜けた場合以外には）放射する物体へ当たるしかほかにありません。
視点から出たレイはいろいろな物質を反射して、最後にライトに交差した時点で計算は終了し、その経路にわたった色がその試行の色として採用されます。
といっても、何度も反射して光にたどり着かないと計算時間を食うだけなので、適当な回数の反射で計算を打ち切ることにします。
もちろん、CG的には、放射しつつ光を反射する材質もあるでしょう。 したがって、放射の割合も確率事象に含める必要があります。
モンテカルロ法は、乱数が性質のよいものならば、試行回数をNとして 1/√N で誤差が消えていき、無限回の試行を行ったとすれば正しい値に収束することが知られています（中央極限定理）。
で、実際にやってみると、待ち時間に耐えられる程度の試行回数ではノイズがおおいんだなこれが…
今回も、質感情報の構造体 Material を変更しました。
ライトのような光を出す物質もプリミティブで表現することにします。そのための、物質がどの程度光を出しているかを表現するためのパラメータ emmisive を追加しました。
また、Blinn-Phong 鏡面反射は使わないので、（主に）拡散反射のための色のパラメータを１つだけ残しておきました。
prim.h
0021: typedef struct{
0022: float COLOR_DIF[4];
0023: float reflection; // 鏡面反射率
0024: float refraction; // 屈折反射率
0025: float diffuse; // 屈折反射率
0026: float emmisive; // 放射率
0027: }Material;
たとえば、次のようなデータになります。拡散や反射、屈折しかしない物体は、emmisive の項に0が追加されるだけです。
render.cpp
0072: OBJ_DATA sphere_data[] = {
0073: {// 小さいのに乗ってるほう
0074: OBJ_TYPE_SPHERE,
0075: {
0076: {0.10f, 0.100f, 0.10f},
0077: 0.2f, 0.8f, 0.0f, 0.0f,
0078: },{{0,0,0},{0,0,0},{0,0,0},},{
0079: {125.5, 0+100,169},
0080: // {185.5, 165+100,169},
0081: 100.0f,
0082: }
0083: },{// 大きいのに乗ってるほう
0084: OBJ_TYPE_SPHERE,
0085: {
0086: {0.05f, 0.050f, 0.05f},
0087: 0.80f, 0.00f, 0.20f, 0.0f,
0088: },{{0,0,0},{0,0,0},{0,0,0},},{
0089: {368.5, 330+100,351},
0090: 100.0f,
0091: }
0092: }};
ライトのプリミティブは、emmisive が1.0で、その他が0になります。 また、色のパラメータを光の強さにも使うことにしました。
render.cpp
0093: OBJ_DATA room_data[14] = {
0094: {
0095: // ライト１
0096: OBJ_TYPE_TRIANGLE,
0097: {
0098: {100.0f*1.00f, 100.0f*0.90f, 100.0f*0.50f},
0099: 0.0f, 0.0f, 0.0f, 1.0f,
0100: },{
0101: {213.0, 548.799f, 227.0},
0102: {213.0, 548.799f, 332.0},
0103: {343.0, 548.799f, 227.0},
0104: }
0105: },{// ライト２
0106: OBJ_TYPE_TRIANGLE,
0107: {
0108: {100.0f*1.00f, 100.0f*0.90f, 100.0f*0.50f},
0109: 0.0f, 0.0f, 0.0f, 1.0f,
0110: },{
0111: {343.0, 548.799f, 227.0},
0112: {213.0, 548.799f, 332.0},
0113: {343.0, 548.799f, 332.0},
0114: }
0115: },{ // 左１
0116: OBJ_TYPE_TRIANGLE,
0117: {
0118: {0.70f, 0.15f, 0.15f},
0119: 0.0f, 0.0f, 1.0f, 0.0f,
0120: },{
0121: {552.8f, 0.0f, 0.0f},
0122: {556.0f, 548.8f, 0.0f},
0123: {549.6f, 0.0f, 559.2f},
0124: }
0125:
さて、レイトレースする部分を見ていきましょう。
まずは、探査する深さが深くなりすぎたときに適当に計算を打ち切るための判定を行います。
今回は、反射、屈折の時には深さは深くしないで、拡散反射のときだけ再帰の深さが深くなるようにしました。
render.cpp
0606: bool GetColor(D3DXVECTOR3 *dest, const D3DXVECTOR4 *x, const D3DXVECTOR4 *v, int depth)
0607: {
0608: D3DXVECTOR3 diffuse_color=D3DXVECTOR3(0,0,0);
0609:
0615: // -----------------------------------------------------------------------
0616: // 再帰しすぎた場合は光が届いていないものと考える
0617: // -----------------------------------------------------------------------
0618: const int DEPTH_MAX = 3;
0619: if(DEPTH_MAX <= depth) { *dest=D3DXVECTOR3(0,0,0); return FALSE;}
そしたら、交差判定です。 ここはいままでも出てきたように計算します。
render.cpp 0621: // ----------------------------------------------------------------------- 0622: // 視線から交点を求める 0623: // ----------------------------------------------------------------------- 0624: float t = CPrimitive::INFINTY_DIST; 0625: CPrimitive *pObj = NULL; 0626: D3DXVECTOR4 p, n;// 交点の位置と法線 0627: 0628: t = pShpereS ->IsAcross(t, &n, &p, &pObj, x, v); 0629: t = pShpereT ->IsAcross(t, &n, &p, &pObj, x, v); 0630: t = pRoom ->IsAcross(t, &n, &p, &pObj, x, v); 0632: t = pBlockTall ->IsAcross(t, &n, &p, &pObj, x, v);
オブジェクトにあたらなかった場合は、背景色を出力します。 （今回は、手前もふさいでしまったので、この効果は実際には働きません）
render.cpp
0634: if( NULL == pObj ){
0635: // 視線はオブジェクトから外れた！
0636: *dest = BG_COLOR;
0637: return TRUE;
0638: }
次に、ロシアンルーレットをして、どのような反射をするか決めます。
反射、屈折、拡散、放射の各確率から、今回の試行がどのプロセスを行うか決定します。
より具体的には、それぞれの確率の和から乱数の上限を決定して、乱数の得られた値からどの種類の反射を行うか決定します。
ちなみに、吸収する確率は、1-cd-cr-cnになりますが、吸収の効果は今回はロシアンルーレットでは考慮せず、最後に色を出力するときに考慮します。
render.cpp 0640: // ----------------------------------------------------------------------- 0641: // 反射方向をロシアンルーレットで決める 0642: // ----------------------------------------------------------------------- 0643: float cd = pObj->m_material.diffuse; // 拡散率 0644: float cr = pObj->m_material.reflection; // 反射率 0645: float cn = pObj->m_material.refraction; // 屈折率 0646: float ce = pObj->m_material.emmisive; // 放射率 0647: float total = cd+cr+cn+ce; 0648: float prob = total * frand(); 0649: 0650: int type = 0;// 拡散 0651: if(prob < cr ) type = 1; else// 反射 0652: if(prob < cr+cn ) type = 2; else// 屈折 0653: if(prob < cr+cn+ce) type = 3; // 放射
拡散する時には、レイはランダムに飛ばします。ただし、物体にめり込むことはないので、新たなレイの方向 dir と、法線ベクトル n の内積を計算して、レイの方向が法線ベクトルと反対向きなら、レイの向きを反対にして、物体の面の方向にレイを飛ばすように調整します。
ランダムな方向にレイを飛ばして、再帰的に得た色の成分は、物体の色と、法線ベクトルとレイの方向の内積を掛けて出力する色とします。
これは、ランバートの余弦則そのままです。
render.cpp
0655: // -----------------------------------------------------------------------
0656: // 交点が、影に隠れているかどうか調べて、拡散光を算出する
0657: // -----------------------------------------------------------------------
0658: if(0==type){
0659: if(++depth<DEPTH_MAX){
0660: // さらに奥を検索
0661: D3DXVECTOR4 dir, pos;
0662:
0663: // 等方的にレイを飛ばす
0664: float theta = D3DX_PI * frand();
0665: float phi = 2.0f*D3DX_PI * frand();
0666: dir.x = sinf(theta) * cosf(phi);
0667: dir.y = sinf(theta) * sinf(phi);
0668: dir.z = cosf(theta);
0669: float dn = D3DXVec3Dot((D3DXVECTOR3 *)&dir, (D3DXVECTOR3 *)&n);
0670: if(dn<0){
0671: // 法線と反対向きに飛ぶつもりなら、方向を反対にする
0672: dn = -dn;
0673: dir.x *= -1;
0674: dir.y *= -1;
0675: dir.z *= -1;
0676: }
0677:
0678: pos = p + 0.01f*(dir);
0679: GetColor(&diffuse_color, &pos, &dir, depth);
0680: // プリミティブの色や拡散光の余弦側を適用
0681: diffuse_color.x *= pObj->m_material.COLOR_DIF[0] * dn * 0.5f*D3DX_PI;
0682: diffuse_color.y *= pObj->m_material.COLOR_DIF[1] * dn * 0.5f*D3DX_PI;
0683: diffuse_color.z *= pObj->m_material.COLOR_DIF[2] * dn * 0.5f*D3DX_PI;
0684: }
0685: }
反射する場合は、反射ベクトルを求めて、さらに探索します。
render.cpp
0687: // -----------------------------------------------------------------------
0688: // 再帰的に光線を飛ばして、より高次の反射を考慮に入れる
0689: // -----------------------------------------------------------------------
0690: // 反射
0691: D3DXVECTOR3 reflect_color=D3DXVECTOR3(0,0,0);
0692: if(1==type){
0693: // 反射ベクトル
0694: D3DXVECTOR4 r = *v - 2.0f*D3DXVec3Dot((D3DXVECTOR3 *)&n, (D3DXVECTOR3 *)v) * n;
0695: D3DXVECTOR4 pos = p + 0.01f*n;// ちょっと浮かす
0696: if(!GetColor( &reflect_color, &pos, &r, depth)){
0697: cr = 0;// 結局、光が届かなかった
0698: }
0699: }
屈折の場合は、屈折ベクトルを求めて、再帰的にレイトレースします。
render.cpp
0701: // -----------------------------------------------------------------------
0702: // 屈折
0703: // -----------------------------------------------------------------------
0704: D3DXVECTOR3 refract_color=D3DXVECTOR3(0,0,0);
0705: if(2==type){
0706: // 屈折ベクトル
0707: D3DXVECTOR4 t, pos;
0708: float VN = D3DXVec3Dot((D3DXVECTOR3 *)v, (D3DXVECTOR3 *)&n);
0709: float eta = 1.5f;
0710: if(VN<0){
0711: // 入射
0712: float D = 1-(1+VN)*(1+VN)/(eta*eta);
0713: if(D<0){cn=0;goto no_refract;}// 全反射
0714: t = (-VN/eta-sqrtf(D))*n+(1/eta)*(*v);
0715: pos = p - 0.01f*n;// ちょっと浮かす
0716: }else{
0717: // 出て行く
0718: float D = 1-(1-VN)*(1-VN)*(eta*eta);
0719: if(D<0){cn=0;goto no_refract;}// 全反射
0720: t = -(VN*eta-sqrtf(D))*n+eta*(*v);
0721: pos = p + 0.01f*n;// ちょっと浮かす
0722: }
0723: D3DXVec3Normalize( (D3DXVECTOR3 *)&t, (D3DXVECTOR3 *)&t );
0724: if(!GetColor( &refract_color, &pos, &t, depth)){
0725: cn = 0;// 結局、光が届かなかった
0726: }
0727: }
0728: no_refract:
放射光は、マテリアルの色そのものになります。
render.cpp 0730: // ----------------------------------------------------------------------- 0731: // 放射 0732: // ----------------------------------------------------------------------- 0733: D3DXVECTOR3 emmisive_color = *(D3DXVECTOR3*)&pObj->m_material.COLOR_DIF;
最後に、今までに計算した色をそれぞれの状況に応じて合成します。
放射光以外は、cd+cr+cn の係数を掛けました。これは、1-cd+cr+cn が光の吸収率になるので、そうでない成分を出力するように調整するためです。
render.cpp
0735: // -----------------------------------------------------------------------
0736: // 色を出力する
0737: // -----------------------------------------------------------------------
0738: switch(type){
0739: case 0:// 拡散
0740: *dest = (cd+cr+cn) * diffuse_color;
0741: break;
0742: case 1:// 反射
0743: *dest = (cd+cr+cn) * reflect_color;
0744: break;
0745: case 2:// 屈折
0746: *dest = (cd+cr+cn) * refract_color;
0747: break;
0748: case 3:// 放射
0749: *dest = ce * emmisive_color;
0750: break;
0751: }
0752:
0753: return TRUE;
0754: }
さて、モンテカルロ法は回数を重ねれば重ねるほど正確な画像に近づいていく手法です。
したがって、何度も何度もレンダリングすることにします。
今回は、１秒間処理したら終了するように Render() 関数を書き換えました。
timeGetTime() 関数を使えば、ミリ秒単位で現在時刻を知ることができるので、関数が呼ばれた時刻を保存しておいて、短い処理を終えるたびに時間を計測し、その時間が１秒を超えていたら関数を終了するようにします。
なお、timeBeginPeriod() で、timeGetTime() の時間分解能を変更できるのですが、時間が正確でもほとんど意味がないので、0.5 秒の時間分解能にしました。
render.cpp
0544: int Render()
0545: {
0546: if(STATE_IDLE==s_state) return 0;
0547:
0548: int i,j,no;
0549: D3DXVECTOR3 col;
0550: int ret = -1;
0551:
0552: timeBeginPeriod( 500 );// タイマー精度の設定
0553:
0554: for(unsigned long t = timeGetTime()
0555: ; timeGetTime()-t < 1000 // 1秒たったか
0556: && STATE_IDLE!=s_state;){ // 終了したか
0557:
0558: switch(s_state){
0559: case STATE_IDLE: // 待機
0560: break;
0561: case STATE_RENDER_BEGIN:// 初期化
0562: s_x = 0;
0563: s_y = 0;
0564: render_cnt+=1.0f;
0565: s_state = STATE_RENDER;
0566: break;
0567: case STATE_RENDER: // 画面の描画
0568: GetColor(&col, ((float)s_x+frand())/(float)RENDER_WIDTH
0569: , ((float)s_y+frand())/(float)RENDER_HEIGHT);
0570: s_total[4*(s_y*RENDER_WIDTH+s_x)+0]+=col.x;// R
0571: s_total[4*(s_y*RENDER_WIDTH+s_x)+1]+=col.y;// G
0572: s_total[4*(s_y*RENDER_WIDTH+s_x)+2]+=col.z;// B
0573:
0574: // 次のピクセルに移動
0575: if(RENDER_WIDTH<=++s_x){
0576: // 行を変える
0577: s_x = 0;
0578: if(RENDER_HEIGHT<=++s_y){
0579: s_state = STATE_RENDER_END;
0580: break;
0581: }
0582: }
0583: break;
0584: case STATE_RENDER_END: // 描画の終了
0585: // 画面の反映
0586: no = 0;
0587: for(j=0; j<RENDER_HEIGHT; j++){
0588: for(i=0; i<RENDER_WIDTH ; i++){
0589: s_data [no+0]=(char)(255.9*min(1,s_total[no+0]/render_cnt));// R
0590: s_data [no+1]=(char)(255.9*min(1,s_total[no+1]/render_cnt));// G
0591: s_data [no+2]=(char)(255.9*min(1,s_total[no+2]/render_cnt));// B
0592: no += 4;
0593: }
0594: }
0595: s_state = STATE_RENDER_BEGIN;// 無限ループ
0596: break;
0597: }
0598: }
0599:
0600: timeEndPeriod( 500 );// タイマー精度を戻す
0601:
0602: return ret;
0603: }
Render() 関数の中では、適当に状態を設定して、レンダリングの手順を整頓しています。
最初に STATE_IDLE の状態にいます。ここで、ボタンが押されるのを待ちます。
ボタンが押されたら、状態は STATE_RENDER_BEGIN になります。STATE_RENDER_BEGIN では、レイとレースするピクセルを左上に設定するなどのことをしています。
そうしたら、STATE_RENDER へ移ります。STATE_RENDER では、ピクセルごとにレイトレースをします。
得られたデータは、ピクセルごとに値を配列 s_total に加算していきます。
レイトレースが終わったら次のピクセルを指定するか、一通り画面をレイトレースしたら次の状態 STATE_RENDER_END へ移ります。
STATE_RENDER_END では、今までの結果を描画回数で割って、フレームバッファに色を出力します。STATE_RENDER_END を通過しない限り画面は更新されないので、今回は１ピクセルごとにレンダリングする姿は見れません。STATE_RENDER_END の次は、再びSTATE_RENDER_BEGIN へ移って次のパスをレンダリングします。
また、今回は何度もレンダリングするのですが、この方法に都合がいいアンチエイリアスも試みました。
それは、レイを飛ばすときに、ピクセルの中をランダムに飛ばす方法です。
いつもはピクセルの中央に対してレイを飛ばすのですが、
ピクセルの中を適当にずらしてレイをサンプリングしたほうがアンチエイリアスがかかりより正確な画像になります。
本当はサンプリングする場所も適当に決めたほうがいいのですが、モンテカルロ法を信じるなら、多くのレイを飛ばしたアンサンブル平均はピクセルの中を一様にレイを飛ばした結果と等しくなるので、今回はレイをランダムに飛ばしてサンプリングします。
render.cpp 0568: GetColor(&col, ((float)s_x+frand())/(float)RENDER_WIDTH 0569: , ((float)s_y+frand())/(float)RENDER_HEIGHT);
Render() 関数自体が１秒ごとに帰ってくる関数になったので、アプリケーションの呼び出し側も変更しましょう。
Render::Begin() 関数で描画を始めるのは同じですが、今回は Render() 関数が帰ってきたときに、イベントの処理や画面の再描画をしています。
また、現在何回目のレンダリングをしているかどうかをウィンドウのタイトルに表示するようにしました。SetWindowText() メソッドでタイトルを変更することができるので、レンダリング回数を所得できる関数を作っておいて、画面を更新するたびにタイトルを書き換えました。
mainDlg.cpp
0210: void CMainDlg::OnButtonRender()
0211: {
0212: bQuit = FALSE;
0213:
0214: // TODO: この位置にコントロール通知ハンドラ用のコードを追加してください
0215: for(Render::Begin(); Render::Render(); this->OnPaint()){
0216: // タイトル変更
0217: char title[256];
0218: sprintf(title, "%d回目", Render::GetRenderCount());
0219: SetWindowText(title);
0220: // メッセージの処理
0221: PumpMessages();
0222: // 終了判定
0223: if(bQuit) return;
0224: }
0225: this->OnPaint();
0226: }
さて、何回レンダリングするとどのような絵ができるのでしょうか？
ちなみに 256x256の大きさで Athron 2000+ で８０分で１００００回程度レンダリングできました。
|
|
１回のレンダリングでは、まばらにしか点がありません。つまり、ほとんどの計算が無駄になっています |
|
|
１０回のレンダリングでは、うっすらと色が見えますが、概要は全然わかりません。 |
|
|
５０回のレンダリングでは、なんとなく形が見えてきました。 集光効果はこの時点でもはっきりわかります |
|
| １００回のレンダリングでは、もう形がわかります。 |
|
| ５００回のレンダリングでは、形がはっきりとわかります。 これからがノイズとの戦いです。 |
|
| １０００回のレンダリングでは、透明な球もはっきりしだしています。。 |
|
| ５０００回のレンダリングでは、だんだんノイズも取れてきました。 |
|
| １３０００回のレンダリングでは、だいぶすべすべです。 でも、ノイズは目立ちますね。 |
|
| ５００００回のレンダリングでは、よりよくなっていますが。 ノイズが取れる気はしません… |
この後、８００００回ほどレンダリングしたら、マシンがフリーズしてしまいました。
比較として、計算打ち切り深度を深くして計算してみました。 １００回、レイを反射、屈折させても光源にたどり着かなかったときに、計算をやめるように変更して、プログラムをまわしてみたら５０００８回まわした時点で次のようになりました。
こちらのほうが間接光が強く出ているカナ？
これでグローバルイルミネーションを取り込むことができました。
といっても間接照明はかなり弱いですね。
このモンテカルロ法はノイズが多すぎで時間もかかります。
このままでは、とてもリアルタイムにはなりそうもありません…