パストレがしたい - ややプログラム紀行

ｽﾞｻ━━━━⊂（ﾟДﾟ⊂⌒｀つ≡≡≡━━━━!!

この記事はISer Advent Calendar 2018 - Adventarの22日目の記事として書かれました。

10月の半ばにPhysically Based Rendering: From Theory to Implementation(通称pbr)という本がネットで無料公開されたので、今回はこれを読んでパストレ(path tracing)をやってみようと思います。*1

1日分の記事にまとめるにはあまりにも書くべき内容が多いことに気づきましたがもはや手遅れなので、内容を希釈してpbrが読みやすくなる程度のあらすじを書いていきたいと思います。pbrに載っている内容以上のことは説明しない(できない)ので、詳しくは本を参照してください。

レイトレーシング

3年生以上はCPU実験でmin-rtでお馴染みと思いますが、レイトレーシングとは光線(レイ)をシミュレートして画像を生成する技術です(雑
pbrで取り上げられるレイトレーシングは特にPhysically Based Ray Tracing(物理ベースレイトレーシング)といって、従来の方法と違い物理法則に忠実に計算を行うことでよりリアルな画像を生成することができます。
物理ベースは本来金属のレンダリングなどで本領を発揮するものですが、今回は尺の都合上省きます。

実装方針

pbr本で紹介されている実装は基本的にpbrt*2で実際に実装されているので参考にしてください。

github.com

pbrtはC++のコードとしても参考になるほど綺麗に書かれているので、今回はクラスの構造などはなるべく保ちつつかなり単純化したものを実装していきたいと思います。

今回重要なクラスは以下の通りです

Shape...図形の形状に関連するクラス。衝突判定に使う。
BRDF...レイの反射を扱うクラス。
Primitive...Shapeと物体の素材(今回は色のみ)をまとめるクラス。
Scene...Primitiveをまとめるクラス。
Integrator...実際に描画を行うクラス。

画像出力

早速実装していこうと思いますが、まずは画像の出力ができる下地を整えておきましょう。
といっても、今回はヘッダーオンリーライブラリのstb_image/stb_image_writeを使わせてもらいます。
github.com
インクルードを通し、謎のものをdefineすれば導入できます。めちゃくちゃ簡単ですね。

#define STB_IMAGE_IMPLEMENTATION
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb_image.h"
#include "stb_image_write.h"

あとはchar型の配列を用意し、いろいろ値をいれたあとにstbi_image_writeを呼び出せば画像が出力されます。

const int width = 256;
const int height = 256;
char pixels[width * height * 3];
stbi_write_png("out.png", width, height, STBI_rgb, pixels, 0);

数学関連のクラス

早速実装していこうと思いますが、まずはベクトル演算やレイのためのクラスを用意する必要があります。
以下のような動作ができるVectorクラスを用意しておいてください(丸投げ*3

Vector<float, 3> v1{1.0,2.0,3.0};
Vector<float, 3> v2{1.0,2.0,3.0};
auto v = v1 + v2;
auto v3 = v - v1;   
std::cout<<"v.dot(v3) = "<<v.dot(v3)<<std::endl;    // 内積
auto n = v.normalize();  // 正規化
auto ts = cross(v, v3);    // 外積

実際にはVectorばかり使うことになるのでVector3fなどとエイリアス宣言しときましょう。

以下のようなRayクラスも準備してください。Rayクラスは内部的に始点を示すベクトルと方向ベクトルを保持し、float値tを与えると始点+方向ベクトル*tを表すベクトルを返すような関数を用意すると楽です。

Ray ray(v, v3);
Vector3f dest = ray(0.5);

最後に座標変換のためのMatrixクラス、変換、逆変換をまとめたTransformクラス及び座標変換用のヘルパー関数も実装しときましょう。
Transformクラスは順変換と逆変換用に4行4列の行列を2つ保持し、VectorおよびRayクラスに作用させられるように関数を用意してください。
4行4列なのは3次元での平行移動、回転、スケール変換を統一的に扱うことができて便利だからですが、詳しくはpbrを見てみてください。(Transformations)

Matrix<float, 2, 3> m1{{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}};                                                                                                                                                        
Matrix<float, 3, 2> m2{{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
auto v4 = (m1 * m2) * v;
auto v5 = (rotateX(0.5) * rotateZ(0.3)).act(v4);

以上でとりあえずのところはなんとかなると思います。~~正直ここら辺はどうでもいい~~

衝突判定

レイトレーシングではカメラからレイを飛ばして、衝突判定したら衝突地点の色やら明るさやらを調べて...という流れなので、まずは衝突判定を実装しましょう。
判定はShapeクラスで行います。

class Shape {
public:
    Shape() = default;
    virtual ~Shape() = default;
    virtual bool intersect(const Ray&, float*) const = 0;
};

このクラスをオブジェクトの形ごとに継承していこうと思います。実際に重要なのはintersect関数で、これの返り値がtrueの時は衝突したことを意味します。第2引数のfloat*には「Rayをtだけ伸ばしたら衝突する」のtの値を格納します。

とりあえず板を実装してみます。intersect関数に渡されるRayはそのオブジェクトのローカル座標上にあると想定するので、奥行きz軸としたときレイとz=0の交点を調べればよいです。

class Plane : public Shape {                                                                                                                                                                                         
public:
    Plane() = default;
    ~Plane() = default;
    bool intersect(const Ray& _ray, float* _t) const {
        Vector3f o = _ray.getOrigin();
        Vector3f d = _ray.getDir();
        if (d.at(2) == 0.0) return false;

        float t = std::abs(o.at(2)) / d.at(2);
        if (t < 0.0000001) return false;
        Vector3f point = _ray(t);
        if (point.at(0) < -0.5 || 0.5 < point.at(0)
                || point.at(1) < -0.5 || 0.5 < point.at(1)) return false;
        *_t = t; 
        return true;
    }
};

球も同様にして2次方程式を解くことで実装できます。

描画担当部分

描画はおそらくレイトレーシングのコア的部分であるIntegratorクラスで行います。この名前は、このクラスのメイン作業が積分処理だからだそうです。

class Integrator {
public:
    virtual ~Integrator() = default;
    void render() {
        std::shared_ptr<Shape> shape(static_cast<Shape*>(new Plane));
        
        const int width = 256;
        const int height = 256;
        char pixels[width * height * 3];
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
                float x = (float)(j - 128) / (float)width * 2.0;
                float y = (float)(i - 128) / (float)height * 2.0;
                Ray ray(Vector3f{x, -y, -1.0f}, Vector3f{0.0f, 0.0f, 1.0f});
                
                int pos = (i * width + j) * 3;
                float t;
                if (!shape->intersect(ray, &t)) {
                    pixels[pos] = 0; 
                    pixels[pos+1] = 0; 
                    pixels[pos+2] = 0; 
                } else {
                    pixels[pos] = 255; 
                    pixels[pos+1] = 255; 
                    pixels[pos+2] = 255; 
                }
            }
        }
        stbi_write_png("out.png", width, height, STBI_rgb, pixels, 0);
    }
};

今回は画像の座標に応じてz軸に平行なレイをとばしてみましたが、実行すると以下のような画像が出力されます。
f:id:dora119:20181222044337p:plain
確かに原点のあたりに板がおいてあるのが分かると思います、レイトレできた😭

シーン管理

上の実装だとIntegratorクラスにオブジェクトを直書きしていますが、このままだと柔軟性に欠けるので、まずShapeクラスをまとめるPrimitiveクラスというものを実装します。
これは本来pbrtではShapeクラスとMaterialクラスを結びつける的な役割があるのですが、今回は尺の都合でMaterialクラスを排除する予定です。

class Primitive {
    std::shared_ptr<Shape> shape;
    Transform worldToLocal;

public:
    Primitive(const std::shared_ptr<Shape> _shape, const Transform& _t) : shape(_shape), worldToLocal(_t) {}
    bool intersect(const Ray& _ray, float* _t) const {
        Ray r = worldToLocal.inv(_ray);

        bool result = false;
        if (shape) result = shape->intersect(r, _t);
        if (result) *_t = worldToLocal.act(*_t);
        return result;
    }
};

これに加え、Primitiveクラスを管理するSceneクラスを作ります。

class Scene {
    std::vector<std::shared_ptr<Primitive>> primitives;

public:
    bool intersect(const Ray& _ray) const {
        float result = false;
        float min_t = std::numeric_limits<float>::max();
        for (auto&& primitive : primitives) {
            float t;
            if (primitive->intersect(_ray, &t)) {                                                                                                                                                                    
                result = true;
                min_t = std::min(min_t, t);
            }
        }
        return result;
    }
    void addPrimitive(const std::shared_ptr<Primitive> _primitive) { primitives.push_back(_primitive); }
};

PrimitiveとSceneを追加したことで次のように書けます。*4

Transform plane_wtl = translate({0.0f, -0.3f, 0.0f}) * rotateX(Pi / 4.0f) * rotateZ(Pi / 4.0f) * scale(2.0f, 2.0f, 2.0f);
std::shared_ptr<Shape> plane(new Plane);
std::shared_ptr<Primitive> p_plane(new Primitive(plane, plane_wtl));

Transform sphere_wtl = scale(0.3f, 0.3f, 0.3f);
std::shared_ptr<Shape> sphere(new Sphere);
std::shared_ptr<Primitive> p_sphere(new Primitive(sphere, sphere_wtl));

Scene scene; 
scene.addPrimitive(p_plane);
scene.addPrimitive(p_sphere);

Integrator integrator;
integrator.render(scene);

f:id:dora119:20181222055555p:plain — 上のコードにより生成された画像

とりあえず衝突判定および根幹になるクラスは作れたので、次は上の画像に色をつけていきましょう。

光学の用語

物理ベースレイトレーシングではいくつかの光に関する簡略化*5が行われていますが、基本的には光というエネルギーのやりとりをレイを通して計算することで画像を生成します。そこでRadiance(放射輝度)という光学の概念が必要になってきます。

正直これは図があったほうがわかりやすいので、ここでは文章で軽く説明して後はpbrを参照してもらうのが速いです。(Radiometry)

まず、Flux(放射束)とは、ある領域を通過する光のエネルギーを時間で微分したもので、式であらわすと $\phi = \frac{\mathrm{d}Q}{\mathrm{d}t}$ です。

さらに、Irradiance(放射照度)、Radiant Exitance(放射発散度)とはFluxを通過領域で微分したもので、式であらわすと $E(p) = \frac{\mathrm{d}\phi}{\mathrm{d}A}$ となります。

最後に、Radiant Exitance(もしくはIrradiance)に方向の区別もつけたものがRadiance(放射輝度)であり、定義的には放射束を光の放射源から微小面積 $\mathrm{d}A$ を投影した $\mathrm{d}A^\perp$ と立体角 $\omega$ で微分したものになります。つまり $L(p, \omega) = \frac{\mathrm{d}^2\phi}{\mathrm{d}A^\perp \mathrm{d}\omega} = \frac{1}{\cos \theta} \frac{\mathrm{d}^2 \phi}{\mathrm{d}A \mathrm{d}\omega}$ です(ここで $\theta$ は地点 $p$ での法線と放射源のなす角)。ここでちょっとpbrの図を引用してみます。

なぜこれを引用したかというと、僕の個人的なちょっとした不満ですが、この図だと微小立体角を考えると自然と投影された面積も微小になるように感じられるくないか？と思ったからです!(T_T)
混乱のもとになるような気がするんですけど僕だけですかね :thinking_face:

ともかく、このRadianceが分かっていればその他IrradianceやRadiant Exitanceも導出できるので、Radianceを中心に話が進んでいきます。

BRDF(理論)

上の議論より入射するRadianceを $L_i(p,\omega_i)$ と書くと $\mathrm{d}E(p, \omega_i) = L_i(p, \omega_i)\cos \theta_i \mathrm{d}\omega_i$ ということが分かりますが、放射するRadianceを $L_o(p, \omega_o)$ と書くとき $\frac{\mathrm{d}L_o}{\mathrm{d}E}$ を $f_r(p, \omega_o, \omega_i)$ という関数としてまとめてしまいます。*6
すると $f_r(p, \omega_o, \omega_i) = \frac{\mathrm{d}L_o}{\mathrm{d}E} = \frac{\mathrm{d}L_o}{L_i(p, \omega_i)\cos \theta_i \mathrm{d}\omega_i}$ なので、放射するRadianceは
$L_o(p, \omega_o) = \int_{H^2(n)}f_r(p, \omega_o, \omega_i) L_i(p, \omega_i) |\cos \theta_i| \mathrm{d}\omega_i$
となることがわかります。ここで $H^2(n)$ は法線 $n$ まわりの半球です。
さらに $f_r$ にはヘルムホルツの相反性とエネルギー保存則を満たしている必要があります。すなわち、

$f_r(p, \omega_o, \omega_i) = f_r(p, \omega_i, \omega_o)$
$\int_{H^2(n)} f_r(p, \omega_o, \omega_i) \cos \theta_i \mathrm{d}\omega_i \leq 1$

が要請されます。
このような $f_r$ をBidirectional Reflectance Distribution Function (BRDF、双方向反射率分布関数)といいます。BRDFをこちら側で用意すれば、後はこの積分を解くことで $L_o$ を求められるので画像の生成もできるはずです！

モンテカルロ法

上の積分を解く手法としてモンテカルロ法があります。まずはおさらいをしましょう。
ある関数 $f(x)$ の積分 $\int_a^b f(x) dx$ を求めたいとき、各 $X_i$ を $(a, b)$ の一様分布からサンプリングして $\frac{b-a}{N} \sum_{i=1}^N f(X_i)$ を計算するというものが基本的なモンテカルロ法です。実際に期待値を計算してみると、
$E\left(\frac{b-a}{N} \sum_{i=1}^N f(X_i)\right) = \frac{b-a}{N} \sum_{i=1}^N E(f(X_i)) = \frac{b-a}{N} \sum_{i=1}^N \int_a^b f(x)\frac{1}{b-a}dx = \int_a^b f(x) dx$
となり確かに期待値が積分値と一致することがわかります。これをもう少し発展させたものとしてimportance sampling(重点的サンプリング)というものがあります。
確率変数 $X_i$ をなるべく被積分関数の定数倍に近い確率密度関数 $p(x)$ を用いてサンプリングしたとき、 $\frac{1}{N}\sum_{i=1}^N \frac{f(X_i)}{p(X_i)}$ は同様に $\int f(x)dx$ と期待値が一致しますが、サンプリングをうまくできればできるほど分散が小さくなります。(実際に分散を計算してみてください)

しかし、上記のレンダリング方程式を解く場合、基本的に $f_r(p, w_o, w_i)L_i(p, w_i)|\cos \theta_i|$ がどのような形になっているか分からないため、そのままではImportance Samplingを行うことができません。*7
そこでさらに発展形としてMultiple Importance Sampling(MIS、多重重点的サンプリング)というものが用いられます。これは名前の通りImportance Samplingの複合形で、基本アイディアは $f(x)g(x)$ の積分値を求めたいときに各 $f(x)$ と $g(x)$ のそれぞれに対しImportance Samplingを適用し、結果の値をいいかんじ*8に合成しようというものです。
ただ、今回は尺の都合上 $f_r(p, w_o, w_i)L_i(p, w_i)$ のうち $L_i(p, w_i)$ 側だけのImportance Samplingのみで済むようにするので、MISは頭の片隅におくくらいにしてください。

BRDF(実装)

長い話も終わったので実際に実装をしていこうと思います。
まずBRDFには $f_r(p, w_o, w_i)$ と、この関数に合わせて $w_i$ をサンプリングするための $sample_f$ という関数が必要ですが、今回はライト側でのみサンプリングを行いたいので、BRDFには $f_r$ 関数があれば充分です。

using RGB = Vector3f;

class BRDF {  
public:
    BRDF() = default;
    virtual ~BRDF() = default;
    virtual RGB f(const Vector3f& _wo, const Vector3f& _wi) const = 0;
};

BRDFは衝突判定後に各点について生成されるので、BRDF::fに地点 $p$ を表す引数は存在しません。

今回は具体的なBRDFとして、もっとも単純なLambertian Reflectcion(ランバート反射)というものを実装しておきましょう。これは入射角と反射角に関係なく定数を返すBRDFです。

class LambertianReflection : public BRDF {
    const RGB color;
  
public:         
    LambertianReflection(const RGB& _color) : color(_color) {}
    ~LambertianReflection() = default;
    RGB f(const Vector3f& _wo, const Vector3f& _wi) const {
        return color / Pi;
    }
};

円周率で割っているのはBRDFのエネルギー保存則を満たすという要請のためです。めちゃくちゃシンプルですが、とりあえずこれでBRDFはできました。

ライト

ライトはLightクラスにより制御されます。BRDFと逆に今回はライト側にあわせてサンプリングを行うので $sampleLi$ 関数を実装する必要があります。

class Light {
public:
    Light() = default;
    virtual ~Light() = default;
    virtual RGB sampleLi(const Vector3f& _p, Vector3f* _wi, float* _pdf) const = 0;                                                                                                                                  
};

sampleLi関数には地点pと光の進む方向サンプリングする入射角に加え、その入射角がサンプリングされる確率を入れておくpdfという引数もあります。というのも、(Multiple) Importance Samplingでは確率密度も計算に用いるからです。

具体的なライトとして点ライトを実装してみましょう。

class PointLight : public Light {
    Vector3f pos;
    RGB color;

public:
    PointLight(const Vector3f& _p, const RGB& _c) : pos(_p), color(_c) {}
    ~PointLight() {}
    RGB sampleLi(const Vector3f& _p, Vector3f* _wi, float* _pdf) const { 
        *_wi = _p - pos;
        *_pdf = 1.0;
        float length = _wi->length();
        return color / (length * length);
    }
};

サンプリングされた入射角はかならず一つに定まるのでpdfには1を代入し、返り値のradianceは距離の二乗分だけ減衰するのでlength * lengthで割っています。*9

Integratorの実装

BRDFとライトが用意できたのでようやくIntegratorの実装していきましょう。まず、積分計算時に衝突地点にのみ集中できるようにSurfaceInteractionクラスを作って、衝突判定時にその後必要そうな情報をそこにまとめておきます。

struct SurfaceInteraction {
    std::unique_ptr<BRDF> brdf;
    Vector3f point;
    Vector3f normal;
    Vector3f wout;

    void actedBy(Transform _t) {
        point = _t.act(point);
        normal = _t.act(normal);
        wout = _t.act(wout);
    }
};

Scene::Intersect実行時にこれらの情報を埋めておきます。さらにSceneクラスのコンストラクタでLightクラスを指定できるようにしてください。*10

パストレーシングをする前にまずは直射光だけを考慮したIntegratorを作ってみたいので、IntegratorにestimateDirect/liという関数を追加しましょう。estimateDirect関数は衝突地点を受け取ったらその地点のRadianceを求める関数で、li関数はScene::intersectとestimateDirectをまとめた関数です。
これによってIntegrator::render関数内ではli関数の結果のみを見ればよいことになります。

RGB li(const Ray& _ray, const Scene& _scene) const {
    SurfaceInteraction si;
    if (_scene.intersect(_ray, &si)) {
        auto light = _scene.getLight();
        if (light) return estimateDirect(si, light, _scene);
    }
    return RGB{0.0, 0.0, 0.0};    
}

次にestimateDirectを実装しましょう。基本的な流れは

sampleLiによってライトからの入射角とRGBをサンプリング
ライトとオブジェクト間に遮蔽物が存在するかチェック
反射角と入射角を衝突地点のローカル座標に変換後BRDFの計算
$f_r(p, w_o, w_i) L_i(p, w_i) | \cos \theta$ の計算

となります。

1、2、4に関しては既に説明済みですが、3の反射角と入射角を衝突地点のローカル座標に変換というのだけ説明がまだでした。
衝突地点のローカル座標というのは、その地点での法線を上方向の軸とする座標のことですが、この変換は法線に加え接線が1つ分かっていればできます。
というのも、法線と接線が分かっていれば外積をとることにより正規直交基底を構成でき、これを列方向にならべたものはローカル座標からグローバル座標への変換行列となるからです。*11
正規直交基底を並べた行列は直交行列なので、その逆行列は元の行列を転置したものになります(線形代数)

以上をまとめると次のような感じになります。(他にも変更はしている)

RGB estimateDirect(const SurfaceInteraction& _si, const std::shared_ptr<Light> _light, const Scene& _scene) const {
    RGB ld{0.0, 0.0, 0.0};

    Ray wi;  
    float lightPdf;
    RGB l = _light->sampleLi(_si.point, &wi, &lightPdf);
    if (lightPdf > 0.0f) {
        float t;
        SurfaceInteraction si;
        if (!(_scene.intersect(wi, &si, &t) && t <= 0.99)) {
            auto l_wout = _si.worldToLocal(_si.wout);
            auto l_wi = _si.worldToLocal((-wi.getDir()).normalize());
            RGB f = _si.brdf->f(l_wout, l_wi) * std::abs(l_wi.dot(Vector3f{0.0, 1.0, 0.0}));
            ld = f * l / lightPdf; 
        }
    }

    return ld;
}