勉強会ネタ。

Share

頭が悪すぎるので,もう一度勉強し直したい。
一人で勉強し直すのは当然やるとして,
他にも勉強し直したい人がいれば,なんかイベントとしてやった方がいいのかね?

以下勉強しなおしたい内容。

・線形代数学
・微積分・極限
・射影幾何
・四元数
・モンテカルロ法
・QEM(Quad Error Metric)
・Progressive Mesh
・交差判定
・IK

執筆活動終わったら,ちまちま勉強し直したい。

レイトレ合宿に行ってきます。

Share

今年もやってまいりました。
そうです,レイトレ合宿です。

今日・明日と2日間参加してきます。
今年は作品と呼べるものは作れなかった。
今年の失敗を忘れないようにメモとして残しておこうと思う。

そもそも,レイトレに関しての知識がゼロに等しいので,
去年まで作ってきたプログラムはバグが多く,今年は捨てて来年に向けたベースコードを作成することにした。

…でも提出3日前で下記のような状態。

何も表示されなくなって,かなり焦る。
原因はわかっているので,そのまま放置して実装を進める。

あと今年の提出2週間前のテストで,起動できない問題にハマった。たぶんDLLが問題。
OpenMP使っているとdllとか必要になるらしいので,脱OpenMP化した。

これで絵がでるはず!

…と思って,試しにレンダリングさせてみる。
何か絵は出たが,どうもおかしい。その時の描画結果が下図。

OpenMPの時はグローバル変数で参照していたのだが,流石にまずいだろうと思ってスレッドごとに別データにするようにしたら,上記のようなバグが出た。
で,上にも書いてある通りに乱数の更新をしないため,ずっと同じ値を参照していたようで変な結果になっていた。
乱数更新を計算後に入れ忘れていたので,更新するようにしたら解決。焦る焦る。

この時点で,まだメッシュは未対応。これが提出前日です。
ヤバい感じが出てきたので,そろそろメッシュを表示しようと思って対応を始める。
メッシュフォーマットを大きく変えるとめんっどちいので,昨年まで使っていたやつに修正を加えて,表示を試みる。

絵が出た!

よくみると,ちゃんと表示されている。
…が,たぶん普通の人には分からないだろうと思って,ネタとして投稿した。
ご覧のように,全く最適なしていないため計算がゲロ重くてレイが全然飛ばない。
結局前日は,BVHを8割ぐらい実装したところで,力尽きて寝た。

提出当日。
実家から自宅に帰る前にBVHを実装。いざ実行!

「…スタックオーバーフローだってさ」

またもや,バグる。
あぁもう!!なんで~

デバッグ時間が取れないまま実家を後にして,自宅に向かう。
移動時間で2時間をロス。
夜飯を買いに行く時間が無いと思ったので,帰り際にコンビニで昼食と夕食を購入し,自宅に到着。
BVHのデバッグを開始。

「うん!さっぱりわからない…」

時間をものすごく使いそうだったので,SAHの実装は諦め中間分割の実装に切り替え。こっちのほうが実装楽だし,時間が無いということで。
で,ようやくBVHが入る。

メッシュの位置を直すのもめんどいので,そのまま。
実装終えたところで,もう時間が5時間を切っている。
ここから,QBVHの実装に入る。
QBVHは一度躓いているので,サクッと実装できた。
で,実行してみると…

「なんか,普通のBVHよりも遅い」

以前実装したときは,BVHよりも確実に速かったのだが,どうも体感的に何も感じない。
で, またバグっていた。衝突判定の実装ミス。
無駄な計算をかなりやっているコードになっていた。焦ってコード書くとろくなことにならないね。皆さん落ち着いて書きましょう。
直したところで,実行。
時間計測するコードを書く時間も惜しいので,レンダリング画像を目視して確認。とりあえず,速くなっているっポイ。
続いて,OBVHの実装に移る。
こちらも,すんなり終わるかと思いきや…。

「なぜかビルドエラーが出る」

なんか,_mm256_cvtss_f32() が無いって怒られる。なぜだ?clangのドキュメントには書いてあるのに。
残り時間も3時間ぐらいになってきたので,ネットの情報を元に自前で関数追加してビルドエラーに対処。
自宅のメインPCがAVXサポートしていないので,サブに使っているノートPCで動作チェック。
とりあえず,絵は出るっポイ。

