Chapter 03

モンテカルロレイトレーシング

mebiusbox
mebiusbox
2021.02.23に更新
このチャプターの目次
  1. 📌 光学
  2. 📌 光の物理量
  3. 📌 モンテカルロレイトレーシング
  4. 📌 モンテカルロ法について
  5. 📌 光の散乱
  6. 📌 重点的サンプリング
    1. 散乱 pdf と重み
    2. 半球上の無作為な方向の生成
    3. 正規直交基底
    4. 光源のサンプリング
    5. 合成PDF
    6. 金属の復活
  7. 📌 この次はどうするか
    1. ロシアンルーレット
    2. パストレーシング
    3. 関与媒質
    4. フォトンマッピング
  8. 📌 参考となる資料
    1. Peter Shirley 著の「Ray Tracing」
    2. smallpt
    3. edupt
    4. memoRANDOM
    5. 「フォトンマッピングー実写に迫るコンピュータグラフィックス」
    6. 「CG Magic」
    7. Robust Monte Carlo Methods for Light Transport Simulation
    8. 「Physically Based Rendering」
    9. レイトレ合宿
  9. 📌 最後に
  10. 📌 参考文献

📌 光学

ここからはモンテカルロ法を使ってレイトレーシングを行う手法を説明します．そのためにまずは少し光学について触れておきます．これまで行ってきた光線を飛ばして光伝達をシミュレーションするような方法をレイトレーシングといいました．これは幾何光学という分野になります．なので，レイトレーシングは幾何光学で説明できる物理現象（例えば反射や屈折）を実現することが出来ます．それとは別に，光は電磁波の一種なので，波動の特性を持っています．光の波動による物理現象（例えば回折や干渉）は波動光学の分野になります．もちろん電磁波なので電磁学にも関連しています．例えば偏光がこの分野になります．光の物理現象をシミュレーションするには，こういった分野がすべて必要になってくるのですが，実際は計算が複雑になったり，よくわからん（切実）ので，基本的に幾何光学だけを考えます．フレネルの方程式のときに偏光が出てきましたが，無視して考えていたのはそういうことです．そのため，写実的なものをレンダリングするためには，まず幾何光学について知っておく必要があります．

📌 光の物理量

コンピュータグラフィックスにおいて，光の最も小さい単位は「光子」です．これを「フォトン（photon）」といいます．光のエネルギーというのはこの光子が集まったもので，放射エネルギーといいます．細かいことをいうと波長が関係してくるのですが，ここでは割愛します．

この放射エネルギーの時間的な変化を表したものを「放射束」といいます．これは放射エネルギーの仕事量のことでもあるので単位は「ワット（W）」です．電球とかで 60W とかありますよね？　あれです．

これまで計算してきた数値は光学においてどのような値（量）なのでしょうか．それは「放射束」です．光伝達というのは放射束を輸送することを指します．

余談ですが，人間の眼が知覚する光の量を測光値といい，この分野を測光学といいます．測光学における放射束を「光束」といいます．

📌 モンテカルロレイトレーシング

古典的なレイトレーシングは，カメラから光線を飛ばして鏡面反射させていました．その後，鏡面反射の方向の周りに向かって複数の光線を飛ばして計算するように改善されました．これを分散レイトレーシングといいます．これから実装していくモンテカルロレイトレーシングはこの分散レイトレーシングをさらに改善させたものです．

光源から放射された放射束は反射を繰り返して物体表面に入射します．物体の表面は拡散反射や鏡面反射などを行ってあらゆる方向に放射束を反射します．そして，物体表面のある点を考えると，半球上のあらゆるところから放射束が入ってくることになります．つまり，半球上から入ってくるすべての放射束を積分する必要があります．これをシミュレーションしようとするのは困難なので，統計的にサンプリングを行って，結果を推定します．そこでモンテカルロ法を使います．積分範囲内で「無作為」に光線を発生させてサンプリングし，平均値を求めて積分の推定値とします．このように表面のある点に入ってくる放射束をモンテカルロ法で求めて，反射や屈折などを計算する方法をモンテカルロレイトレーシングといいます．

📌 モンテカルロ法について

モンテカルロ法，またそれに関する確率と統計については「CGのための確率・統計入門」を参照してください．また，立体角といった新しい用語が出てきます．それらについては「基礎からはじめる物理ベースレンダリング」で説明していますので，分からない場合は参照してください．これからはモンテカルロ法に基づいて，各種の物理現象をシミュレーションしていきます．前回の最後でコーネルボックスをレンダリングしました．そのコーネルボックスを使ってモンテカルロレイトレーシングを実装していきます．

tuto-raytracing-cornelbox-s2000-output

📌 光の散乱

材質のところで散乱について説明しました．ここでもう一度定義すると，散乱とは物体表面に光線が当たったときに，光線の向きが変わる現象のことです．また，物体表面で吸収されずに，反射される光の割合を反射率（アルベド）といいます．ここで，光とは放射束のことです．

散乱したときの方向分布は，立体角の確率密度関数（以降，pdfと記述します）によって表され，これを s(direction) と表します．この pdf は入射角に依存していて，グレージング角に近づくほど小さくなっていきます．

ある表面上の点に入射してくる全放射束は

color = \int_{\Omega} albedo \cdot s(direction) \cdot color(direction).

ここで \Omega は半球を表します．上記のように半球上の全ての方向から入射した放射束を積分したものを放射照度（Irradiance）といいます．これをモンテカルロ法で解くと

color = \frac{1}{N} \sum_{i=1}^{N} \frac{albedo \cdot s(direction) \cdot color(direction)}{p(direction)}.

color(direction) は direction 方向から入射してくる放射束で，これを放射輝度（Radiance）といいます．p(direction) は無作為に生成した方向の pdf です．

ランバート面では s(direction) は \cos(\theta) に比例します．pdf は積分すると 1 になるので，\cos(\theta) < 1 を満たします．また，角度が大きくなる（グレージング角に近づく）と値は減少していき，90^{\circ} のときは s(direction) = 0 となります．\cos(\theta) を半球積分すると \pi になりますので，

s(direction) = \frac{\cos(\theta)}{\pi}.

s(direction) と p(direction) が同じなら

color = albedo \cdot color(direction),

となり，今までの計算と変わらないことになります．

余談ですが，物理ベースレンダリングとかでよくでてくる双方向反射率分布関数（Bidirectional reflectance distribution function: BRDF）との関係は

BRDF = \frac{albedo \cdot s(direction)}{\cos\theta},

であり，ランバート面の s(direction) = \frac{\cos\theta}{\pi} を上の式に代入すると

