Note
//
// expanded smallppm (code is exactly the same as smallppm.cpp but with more comments)
//
// Japanese comments are added and formatted by Gregory Igehy
//
// smallppm, Progressive Photon Mapping by T. Hachisuka
// originally smallpt, a path tracer by Kevin Beason, 2008
// Usage: ./smallppm 100000 && xv image.ppm
// ^^^^^^:number of photons emitted
#define PI ( ( double ) 3.14159265358979 )
// PPM の α のパラメータ
// 繰り返しごとに追加するフォトン数をスケールする
#define ALPHA ( ( double ) 0.7 )
// ハルトン・シーケーンス (with reverse permutation)
int primes[ 61 ] =
{
2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,
83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,
191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283
};
// 便利関数
inline int rev( const int i, const int p )
{
if ( i == 0 )
{
return i;
}
else
{
return p - i;
}
}
// ハルトン・シーケンスの関数
double hal( const int b, int j )
{
const int p = primes[ b ];
double h = 0.0,
f = 1.0 / ( double ) p,
fct = f;
while ( j > 0 )
{
h += rev( j % p, p ) * fct;
j /= p;
fct *= f;
}
return h;
}
// 3 次元のベクトルの構造体
struct Vec
{
double x, y, z;
Vec(double x_ = 0, double y_ = 0, double z_ = 0) {x = x_; y = y_; z = z_;}
inline Vec operator+(const Vec &b) const {return Vec(x+b.x, y+b.y, z+b.z);}
inline Vec operator-(const Vec &b) const {return Vec(x-b.x, y-b.y, z-b.z);}
inline Vec operator+(double b) const {return Vec(x + b, y + b, z + b);}
inline Vec operator-(double b) const {return Vec(x - b, y - b, z - b);}
inline Vec operator*(double b) const {return Vec(x * b, y * b, z * b);}
inline Vec mul(const Vec &b) const {return Vec(x * b.x, y * b.y , z * b.z);}
inline Vec norm() {return (*this) * (1.0 / sqrt(x*x+y*y+z*z));}
inline 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);}
};
// 最大値を返すマクロ
#define MAX( x, y ) ( ( x > y ) ? x : y )
// AABB の構造体
struct AABB
{
Vec min, max;
// 位置座標 p を考慮してフィットする
inline void fit( const Vec& p )
{
if ( p.x < min.x )
{
min.x = p.x;
}
if ( p.y < min.y )
{
min.y = p.y;
}
if ( p.z < min.z )
{
min.z = p.z;
}
max.x = MAX( p.x, max.x );
max.y = MAX( p.y, max.y );
max.z = MAX( p.z, max.z );
}
// 初期化
inline void reset()
{
min = Vec( 1e20, 1e20, 1e20 );
max = Vec( - 1e20, - 1e20, - 1e20 );
}
};
// カメラレイからヒットした点の構造体
// L(S)*D
struct HitPoint
{
// 重み (BRDF とピクセルのフィルタリングの値 )
Vec scale_factor,
// 位置座標
pos,
// 法線
nrm,
// 光束
flux;
// 密度推定用の半径^2
double radius2;
// n = ( N / ALPHA )
// N : 追加するフォトンの数
unsigned int photon_n;
// ヒット点が該当するピクセルのインデックス
int pixel_index;
};
// リストの要素
struct ListElement
{
// ヒット点の構造体 i へのポインタ
HitPoint* hp;
// 次のリスト要素へのポインタ
ListElement* next;
};
// 新規のリストの要素を確保して, 初期化する
ListElement* CreateListElement( HitPoint* hp,
// 次のリストの要素へのポインタ
ListElement* list_element )
{
ListElement* p = new ListElement;
p->hp = hp;
p->next = list_element;
return p;
}
// ハッシュテーブルの最大数を設定
unsigned int num_hash,
// カメラレイを飛ばすときのピクセルのインデックス
pixel_index,
num_photon;
// グリッドセルのサイズの逆数
double hash_s;
// ヒット点のハッシュのグリッド
// ( フォトンレイが 1 つのグリッドセルに対応し,
// それが影響を与えるヒット点のリストを作っておく )
ListElement** hitpoints_hash_grid;
// ヒット点のリストの先頭へのポインタ
ListElement* hitpoints_head = NULL;
// ヒット点の AABB
AABB hit_point_aabb;
// 空間のインデックスから,ハッシュ値を求める関数
inline unsigned int hash(const int ix, const int iy, const int iz)
{
return ( unsigned int )
( ( ( ix * 73856093 ) ^ ( iy * 19349663 ) ^ ( iz * 83492791 ) )
%
num_hash );
}
// ハッシュグリッドの作成
// - それぞれのグリッドセルに, 影響が及ぶヒット点を登録する
void build_hitpoints_hash_grid( const int w, const int h )
{
// AABB のサイズを, 全てのヒット点から求める
hit_point_aabb.reset();
ListElement* lst = hitpoints_head;
while ( lst != NULL )
{
// リストのループに使う変数
HitPoint* hp = lst->hp;
lst = lst->next;
hit_point_aabb.fit( hp->pos );
}
// 半径の初期値を推定
Vec ssize = hit_point_aabb.max - hit_point_aabb.min;
double irad = ( ( ssize.x + ssize.y + ssize.z) / 3.0 )
/
( ( w + h ) / 2.0 ) * 2.0;
// AABB をリセットする
hit_point_aabb.reset();
// ハッシュテーブルの大きさを求める
lst = hitpoints_head;
int vphoton = 0;
while ( lst != NULL )
{
// リストのループに使う変数
HitPoint *hp = lst->hp;
lst = lst->next;
// 値を初期化しておく
hp->radius2 = irad * irad;
hp->photon_n = 0;
hp->flux = Vec();
vphoton++;
// AABB を作り直す
hit_point_aabb.fit( hp->pos - irad );
hit_point_aabb.fit( hp->pos + irad );
}
// グリッドセルの大きさは, 初期半径の 2 倍にする
// ( make each grid cell two times larger than the initial radius )
hash_s = 1.0 / ( irad * 2.0 );
// ハッシュテーブルの最大数を設定
num_hash = vphoton;
// ハッシュテーブルの構築
hitpoints_hash_grid = new ListElement* [ num_hash ];
// nullptr で初期化しておく
for ( unsigned int i = 0; i < num_hash; i++ )
{
hitpoints_hash_grid[ i ] = NULL;
}
// ヒット点自身が,周囲から影響を受ける全てのグリッドセルに自分を登録する
lst = hitpoints_head;
while ( lst != NULL )
{
// ヒット点を順に,ループするための変数
HitPoint *hp = lst->hp;
lst = lst->next;
// (ヒット点が影響を受ける) グリッドセルの範囲のインデックス
Vec BMin = ( ( hp->pos - irad ) - hit_point_aabb.min ) * hash_s;
Vec BMax = ( ( hp->pos + irad ) - hit_point_aabb.min ) * hash_s;
for ( int iz = abs( int( BMin.z ) ); iz <= abs( int( BMax.z ) ); iz++ )
{
for ( int iy = abs( int( BMin.y ) ); iy <= abs( int( BMax.y ) ); iy++ )
{
for ( int ix = abs( int( BMin.x ) ); ix <= abs( int( BMax.x ) ); ix++ )
{
// グリッドセルのインデックスから,ハッシュ値を取得
int hv = hash( ix, iy, iz );
// グリッドセルに, ヒット点である自分をリストに登録する
hitpoints_hash_grid[ hv ] = CreateListElement( hp, hitpoints_hash_grid[ hv ] );
}
}
}
}
}
// レイの構造体
struct Ray
{
Vec o, d;
Ray()
{
};
Ray( const Vec& o_, const Vec& d_) :
o( o_ ), d( d_ )
{
}
};
// 反射モデル
enum Refl_t
{
DIFF,
SPEC,
REFR
};
// 球の構造体
struct Sphere
{
// 半径
double rad;
// 中心の位置座標と, ベースカラー
Vec p, c;
// 反射モデルの種類
Refl_t refl;
Sphere( const double r_, const Vec& p_, const Vec& c_, const Refl_t re_ ) :
rad( r_ ), p( p_ ),
c( c_ ), refl( re_ )
{
}
// 球とレイとのヒット判定をする関数
inline double intersect( const Ray& r ) const
{
Vec op = p - r.o;
double t, b = op.dot(r.d),
det=b*b-op.dot(op)+rad*rad;
if ( det < 0 )
{
return 1e20;
}
else
{
det = sqrt(det);
}
return ( ( t = b - det ) > 1e-4 ) ?
t :
( ( ( t = b + det ) > 1e-4 ) ?
t : 1e20 );
}
};
// シーン中の球の配列
Sphere sph[] =
{
// Scene: radius, position, color, material
Sphere( 1e5, Vec( 1e5+1 , 40.8, 81.6 ), Vec( .75,.25,.25), DIFF ), // Left
Sphere( 1e5, Vec( -1e5+99, 40.8, 81.6 ), Vec( .25,.25,.75), DIFF ), // Right
Sphere( 1e5, Vec( 50, 40.8, 1e5 ), Vec( .75,.75,.75), DIFF ), // Back
Sphere( 1e5, Vec( 50, 40.8, -1e5+170 ), Vec( 0, 0, 0 ), DIFF ), // Front
Sphere( 1e5, Vec( 50, 1e5, 81.6 ), Vec( .75,.75,.75), DIFF ), // Bottomm
Sphere( 1e5, Vec( 50, -1e5+81.6, 81.6 ), Vec( .75,.75,.75), DIFF ), // Top
//
Sphere( 16.5, Vec( 27, 16.5, 47), Vec( 1, 1, 1)*.999, SPEC ), // Mirror
Sphere( 16.5, Vec( 73, 16.5, 88), Vec( 1, 1, 1)*.999, REFR ), // Glass
Sphere( 8.5, Vec( 50, 8.5, 60), Vec( 1, 1, 1)*.999, DIFF ) // Middle
};
// トーンマッピングとガンマ補正
int toInt( double x )
{
return int( pow( 1 - exp( - x ), 1 / 2.2 ) * 255 + .5);
}
// 最初のヒット点までのレイのパラメータ t を求める関数
inline bool intersect( const Ray& r, double& t, int& id )
{
const s32 n = sizeof( sph ) / sizeof( Sphere );
double d,
inf = 1e20;
t = inf;
for( int i = 0; i < n; i++ )
{
d = sph[ i ].intersect( r );
if( d < t )
{
t = d;
id = i;
}
}
return ( t < inf );
}
// ポイントライトから, QMC を利用して フォトンレイを生成
void genp( Ray* photon_ray, Vec* flux, const s32 hal_i )
{
( * flux ) = Vec( 2500, 2500, 2500 ) * ( 4.0 * PI );
// QMC
double p = 2. * PI * hal( 0, hal_i ),
t = 2. * acos( sqrt( 1. - hal( 1 , hal_i ) ) );
double st = sin( t );
// フォトンレイ
photon_ray->d = Vec( cos( p ) * st, cos( t ), sin( p ) * st );
photon_ray->o = Vec( 50, 60, 85 );
}
void trace( const Ray& ray,
int depth,
const bool is_camera_ray,
const Vec& flux,
// BRDF とピクセルのフィルタリングの値
const Vec& scale_factor,
int hal_i )
{
double t;
int id;
// デプス
depth++;
// レイがシーンにヒットしない場合
if ( ! intersect( ray, t, id ) ||
( depth >= 20 ) )
{
return;
}
const s32 depth_3 = 3 * depth;
// ヒット点の情報
const Sphere& obj = sph[ id ];
const Vec x = ray.o + ray.d * t,
n = ( x - obj.p ).norm(),
material_color = obj.c;
const Vec nl = ( n.dot( ray.d ) < 0 ) ? n : ( n * - 1 );
// RGB の反射率の最大値
double p =
( material_color.x > material_color.y &&
material_color.x > material_color.z ) ?
material_color.x :
( ( material_color.y > material_color.z ) ? material_color.y : material_color.z );
// ランバート面の場合
if ( obj.refl == DIFF )
{
// 乱数の生成 (QMC を利用)
const double r1 = 2. * PI * hal( depth_3 - 1 , hal_i ),
r2 = hal( depth_3 + 0 , hal_i );
const double r2s = sqrt( r2 );
// ディフューズ用のフォトンレイの作成
Vec w = nl,
u = ( ( fabs( w.x ) > .1 ? Vec( 0, 1 ) : Vec( 1 ) ) % w ).norm();
Vec v = w % u,
d = ( u * cos( r1 )* r2s + v * sin( r1 ) * r2s + w * sqrt( 1 - r2 ) ).norm();
// カメラレイの場合
if ( is_camera_ray )
{
// ヒット点の情報の構造体を新規に確保
HitPoint* hp = new HitPoint;
// ヒット点の
hp->scale_factor = material_color.mul( scale_factor );
// ヒット点の位置座標
hp->pos = x;
// ヒット点の法線
hp->nrm = n;
// ピクセルのインデックス(グローバル変数)を設定
hp->pixel_index = pixel_index;
// ヒット点の情報の構造体を, リストの先頭に追加
hitpoints_head = CreateListElement( hp, hitpoints_head );
}
// フォトンレイの場合
else
{
// 今のヒット点のグリッドのインデックスを求める
// find neighboring measurement points and accumulate flux via progressive density estimation
Vec hh = ( x - hit_point_aabb.min ) * hash_s;
int ix = abs( int( hh.x ) ),
iy = abs( int( hh.y ) ),
iz = abs( int( hh.z ) );
// 本当は並列化した方が速いけど, 今回はしない
// strictly speaking, we should use #pragma omp critical here.
// it usually works without an artifact due to the fact that photons are
// rarely accumulated to the same measurement points at the same time (especially with QMC).
// it is also significantly faster.
{
ListElement* hp = hitpoints_hash_grid[ hash( ix, iy, iz ) ];
// ヒット点に影響を与えるフォトンをループする
while ( hp != NULL )
{
HitPoint* hitpoint = hp->hp;
hp = hp->next;
Vec v = hitpoint->pos - x;
if ( // 法線間の角度が 90° 未満であるか ?
// check normals to be closer than 90 degree (avoids some edge brightning)
( hitpoint->nrm.dot( n ) > 1e-3 ) &&
// 影響範囲内にあるか ?
( v.dot(v) <= hitpoint->radius2 ) )
{
// ヒット点の密度推定用の半径の更新
// unlike N in the paper, hitpoint->n stores "N / ALPHA" to make it an integer value
const double g = ( hitpoint->photon_n * ALPHA + ALPHA ) /
( hitpoint->photon_n * ALPHA + 1.0 );
hitpoint->radius2 = hitpoint->radius2 * g;
// フォトンの数の更新
hitpoint->photon_n++;
// flux の更新
hitpoint->flux = (
hitpoint->flux +
hitpoint->scale_factor.mul( flux ) * ( 1. / PI )
) * g;
}
}
}
// ロシアンルーレット
if ( hal( depth_3 + 1, hal_i ) < p )
{
// モンテカルロ法
const Vec next_flux = material_color.mul( flux ) * ( 1. / p );
trace( Ray( x, d ),
depth,
is_camera_ray,
next_flux,
scale_factor,
hal_i );
}
}
}
// 完全なスペキュラ反射の場合
else if ( obj.refl == SPEC )
{
const Vec next_flux = material_color.mul( flux );
const Vec next_scaling_factor = material_color.mul( scale_factor );
trace( Ray( x, ray.d - n * 2.0 * n.dot( ray.d ) ),
depth,
is_camera_ray,
next_flux,
next_scaling_factor,
hal_i );
}
// 完全に滑らかなガラスの場合
else
{
const Ray reflect_ray( x, ray.d - n * 2.0 * n.dot( ray.d ) );
const bool into = ( n.dot( nl ) > 0.0 );
double nc = 1.0,
nt = 1.5,
nnt = into ? ( nc / nt ) : ( nt / nc ),
ddn = ray.d.dot( nl ),
cos2t;
// 全反射の場合
if ( ( cos2t = 1 - nnt * nnt * ( 1 - ddn * ddn ) ) < 0 )
{
return trace( reflect_ray, depth, is_camera_ray,
flux, scale_factor, hal_i );
}
Vec td = ( ray.d * nnt - n * ( ( into ? 1 : -1 ) *
( ddn * nnt + sqrt( cos2t ) ) ) ).norm();
double a = nt - nc,
b = nt + nc,
R0 = a * a / ( b * b ),
c = 1 - ( into ? - ddn : td.dot( n ) );
double Re = R0 + ( 1 - R0 ) * c * c * c * c * c,
P = Re;
const Ray refract_ray( x, td );
const Vec fa = material_color.mul( scale_factor );
// カメラレイの場合
if ( is_camera_ray )
{
// 反射レイ
const Vec reflect_scale_factor = fa * Re;
trace( reflect_ray, depth, is_camera_ray, flux,
reflect_scale_factor, hal_i );
// 屈折レイ
const Vec refract_scale_factor = fa * ( 1.0 - Re );
trace( refract_ray, depth, is_camera_ray, flux,
refract_scale_factor, hal_i );
}
// フォトンレイの場合
else
{
// ロシアンルーレット
( hal( depth_3 - 1, hal_i ) < P ) ?
trace( reflect_ray, depth, is_camera_ray, flux, fa, hal_i ) :
trace( refract_ray, depth, is_camera_ray, flux, fa, hal_i );
}
}
}
// エントリポイント
int main( int argc, char *argv[] )
{
const u32 default_samps_num = 1000;
// const u32 default_samps_num = 5000;
// フォトンパスの合計は サンプル数 * 1000
const int w = 1024 / 2,
h = 768 / 2,
samps = ( ( argc == 2 ) ?
MAX( atoi( argv[ 1 ] ) / 1000, 1 ) :
default_samps_num );
// カメラが向いているベクトル
Ray cam( Vec( 50, 48, 295.6 ),
Vec( 0, - 0.042612, - 1 ).norm() );
// カメラの座標系の軸
Vec cx = Vec( w*.5135 / h ),
cy = ( cx % cam.d ).norm()*.5135;
// レンダリング結果を保存する配列
Vec* c = new Vec[ w * h ],
vw;
// カメラからのレイトレースのパス
for ( int y = 0; y < h; y++ )
{
fprintf( stderr, "\rCameraRay pass %5.2f%%",
100.0 * y / ( h - 1 ) );
for ( int x = 0; x < w; x++ )
{
// ピクセルのインデックス (グローバル変数)
pixel_index = x + y * w;
Vec d = cx * ( ( x + 0.5 ) / w - 0.5 ) +
cy * ( - ( y + 0.5 ) / h + 0.5 ) + cam.d;
const s32 depth = 0;
const bool is_camera_ray = true;
const Vec flux = Vec( 0, 0, 0 );
const Vec scale_factor = Vec( 1, 1, 1 );
const s32 hal_i = 0;
// カメラレイごとの,ヒット点の情報を作る
trace( Ray( cam.o + d * 140, d.norm() ),
depth,
is_camera_ray,
flux,
scale_factor,
hal_i );
}
}
fprintf( stderr, "\n" );
// ハッシュグリッドの作成
{
fprintf( stderr, "\rBuilding HashGrid" );
build_hitpoints_hash_grid( w, h );
fprintf( stderr, "\n" );
}
num_photon = samps;
vw = Vec( 1, 1, 1 );
// 全てのフォトンを追跡するパス
// OpenMP で並列化
#pragma omp parallel for schedule(dynamic, 1)
// サンプル数のループ
for( s32 i = 0; i < ( s32 ) num_photon; i++ )
{
// 進捗のデバッグ表示
const double p = 100. * ( i + 1 ) / num_photon;
fprintf( stderr, "\rPhotonRay pass %5.2f%%", p );
const s32 m = 1000 * i;
Ray photon_ray;
Vec flux;
// フォトンレイのループ (1000 個)
for( int j = 0; j < 1000; j++ )
{
const s32 hal_i = m + j;
// フォトンレイの生成
genp( ( & photon_ray ),
( & flux ),
hal_i );
const s32 depth = 0;
const bool is_camera_ray = false;
trace( photon_ray,
depth,
is_camera_ray,
flux,
vw,
hal_i );
}
}
// ピクセルごとの,フォトン密度の推定
ListElement* lst = hitpoints_head;
{
fprintf( stderr, "\rRadianceEstimate pass" );
while ( lst != NULL )
{
// ループするための変数
HitPoint* hp = lst->hp;
lst = lst->next;
// ヒット点が該当するピクセルのインデックス
const s32 pixel_index = hp->pixel_index;
// 放射輝度を加算(フォトンマッピングの公式より)
c[ pixel_index ] = c[ pixel_index ] + hp->flux * ( 1.0 / ( PI * hp->radius2 * num_photon * 1000.0 ) );
}
fprintf( stderr, "\n" );
}
// 結果を PPM として保存
{
FILE* f;
fopen_s( ( &f ), "image.ppm", "w" );
fprintf( f, "P3\n%d %d\n%d\n", w, h, 255 );
for ( int i = 0; i < w * h; i++ )
{
fprintf( f, "%d %d %d ",
toInt( c[ i ].x ), toInt( c[ i ].y ), toInt( c[ i ].z ) );
}
}
}