ここで,残り時間が2時間。
流石に,IBLをMISを入れるのは無理と判断。
最終レンダリング用のアセットをここから探し始める。一応,用意していたアセットはあるんだけども,あんまり見栄えしないのでボツにした。
で,もう時間もなく探している時間もないので適当に基本形状を配置することに。

プログラムが起動できない問題が解決できていないので,ようやくチェックすることに。

Visual Studioでexe起動時にログ出力されるdllを取り合えず片っ端から,batでコピーして実行ファイルに同梱してみたのだが,上の結果。
とりあえず,エラーを再現するために,Virtual BoxにWindows 7をインストール。
VirtualBoxで動作させてみたら,案の定同じエラーが出た。
で,結局DLLが問題っぽいので,DLLを使わないようにVCのビルド設定で「マルチスレッドDLL」になっている設定を「マルチスレッド」に変更して,軌道を試みる。

「お、出た」

これで,めでたく起動した!
この時点で,残り時間が既に1時間を切っている。

流石にそろそろzip用意しないとまずいので,提出用スライドの作成に入る。
前日に力尽きて寝る前に,色々とページ作っていたけど時間が無いし,今年は負け戦なのでアピールするべきことも無いのでバッサリ省略。
スライド作成中にレンダリングをさせておいて,最新のレンダリング画像を添付。

exeの起動を確認して,急いでzipに固めて提出。
残り時間20分。

レンダリング画像を見ると,やっぱり背景がどうも気に入らないので
HDR画像を差し替えて,シーンファイルも更新。
起動して,1枚目の画像が出力されるところまで確認して,大慌てでzip化して再び提出。

残り4分だった。

来年はこういうことが無いように,ちゃんと準備したい。

上記で書いた以外にも,画像サイズを変えると描画結果がおかしくなるバグとか色々とあって,
そのあたり直したので,来年はもう少しまともなレンダリングで出来るはず

…であってほしいなぁ。

Mathjaxで式番号をフォーマットするやりかた

Share

かなり強引だけど,式番号をフォーマットするやり方。

<script type=”text/x-mathjax-config”>
MathJax.Hub.Config({
   TeX: {
       equationNumbers: {
           autoNumber: “AMS”,
           formatNumber: function (n) {return ‘8.’+n}
       }
   }
});
MathJax.Hub.Queue([“resetEquationNumbers”,MathJax.InputJax.TeX]);
</script>

これで式(8.1)みたいな感じにできる。

参考プログラム

Share

簡単なツールを作るのに,ImGuiがかなり重宝してきたので,
これを機にImGuiを使ったライブラリをちょっと作ろうかと思い始めました。

そんなわけで,今すぐ動けるわけではないので,後で実装するために参考プログラムを一覧を書いておこうと思います。

ImGuizmo

ImGuiを使ったマニュピュレーター(ギズモ)のライブラリ。
https://github.com/CedricGuillemet/ImGuizmo

カラーピッカー

https://github.com/ocornut/imgui/issues/346

ノードベースエディタ

https://github.com/ocornut/imgui/issues/306

ドッキングパネル

https://github.com/ocornut/imgui/issues/351
https://github.com/nem0/lumixengine

タイムライン

https://github.com/nem0/LumixEngine/blob/timeline_gui/external/imgui/imgui_user.inl

プロファイラ

https://github.com/mikesart/gpuvis

タブコントロール

https://github.com/ocornut/imgui/issues/261
https://github.com/ocornut/imgui/issues/1083

シリアライズ・デシリアライズ

http://qiita.com/Ushio/items/827cf026dcf74328efb7

Undo/Redo

まともに実装しようとするとコマンドシステムを自前で構築するしかなさそう。

ImGuiのLuaバインディング

https://github.com/patrickriordan/imgui_lua_bindings

ドラッグ&ドロップ

https://github.com/ocornut/imgui/issues/143

GGXの語源って?

Share

“Ground Glass” and solving for the “X” in the BSDF. を略してGGXらしい…。

GGXという単語が初出であるWalterの論文に曇りガラスが出ているので,なんとなく納得。

レイトレ再入門

Share

この記事はレイトレ合宿5!?のアドベントカレンダーの4週目の記事です。

こんにちわ,Pocolです。
皆さんレイを飛ばしていますか?

