古典的レイトレーサー&簡単なパストレーサーの実装(第4回)

今日やること

色を指定できるようにする
マテリアルを指定できるようにする
ミラーの実装
ガンマ補正
はじめてのPath Tracing
光の伝わり方を数式で表す
モンテカルロ積分

色を指定できるようにする

複数の球を描画できるようになったので、球に個別に色を指定できるようにしてみましょう。

Sphere自体が色の情報を持つようし、rayがどのSphereと当たったのか区別できるようになれば色をつけることができそうです。

まずはSphereクラスに色を表すメンバ変数を追加しましょう。

class Sphere {
  public:
    Vec3 center;
    double radius;
    Vec3 color //追加
    ...
}

こうすれば球を作るときに色を指定できるようになります。

Sphere s = Sphere(Vec3(0, 0, 0), 1.0, Vec3(0, 1, 0)); //緑色の球

あとはrayがどの球体と当たったのか区別できるようにする必要があります。一つの方法として、衝突情報(Hitクラス)に当たった球体へのポインタを保持させる方法があります。

class Hit {
  public:
    double t;
    Vec3 hitPos;
    Vec3 hitNormal;
    const Sphere* hitSphere; //追加
}

これを利用するとrayが衝突した物体の色は次のように取得できます。

Hit hit;
accel.intersect(ray, hit)
Vec3 color = hit.hitSphere->color; //Sphereのcolorメンバ変数を参照する

演習

上の内容を参考にして、個別に色を指定できるようにしてみましょう。

マテリアルを指定できるようにする

これも個別に色を指定するのと同じようにできます。Sphereクラスにマテリアルの種類を表すメンバ変数を追加してあげればよいのです。

class Sphere {
  public:
    Vec3 center;
    double radius;
    Vec3 color;
    int material; //追加
}

例えば0だったらDiffuse、1だったらMirrorを表すようにできます。

Sphere s = Sphere(Vec3(0, 0, 0), 1.0, Vec3(0, 1, 0), 0); //緑色のDiffuse球

演習

マテリアルを指定できるようにしてみよう

ミラーの実装

前に宿題でレイが反射した方向を返す関数reflect()を作りましたね。それを利用してミラーを実装してみましょう。

鏡に映るのはrayが反射した先の色です。つまり、鏡に当たったときに新たに反射方向へのrayを生成して衝突計算を行い、その結果の色を画素に書き込むということを行います。

プログラムにすると次のようになります。

Ray ray = cam.getRay(u, v);
Hit hit;
if(accel.intersect(ray, hit)) {
  //Diffuse
  if(hit.hitSphere->material == 0) {
    //Diffuse面に対する陰影計算
    ...
  }
  //Mirror
  if(hit.hitSphere->material == 1) {
    //反射したレイを生成
    Ray nextRay(hit.hitPos + 0.001*hit.hitNormal, reflect(ray.direction, hit.hitNormal));
    //反射した先で陰影計算
    ...
  }
}

色を求めるには、反射した先で再び陰影計算を行う必要があります。

これを行うために陰影計算を行う処理を次のような関数にまとめると便利です。

Vec3 getColor(const Ray& ray, int depth = 0);

getColor()は受け取ったrayの方向から来る光の強さを計算して返してくれる関数です。depthについてはあとで説明します。

今まで書いてきた陰影処理をこのgetColor()の中に移しましょう。

Accel accel; //グローバルに定義する必要がある
Vec3 lightDir = normalize(Vec3(1, 1, -1)); //光源の方向

Vec3 getColor(const Ray& ray, int depth = 0) {
  Hit hit;
  if(accel.intersect(ray, hit )) {
    //Diffuse
    if(hit.hitSphere->material == 0) {
      //Shadow Rayを生成
      Ray shadowRay(hit.hitPos + 0.001*hit.hitNormal, lightDir);
      Hit shadow_hit;
      //Shadow Rayが当たったら影をつける
      if(!accel.intersect(shadowRay, shadow_hit)) {
        double I = std::max(dot(lightDir, hit.hitNormal), 0.0);
        return I * hit.hitSphere->color;
      }
      else {
        return Vec3(0, 0, 0);
      }
    }
    //Mirror
    else if(hit.hitSphere->material == 1) {
      //反射した先で陰影計算
      Ray nextRay(hit.hitPos + 0.001*hit.hitNormal, reflect(ray.direction, hit.hitNormal));
      ...
    }
  }
  else {
    return Vec3(0, 0, 0);
  }
}