BRDF = \frac{albedo}{\pi}.

これはランバートの BRDF と一致します．
また，関与媒質において s(direction) は位相関数（phase function）に相当します．

📌 重点的サンプリング

Importance sampling.

散乱 pdf と重み

材質クラスに散乱pdfを返す関数 scattering_pdf を追加します．デフォルトで 0 を返します．

virtual float scattering_pdf(const Ray& r, const HitRec& hrec) const { return 0; }

ランバート材質の scattering_pdf は \cos(\theta)/\pi を返します．また，裏面に当たったときは散乱しないよう（つまり 0）にしています．

virtual float scattering_pdf(const Ray& r, const HitRec& hrec) const override {
    return std::max(dot(hrec.n, normalize(r.direction())), 0.0f) / PI;
}

また，モンテカルロ法で積分するために p(direction) が必要です．この p(direction) は確率の重みを表しているともいえます．ScatterRec に pdf_value を追加し，scatter 関数で値を計算します．

class ScatterRec {
public:
    Ray	ray;
    vec3 albedo;
    float pdf_value;
};

virtual bool scatter(const Ray& r, const HitRec& hrec, ScatterRec& srec) const override {
    vec3 target = hrec.p + hrec.n + random_in_unit_sphere();
    srec.ray = Ray(hrec.p, normalize(target - hrec.p));
    srec.albedo = m_albedo->value(hrec.u, hrec.v, hrec.p);
    srec.pdf_value = dot(hrec.n, srec.ray.direction()) / PI;
    return true;
};

color 関数で，散乱 pdf の値とアルベド，放射輝度 (color関数の返り値）の積を求めて，それに ScatterRec.pdf_value で除算します．

vec3 color(const rayt::Ray& r, const Shape* world, int depth) {
    HitRec hrec;
    if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
        vec3 emitted = hrec.mat->emitted(r, hrec);
        ScatterRec srec;
        if (depth < MAX_DEPTH && hrec.mat->scatter(r, hrec, srec) && srec.pdf_value > 0) {
            float spdf_value = hrec.mat->scattering_pdf(srec.ray, hrec);
            vec3 albedo = srec.albedo * spdf_value;
            return emitted + mulPerElem(albedo, color(srec.ray, world, depth + 1)) / srec.pdf_value;
        }
        else {
            return emitted;
        }
    }
    return background(r.direction());
}

これを実行すると次のようになります．（rayt301.cpp)

tuto-raytracing-monte-spdf-output

ノイズが多くありますね．これからこのノイズを無くしていきます．

半球上の無作為な方向の生成

重点的サンプリングを行うには表面の法線を中心に半球上の無作為な方向を生成する必要があります．まずは z 軸が法線方向，\theta が法線方向からの角度を表しているとします．方向の作成には方位角 \phi の確率変数 r_1 と，極角の確率変数 r_2 が必要です．最初にどちらも一様分布に従う確率変数を求めて，それから確率分布が \cos になるような極角の確率変数 r_2 を求めます．

方位角 \phi の区間は[0,2\pi]なので一様分布に従う pdf は

p(\phi) = \frac{1}{2\pi}.

確率変数 r_1 は

r_1 = \int_{0}^{\phi}\frac{1}{2\pi} = \frac{\phi}{2\pi}.

\phi について解くと

\phi = 2\pi r_1

極角 \theta の pdf は p(direction) = f(\theta)，立体角 \sin\theta d\theta d\phi なので，確率変数 r_2 は

\begin{aligned} r_2 &= \int_{0}^{2\pi} \int_{0}^{\theta} f(\theta) \sin\theta d\theta d\phi \\[2ex] &= 2\pi \int_{0}^{\theta} f(\theta) \sin\theta d\theta \\[2ex] \end{aligned}

単位球の面積は 4\pi なので，f(\theta) は \frac{1}{4\pi} となります．

\begin{aligned} r_2 &= 2\pi \int_{0}^{\theta} \frac{1}{4\pi} \sin(t) dt \\[2ex] &= \frac{1}{2} \int_{0}^{\theta}\sin(t) dt\\[2ex] &= \frac{1}{2} [-\cos]_{0}^{\theta} \\[2ex] &= -\frac{\cos(\theta)}{2} - \left(-\frac{\cos(0)}{2}\right) \\[2ex] &= -\frac{\cos(\theta)}{2} + \frac{1}{2} \\[2ex] &= \frac{(1-\cos(\theta))}{2}. \\[2ex] \end{aligned}

\cos(\theta) について解くと

\cos(\theta) = 1-2r_2.

極座標（\phi,\theta）からデカルト座標（x,y,z）に変換するときは次の式を使います．

\begin{aligned} x &= \cos(\phi)\sin(\theta) \\[2ex] y &= \sin(\phi)\sin(\theta) \\[2ex] z &= \cos(\theta). \end{aligned}

\cos^2 + \sin^2 + 1 を使って変形し，代入すると

\begin{aligned} x &= \cos(2\pi r_1)\sqrt{1-(1-2r_2)^2} \\[2ex] y &= \sin(2\pi r_2)\sqrt{1-(1-2r_2)^2} \\[2ex] z &= 1-2r_2. \end{aligned}

平方根の項を整理すると

\begin{aligned} \sqrt{1-(1-2r_2)^2} &= \sqrt{1-(1-4r_2+4_r^2)} \\[2ex] &= \sqrt{4r_2 - 4r_2^2} \\[2ex] &= \sqrt{4(r_2 - r_2^2)} \\[2ex] &= 2\sqrt{r_2 - r_2^2} \\[2ex] &= 2\sqrt{r_2(1-r_2)}. \end{aligned}

なので，

\begin{aligned} x &= \cos(2\pi r_1)2\sqrt{r_2(1-r_2)} \\[2ex] y &= \sin(2\pi r_1)2\sqrt{r_2(1-r_2)} \\[2ex] z &= 1-2r_2. \end{aligned}

これは単位球上の無作為な方向なので，半球上に制限します．半球上で一様分布な pdf は p(direction) = \frac{1}{2\pi} なので，

\begin{aligned} r_2 &= 2\pi\int_{0}^{\theta}\frac{1}{2\pi}\sin(t)dt \\[2ex] &= \int_{0}^{\theta}\sin(t)dt \\[2ex] &= \left[-cos(t)\right]_0^{\theta} \\[2ex] &= -\cos(\theta)-(-\cos(0)) \\[2ex] &= 1-\cos(\theta). \end{aligned}

\cos(\theta) について解くと

\cos(\theta) = 1-r_2.