さて,昨年のレイトレ合宿ではビリという結果を残してしまったので,もう一度一から出直そうと思いました。…ということで,レイトレ再入門と題して,勉強し直してみようと思います。

レンダーターゲットの用意

まずは,描画先が必要になるのでレンダーターゲットを用意します。
HDRを表現できるように倍精度浮動小数型のVector3クラスを使いました。

    // レンダーターゲット生成.
    std::vector<Vector3> image;
    image.resize(width * height);

レンダーターゲットを作成したら,忘れずにクリアしておきます。

    // レンダーターゲットをクリア.
    for (size_t i = 0; i < image.size(); ++i)
    { image[i] = g_back_ground; }

今までは,自分でビットマップクラスを用意していたのですが,さすがにダルくなってきたので,
stb_image_writeを使いました。
ビットマップ出力は次のような感じです。

//-------------------------------------------------------------------------------------------------
//      BMPファイルに保存します.
//-------------------------------------------------------------------------------------------------
void save_to_bmp(const char* filename, int width, int height, const double* pixels)
{
    std::vector<uint8_t> images;
    images.resize(width * height * 3);

    const double inv_gamma = 1.0 / 2.2;

    for(auto i=0; i<width * height * 3; i+=3)
    {
        auto r = pow(pixels[i + 0], inv_gamma);
        auto g = pow(pixels[i + 1], inv_gamma);
        auto b = pow(pixels[i + 2], inv_gamma);

        r = saturate(r);
        g = saturate(g);
        b = saturate(b);

        images[i + 0] = static_cast<uint8_t>( r * 255.0 + 0.5 );
        images[i + 1] = static_cast<uint8_t>( g * 255.0 + 0.5 );
        images[i + 2] = static_cast<uint8_t>( b * 255.0 + 0.5 );
    }

    stbi_write_bmp(filename, width, height, 3, images.data());
}

シーンを用意する

次に簡単なシーンを用意します。球が最もやりやすいので,eduptと同じように大きな球を用意してコーネルボックスを表現します。
データは次の通りです。

const Sphere  g_spheres[] = {
    Sphere(1e5,     Vector3( 1e5 + 1.0,    40.8,          81.6), Vector3(0.25,  0.75,  0.25)),
    Sphere(1e5,     Vector3(-1e5 + 99.0,   40.8,          81.6), Vector3(0.25,  0.25,  0.75)),
    Sphere(1e5,     Vector3(50.0,          40.8,           1e5), Vector3(0.75,  0.75,  0.75)),
    Sphere(1e5,     Vector3(50.0,          40.8,  -1e5 + 170.0), Vector3()                  ),
    Sphere(1e5,     Vector3(50.0,           1e5,          81.6), Vector3(0.75,  0.75,  0.75)),
    Sphere(1e5,     Vector3(50.0,   -1e5 + 81.6,          81.6), Vector3(0.75,  0.75,  0.75)),
    Sphere(16.5,    Vector3(27.0,          16.5,          47.0), Vector3(0.75,  0.25,  0.25)),
    Sphere(16.5,    Vector3(73.0,          16.5,          78.0), Vector3(0.99,  0.99,  0.99))
};

上記で定義している球は次のようにコーディングしています。

///////////////////////////////////////////////////////////////////////////////////////////////////
// Sphere sturcture
///////////////////////////////////////////////////////////////////////////////////////////////////
struct Sphere
{
    double          radius;     //!< 半径です.
    Vector3         pos;        //!< 位置座標です.
    Vector3         color;      //!< 色です.

    Sphere
    (
        double          r,
        const Vector3&  p,
        const Vector3&  c
    )
    : radius    (r)
    , pos       (p)
    , color     (c)
    { /* DO_NOTHING*/ }

    inline double intersect(const Ray& ray) const
    {
        auto p = pos - ray.pos;
        auto b = dot(p, ray.dir);
        auto det = b * b - dot(p, p) + radius * radius;
        if (det >= 0.0)
        {
            auto sqrt_det = sqrt(det);
            auto t1 = b - sqrt_det;
            auto t2 = b + sqrt_det;
            if (t1 > D_HIT_MIN)
            { return t1; }
            else if (t2 > D_HIT_MIN)
            { return t2; }
        }
 
        return D_HIT_MAX;
    }
};