では鏡面反射部分を実装していきましょう。反射した先の色を計算させればいいので、反射rayを引数として再びgetColor()を呼び出せば良いです。

else if(hit.hitSphere->material == 1) {
  //反射Rayを生成
  Ray nextRay(hit.hitPos + 0.001*hit.hitNormal, reflect(ray.direction, hit.hitNormal));
  //単純に反射先の色を返すだけ
  return getColor(nextRay, depth + 1);
}

ここでdepthに+1した値を引数として与えます。無限回の反射を計算させようとすると Stack Overflow というエラーが出てしまうので、100回以上の反射は追跡しないようにします。

Vec3 getColor(const Ray& ray, int depth = 0) {
  if(depth > 100) return Vec3(0, 0, 0);
  ...

演習

以上を踏まえてミラーを自分で実装してみよう。

ガンマ補正

コンピューターのディスプレイで表示される色というのは実際に表示されるべき色とはかなり異なる色になっています。そこで表示されるべき色に対して補正をかけてあげることで、ディスプレイに表示される色をより自然な色に近づけることができます。

次のような補正をかけてあげます。

I_{display} = I_{original}^{\frac{1}{2.2}}

これを行うためにImageクラスにgamma_correction()という関数を実装しましょう。

void gamma_correction() {
    for(int i = 0; i < width; i++) {
        for(int j = 0; j < height; j++) {
            Vec3 col = this->getPixel(i, j);
            col.x = std::pow(col.x, 1.0/2.2);
            col.y = std::pow(col.y, 1.0/2.2);
            col.z = std::pow(col.z, 1.0/2.2);
            this->setPixel(i, j, col);
        }
    }
};

画素データ一つ一つに対して補正をかけてあげるだけです。

実際に使うときにはppm_output()を行う前に使うようにします。

img.gamma_correction();
img.ppm_output();

はじめてのPath Tracing

ミラーを実装することで再帰的なRayの追跡を行うことが可能になりました。Diffuse面に対しても再帰的な追跡を行うようにしてみましょう。

Diffuse面では光は半球面の全ての方向に等確率でランダムに光が反射されます。

なのでまずは半球面の全ての方向を等確率でランダムに生成する関数randomHemisphere(const Vec3& n)を作る必要があります。 $\vec{n}$ は上方向を表すベクトルです。

半球面でのランダムな方向の生成

$[0, 1]$ の一様乱数を $u, v$ としたとき、次のように $x, y, z$ を計算すると $(x, y, z)$ は半球面の全方向の一様サンプリングになります。

\begin{aligned} y &= u \\ x &= \sqrt{1 - u^2}\cos{2\pi v} \\ z &= \sqrt{1 - u^2}\sin{2\pi v} \end{aligned}

Vec3 randomHemisphere(const Vec3& n) {
  double u = rnd();
  double v = rnd();
  
  double y = u;
  double x = std::sqrt(1 - u*u)*std::cos(2*M_PI*v);
  double z = std::sqrt(1 - u*u)*std::sin(2*M_PI*v);
  
  return Vec3(x, y, z);
}

これで半球面でのランダムな方向を手に入れることができました。ただこのままでは必ず上を向いた半球でのランダムな方向が返ってくることになります。あらゆる方向を向いた半球に対してランダムな方向を生成できるようにしたいですね。

これを実現するには、受けとった $\vec{n}$ を $y$ 成分の基底として使えば良いです。同様に $x, z$ 成分の基底が求まれば、任意の方向を向いた半球に対してランダムな方向を生成できるようになります。

ここで必要になるのが、

\vec{n}

のみから正規直交基底を生成してくれる関数です。これを行うorthonormalBasis(const Vec3& n, Vec3& x, Vec3& z)を作ってみましょう。

void orthonormalBasis(const Vec3& n, Vec3& x, Vec3& z) {
  if(n.x > 0.9) x = Vec3(0, 1, 0);
  else x = Vec3(1, 0, 0);
  x = x - dot(x, n)*n;
  x = normalize(x);
  z = normalize(cross(n, x));
}

まず最初に適当に横方向のベクトルを決めて、それが $\vec{n}$ と直角となるように $\vec{x}$ をとり、あとは $\vec{n}$ と $\vec{x}$ に直交するように外積で $\vec{z}$ を定義してあげます。

これを使えばrandomHemisphere(const Vec3& n)は次のように書き直せます。

Vec3 randomHemisphere(const Vec3& n) {
  double u = rnd();
  double v = rnd();
  
  double y = u;
  double x = std::sqrt(1 - u*u)*std::cos(2*M_PI*v);
  double z = std::sqrt(1 - u*u)*std::sin(2*M_PI*v);
  
  Vec3 xv, zv;
  orthonormalBasis(n, xv, zv);
  
  return x*xv + y*n + z*zv;
}

Diffuse面に対する再帰的な追跡

あとはgetColor()のDiffuseに対する実装を変えるだけです。反射した先の色を光源と考えて以前と同様の計算を行えば良いです。

Vec3 getColor(const Ray& ray, int depth = 0) {
  //100回以上の反射は追跡しない
  if(depth > 100) return Vec3(0, 0, 0);
  
  Hit hit;
  if(accel.intersect(ray, hit )) {
    //Diffuse
    if(hit.hitSphere->material == 0) {
      Ray nextRay(hit.hitPos + 0.001*hit.hitNormal, randomHemisphere(hit.hitNormal));
      double cos_term = std::max(dot(nextRay.direction, hit.hitNormal), 0.0);
      return cos_term * hit.hitSphere->color * getColor(nextRay, depth + 1);
    }
    ...

反射した先の色を掛け合わせているので、緑の壁の色が青の床にも染み渡るようになります。白い壁に蛍光色のものを近づけるとその色が壁に映り込むのと同じです。

Ray空に飛んでいった場合には白色を返すようにしてみましょう。この場合、空全体が光源になります。

if(accel.intersect(ray, hit)) {
  ...
}
else {
  return Vec3(1, 1, 1);
}

リアルな見た目になってきましたね。下の床に注目すると若干緑が映り込んでいるのが分かります。

演習

Diffuseも再帰的な追跡を行うようにしてみよう。

光の伝わり方を数式で表す

$L_i(x, \vec{\omega})$ が点 $x$ に $\vec{\omega}$ から入ってくる光の強さを表し、 $L_o(x, \vec{\omega})$ が点 $x$ から $\vec{\omega}$ に出ていく光の強さを表すとします。

$f(x, \vec{\omega_i}, \vec{\omega_o})$ は物体の反射特性を表す関数で、 $\vec{\omega_i}$ から入射した光がどれだけ $\vec{\omega_o}$ に反射されるかを表します。この関数を 双方向反射率分布関数(BRDF) といいます。

$\vec{n}$ を物体の法線とすると、ある特定の方向 $\vec{\omega_i}$ から入射してくる光が $\vec{\omega_o}$ にどれだけ反射されるかは

L_o(x, \vec{\omega_o}) = f(x, \vec{\omega_i}, \vec{\omega_o})L_i(x, \vec{\omega_i})(\vec{\omega_i} \cdot \vec{n})

で計算することができます。

入射光が半球全体 $\Omega$ から入ってくるとすると、 $\vec{\omega_o}$ に反射される光は積分を取ることで

L_o(x, \vec{\omega_o}) = \int_{\Omega} f(x, \vec{\omega_i}, \vec{\omega_o})L_i(x, \vec{\omega_i})(\vec{\omega_i} \cdot \vec{n}) d\vec{\omega_i}

で計算することができます。

上のプログラムではgetColor()が $L_i(x, \vec{\omega_i})$ に対応します。この積分は解析的に評価することができないので、コンピューターを用いて数値積分してあげる必要があるのですが、その一手法として、乱数を用いて積分を評価するモンテカルロ積分という方法があります。

モンテカルロ積分

下のような積分をコンピューターを用いて数値積分することを考えてみましょう。

I = \int_0^1 f(x)dx

一般的には $x$ 軸を細かい区間に分けて関数値を足し合わせる 区分求積法 が有名ですが、 $x$ の値として乱数を用いて関数値を足し合わせる モンテカルロ積分(Monte Carlo Integration) と呼ばれる方法があります。

0から1の範囲にある乱数列 $\{ x_1, x_2, \dots, x_N \}$ を確率密度関数 $p(x)$ を持つある分布から発生させるとすると

I \approx \frac{1}{N}\sum_{i = 1}^N \frac{f(x_i)}{p(x_i)}

で計算することができます。

モンテカルロ積分の利点は高次元積分の場合でも乱数を用いて評価する点を生成できることです。区分求積法などの方法では軸を等間隔に分割する必要がありますが、高次元だとそもそも軸を分割する計算に時間がかかり、次元の呪いという問題も発生します。

モンテカルロ積分の収束のオーダーは $O(\frac{1}{\sqrt{N}})$ です。つまり、精度を10倍にするには100倍のサンプリング数が必要になります。収束が遅いのがモンテカルロ積分の難点です。

収束スピードを改善する方法として 準モンテカルロ法(Quasi Monte Carlo Method) や 重点的サンプリング(Importance Sampling) といった方法があります。これについては次回以降見ていきます。

例1

0から1の範囲の一様分布から乱数を生成するとすると $p(x) = 1$ で

I \approx \frac{1}{N}\sum_{i = 1}^N f(x_i)

で計算することができます。単純に乱数で評価した値の平均をとっているだけですね。

例2

$[0, 2]\times[0, 2]$ の正方形から一様乱数を作ると $p(x) = \frac{1}{2\times 2} = \frac{1}{4}$ で

I \approx \frac{4}{N}\sum_{i = 1}^N f(x_i)

となります。前回の宿題はこれを行っていました。

モンテカルロ積分による反射光の評価

L_o(x, \vec{\omega_o}) = \int_{\Omega} f(x, \vec{\omega_i}, \vec{\omega_o})L_i(x, \vec{\omega_i})(\vec{\omega_i} \cdot \vec{n}) d\vec{\omega_i}

をモンテカルロ積分を使って計算してみましょう。まずは半球面の方向を一様に発生させ、これを $\{ \vec{\omega_1}, \vec{\omega_2}, \dots, \vec{\omega_N} \}$ とします。確率密度関数は $p(\vec{\omega}) = \frac{1}{2\pi}$ です。(半球の表面積は $2\pi$ なので)

すると

L_o(x, \vec{\omega_o}) \approx \frac{1}{N}\sum_{i = 1}^N \frac{f(x, \vec{\omega_i}, \vec{\omega_o})L(x, \vec{\omega_i})(\vec{\omega_i}\cdot \vec{n})}{p(\vec{\omega_i})} = \frac{1}{N}\sum_{i = 1}^N 2\pi f(x, \vec{\omega_i}, \vec{\omega_o})L(x, \vec{\omega_i})(\vec{\omega_i}\cdot \vec{n})

となります。

実際に計算してみる

上のプログラムを修正して、Diffuse面の計算をモンテカルロ積分に置き換えてみましょう。

Diffuse面の場合は $f(x, \vec{\omega_i}, \vec{\omega_o}) = \frac{\rho}{\pi}$ です。ここで $\rho$ は反射率です。プログラム中では物体の色に対応します。

変更するのはreturnの部分だけです。

//Diffuse
if(hit.hitSphere->material == 0) {
  Ray nextRay(hit.hitPos + 0.001*hit.hitNormal, randomHemisphere(hit.hitNormal));
  double cos_term = std::max(dot(nextRay.direction, hit.hitNormal), 0.0);
  return 2*M_PI * hit.hitSphere->color/M_PI * cos_term * getColor(nextRay, depth + 1);
}

先ほどより明るい見た目になりましたね。これは上記の積分を正しく近似しているので、サンプル数を増やすと物理的に正しい光の強さに収束します。

下をミラーにするとこんな感じです。

とても綺麗ですね!!

練習問題

問1

モンテカルロ積分を用いて以下の積分を評価せよ。

\int_D \frac{xy}{(x + y + 1)^2} dxdy

$D$ は $[0, 1]\times[0, 1]$ の正方形とする。

問2

モンテカルロ積分を用いて以下の積分を評価せよ。

\int_D f(x, y, z)dxdydz

$f(x, y, z) = \begin{cases} 1 & (x^2 + y^2 + z^2 < 1) \\ 0 & (otherwise) \end{cases}$ とし、 $D$ は $[-1, 1]\times[-1, 1]\times[-1, 1]$ の立方体とする。

問3

球をたくさん追加したり、色を変えたりして面白いレンダリング画像を作成せよ。