デカルト座標に変換すると

\begin{aligned} x &= \cos(2\pi r_1)\sqrt{r_2(2-r_2)} \\[2ex] y &= \sin(2\pi r_1)\sqrt{r_2(2-r_2)} \\[2ex] z &= 1-r_2. \end{aligned}

次に pdf が p(direction) = \frac{\cos(\theta)}{\pi} の方向を生成します．

\begin{aligned} r_2 &= 2\pi \int_{0}^{\theta} \left(\frac{\cos(t)}{\pi}\right)\sin(t) dt \\[2ex] &= 2 \int_{0}^{\theta} \cos(t) \sin(t) dt \\[2ex] \end{aligned}
t = \cos(\theta), \frac{dt}{d\theta} = -\sin(\theta), dt = -\sin\theta d\theta

で置換積分します．

\begin{aligned} r_2 &= -2 \int_{0}^{\cos\theta} t dt \\[2ex] &= \left[\frac{1}{2}t^2\right]_0^{\cos\theta} \\[2ex] &= -2\left\{\left(\frac{1}{2}\cos^2(\theta)\right)- \left(\frac{1}{2}\cos^2(0)\right)\right\} \\[2ex] &= -(\cos^2(\theta)-1) \\[2ex] &= 1-\cos^2(\theta) \end{aligned}

よって

\cos(\theta) = \sqrt{1-r_2}.

デカルト座標で表すと

\begin{aligned} z &= \cos(\theta) = \sqrt{1-r_2} \\[2ex] x &= \cos(\phi)\sin(\theta) = \cos(2\pi r_1)\sqrt{1-z^2} = \cos(2\pi r_1)\sqrt{r_2} \\[2ex] y &= \sin(\phi)\sin(\theta) = \cos(2\pi r_1)\sqrt{1-z^2} = \sin(2\pi r_1)\sqrt{r_2}. \end{aligned}

これをコードに書き起こしたものが以下になります．

// rayt.h
inline vec3 random_cosine_direction() {
    float r1 = drand48();
    float r2 = drand48();
    float z = sqrt(1.f - r2);
    float phi = PI2*r1;
    float x = cos(phi) * sqrt(r2);
    float y = sin(phi) * sqrt(r2);
    return vec3(x, y, z);
}

この関数がどのようなベクトルを作成するのか統計的に調べてみます．Python+Numpy+Matplotlib で，上記のベクトルを適当な個数作成してプロットしてみました．プログラムは次のようになっています．

import numpy as np
import matplotlib.pyplot as plt

t = np.random.random((2, 1000))
x = np.cos(t[0]*2*np.pi) * np.sqrt(t[1])
y = np.sin(t[0]*2*np.pi) * np.sqrt(t[1])
z = np.sqrt(1-t[1])
fig = plt.figure()
ax = fig.add_subplot(131, projection='3d')
ax.scatter3D(x, y, z)
ax = fig.add_subplot(132)
ax.scatter(x, z, s=2)
ax = fig.add_subplot(133)
ax.scatter(x, y, s=2)
plt.show()

fig = plt.figure()
ax = fig.add_subplot(131)
ax.hist(x)
ax = fig.add_subplot(132)
ax.hist(y)
ax = fig.add_subplot(133)
ax.hist(z)
plt.show()

結果は次の通りです．

tuto-raytracing-random_cosine_direction_scat

tuto-raytracing-random_cosine_direction_hist

一段目の一番左を見ると半球上に点が生成されているのがわかります．一段目の中央，XZ のプロットを見てみるとZ値が高くなると点の数も増えているように見えませんか？　ちょっとわかりづらいかもしれません．そんなときは２段目のヒストグラムを見てみます．一番右が Z値ですが，段階的になっていることがわかります．つまり，Z値が高くなるほど多くなっています．Zと比べてXYはどうでしょうか．１段目の一番右を見てみると円の中で一様分布しているように見えます．先ほどの関数はこのような分布をするベクトルを無作為に生成します．

正規直交基底

Orhto-normal bases.

z方向を中心とした方向を無作為に作成できるようになったので，任意の方向（例えば表面の法線）を中心にできるようにします．そのため，互いに直交する３つの単位ベクトルを決定し，基底変換を行います．基底変換はカメラのときにもやったので，ここでは詳しく説明しません．今回は法線から正規直交基底を作成します．カメラのときには vup を指定していましたが，ここでは基本的にY方向のベクトルを使います．ただし，法線と平行の場合に正しく計算できないので，その時はX方向のベクトルを使います．正規直交なので，基底ベクトルは正規化します．

// rayt.h
class ONB {
public:
    ONB() {}
    inline vec3& operator[](int i) {
        return m_axis[i];
    }
    inline const vec3& operator[](int i) const {
        return m_axis[i];
    }
    const vec3& u() const {
        return m_axis[0];
    }
    const vec3& v() const {
        return m_axis[1];
    }
    const vec3& w() const {
        return m_axis[2];
    }
    vec3 local(float a, float b, float c) const {
        return a*u() + b*v() + c*w();
    }
    vec3 local(const vec3& a) const {
        return a.getX()*u() + a.getY()*v() + a.getZ()*w();
    }
    void build_from_w(const vec3& n) {
        m_axis[2] = normalize(n);
        vec3 a;
        if (fabs(w().getX()) > 0.9) {
            a = vec3(0, 1, 0);
        }
        else {
            a = vec3(1, 0, 0);
        }
        m_axis[1] = normalize(cross(w(), a));
        m_axis[0] = cross(w(), v());
    }

private:
    vec3 m_axis[3];
};

これを使って無作為に作成した方向ベクトルを，物体表面の法線方向を中心として基底変換します．

virtual bool scatter(const Ray& r, const HitRec& hrec, ScatterRec& srec) const override {
    ONB onb; onb.build_from_w(hrec.n);
    vec3 direction = onb.local(random_cosine_direction());
    srec.ray = Ray(hrec.p, normalize(direction));
    srec.albedo = m_albedo->value(hrec.u, hrec.v, hrec.p);
    srec.pdf_value = dot(onb.w(), srec.ray.direction()) / PI;
    return true;
};

レンダリングすると次のようになります．(rayt302.cpp, rayt302.h)

tuto-raytracing-monte-cosine

光源のサンプリング

レイトレーシングでは，カメラから光線を飛ばして，最終的に光源まで辿りつかないと入射する放射束の精度が大きく下がります．光線の反射は基本的にランダムな方向なのですが，意図的に光源方向に光線を飛ばしてサンプリングすれば光源にたどり着く可能性が上がります（光線の開始点と光源との間に物体があると，光源に辿り着かないことがあります）