交差判定はintersect()メソッドで行い,実装は2次方程式の判別解\(\frac{D}{4}\)を用いて,当たったかどうかを計算しています。
ここまででシーンデータが準備出来ました。

レイを飛ばす

いよいよレイを飛ばします。レイはカメラから飛ばすので,最初にカメラクラスを実装しておきます。

///////////////////////////////////////////////////////////////////////////////////////////////////
// Camera class
///////////////////////////////////////////////////////////////////////////////////////////////////
class Camera
{
public:
    Camera
    (
        const Vector3&  position,
        const Vector3&  dir,
        const Vector3&  upward,
        double          fov,
        double          aspect,
        double          znear
    )
    {
        pos         = position;
        axis_x      = normalize(cross(dir, upward)) * fov * aspect;
        axis_y      = normalize(cross(dir, axis_x)) * fov;
        axis_z      = dir;
        near_clip   = znear;
    }

    inline Ray emit(double x, double y) const
    {
        auto d = axis_x * x + axis_y * y + axis_z;
        auto p = pos + d * near_clip;
        return Ray(p, normalize(d));
    }

private:
    Vector3 pos;        //!< 位置座標です.
    Vector3 axis_x;     //!< 基底ベクトル(X軸)
    Vector3 axis_y;     //!< 基底ベクトル(Y軸)
    Vector3 axis_z;     //!< 基底ベクトル(Z軸).
    double  near_clip;  //!< ニア平面までの距離.
};

やっていることはnear_clipをスクリーンの位置と見立てて,posを中心としたビュー空間を構成する基底ベクトルを求めます。
正規直交基底ベクトルを用いるのが普通なのですが,レイを飛ばすたびに同じ計算を行うと処理効率が悪いので,毎回同じ計算をする箇所は最初に1度計算してしまって,値を保存して使いまわすことにしています(axis_xとaxis_yのところが該当箇所です)。
ビュー空間を構成するベクトルが定まったら,emit()メソッドを使うことでレイを発射することができます。emit()メソッドでやっていることは,視点位置から,視線ベクトル方向に一旦移動し,xとyで指定されるスクリーン位置までaxis_xとaxis_yを使って移動します。このベクトルがレイを発射する方向となるので,正規化することでレイの方向ベクトルが定まります。レイの位置座標はカメラ位置から,方向ベクトルにnear_clip分だけ進んだところがスクリーンのヒット位置になります。この点を起点としてレイを発射します。

次に,レイと物体との交差判定を行います。すでにSphereクラスに交差判定用の処理があるので,これをシーンに配置されている球の数分for分で回して,一番交差距離が短いものを衝突物体として採用します。

//-------------------------------------------------------------------------------------------------
//      シーンとの交差判定を行います.
//-------------------------------------------------------------------------------------------------
inline bool intersect_scene(const Ray& ray, double* t, int* id)
{
    auto n = static_cast<int>(sizeof(g_spheres) / sizeof(g_spheres[0]));

    *t  = D_MAX;
    *id = -1;

    for (auto i = 0; i < n; ++i)
    {
        auto d = g_spheres[i].intersect(ray);
        if (d > D_HIT_MIN && d < *t)
        {
            *t  = d;
            *id = i;
        }
    }

    return (*t < D_HIT_MAX);
}

交差判定をして,ヒットしたら,そのピクセルに色を塗ります。

//-------------------------------------------------------------------------------------------------
//      交差物体の色を求めます.
//-------------------------------------------------------------------------------------------------
Vector3 shade(const Ray& ray)
{
    double t;
    int   id;

    // シーンとの交差判定.
    if (!intersect_scene(ray, &t, &id))
    { return g_back_ground; }

    // 交差物体の色を返却.
    return g_spheres[id].color;
}

これで,ヒット判定も実装できたので,レイを飛ばしてみます。
実装は次のようになります。


