Gregory Igehy

Dancing at hemisphere coordinate

Adding Japanese comments for smallppm

About

  • smallppm に日本語のコメントを追加してみた

Reference URL

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 ) );
        }
    }
}