コーネルボックスでは光源は四角形です．この領域上をサンプリングするために適切な確率密度関数を決定する必要があります．

表面上のある点 p から光源上の点 q の関係は次の図のようになっています．

tuto-raytracing-solidangle-relationship

点 q の面積は dA です．この領域の確率はどれくらいでしょうか．四角形の面積は A なので，一様分布ならば \frac{dA}{A} です．d\omega は立体角で球の半径を r，立体角当たりの面積を da とすると d\omega = \frac{da}{r^2} という関係があります．光源上の微小面積 dA の投影面積は dA\cos(\theta)なので da = dA\cos(\theta)，なので以下の関係式が成り立ちます．

d\omega = \frac{dA\cos(\theta)}{distance(\vec{p},\vec{q})^2}.

球上の点をサンプリングする確率は p(direction) なので，確率変数は p(direction)\cdot d\omega です．今，球上の d\omega と光源上の dA の確率を一致させたいので，

p(direction)\cdot d\omega = p(direction)\cdot\frac{dA\cos(\theta)}{distance(\vec{p},\vec{q})^2} = \frac{dA}{A}.

上の式から p(direction) を求めると

p(direction) = \frac{distance(\vec{p},\vec{q})^2}{\cos(\theta)\cdot A}.

この式を使って，color 関数の中で光源を直接サンプリングするようにしてみます．

vec3 color(const rayt::Ray& r, const Shape* world, int depth) {
    HitRec hrec;
    if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
        vec3 emitted = hrec.mat->emitted(r, hrec);
        ScatterRec srec;
        if (depth < MAX_DEPTH && hrec.mat->scatter(r, hrec, srec)) {

            vec3 on_light = vec3(213 + drand48()*(343 - 213), 554, 227 + drand48()*(332 - 227));
            vec3 to_light = on_light - hrec.p;
            float distance_squared = lengthSqr(to_light);
            to_light = normalize(to_light);
            if (dot(to_light, hrec.n) < 0) {
                return emitted;
            }
            float light_area = (343 - 213)*(332 - 227);
            float light_cosine = fabs(to_light.getY());
            if (light_cosine < 0.000001) {
                return emitted;
            }

            srec.pdf_value = distance_squared / (light_cosine * light_area);
            srec.ray = Ray(hrec.p, to_light);

            float spdf_value = hrec.mat->scattering_pdf(srec.ray, hrec);
            vec3 albedo = srec.albedo * spdf_value;
            return emitted + mulPerElem(albedo, color(srec.ray, world, depth + 1)) / srec.pdf_value;
        }
        else {
            return emitted;
        }
    }
    return background(r.direction());
}

これを実行すると次のようになります．(rayt303.cpp)

tuto-raytracing-direct-light-output

光源の周りにつぶつぶがありますね．これは光源と天井の間に隙間があるためです．今の光源は両面から当たっても発光するので，片面のみにします．そのため，DiffuseLight::emitted を変更します．

virtual vec3 emitted(const Ray& r, const HitRec& hrec) const override {
    if (dot(hrec.n, r.direction()) < 0) {
        return m_emit->value(hrec.u, hrec.v, hrec.p);
    }
    else {
        return vec3(0);
    }
}

また，コーネルボックスの光源を反転させます．

world->add(builder.rectXZ(213, 343, 227, 332, 554, light).flip().get());

レンダリングすると次のようになります．(rayt304.cpp)

tuto-raytracing-direct-light-flip-output

改善されているのがわかると思います．

合成PDF

これまでに，光源をサンプリングするための pdf と，\cos(\theta) に従った 余弦 pdf を作成しました．これらは合成することができます．pdf は区間[0,1]であればいいので，複数の pdf を重み加算します．例えば，光源 pdf と 反射 pdf をそれぞれ 0.5 の重みとすると次のような合成 pdf が作れます．

\text{mix\_pdf(direction)} = \frac{1}{2} \text{light\_pdf(direction)} + \frac{1}{2} \text{reflection\_pdf(direction)}.

合成PDFを作るために，pdf をクラスにします．

class Pdf {
public:
    virtual float value(const HitRec& hrec, const vec3& direction) const = 0;
    virtual vec3 generate(const HitRec& hrec) const = 0;
};

方向ベクトルから確率密度を求める value 関数と，確率分布に従った方向ベクトルを生成する generate 関数があります．hrec は生成に必要な情報が含まれていることがあるので渡しています．

まずは 余弦 pdf を実装すると

class CosinePdf : public Pdf {
public:
    CosinePdf() { }

    virtual float value(const HitRec& hrec, const vec3& direction) const override {
        float cosine = dot(normalize(direction), hrec.n);
        if (cosine > 0) {
            return cosine / PI;
        }
        else {
            return 0;
        }
    }

    virtual vec3 generate(const HitRec& hrec) const override {
        ONB uvw; uvw.build_from_w(hrec.n);
        vec3 v = uvw.local(random_cosine_direction());
        return v;
    }
};

この余弦 pdf を試してみます．color 関数を次のようにしてみます．

vec3 color(const rayt::Ray& r, const Shape* world, int depth) {
    HitRec hrec;
    if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
        vec3 emitted = hrec.mat->emitted(r, hrec);
        ScatterRec srec;
        if (depth < MAX_DEPTH && hrec.mat->scatter(r, hrec, srec)) {
            CosinePdf cpdf;
            srec.ray = Ray(hrec.p, cpdf.generate(hrec));
            srec.pdf_value = cpdf.value(hrec, srec.ray.direction());
            if (srec.pdf_value > 0) {
                float spdf_value = hrec.mat->scattering_pdf(srec.ray, hrec);
                vec3 albedo = srec.albedo * spdf_value;
                return emitted + mulPerElem(albedo, color(srec.ray, world, depth + 1)) / srec.pdf_value;
            }
            else {
                return emitted;
            }
        }
        else {
            return emitted;
        }
    }
    return background(r.direction());
}

この結果は次のようになります．(rayt305.cpp)

tuto-raytracing-cosine-pdf-output

次に物体そのものを pdf として扱えるようにします．これで，四角形を光源としてサンプリングできるようになります．最初に shape クラスに pdf 関連の関数を追加します．

class Shape {
public:
    virtual bool hit(const Ray& r, float t0, float t1, HitRec& hrec) const = 0;
    virtual float pdf_value(const vec3& o, const vec3& v) const { return 0; }
    virtual vec3 random(const vec3& o) const { return vec3(1, 0, 0); }
};

Rect クラスでオーバーライドします．