//-------------------------------------------------------------------------------------------------
//      メインエントリーポイントです.
//-------------------------------------------------------------------------------------------------
int main(int argc, char** argv)
{
    // レンダーターゲットのサイズ.
    int width  = 640;
    int height = 480;

    // カメラ用意.
    Camera camera(
        Vector3(50.0, 52.0, 295.6),                 // カメラ位置.
        normalize(Vector3(0.0, -0.042612, -1.0)),   // 視線ベクトル.
        Vector3(0.0, 1.0, 0.0),                     // 注視点.
        0.5135,                                     // 垂直画角(rad)
        double(width) / double(height),             // アスペクト比.
        130.0                                       // スクリーンまでの距離.
    );

    // レンダーターゲット生成.
    std::vector<Vector3> image;
    image.resize(width * height);

    // レンダーターゲットをクリア.
    for (size_t i = 0; i < image.size(); ++i)
    { image[i] = g_back_ground; }

    for (auto y = 0; y < height; ++y)
    {
        for (auto x = 0; x < width; ++x)
        {
            auto idx = y * width + x;
            auto fx = double(x) / double(width)  - 0.5;
            auto fy = double(y) / double(height) - 0.5;

            // Let's レイトレ!
            image[idx] += shade(camera.emit(fx, fy));
        }
    }

    // レンダーターゲットの内容をファイルに保存.
    save_to_bmp("image.bmp", width, height, &image.data()->x);

    // レンダーターゲットクリア.
    image.clear();

    return 0;
}

実装出来たら,実行してみましょう。
キチンとヒットしていれば下記のように色がつくはずです。

ここまでの,実装プログラムをGithubにアップロードしておきました。
https://github.com/ProjectAsura/sample_hit

これでレイが飛ばせるようになりました。

古典的レイトレーシング

続いて,古典的レイトレーシングを実装してみます。マテリアルの概念を導入し,マテリアルに沿って反射レイの飛ばし方を変えてみます。
まず,反射タイプを追加します。

///////////////////////////////////////////////////////////////////////////////////////////////////
// ReflectionType enum
///////////////////////////////////////////////////////////////////////////////////////////////////
enum ReflectionType
{
    Diffuse             = 0,    //!< 完全拡散反射.
    PerfectSpecular     = 1,    //!< 完全鏡面反射.
    Refraction          = 2,    //!< 屈折.
};

反射タイプを追加したら,shade()メソッドをradiance()メソッドにリネームして,下記のような実装を行います。

//-------------------------------------------------------------------------------------------------
//      放射輝度を求めます.
//-------------------------------------------------------------------------------------------------
Vector3 radiance(const Ray& ray, int depth)
{
    double t;
    int   id;

    // シーンとの交差判定.
    if (!intersect_scene(ray, &t, &id))
    { return g_back_ground; }

    // 交差物体.
    const auto& obj = g_spheres[id];

    // 交差位置.
    const auto hit_pos = ray.pos + ray.dir * t;

    // 法線ベクトル.
    const auto normal  = normalize(hit_pos - obj.pos);

    // 物体からのレイの入出を考慮した法線ベクトル.
    const auto orienting_normal = (dot(normal, ray.dir) < 0.0) ? normal : -normal;

    // 打ち切り深度に達したら終わり.
    if(depth > g_max_depth)
    { return g_back_ground; }

    switch (obj.type)
    {
    case ReflectionType::Diffuse:
        {
            double t_;
            int    id_;

            // ライトベクトル.
            auto light_dir  = g_light_pos - hit_pos;

            // ライトまでの距離.
            auto light_dist = length(light_dir);

            // ライトベクトルを正規化.
            light_dir /= light_dist;

            // ライトとの間に遮蔽物がないことを確認.
            intersect_scene(Ray(hit_pos, light_dir), &t_, &id_);

            // 遮蔽物がない場合.
            if (t_ >= light_dist)
            {
                auto diffuse = obj.color * max(dot(orienting_normal, light_dir), 0.0) / (light_dist * light_dist);
                return g_light_color * diffuse;
            }
            else
            {
                // 遮蔽物がある.
                return g_shadow_color;
            }
        }
        break;

    case ReflectionType::PerfectSpecular:
        {
            // 反射させる.
            return obj.color * radiance(Ray(hit_pos, reflect(ray.dir, normal)), depth + 1);
        }
        break;

    case ReflectionType::Refraction:
        {
            // 反射レイ
            auto reflect_ray = Ray(hit_pos, reflect(ray.dir, normal));

            // 内部侵入するか?
            auto into = dot(normal, orienting_normal) > 0.0;

            // 空気の屈折率
            const auto nc = 1.0;

            // 物体の屈折率
            const auto nt = 1.5;

            // Snellの法則.
            const auto nnt = (into) ? (nc / nt) : (nt / nc);
            const auto vn  = dot(ray.dir, orienting_normal);
            const auto cos2t = 1.0 - nnt * nnt * (1.0 - vn * vn);

            // 全反射かどうかチェック.
            if (cos2t < 0.0)
            { return obj.color * radiance(reflect_ray, depth + 1); }

            // 屈折ベクトル.
            auto refract = normalize(ray.dir * nnt - normal * ((into) ? 1.0 : -1.0) * (vn * nnt + sqrt(cos2t)) );

            // Schlickによる Fresnel の反射係数の近似.
            const auto a  = nt - nc;
            const auto b  = nt + nc;
            const auto R0 = (a * a) / (b * b);

            const auto c  = 1.0 - ((into) ? -vn : dot(refract, normal));
            const auto Re = R0 + (1.0 - R0) * pow(c, 5.0);

            const auto nnt2 = pow((into) ? (nc / nt) : (nt /nc), 2.0);
            const auto Tr = (1.0 - Re) * nnt2;
            const auto p  = 0.25 + 0.5 * Re;

            // 屈性レイ
            Ray refract_ray(hit_pos, refract);

            const auto reflect_result = radiance(reflect_ray, depth + 1) * Re;
            const auto refract_result = radiance(refract_ray, depth + 1) * Tr;

            return obj.color * (reflect_result + refract_result);
        }
        break;
    }

    // どれにもヒットしなかった.
    return g_back_ground;
}

