日本語のコメントをつけてみて, 整形した smallpt.cpp
// 権利表記
// smallpt, a Path Tracer by Kevin Beason, 2008
//
// Japanese commments are added by Gregory Igehy, 2014
// ビルド方法について
// Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
// Remove "-fopenmp" for g++ version < 4.2
// 実行方法(1 ピクセルあたり 5000 個のサンプル )と画像のプレビュー例
// Usage: time ./smallpt 5000 && xv image.ppm
// 標準 C ライブラリのインクルードファイル
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
// 3 次元ベクトル (x, y, z) の構造体
struct Vec
{
// メンバ変数
double x, y, z;
// コンストラクタ
Vec( double x_ = 0, double y_ = 0, double z_ = 0)
{
x = x_; y = y_; z = z_;
}
// ベクトルの加算
Vec operator+(const Vec &b) const
{
return Vec(x+b.x,y+b.y,z+b.z);
}
// ベクトルの減算
Vec operator-(const Vec &b) const
{
return Vec(x-b.x,y-b.y,z-b.z);
}
// ベクトルにスカラーを乗算
Vec operator*(double b) const
{
return Vec(x*b,y*b,z*b);
}
// ベクトルに他のベクトルの要素を乗算
Vec mult(const Vec &b) const
{
return Vec(x*b.x,y*b.y,z*b.z);
}
// 自分自身を正規化
Vec& norm()
{
return ( * this ) = (* this) * (1 / sqrt( x * x + y * y + z * z ) );
}
// ベクトルの内積
double dot(const Vec &b) const
{
return x * b.x + y * b.y + z * b.z ;
}
// ベクトルの外積
Vec operator%(Vec&b)
{
return Vec( y * b.z - z * b.y,
z * b.x - x * b.z,
x * b.y - y * b.x );
}
};
// レイの構造体
struct Ray
{
// メンバ変数 : レイの開始位置座標 o と 方向ベクトル d
Vec o, d;
// コンストラクタ
Ray(Vec o_, Vec d_) :
o(o_), d(d_)
{
}
};
// マテリアル用のレイの反射の種類の列挙定数
enum Refl_t
{
// 完全なディフューズ反射のみ,
DIFF,
// 完全に滑らかなスペキュラ反射のみ
SPEC,
// 完全な屈折も考慮した SPEC
REFR
};
// 球の構造体
struct Sphere
{
// 半径
double rad;
// 中心の位置座標, エミッション, ベースカラー
Vec p, e, c;
// 反射モデルの種類
Refl_t refl;
// コンストラクタ
Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_) :
rad(rad_), p(p_), e(e_), c(c_), refl(refl_)
{
}
// レイとのヒット判定をするメソッド
// もしヒットした場合はレイの開始点からヒット点までの距離(t)を返す
// ヒットしない場合は 0 を返す
double intersect(const Ray &r) const
{
// 2 次の判別式を使ったヒット判定を行うために
// t^2 * d.d + 2 * t * ( o - p ).d + ( o - p ).( o - p ) - R^2 = 0
// を解く
// op は レイの開始点 o から 球の中心の位置座標 p に向かうベクトル
Vec op = p - r.o;
double t,
// エプシロン
eps = 1e-4,
// 2 次の判別式 det
b = op.dot( r.d ),
det = b * b - op.dot( op )+ rad * rad;
// 解がない場合
if ( det < 0 )
{
return 0;
}
// 解がある場合
else
{
det = sqrt( det );
// 解は レイの進行方向 ( t > 0 )かつ
// 解の小さい方を返す
return ( ( t = b - det ) > eps ) ?
t :
( ( t = b + det ) > eps ? t : 0 );
}
}
};
////////////////////////////////////////////////////////////////////////////////
// 部屋のジオメトリである "球の構造体の配列"
//
// 順に: 球の半径, 球の中心の位置座標(ワールド座標系),
エミッション, ベースカラー, マテリアルの反射の種類
Sphere spheres[] =
{
// 左の壁
Sphere( 1e5 ,Vec( 1e5+1,40.8,81.6) ,
Vec(), Vec(.75,.25,.25) , DIFF),
// 右の壁
Sphere( 1e5 , Vec(-1e5+99,40.8, 81.6),
Vec(), Vec(.25,.25,.75) ,DIFF),
// 背面の壁
Sphere( 1e5, Vec(50,40.8, 1e5),
Vec(), Vec(.75,.75,.75) , DIFF),
// 正面の壁
Sphere( 1e5, Vec(50,40.8,-1e5+170),
Vec(), Vec() , DIFF),
// 床面
Sphere( 1e5, Vec(50, 1e5, 81.6),
Vec(), Vec(.75,.75,.75) , DIFF),
// 天井
Sphere( 1e5, Vec(50,-1e5+81.6,81.6),
Vec(), Vec(.75,.75,.75) , DIFF),
// ミラーの球
Sphere( 16.5, Vec(27,16.5,47),
Vec(), Vec(1,1,1)*.999 , SPEC),
// ガラスの球
Sphere( 16.5, Vec(73,16.5,78),
Vec(), Vec(1,1,1)*.999 , REFR),
// 天井のライト(エミッションしている球)
Sphere( 600, Vec(50,681.6-.27,81.6), Vec(12,12,12),
Vec() , DIFF)
};
////////////////////////////////////////////////////////////////////////////////
// 便利関数 : double を 0 ～ 1 にクランプ
inline double clamp(double x)
{
return ( x < 0 ) ? ( 0 ) : ( ( x > 1 ) ? 1 : x );
}
// 便利関数 : double を 整数に変換
// ( 具体的には double を 0～1 にクランプしてから, ガンマ補正済みの整数(0～255)に変換 )
inline int toInt(double x)
{
return int(
// クランプ後にガンマ補正
pow( clamp(x), 1 / 2.2 )
* 255 + .5
);
}
// レイがシーン中のジオメトリとヒットしたか ? を判定する関数
inline bool intersect( const Ray& r, double& t, int& id )
{
double n = sizeof( spheres ) / sizeof( Sphere ),
d,
// 大きな値で初期化しておく
inf = t =1e20;
// 球の構造体を走査するループ
for( int i = int( n ); i-- ; )
{
// レイが球にヒットしている場合 かつ
if ( ( d = spheres[i].intersect( r ) ) &&
// 今までのよりもヒット時のレイの距離 t が小さい場合
( d < t ) )
{
// ヒット時のレイの開始位置とヒット点間の距離 t を保存
t = d;
// ヒットした球のインデックスを保存
id = i;
}
}
// t が初期値よりも小さかったら, ヒットしていると判定
return ( t < inf );
}
////////////////////////////////////////////////////////////////////////////////
// 指定したレイで, シーンをパストレースした結果の 放射輝度(radiance) を返す再帰関数
Vec radiance( const Ray& r,
// デプスはこのレイの反射や屈折の回数
int depth,
// 乱数を取得するための種
unsigned short* Xi )
{
// レイがヒットした位置とレイの開始地点間の距離
double t;
// レイがヒットした球の配列のインデックス
int id = 0;
// レイがジオメトリとヒットしないときは, 背景色の黒を返す
if ( ! intersect( r, t, id ) )
{
return Vec();
}
// レイがヒットした球
const Sphere &obj = spheres[ id ];
// レイがヒットした点の位置座標(ワールド座標系)
Vec x = r.o + r.d * t,
// レイがヒットした点の法線(ワールド座標系)
n = ( x - obj.p ).norm(),
// ヒットした点を中心とした半球座標系を作るための法線の向き
// レイが球の内側からヒットした場合だけ, 法線を反転する
nl = ( n.dot( r.d ) < 0 ) ? n : ( n * ( - 1 ) ),
// マテリアルのベースカラー RGB
f = obj.c;
// マテリアルのベースカラー RGB のそれぞれの要素のうち, 最大値を求める
double p = ( ( f.x > f.y ) && ( f.x > f.z ) ) ?
f.x :
( ( f.y > f.z ) ? f.y : f.z );
// インクリメント後のデプスが 5 より大きい場合
if ( ++depth > 5 )
{
// 乱数とベースカラーの最大値でロシアンルーレット
// erand48 は [ 0.0, 1.0 ) の乱数を返す
// ( 詳しくは http://kaworu.jpn.org/doc/FreeBSD/jman/man3/erand48.3.php を参照 )
if ( erand48( Xi ) < p )
{
// モンテカルロ法
// (洋書 "Realistic Ray Tracing" の p149-150 を参照 )
// ∫ ( f * p ) du ≒ ( 1 / N ) * Σ f [0]
//
// f * p = f' とする [1]
//
// [0] に [1]を代入して f を消去すると
// ∫ f' du ≒ ( 1 / N ) * Σ ( f' / p ) [2]
//
// 最後に [2] の f' を f で置換すると
// ∫ f du ≒ ( 1 / N ) * Σ ( f / p ) [3]
//
f = f * ( 1 / p );
}
else
{
// 反射しないので, エミッションの放射輝度を返して終了
return obj.e;
}
}
// ディフューズシェーディングのみの場合
if ( obj.refl == DIFF )
{
// 半径が 1 の円板をランダムサンプリング
//
// r1 : 円板上の角度の乱数
double r1 = 2 * M_PI * erand48( Xi ),
// r2s : 円板の中心の位置座標からの半径
r2 = erand48( Xi ),
r2s = sqrt( r2 );
// 半球座標系 ( u, v, w )
// ヒット点の位置座標を中心として, w が天頂を向いている半球
Vec w = nl,
u = ( ( fabs( w.x ) > 0.1 ? Vec( 0, 1 ) : Vec( 1 ) )
% w
).norm(),
v = w % u;
// 反射レイの方向のサンプリング
// 円板のサンプリング結果と 半球座標系 を利用
// (洋書 "Realistic Ray Tracing" の p153-154 を参照 )
Vec d = ( u * cos( r1 ) * r2s + v * sin( r1 ) * r2s + w * sqrt( 1 - r2 ) ).norm();
// ヒット点のマテリアルのエミッションの放射輝度 と
return obj.e +
// 反射レイの放射輝度 の両方
f.mult( radiance( Ray( x , d ), depth, Xi ) );
}
// 完全なスペキュラ反射
else if ( obj.refl == SPEC )
{
// ヒット点のマテリアルのエミッションの放射輝度 と
return obj.e +
// 完全スペキュラ反射の反射レイによる放射輝度
f.mult(
radiance( // 完全スペキュラ反射の反射レイ
Ray( x, r.d - n * 2 * n.dot( r.d ) ),
depth, Xi )
);
}
// 完全な屈折をする際の屈折レイ
Ray reflRay( x, r.d - n * 2 * n.dot( r.d ) );
// レイが ジオメトリの外側から内側に向かって進入しているときに true
bool into = n.dot( nl ) > 0;
// 真空の屈折率
double nc = 1,
// ガラスの屈折率 ( http://ww1.tiki.ne.jp/~uri-works/tmp/ )
nt = 1.5,
// 相対屈折率
nnt = into ? ( nc / nt ) : ( nt / nc ),
ddn = ( r.d ).dot( nl ),
cos2t;
// 完全な全反射の場合
if ( ( cos2t = 1 - nnt * nnt * ( 1 - ddn * ddn ) ) < 0 )
{
// ヒット点のマテリアルのエミッションによる放射輝度 と
return obj.e +
// 全反射した反射レイによる放射輝度 の両方
f.mult( radiance( reflRay, depth, Xi ) );
}
// 屈折レイの方向
Vec tdir = ( r.d * nnt -
n * ( ( into ? 1 : ( - 1 ) ) *
( ddn * nnt + sqrt( cos2t ) ) ) ).norm();
// R0 : 垂直入射時の方向のフレネル反射率 F(0°)
double a = nt - nc,
b = nt + nc,
R0 = a * a / ( b * b ),
// c : 1 - cos θ
c = 1 - ( into ? ( - ddn ) : tdir.dot( n ) );
// Re : フレネル反射率 F
double Re = R0 + ( 1 - R0 ) * c * c * c * c * c,
// Tr : 透過率
Tr = 1 - Re,
// P : 反射する確率
// P = Re にすると, ( depth > 2 ) のときに
// Re が小さいときや反射レイの放射輝度が高いときに問題が起きる
// 従って, 妥協策として
// P = k / 2 + ( 1 - k ) * Re
// として, k をヒューリスティックに 0.5 とした
// (洋書 "Realistic Ray Tracing" の p177-178 を参照 )
//
P = 0.25 + 0.5 * Re,
// 下でのロシアンルーレットを考慮して, モンテカルロ法として確率で除算した値
RP = Re / P,
TP = Tr / ( 1 - P );
// ヒット点のマテリアルのエミッションによる放射輝度 と
return obj.e +
f.mult( // デプスによる条件分岐
( depth > 2 )
?
// デプスが 2 より大きい場合はロシアンルーレット
( ( erand48( Xi ) < P ) ?
// 反射レイの場合の放射輝度 あるいは
radiance( reflRay , depth, Xi) * RP :
// 屈折レイの場合の放射輝度 のどちらか
radiance( Ray( x, tdir ), depth, Xi) * TP )
:
// デプスが 0,1,2 の場合
// 反射レイの放射輝度 と
radiance( reflRay , depth, Xi ) * Re +
// 屈折レイの放射輝度 の両方
radiance( Ray ( x, tdir ), depth, Xi ) * Tr );
}
////////////////////////////////////////////////////////////////////////////////
// エントリポイント
int main(int argc, char *argv[])
{
// レイトレースで作る画像の大きさは 1024 x 768 ピクセル
int w = 1024, h = 768,
// ピクセルのサンプル数が argv[ 1 ] として指定された場合は使う
//
// 但し, サブピクセルが 2x2 = 4 なので
// サブピクセルのサンプル数 = ピクセルのサンプル数 / 4 とする
samps = ( argc == 2 ) ? ( atoi( argv[ 1 ] ) / 4 ) : 1;
// カメラの座標系 ( cam, cx, cy の 3 軸)
/
// カメラのレイを
// カメラの位置座標 (ワールド座標系 ) と
Ray cam( Vec( 50, 52, 295.6 ),
// カメラが向いている向き で初期化
Vec( 0, - 0.042612, - 1 ).norm() );
// cx : カメラの真右方向のベクトル
// ( 0.5135 は fovy )
Vec cx = Vec( w * .5135 / h),
// cy : カメラの真上方向のベクトル
cy = ( cx % cam.d ).norm() * .5135,
// r : サブピクセルごとの放射輝度を保存しておく一時変数
r,
// レイトレースの結果を保存するためのピクセル数分の RGB の配列
*c = new Vec[ w * h ];
// OpenMP の指示文
// omp parallel for で次の for ループを実行スレッドの数だけ並列化させる
//
// まず ここでのチャンクは(1 つのスレッドが処理を行う最小単位である)ループの反復回数
//
// dynamic の場合は処理を終了した順に(遊んでいる)スレッドに対して, 次のチャンクの割り当てる
// また チャンクはループ 1 回分と明示的に指定している
//
// 変数 r はスレッドごとにプライベートになるように指定
//
// ( 詳しくは "C 言語による OpenMP 入門" 黒田久泰
http://www.cc.u-tokyo.ac.jp/support/kosyu/03/kosyu-openmp_c.pdf を参照 )
#pragma omp parallel for schedule( dynamic, 1 ) private( r )
// y = 0～767 ピクセルの縦方向のループ
for (int y = 0; y < h; y++ )
{
// 進捗度合いを標準エラーに出力
fprintf( stderr, "\rRendering (%d spp) %5.2f%%",
samps*4, 100. * y / ( h - 1 ) );
// x = 0 ～ 1023 ピクセルの横方向のループ
for ( unsigned short x = 0,
// 乱数用の種
Xi[ 3 ]={ 0, 0 ,y * y * y };
x < w; x++ )
{
// ピクセルを 2x2 のサブピクセルの分割
//
// サブピクセル sy = 0～1の縦方向のループ
for ( int sy = 0,
// 結果の書き込み先の配列のインデックス
i = ( h - y - 1 ) * w + x;
sy < 2; sy++ )
{
// サブピクセル sx = 0～1 の横方向のループ
for ( int sx = 0; sx < 2;
sx++,
// サブピクセルごとの放射輝度を保存する一時変数 r を
// (0.0, 0.0, 0.0) で初期化しておく
r = Vec() )
{
// サブピクセルごとのサンプル数 s = 0～ ( samps - 1 )のループ
for (int s = 0; s < samps; s++ )
{
// サンプリング
// サンプリングに相関があるとエイリアシングの問題が発生
// そこで, ジッタリングや層化抽出法を利用
//
// テントフィルタ
// p( x ) = - x + 1 ( if 0 <= x <= 1 )
// x + 1 ( if -1 <= x < 0 )
// 0 ( if x > 1 or x < -1 )
// ピクセルの中心に近いサンプルの重要度を上げる非一様サンプリング
// 重みを可変にする代わりに, サンプリングの分布を非一様にしている
// (洋書 "Realistic Ray Tracing" の p54-56 と p60-61 を参照 )
//
double r1 = 2 * erand48( Xi ),
dx = ( r1 < 1 ) ? ( sqrt( r1 ) - 1 ) : ( 1 - sqrt( 2 - r1 ) );
//
double r2 = 2 * erand48( Xi ),
dy = ( r2 < 1 ) ? ( sqrt( r2 ) - 1 ) : ( 1 - sqrt( 2 - r2 ) );
// サンプルを貫くカメラのレイを求める
//
// cx : カメラの真右方向のベクトル
// cy : カメラの真上方向のベクトル
// cam.d : カメラが向いている向き
//
// ピクセルの中心はスクリーン座標系では整数
// ピクセルの範囲はそれの縦横 -0.5～0.5
Vec d = cx * ( ( ( sx + .5 + dx ) / 2 + x ) / w - .5 ) +
cy * ( ( ( sy + .5 + dy ) / 2 + y ) / h - .5 ) +
cam.d;
// カメラのレイの方向の放射輝度を求める
// カメラのレイの開始位置を少し前にずらして, 部屋の内部に入れる
r = r + radiance( Ray ( cam.o + d * 140,
d.norm() ),
0, Xi) *
// モンテカルロ法
// 放射輝度の合計値を, サブピクセルごとのサンプル数で除算
// (洋書 "Realistic Ray Tracing" の p149-150 を参照 )
( 1. / samps );
}
// ピクセルごとの放射輝度に, サブピクセルごとの放射輝度を加算
c[i] = c[i] +
// モンテカルロ法
// 放射輝度の合計値を, ピクセルごとのサブピクセル数 ( 4 個 )の逆数で乗算
Vec( clamp( r.x ), clamp( r.y ), clamp( r.z ) ) * .25;
} // for sx
} // for sy
} // for x
} // for y
//
// 以下は計算した結果を, ppm 画像としてファイル出力する処理
//
// image.ppm のファイルをオープン
FILE* f = fopen("image.ppm", "w");
// PPM (P3) のヘッダの記述
fprintf( f, "P3\n%d %d\n%d\n", w, h, 255);
// 各ピクセルのループ
for ( int i = 0; i < w * h; i++ )
{
// double 型の RGB を, 整数型に変換
fprintf( f, "%d %d %d ",
toInt( c[i].x ), toInt( c[i].y ), toInt( c[i].z ) );
}
}