virtual float pdf_value(const vec3& o, const vec3& v) const override {
    if (m_axis != kXZ) return 0;
    HitRec hrec;
    if (this->hit(Ray(o, v), 0.001f, FLT_MAX, hrec)) {
        float area = (m_x1 - m_x0) * (m_y1 - m_y0);
        float distance_squared = pow2(hrec.t) * lengthSqr(v);
        float cosine = fabs(dot(v, hrec.n)) / length(v);
        return distance_squared / (cosine * area);
    }
    else {
        return 0;
    }
}

virtual vec3 random(const vec3& o) const override {
    if (m_axis != kXZ) return vec3(1, 0, 0);
    float x = m_x0 + drand48()*(m_x1 - m_x0);
    float y = m_y0 + drand48()*(m_y1 - m_y0);
    vec3 random_point;
    switch (m_axis) {
    case kXY:
        random_point = vec3(x, y, m_k);
        break;
    case kXZ:
        random_point = vec3(x, m_k, y);
        break;
    case kYZ:
        random_point = vec3(m_k, x, y);
        break;
    }
    vec3 v = random_point - o;
    return v;
}

color 関数も書き換えます．

vec3 color(const rayt::Ray& r, const Shape* world, int depth) {
    HitRec hrec;
    if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
        vec3 emitted = hrec.mat->emitted(r, hrec);
        ScatterRec srec;
        if (depth < MAX_DEPTH && hrec.mat->scatter(r, hrec, srec)) {
            ShapePtr p(std::make_shared<Rect>(
                213, 343, 227, 332, 554, Rect::kXZ, MaterialPtr()));
            ShapePdf cpdf(p, hrec.p);
            srec.ray = Ray(hrec.p, cpdf.generate(hrec));
            srec.pdf_value = cpdf.value(hrec, srec.ray.direction());

            float spdf_value = hrec.mat->scattering_pdf(srec.ray, hrec);
            vec3 albedo = srec.albedo * spdf_value;
            return emitted + mulPerElem(albedo, color(srec.ray, world, depth + 1)) / srec.pdf_value;
        }
        else {
            return emitted;
        }
    }
    return background(r.direction());
}

この結果は次のようになります．(rayt306.cpp)

tuto-raytracing-shape-pdf-output

それでは合成 pdf を実装します．今回は２つの pdf を合成するクラス MixturePdf とします．重みは 0.5 です．無作為なベクトルはそれぞれ 50% の確率でどちらかの pdf を使用して生成するようにしています．

class MixturePdf : public Pdf {
public:
    MixturePdf(const Pdf* p0, const Pdf* p1) { m_pdfs[0] = p0; m_pdfs[1] = p1; }

    virtual float value(const HitRec& hrec, const vec3& direction) const override {
        float pdf0_value = m_pdfs[0]->value(hrec, direction);
        float pdf1_value = m_pdfs[1]->value(hrec, direction);
        return 0.5f*pdf0_value + 0.5f*pdf1_value;
    }

    virtual vec3 generate(const HitRec& hrec) const override {
        if (drand48()<0.5f) {
            return m_pdfs[0]->generate(hrec);
        }
        else {
            return m_pdfs[1]->generate(hrec);
        }
    }

private:
    const Pdf* m_pdfs[2];
};

これを使って，光源 pdf と余弦 pdf を合成してみます．

vec3 color(const rayt::Ray& r, const Shape* world, int depth) {
    HitRec hrec;
    if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
        vec3 emitted = hrec.mat->emitted(r, hrec);
        ScatterRec srec;
        if (depth < MAX_DEPTH && hrec.mat->scatter(r, hrec, srec)) {
            ShapePtr p(std::make_shared<Rect>(
                213, 343, 227, 332, 554, Rect::kXZ, MaterialPtr()));
            ShapePdf shapePdf(p, hrec.p);
            CosinePdf cosPdf;
            MixturePdf mixPdf(&shapePdf, &cosPdf);
            srec.ray = Ray(hrec.p, mixPdf.generate(hrec));
            srec.pdf_value = mixPdf.value(hrec, srec.ray.direction());
            if (srec.pdf_value > 0) {
                float spdf_value = hrec.mat->scattering_pdf(srec.ray, hrec);
                vec3 albedo = srec.albedo * spdf_value;
                return emitted + mulPerElem(albedo, color(srec.ray, world, depth + 1)) / srec.pdf_value;
            }
            else {
                return emitted;
            }
        }
        else {
            return emitted;
        }
    }
    return background(r.direction());
}

この結果は次のようになります．(rayt307.cpp)

tuto-raytracing-mixture-pdf-output

これまでは color 関数で pdf を作成していましたが，本来は材質によって pdf は変わります．この仕組みを入れましょう．まずは，ScatterRec の pdf_value を Pdf のポインタに変更します．

class ScatterRec {
public:
    Ray	ray;
    vec3 albedo;
    const Pdf* pdf;
};

ランバート材質の pdf を設定します．

class Lambertian : public Material {
public:
    Lambertian(const TexturePtr& a)
        : m_albedo(a) {
    }

    virtual bool scatter(const Ray& r, const HitRec& hrec, ScatterRec& srec) const override {
        srec.albedo = m_albedo->value(hrec.u, hrec.v, hrec.p);
        srec.pdf = &m_pdf;
        return true;
    };

    virtual float scattering_pdf(const Ray& r, const HitRec& hrec) const override {
        return std::max(dot(hrec.n, normalize(r.direction())), 0.0f) / PI;
    }

private:
    TexturePtr m_albedo;
    CosinePdf  m_pdf;
};

次に color 関数を書き換えます．

vec3 color(const rayt::Ray& r, const Shape* world, int depth) {
    HitRec hrec;
    if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
        vec3 emitted = hrec.mat->emitted(r, hrec);
        ScatterRec srec;
        if (depth < MAX_DEPTH && hrec.mat->scatter(r, hrec, srec)) {
            ShapePtr p(std::make_shared<Rect>(
                213, 343, 227, 332, 554, Rect::kXZ, MaterialPtr()));
            ShapePdf shapePdf(p, hrec.p);
            MixturePdf mixPdf(&shapePdf, srec.pdf);
            srec.ray = Ray(hrec.p, mixPdf.generate(hrec));
            float pdf_value = mixPdf.value(hrec, srec.ray.direction());
            if (pdf_value > 0) {
                float spdf_value = hrec.mat->scattering_pdf(srec.ray, hrec);
                vec3 albedo = srec.albedo * spdf_value;
                return emitted + mulPerElem(albedo, color(srec.ray, world, depth + 1)) / pdf_value;
            }
            else {
                return emitted;
            }
        }
        else {
            return emitted;
        }
    }
    return background(r.direction());
}