リアルタイムレンダリングのようにLambertのBRDFを計算しているのと,反射・屈折の処理が追加されています。
反射・屈折処理についてはeduptとほぼ同じなので説明はeduptのスライドを参照してください。

名前をradiance()に変えたのでmain()関数側も忘れずに修正しておきます。

    for (auto y = 0; y < height; ++y)
    {
        for (auto x = 0; x < width; ++x)
        {
            auto idx = y * width + x;
            auto fx = double(x) / double(width)  - 0.5;
            auto fy = double(y) / double(height) - 0.5;

            // Let's レイトレ!
            image[idx] += radiance(camera.emit(fx, fy), 0);
        }
    }

修正したら実行してみます。
きちんと実装されていれば次のようになるはずです。

この程度の処理なら,現代的なPCを使っていれば1秒かからないうちに終わります。
ここまでのサンプルプログラムをGithubにアップロードしておきました。
https://github.com/ProjectAsura/sample_rt

Path Tracing

 さて,続いてパストレです。radiance()メソッドでは,ここまで再帰を使っていましたが,あんまり呼び出すとスタックオーバーフローになる可能性があるので,再帰を使わずにループ文に書き直します。実装は次のような感じです。

//-------------------------------------------------------------------------------------------------
//      放射輝度を求めます.
//-------------------------------------------------------------------------------------------------
Vector3 radiance(const Ray& input_ray, int depth, Random* random)
{
    Vector3 L(0, 0, 0);
    Vector3 W(1, 1, 1);
    Ray ray(input_ray.pos, input_ray.dir);

    while(true)
    {
        double t;
        int   id;

        // シーンとの交差判定.
        if (!intersect_scene(ray, &t, &id))
        { break; }

        // 交差物体.
        const auto& obj = g_spheres[id];

        // 交差位置.
        const auto hit_pos = ray.pos + ray.dir * t;

        // 法線ベクトル.
        const auto normal  = normalize(hit_pos - obj.pos);

        // 物体からのレイの入出を考慮した法線ベクトル.
        const auto orienting_normal = (dot(normal, ray.dir) < 0.0) ? normal : -normal;

        auto p = max(obj.color.x, max(obj.color.y, obj.color.z));

        L += W * obj.emission;

        // 打ち切り深度に達したら終わり.
        if(depth > g_max_depth)
        {
            if (random->get_as_double() >= p)
            { break; }
        }
        else
        {
            p = 1.0;
        }

        switch (obj.type)
        {
        case ReflectionType::Diffuse:
            {
                // 基底ベクトル.
                Vector3 u, v, w;

                w = orienting_normal;
                if (abs(w.x) > 0.1)
                { u = normalize(cross(Vector3(0, 1, 0), w)); }
                else
                { u = normalize(cross(Vector3(1, 0, 0), w)); }
                v = cross(w, u);

                const auto r1 = D_2PI * random->get_as_double();
                const auto r2 = random->get_as_double();
                const auto r2s = sqrt(r2);

                auto dir = normalize(u * cos(r1) * r2s + v * sin(r1) * r2s + w * sqrt(1.0 - r2));

                ray = Ray(hit_pos, dir);
                W *= (obj.color / p);
            }
            break;

        case ReflectionType::PerfectSpecular:
            {
                ray = Ray(hit_pos, reflect(ray.dir, normal));
                W *= (obj.color / p);
            }
            break;

        case ReflectionType::Refraction:
            {
                Ray reflect_ray = Ray(hit_pos, reflect(ray.dir, normal));
                auto into = dot(normal, orienting_normal) > 0.0;

                const auto nc = 1.0;
                const auto nt = 1.5;
                const auto nnt = (into) ? (nc / nt) : (nt / nc);
                const auto ddn = dot(ray.dir, orienting_normal);
                const auto cos2t = 1.0 - nnt * nnt * (1.0 - ddn * ddn);

                if (cos2t < 0.0)
                {
                    ray = reflect_ray;
                    W *= (obj.color / p);
                    break;
                }

                auto dir = normalize(ray.dir * nnt - normal * ((into) ? 1.0 : -1.0) * (ddn * nnt + sqrt(cos2t)));

                const auto a = nt - nc;
                const auto b = nt + nc;
                const auto R0 = (a * a) / (b * b);
                const auto c = 1.0 - ((into) ? -ddn : dot(dir, normal));
                const auto Re = R0 + (1.0 - R0) * pow(c, 5.0);
                const auto Tr = 1.0 - Re;
                const auto prob = 0.25 + 0.5 * Re;

                if (random->get_as_double() < prob)
                {
                    ray = reflect_ray;
                    W *= (obj.color * Re / prob) / p; 
                }
                else
                {
                    ray = Ray(hit_pos, dir);
                    W *= (obj.color * Tr / (1.0 - prob)) / p;
                }
            }
            break;
        }

        depth++;
    }

    return L;
}