この結果は前回と変わらないはずです．(rayt308.cpp)

金属の復活

コーネルボックスの材質はランバートのみでしたが，金属材質も現状に合わせてみます．今回の鏡面反射では反射方向が決まっているので，pdf は必要ないと考えます．そのため，材質が金属の場合は単純に発光と反射方向からの放射束のみとします．それには，鏡面反射かどうか知る必要があるので，ScatterRec にその情報を追加します．

class ScatterRec {
public:
    Ray	ray;
    vec3 albedo;
    const Pdf* pdf;
    bool is_specular;
};

金属材質の scatter 関数を書き換えます．また，dielectric の方にも同様に追加します．

virtual bool scatter(const Ray& r, const HitRec& hrec, ScatterRec& srec) const override {
    vec3 reflected = reflect(normalize(r.direction()), hrec.n);
    reflected += m_fuzz*random_in_unit_sphere();
    srec.ray = Ray(hrec.p, reflected);
    srec.albedo = m_albedo->value(hrec.u, hrec.v, hrec.p);
    srec.pdf = nullptr;
    srec.is_specular = true;
    return dot(srec.ray.direction(), hrec.n) > 0;
}

ランバート材質の scatter 関数で is_specular を false に設定します．

virtual bool scatter(const Ray& r, const HitRec& hrec, ScatterRec& srec) const override {
    srec.albedo = m_albedo->value(hrec.u, hrec.v, hrec.p);
    srec.pdf = &m_pdf;
    srec.is_specular = false;
    return true;
};

そして，color 関数で，鏡面反射のときの処理を追加します．

vec3 color(const rayt::Ray& r, const Shape* world, int depth) {
    HitRec hrec;
    if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
        vec3 emitted = hrec.mat->emitted(r, hrec);
        ScatterRec srec;
        if (depth < MAX_DEPTH && hrec.mat->scatter(r, hrec, srec)) {
            if (srec.is_specular) {
                return emitted + mulPerElem(srec.albedo, color(srec.ray, world, depth + 1));
            }
            else {
                ShapePtr p(std::make_shared<Rect>(
                    213, 343, 227, 332, 554, Rect::kXZ, MaterialPtr()));
                ShapePdf shapePdf(p, hrec.p);
                MixturePdf mixPdf(&shapePdf, srec.pdf);
                srec.ray = Ray(hrec.p, mixPdf.generate(hrec));
                float pdf_value = mixPdf.value(hrec, srec.ray.direction());
                if (pdf_value > 0) {
                    float spdf_value = hrec.mat->scattering_pdf(srec.ray, hrec);
                    vec3 albedo = srec.albedo * spdf_value;
                    return emitted + mulPerElem(albedo, color(srec.ray, world, depth + 1)) / pdf_value;
                }
                else {
                    return emitted;
                }
            }
        }
        else {
            return emitted;
        }
    }
    return background(r.direction());
}

あとはコーネルボックスの箱を１つをアルミニウム（金属）にしてみます．

MaterialPtr aluminum = std::make_shared<Metal>(
    std::make_shared<ColorTexture>(vec3(0.8f, 0.85f, 0.88f)), 0.0f);

ShapeList* world = new ShapeList();
ShapeBuilder builder;
world->add(builder.rectYZ(0, 555, 0, 555, 555, green).flip().get());
world->add(builder.rectYZ(0, 555, 0, 555, 0, red).get());
world->add(builder.rectXZ(213, 343, 227, 332, 554, light).flip().get());
world->add(builder.rectXZ(0, 555, 0, 555, 555, white).flip().get());
world->add(builder.rectXZ(0, 555, 0, 555, 0, white).get());
world->add(builder.rectXY(0, 555, 0, 555, 555, white).flip().get());
world->add(builder.box(vec3(0), vec3(165), white)
    .rotate(vec3::yAxis(), -18)
    .translate(vec3(130, 0, 65))
    .get());
world->add(builder.box(vec3(0), vec3(165, 330, 165), aluminum)
    .rotate(vec3::yAxis(), 15)
    .translate(vec3(265, 0, 295))
    .get());
m_world.reset(world);

これを実行すると次のような画像になります．(rayt309.cpp)

tuto-raytracing-metal-pdf-output

左上のところにノイズが強く発生してしまっています．これはアルミニウム箱へのサンプリングが減ってしまったためです．そのため，その箱に重点的サンプリングすると改善できるはずです．ただし，今回は箱ではなく球にして重点的サンプリングします．その理由は箱より簡単だからです．

ある点から球に対して重点的サンプリングする場合，球の立体角を一様分布でサンプリングすればいいと考えられます．一様分布に従う球のベクトルについてはすでにやりました．それによると，

r_2 = 2\pi \int_{0}^{\theta}f(t)\sin(t) dt.

今，f(t) が未知なので，C=f(t) とすると

r_2 = 2\pi \int_{0}^{\theta}C\sin(t) dt.

この C を求めると

\begin{aligned} r_2 &= 2\pi\int_{0}^{\theta}C\sin(t) dt \\[2ex] &= 2\pi C\left[-cos(t)\right]_{0}^{\theta} \\[2ex] &= 2\pi C\left\{(-\cos(\theta))-(-\cos(\theta))\right\} \\[2ex] &= 2\pi C(-\cos(\theta) + 1) \\[2ex] &= 2\pi C(1-\cos(\theta)). \end{aligned}

よって，

C = -\frac{r_2-1}{2\pi\cos(\theta)}.

また，

\cos(\theta) = 1-\frac{r_2}{2\pi C}.

r_2=1 のとき，\theta_{max} が得られるとすると，C は

\begin{aligned} 2\pi C(1-\cos(\theta)) &= r_2 \\[2ex] 2\pi C(1-\cos(\theta_{max})) &= 1 \\[2ex] C(1-\cos(\theta_{max})) &= \frac{1}{2\pi} \\[2ex] C &= \frac{1}{2\pi(1-\cos(\theta_{max}))}. \end{aligned}

これを先ほどの式に代入すると

\begin{aligned} \cos(\theta) &= 1-\frac{r_2}{2\pi C} \\[2ex] &= 1-\frac{r_2}{2\pi}\cdot 2\pi(1-\cos(\theta_{max})) \\[2ex] &= 1-r_2(1-\cos(\theta_{max})). \end{aligned}

これを使ってベクトルを生成できます．\phi は前と同じです．デカルト座標で表すと

\begin{aligned} z &= \cos(\theta) = 1 - r_2(\cos(1-\theta_{max})) \\[2ex] x &= \cos(\phi)sin(\theta) = \cos(2\pi r_1)\sqrt{1-z^2} \\[2ex] y &= \sin(\phi)sin(\theta) = \sin(2\pi r_1)\sqrt{1-z^2}. \end{aligned}

ここで，\theta_{max} は図のような関係になっています．

tuto-raytracing-theta-max

この図から \theta_{max} を求めると

\sin(\theta_{max}) = \frac{R}{distance(\vec{c}-\vec{p})}.

sin^2\theta+cos^2\theta=1 を使って

\cos(\theta_{max}) = \sqrt{1-\frac{R^2}{distance(\vec{c}-\vec{p})^2}}.

次に方向分布を表す pdf を算出する必要があります．球の全立体角は \omega で，ある点から球に向かう方向の微小立体角は \frac{1}{\omega} です．今は単位球上の領域を考えているので

\omega = \int_{0}^{2\pi} \int_{0}^{\theta_{max}} \sin(\theta) d\theta d\phi = 2\pi(1-\cos(\theta_{max})).

これらを使って sphere クラスの pdf_value 関数と random 関数を実装します．

virtual float pdf_value(const vec3& o, const vec3& v) const override {
    HitRec hrec;
    if (this->hit(Ray(o, v), 0.001f, FLT_MAX, hrec)) {
        float dd = lengthSqr(m_center - o);
        float rr = std::min(pow2(m_radius), dd);
        float cos_theta_max = sqrtf(1.0f - rr*recip(dd));
        float solid_angle = PI2*(1.0f - cos_theta_max);
        return recip(solid_angle);
    }
    else {
        return 0;
    }
}

virtual vec3 random(const vec3& o) const override {
    vec3 direction = m_center - o;
    float distance_squared = lengthSqr(direction);
    ONB uvw; uvw.build_from_w(direction);
    vec3 v = uvw.local(random_to_sphere(m_radius, distance_squared));
    return v;
}

ここで，random_to_sphere はある点から球方向に向かう無作為なベクトルを生成します．

inline vec3 random_to_sphere(float radius, float distance_squared) {
    float r1 = drand48();
    float r2 = drand48();
    float rr = std::min(pow2(radius), distance_squared);
    float cos_theta_max = sqrtf(1.f - rr * recip(distance_squared));
    float z = 1.0f - r2*(1.0f - cos_theta_max);
    float sqrtz = sqrtf(1.f - pow2(z));
    float phi = PI2*r1;
    float x = cosf(phi) * sqrtz;
    float y = sinf(phi) * sqrtz;
    return vec3(x, y, z);
}

これで，球をサンプリングできるようになります．ここで，少し Scene クラスを変更して，重点的サンプリング用の物体を渡せるようにします．新たに m_light メンバを追加して，build 関数で生成します．

std::unique_ptr<Shape> m_light;

color 関数の引数に追加して，処理を行います．

vec3 color(const rayt::Ray& r, const Shape* world, const Shape* light, int depth) {
    HitRec hrec;
    if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
        vec3 emitted = hrec.mat->emitted(r, hrec);
        ScatterRec srec;
        if (depth < MAX_DEPTH && hrec.mat->scatter(r, hrec, srec)) {
            if (srec.is_specular) {
                return emitted + mulPerElem(srec.albedo, color(srec.ray, world, light, depth + 1));
            }
            else {
                ShapePdf shapePdf(light, hrec.p);
                MixturePdf mixPdf(&shapePdf, srec.pdf);
                srec.ray = Ray(hrec.p, mixPdf.generate(hrec));
                float pdf_value = mixPdf.value(hrec, srec.ray.direction());
                if (pdf_value > 0) {
                    float spdf_value = hrec.mat->scattering_pdf(srec.ray, hrec);
                    vec3 albedo = srec.albedo * spdf_value;
                    return emitted + mulPerElem(albedo, color(srec.ray, world, light, depth + 1)) / pdf_value;
                }
                else {
                    return emitted;
                }
            }
        }
        else {
            return emitted;
        }
    }
    return background(r.direction());
}

render 関数の中で color 関数を呼んでいるところも変更します．

c += color(r, m_world.get(), m_light.get(), 0);

あとは，コーネルボックスの箱を１つ球に変更して，材質も Metal から Dielectric に変えます．

MaterialPtr metal = std::make_shared<Dielectric>(1.5f);

ShapeList* world = new ShapeList();
ShapeBuilder builder;
world->add(builder.rectYZ(0, 555, 0, 555, 555, green).flip().get());
world->add(builder.rectYZ(0, 555, 0, 555, 0, red).get());
world->add(builder.rectXZ(213, 343, 227, 332, 554, light).flip().get());
world->add(builder.rectXZ(0, 555, 0, 555, 555, white).flip().get());
world->add(builder.rectXZ(0, 555, 0, 555, 0, white).get());
world->add(builder.rectXY(0, 555, 0, 555, 555, white).flip().get());
world->add(builder.sphere(vec3(190, 90, 190), 90, metal).get());
world->add(builder.box(vec3(0), vec3(165, 330, 165), white)
    .rotate(vec3::yAxis(), 15)
    .translate(vec3(265, 0, 295))
    .get());
m_world.reset(world);

m_light.reset(new Sphere(
    vec3(190, 90, 190), 90, MaterialPtr()));

レンダリングすると次のようになります．(rayt310.cpp)

tuto-raytracing-sphere-pdf-output

ノイズがひどいですね．これは重点的サンプリング用の物体が球だけだからです．ShapeList に pdf の関数を追加して，複数の物体に対して重点的サンプリング出来るように対応します．重みは均等とし，ベクトルの生成もリストの中から無作為に選んだ物体を使用します．

virtual float pdf_value(const vec3& o, const vec3& v) const override {
    float weight = 1.0f / m_list.size();
    float sum = 0;
    for (auto& p : m_list) {
        sum += weight * p->pdf_value(o, v);
    }
    return sum;
}

virtual vec3 random(const vec3& o) const override {
    size_t n = m_list.size();
    size_t index = size_t(drand48() * n);
    if (n > 0 && index >= n) {
        index = n - 1;
    }
    return m_list[index]->random(o);
}

これを使って，重点的サンプリング用の物体が光源（四角形）と誘電体（球）で構成されるようにします．

ShapeList* l = new ShapeList();
l->add(builder.rectXZ(213,343,227,332,554,MaterialPtr()).get());
l->add(builder.sphere(vec3(190,90,190),90,MaterialPtr()).get());
m_light.reset(l);