次に,複数サンプルとれるようにmain()関数内の処理を次のように変更します。

//-------------------------------------------------------------------------------------------------
//      メインエントリーポイントです.
//-------------------------------------------------------------------------------------------------
int main(int argc, char** argv)
{
    // レンダーターゲットのサイズ.
    int width   = 640;
    int height  = 480;
    int samples = 512;

    // カメラ用意.
    Camera camera(
        Vector3(50.0, 52.0, 295.6),
        normalize(Vector3(0.0, -0.042612, -1.0)),
        Vector3(0.0, 1.0, 0.0),
        0.5135,
        double(width) / double(height),
        130.0
    );

    // レンダーターゲット生成.
    std::vector<Vector3> image;
    image.resize(width * height);

    Random random(123456);

    // レンダーターゲットをクリア.
    for (size_t i = 0; i < image.size(); ++i)
    { image[i] = g_back_ground; }

    for(auto s = 0; s < samples; ++s)
    {
        printf_s("%.2lf%% complete\r", (double(s)/double(samples) * 100.0));

        for (auto y = 0; y < height; ++y)
        {
            for (auto x = 0; x < width; ++x)
            {   
                auto idx = y * width + x;

                auto fx = double(x) / double(width)  - 0.5;
                auto fy = double(y) / double(height) - 0.5;

                // Let's レイトレ!
                image[idx] += radiance(camera.emit(fx, fy), 0, &random) / samples;
            }
        }
    }

    // レンダーターゲットの内容をファイルに保存.
    save_to_bmp("image.bmp", width, height, &image.data()->x);

    // レンダーターゲットクリア.
    image.clear();

    return 0;
}

見ると分かるように,乱数が導入されているのと,サンプル数分ループする処理が追加されています。
このプログラムを実行すると次のような結果が得られます。

サンプル数が少ないため,ノイジーですが,ガラス玉付近に集光現象が見られたり,柔らかい影が表現されていたりなど,ライティング結果の向上が見て取れます。
ここまでのプログラムをGithubにアップロードしておきました。
https://github.com/ProjectAsura/sample_pt

これでパストレーシングも実装しました。

Next Event Estimation

次に速度向上のために直接光ライティングを導入してみます。
Next Event EstimationについてはSchokerさんが詳しく説明されています。Schokerさんのページを参照しましょう。
http://rayspace.xyz/CG/contents/path_tracing_implementation/

物体と衝突したら,光源方向にシャドウレイを飛ばします。シャドウレイを飛ばした結果,ライト以外の遮蔽物と交差していなければライトの寄与が取れます。遮蔽されている場合は何もしません。今回のサンプルではライトが一つなので,1つのライトを選択して,乱数により適当なライト表面上の位置を決定します。この点に向かってシャドウレイを飛ばして,交差判定を行います。遮蔽が無い場合は,BRDFと,G項とライトの確率密度関数から寄与を計算します。G項については下記のスライドに記載があります。

このNext Event Estimationの処理をDiffuse計算の先頭に追加します。

        case ReflectionType::Diffuse:
            {
                #if 1
                // Next Event Estimation
                {
                    const auto& light = g_spheres[g_lightId];

                    const auto r1 = D_2PI * random->get_as_double();
                    const auto r2 = 1.0 - 2.0 * random->get_as_double();
                    const auto light_pos = light.pos + (light.radius + D_HIT_MIN) * Vector3(sqrt(1.0 - r2 * r2) * cos(r1), sqrt(1.0 - r2 * r2) * sin(r1), r2);

                    // ライトベクトル.
                    auto light_dir   = light_pos - hit_pos;

                    // ライトへの距離の2乗
                    auto light_dist2 = dot(light_dir, light_dir);

                    // 正規化.
                    light_dir = normalize(light_dir);

                    // ライトの法線ベクトル.
                    auto light_normal = normalize(light_pos - light.pos);

                    auto dot0 = dot(orienting_normal, light_dir);
                    auto dot1 = dot(light_normal, -light_dir);
                    auto rad2 = light.radius * light.radius;

                    // 寄与が取れる場合.
                    if (dot0 >= 0 && dot1 >= 0 && light_dist2 >= rad2)
                    {
                        double shadow_t;
                        int    shadow_id;
                        Ray    shadow_ray(hit_pos, light_dir);

                        // シャドウレイを発射.
                        auto hit = intersect_scene(shadow_ray, &shadow_t, &shadow_id);

                        // ライトのみと衝突した場合のみ寄与を取る.
                        if (hit && shadow_id == g_lightId)
                        {
                            auto G = dot0 * dot1 / light_dist2;
                            auto pdf = 1.0 / (4.0 * D_PI * rad2);

                            L += W * light.emission * (obj.color / D_PI) * G / pdf;
                        }
                    }
                }
                #endif

この追加を加えたプログラムを実行すると,次のような結果が得られました。

Next Event Estimationが入っていないプログラムに比べると良くなっているのが見て取れます。
ここまでのプログラムをGithubにアップロードしておきました。
https://github.com/ProjectAsura/sample_pt2

若干,プログラムが怪しいですが…
一応Next Event Estimationが実装できたということにしておきたいと思います。

超雑訳 Robust Monte Carlo Methods For Light Transport Simulation (1)

Share

こんにちわ。
Pocolです。

ようやくレンダリングの仕事にありつけたので,きちんとしたレンダリストになるためには,PBRTとVeachの論文を読んでおくのは当たり前だそうなので…,きちんとしたレンダリストになるために読んでおこうかなと思います。ページ数が英文で約390ページ程度あるのですが,一気に訳す力は無いので,コツコツと読んでいこうと思います。
原文はhttp://graphics.stanford.edu/papers/veach_thesis/にあります。

とりあえず,今回はAbstractだけ読んでみます。

 

(さらに…)

GPU Pro 8は出ないらしい。

Share

そういえば,この間Twitterでbgfxの作者の方に教えてもらったのですが,
“GPU Pro 8″という名前の本は出ないらしいです。
どうもタイトル名が変わって”GPU Zen”という名前で出るみたいです。
この本の目次と表紙カバーが編集者を務めるWolfgang Engelが開設しているhttp://gpuzen.blogspot.jp/にて公開されています。

内容としては…

  • Rendering stable indirect illumination computed from reflective shadow maps
  • Real-Time Linear-Light Shading with Linearly Transformed Cosines
  • Scalable Adaptive SSAO

あたりが気になっています。
まだ,いつ出版か決まっていないようですが,はやく出ないもんですかね…。

ちなみに2番目のやつは,作者が
https://labs.unity.com/article/linear-light-shading-linearly-transformed-cosinesにてPDFを公開してくれています。