これをレンダリングすると次のようになります．(rayt311.cpp)

tuto-raytracing-shape-list-pdf-output

これが目標の画像になります．

さて，とうとう最後になりましたが，どこかの計算でゼロ除算があったり，sqrt 関数に負の値を入れたりとか，初期化していない変数にアクセスしたりとかで，値がおかしくなる場合があるかもしれません．計算結果がそのような値になったら黒になるようにフィルターを追加しておきましょう．今さらな感じがしますけどね．

// rayt.h
class DenanFilter : public ImageFilter {
public:
    DenanFilter() {}
    virtual vec3 filter(const vec3& c) const override {
        return vec3(
            !(c[0] == c[0]) ? 0 : c[0],
            !(c[1] == c[1]) ? 0 : c[1],
            !(c[2] == c[2]) ? 0 : c[2]);
    }
};

これをフィルターリストに追加します．(rayt312.h)

Image(int w, int h) {
    m_width = w;
    m_height = h;
    m_pixels.reset(new rgb[m_width*m_height]);
    m_filters.push_back(std::make_unique<DenanFilter>());
    m_filters.push_back(std::make_unique<GammaFilter>(GAMMA_FACTOR));
    m_filters.push_back(std::make_unique<TonemapFilter>());
}

黒くしていますが，赤とかにしておくと NaN かどうかがわかりやすくなると思います．

以上でこの入門は終了です．お疲れ様でした！

📌 この次はどうするか

ロシアンルーレット

レイトレーシングでは光線が物体の表面に当たるたびに，反射した光線をさらに追跡するといった再帰的な処理になっています．この再帰処理はどこかで打ち切らなければなりません．今回は特定の回数に達したときに打ち切るようになっていますが，その光線に適した再帰回数を決められるのが理想です．そこでロシアンルーレットと呼ばれる手法があります．これは重点的サンプリングのように，どれくらい再帰処理をすればよいかを確率によって決定します．

パストレーシング

コンピュータグラフィックスでは物理ベースレンダリングが主流になっています．この根源となる数式がレンダリング方程式とよばれるもので，モンテカルロレイトレーシングを使ってレンダリング方程式を計算することをパストレーシングといいます．

関与媒質

埃や煙，雲といったものを関与媒質ということがあります．これは微粒子が集まっているものと考えられ，光がそれを通過するときに散乱や吸収が行われます．

フォトンマッピング

光の最も小さい単位は「光子（フォトン）」です．このフォトンが光源からどのようにばらまかれているかを表した分布データ（フォトンマップ）を作成し，それを追跡することでレンダリングを行う手法です．

📌 参考となる資料

Peter Shirley 著の「Ray Tracing」

今回の内容のベースとなっている Peter Shirley 著の３冊をお勧めします．

この入門では取り扱っていない関与媒質，BVH，ノイズ，モーションブラー，被写界深度についても書かれています．一部名前を変えていますが，ほとんどが同じ名前になっているので，わかりやすくなっているのではないかと思います．冒頭でもお伝えした通り，この本の内容は現在インターネットから無料で手に入ります．

https://raytracing.github.io/

smallpt

http://www.kevinbeason.com/smallpt/

とても小さいコードのパストレーシングです．次にパストレーシングに挑戦するなら参考になると思います．
smallpt の解説もありますのでこちらも参考にしてみてください．

Presentation slides about smallpt by David Cline

edupt

http://kagamin.net/hole/edupt/

上記の smallpt を c++ で書き直し，ソースコードには日本語のコメントが沢山付いている物理ベースレンダラです．とりあえず，次のことに迷ったらこちらのソースコードを見ることをお勧めします．edupt の解説もありますので参考にしてみてください．（上記リンク先から参照できます）

また，作者の方は「Computer Graphics Gems JP 2015　コンピュータグラフィックス技術の最前線」という本で，「パストレーシングで始めるオフライン大域照明レンダリング入門」を執筆されているので，そちらもお勧めします．

memoRANDOM

http://rayspace.xyz/CG/

レイトレーシングに関する情報が豊富にある日本語のサイトです．とても重宝しています．

「フォトンマッピングー実写に迫るコンピュータグラフィックス」

amazon: フォトンマッピングー実写に迫るコンピュータグラフィックス

フォトンマッピングについて書かれている翻訳書です．フォトンマッピングを実装するときや，レイトレーシングをするうえでも必読の本だと言えます．この本を読んでわからないところがあればそれを調べていくというやり方でも実力は付いていくと思います．

「CG Magic」

amazon: CG Magic

レンダリングに関する技術的な解説をしている邦書です．とりあえず，どんな技術があるのか確認したい場合にお勧めします．

Robust Monte Carlo Methods for Light Transport Simulation

Robust Monte Carlo Methods for Light Transport Simulation

モンテカルロ法について詳しく書かれている文書です．400ページ以上あるのに，無料で手に入るので英語ができるならぜひ読んでみてください．あわせて， Global Illumination compendium もおすすめします．

「Physically Based Rendering」

http://www.pbrt.org/

1000ページ以上あるトンデモナイ本ですが，内容はとても充実しているそうです．個人的にはさらっと見ただけで，ちゃんと読んではいません．洋書なので，知識と根気があるなら読んでみてはいかがでしょうか．

レイトレ合宿

レイトレーシングの作品を出して技術を競い合っている人たちがいます．また，技術力を高めるためにソースコードや資料の公開をしてくれています．大変勉強になりますので，興味のある方は見てみてください．

📌 最後に

最近勉強したので，そのまとめとして書きました．説明不足や理解不足があるかもしれません．少しでも参考になれば幸いです．

📌 参考文献

  • Peter Shirley, 「Ray Tracing in One Weekend」, 2016
  • Peter Shirley, 「Ray Tracing: the Next Week」, 2016
  • Peter Shirley, 「Ray Tracing: The Rest Of Your Life」, 2016
  • Eric Veach, 「Robust Monte Carlo Methods for Light Transport Simulation」, 1997
  • Henrik Wann Jenson著，苗村健訳「フォトンマッピングー実写に迫るコンピュータグラフィックス」オーム社，2002
  • 倉地紀子，「CG Magic：レンダリング」オーム社，2007
  • Eric Lengyel著，狩野智英訳「ゲームプログラミングのための３Dグラフィックス数学」ボーンデジタル社，2006
  • 鈴木健太郎，「パストレーシングで始めるオフライン大域照明レンダリング入門」『Computer Graphics Gems JP 2015 コンピュータグラフィックス技術の最前線』ボーンデジタル，2015, p. 3