Help us understand the problem. What is going on with this article?

レイトレーシング入門1「光線の基本と反射」

More than 1 year has passed since last update.

はじめに

コンピュータグラフィックをレンダリングするとき,リアルタイムレンダリングかどうかでレンダリングする手法は基本的に変わります.非リアルタイムレンダリングではレイトレーシングが現在でも主流で,様々な手法が開発されたり改良されています.レイトレーシング法の中でもモンテカルロレイトレーシングはランダムにレイを反射させて光伝達を追跡し,それ以外の光の物理現象を正確に表現することができます.ただし,サンプリング数が多くないと計算結果の精度が良くないという問題を抱えています.それについても様々な改良方法が研究されています.また,リアルタイムにおいてもレイマーチングという手法でレンダリングが行われることも活発になってきており,ゲームなどのリアルタイムレンダリングでも活用する機会が増えてきています.

本記事では,モンテカルロレイトレーシングについて解説します.内容としては Peter Shirley 著の3冊

  • Ray Tracing in One Weekend
  • Ray Tracing: the Next Week
  • Ray Tracing: The Rest Of Your Life

を基本としています.ただし,全ての内容は含まれておらず,またコードの設計などには若干手を入れています.関与媒質,BVH,ノイズ,モーションブラー,被写界深度といった内容はここでは扱いませんので,興味があれば読んでみてください.Kindle版のみですが,Kindle Unlimited なら無料,また,どれも 350 円程度です.

すでに多くのレイトレーシングレンダラーのソースコードが公開されていたり,資料が公開されています.最後にいくつか紹介しますので,そちらも興味があれば参照してみてください.

コードは C++ で書いています.STL や c++11(スマートポインタ)など使いますので,ある程度言語について知っていることが前提となります.といってもそんなに多くはないので,わからない機能がありましたら,ネットで検索すれば情報は出てくるはずです.名前空間は特別なものを除いて明示的に記載します.

今回の目標は次のような画像をレンダリングすることです.

tuto-raytracing-shape-list-pdf-output.png

ソースコード

本記事のソースコードは以下の場所にあります.

rayt

Windowsで,VisualStudio 2015 Community を使用しています.また,プログラミング言語は C++ です.

rayt.cpp ファイルがメインファイルになっています.ソースコードには rayt101.cpp という様に連番がついたファイルがあります.これはその内容を全て rayt.cpp に上書きするようになっていますので,試すときはそのようにしてください.レンダリング画像を掲示している場所には,その画像を生成したソースコードファイル名を一緒に明記するようにしています.

画像出力

まずは,レンダリングした画像をファイルに保存できるような機能を作ります.この機能がないと何も始まりません.

コードは rayt.hrayt.cpp に書いていきます.通常はクラスごとにファイルを分割した方がよいのですが,今回はこの2ファイルに集約しました.ユーティリティ関数や一度実装した後に変更のないクラスは rayt.h に記載するようにします.本文中でコードのところに // rayt.h と記載があれば rayt.h に実装し,記載がない場合は rayt.cpp に実装することを表しています.あと,基本的に名前空間 rayt 内に定義するようにしています.

まずは Image クラスを作成します.これはスクリーンのカラー値を格納する領域を管理し,そこに書き込む機能を持っています.

// rayt.h
#include <memory>
#include <iostream>

namespace rayt {

    class Image {
    public:
        struct rgb {
            unsigned char r;
            unsigned char g;
            unsigned char b;
        };

        Image() : m_pixels(nullptr) { }
        Image(int w, int h) {
            m_width = w;
            m_height = h;
            m_pixels.reset(new rgb[m_width*m_height]);
        }

        int width() const { return m_width; }
        int height() const { return m_height; }
        void* pixels() const { return m_pixels.get(); }

        void write(int x, int y, float r, float g, float b) {
            int index = m_width*y + x;
            m_pixels[index].r = static_cast<unsigned char>(r*255.99f);
            m_pixels[index].g = static_cast<unsigned char>(g*255.99f);
            m_pixels[index].b = static_cast<unsigned char>(b*255.99f);
        }

    private:
        int m_width;
        int m_height;
        std::unique_ptr<rgb[]> m_pixels;
    };
} // namespace rayt

スマートポインタ std::unique_ptr を使用するためには memory ヘッダーファイルをインクルードする必要があります.また,コンソール出力を使いますので,iostream もインクルードします.

次にメイン関数を実装します.今回は適当にグラデーションの画像を作成します.
Image クラスのインスタンスを作成して,書き込みます.

#include "rayt.h"

int main()
{
    int nx = 200;
    int ny = 100;
    std::unique_ptr<rayt::Image> image(std::make_unique<rayt::Image>(nx, ny));

    for (int j = 0; j<ny; ++j) {
        std::cerr << "Rendering (y = " << j << ") " << (100.0 * j / (ny - 1)) << "%" << std::endl;
        for (int i = 0; i<nx; ++i) {
            float r = float(i) / float(nx);
            float g = float(j) / float(ny);
            float b = 0.5f;
            image->write(i, j, r, g, b);
        }
    }

    return 0;
}

次に,Image クラスに書き込んだカラー値を画像に出力します.
出力には stb ライブラリを使用します.stb ライブラリはパブリックドメインライブラリです.とても手軽に扱えて便利です.

http://nothings.org/stb/

この中の stb_image.hstb_image_write.h を使います.

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

STB_IMAGE_IMPLEMENTATIONSTB_IMAGE_WRITE_IMPLEMENTATIONstb_image.hstb_image_write.h の前で定義する必要があります.

実際にファイルに出力するには以下のコードを使います.

stbi_write_bmp("render.bmp", nx, ny, sizeof(rayt::Image::rgb), image->pixels());

これはビットマップ形式で render.bmp というファイル名で出力します.これを実行すると以下の画像が得られます.(rayt101.cpp, rayt101.h)

tuto-raytracing-output.png

レイトレーシングの原理

レイトレーシングはまずカメラから光線を飛ばします.飛んでいる光線が何かに当たったら光線は反射します.反射する方向は当たった物体の表面の材質によって変わります.光線は何度か反射を繰り返し,最終的に光源に辿り着きます.この光源から放射された光をカメラまで輸送するシミュレーションを行い,人間の眼で正しく認識できるように変換して画像ファイルとして出力します.

これから実装するために,いくつか定義やユーティリティ関数,ベクトルクラス,基本クラスなど追加していきます.

Vectormath

ベクトル計算に Vectormath ライブラリを使用します.

https://github.com/erwincoumans/sce_vectormath

Vectormath ライブラリは Sony が開発したライブラリでソースが公開され,物理エンジン Bullet で使われています.3Dグラフィックス用のベクトル計算ライブラリということもあって,ベクトルや行列,クォータニオンを扱う十分な機能が含まれています.

名前空間 Vectormath に定義されていて長いので,using マクロで宣言し,よく使うベクトルクラスは typedef で短い名前で再定義しています.

インクルードするディレクトリに vectormath/include を追加する必要があります.

// rayt.h
#include <vectormath/scalar/cpp/vectormath_aos.h>
using namespace Vectormath::Aos;
typedef Vector3 vec3;
typedef Vector3 col3;

C++ のクラスなので,四則演算等は演算子が定義されているので直感的に操作できます.ただし,内積や外積などはメンバ関数ではなく,グローバルな関数として用意されています.よく使う関数を解説します.

関数名 説明
float dot(a,b) 内積を求めます 
vec3 cross(a,b) 外積を求めます
float length(v) ベクトルの長さを求めます
float lengthSqr(v) ベクトルの長さの二乗を求めます
vec3 normalize(v) ベクトルを正規化します
vec3 mulPerElem(a,b) ベクトルの各要素ごとに乗算します
vec3 maxPerElem(a,b) 2つのベクトルから各要素ごとに最大であるベクトルを返します
vec3 minPerElem(a,b) 2つのベクトルから各要素ごとに最小であるベクトルを返します
vec3 lerp(t,a,b) 2つのベクトルをパラメータ t で線形補間したベクトルを返します

数学の定数

数学の定数をいくつか定義しています.

// rayt.h
#include <float.h>  // FLT_MIN, FLT_MAX
#define PI 3.14159265359f
#define PI2 6.28318530718f
#define RECIP_PI 0.31830988618f
#define RECIP_PI2 0.15915494f
#define LOG2 1.442695f
#define EPSILON 1e-6f
#define GAMMA_FACTOR 2.2f

ユーティリティ関数

また,ユーティリティ関数をいくつか定義しています.

// rayt.h
#include <random>
inline float drand48() {
    return float(((double)(rand()) / (RAND_MAX))); /* RAND_MAX = 32767 */
}

inline float pow2(float x) { return x*x; }
inline float pow3(float x) { return x*x*x; }
inline float pow4(float x) { return x*x*x*x; }
inline float pow5(float x) { return x*x*x*x*x; }
inline float clamp(float x, float a, float b) { return x < a ? a : x > b ? b : x; }
inline float saturate(float x) { return x < 0.f ? 0.f : x > 1.f ? 1.f : x; }
inline float recip(float x) { return 1.f / x; }
inline float mix(float a, float b, float t) { return a*(1.f - t) + b*t; /* return a + (b-a) * t; */ }
inline float step(float edge, float x) { return (x < edge) ? 0.f : 1.f; }
inline float smoothstep(float a, float b, float t) { if (a >= b) return 0.f; float x = saturate((t - a) / (b - a)); return x*x*(3.f - 2.f * t); }
inline float radians(float deg) { return (deg / 180.f)*PI; }
inline float degrees(float rad) { return (rad / PI) * 180.f; }

以下にリファレンスを記載します.

関数名 説明
float drand48() [0,1]の一様乱数を返します 
float pow2(x) xの自乗を求めます
float pow3(x) xの3乗を求めます
float pow4(x) xの4乗を求めます
float pow5(x) xの5乗を求めます
float clamp(x,a,b) $a \leq x \leq b$ を満たす x を返します
float saturate(x) $0 \leq x \leq 1$ を満たす x を返します
float recip(x) $\frac{1}{x}$ を返します
float mix(a,b,t) aとbをパラメータ t で線形補間した値を返します
float step(edge,x) edge 以下のとき 0 を,それ以外のときは 1 を返します
float smoothstep(a,b,t) aとbをパラメータ t でスプライン補間した値を返します
float radians(deg) 度数からラジアンに変換します
float degrees(rad) ラジアンから度数に変換します

基本クラス

光線

光線は始点と向きを持っています.コードは次の通りです.

// rayt.h
class Ray {
    public:
        Ray() {}
        Ray(const vec3& o, const vec3& dir)
            : m_origin(o)
            , m_direction(dir) {
        }

        const vec3& origin() const { return m_origin; }
        const vec3& direction() const { return m_direction; }
        vec3 at(float t) const { return m_origin + t*m_direction; }

    private:
        vec3 m_origin;    // 始点
        vec3 m_direction; // 方向(非正規化)
};

at 関数は,始点 $\vec{o}$ から方向 $\vec{d}$ に向かう直線上の任意の点 $\vec{p}$ を,パラメータ $t$ を指定して求めます.$t=0$ のときは始点,$t=1$ のときは始点から方向ベクトル $\vec{d}$ 進んだ位置,$t>1$ なら方向ベクトルより先の位置,$t < 0$ なら,始点から反対方向に向かう直線上の位置を取得できます.式で表すと

p(t) = \vec{o} + t\vec{d}.

tuto-raytracing-ray.png

カメラ

レイトレーシングでは投影するスクリーン上のピクセルごとに光線を飛ばします.光線はカメラの位置から発生し,カメラの向きに飛んでいきます.投影するスクリーンの大きさを $sx,sy$ とし,各ピクセルへの光線はスクリーン上の位置 $u,v$ から算出します.

tuto-raytracing-screen.png

カメラは投影するスクリーンの $X$ 軸と $Y$ 軸,そしてカメラの向きを $Z$ 軸とする直交基底ベクトルを持っています.この基底ベクトルを使って $u,v$ を基底変換します.

光線を飛ばすピクセルの位置が $si,sj$ とすると,

u = \frac{si}{sx}, \quad v = \frac{sj}{sy} \qquad (0 \leq u \leq 1, 0 \leq v \leq 1).

基底ベクトル $\vec{u},\vec{v},\vec{w}$ を使ってこれを基底変換すると,投影するスクリーン上の位置 $\vec{p}$ は

\vec{p} = \vec{u}\cdot u + \vec{v}\cdot v + \vec{w}.

カメラの位置を原点に置き,-z 方向を向いているとします.スクリーンの右上が $(2,1,-1)$ となるような基底ベクトルは次のようになります.

\vec{u} = (4,0,0), \quad \vec{v} = (0,2,0), \quad \vec{w} = (-2,-1,-1).

下図はこの関係を表したものです.

tuto-raytracing-camera.png

カメラの指定には便利な LookAt 方式があります.カメラの位置 lookfrom と視線対象の位置 lookat,カメラの上方向を表す vup から基底ベクトルを計算することができます.

tuto-raytracing-camera-lookat.png

まず,カメラにおける正規直交基底ベクトル $XYZ$ を計算します.$Z$ は $lookat - lookfrom$ です.次に $X$ は上方向を表す vup と求めた $Z$ の外積です.そして,$Y$ は $X$ と $Z$ の外積となります.

この正規直交基底ベクトルから,求めたい基底ベクトルを導出します.まず,スクリーン上の高さ $h$ はカメラからスクリーンまでの距離を 1 とすると $h = \tan(\theta/2)$ です.幅 $w$ はスクリーンのアスペクト比($aspect = sx/sh$)を使って $w = aspect \cdot h$ です.

基底ベクトル $\vec{u}$ と $\vec{v}$ はそれぞれ正規直交ベクトルの $XY$ に対応しているので

\vec{u} = 2wX, \quad \vec{v} = 2hY.

次に $\vec{w}$ です.スクリーン上の位置 $\vec{p}$ を基底ベクトルから求めるときは以下の式を使いました.

\vec{p} = \vec{u}\cdot u + \vec{v}\cdot v + \vec{w}.

この式から $\vec{w}$ を求めると

\vec{w} = \vec{p} - \vec{u}\cdot u - \vec{v}\cdot v.

今,$\vec{p} = \vec{o} - Z$($\vec{o}$ はカメラの位置)なので

\vec{w} = \vec{o} - wX - hY - Z.

getRay 関数では,$u,v$ から光線を生成します.基底変換で得られた位置からカメラの位置との差が方向ベクトルとなります.

以下がカメラクラスのコードです.コンストラクタは基底ベクトルを直接指定するのと,LookAt 方式で指定するものがあります.

// rayt.h
class Camera {
public:
    Camera() {}
    Camera(const vec3& u, const vec3& v, const vec3& w) {
        m_origin = vec3(0);
        m_uvw[0] = u;
        m_uvw[1] = v;
        m_uvw[2] = w;
    }
    Camera(const vec3& lookfrom, const vec3& lookat, const vec3& vup, float vfov, float aspect) {
        vec3 u, v, w;
        float halfH = tanf(radians(vfov)/2.0f);
        float halfW = aspect * halfH;
        m_origin = lookfrom;
        w = normalize(lookfrom - lookat);
        u = normalize(cross(vup, w));
        v = cross(w, u);
        m_uvw[2] = m_origin - halfW*u - halfH*v - w;
        m_uvw[0] = 2.0f * halfW * u;
        m_uvw[1] = 2.0f * halfH * v;
    }

    Ray getRay(float u, float v) const {
        return Ray(m_origin, m_uvw[2] + m_uvw[0]*u + m_uvw[1]*v - m_origin);
    }

private:
    vec3 m_origin;  // 位置
    vec3 m_uvw[3];  // 直交基底ベクトル
};

このカメラを使ってみます.カメラの基底ベクトルを指定し,各ピクセルに向かう光線を作成します.その光線の方向ベクトルを正規化して,そのベクトルの $y$ を使って2色のグラデーションを描きます.

vec3 color(const rayt::Ray& r) {
    vec3 d = normalize(r.direction());
    float t = 0.5f*(r.direction().getY() + 1.0f);
    return lerp(t, vec3(0.5f, 0.7f, 1.0f), vec3(1));
}

int main()
{
    int nx = 200;
    int ny = 100;
    std::unique_ptr<rayt::Image> image(std::make_unique<rayt::Image>(nx, ny));

    vec3 x(4.0f, 0.0f, 0.0f);
    vec3 y(0.0f, 2.0f, 0.0f);
    vec3 z(-2.0f, -1.0f, -1.0f);
    std::unique_ptr<rayt::Camera> camera(std::make_unique<rayt::Camera>(x, y, z));

    for (int j = 0; j<ny; ++j) {
        std::cerr << "Rendering (y = " << j << ") " << (100.0 * j / (ny - 1)) << "%" << std::endl;
        for (int i = 0; i<nx; ++i) {
            float u = float(i) / float(nx);
            float v = float(j) / float(ny);
            rayt::Ray r = camera->getRay(u, v);
            vec3 c = color(r);
            image->write(i, j, c.getX(), c.getY(), c.getZ());
        }
    }

    stbi_write_bmp("render.bmp", nx, ny, sizeof(rayt::Image::rgb), image->pixels());
    return 0;
}

実行すると次のような画像になります.(rayt102.cpp, rayt102.h)

tuto-raytracing-camera-output.png

球の追加

これでは寂しいので物体を追加しましょう.ここでは「球」を追加します.球はとても扱いやすい形状なので,いろんなところで使われます.

球は中心を原点とすると半径 $r$ を使って

x^2+y^2+z^2 = r^2,

という式で表されます.中心位置を $cx,cy,cz$ とすると

(x-cx)^2+(y-cy)^2+(z-cz)^2 = r^2.

ベクトルを使うと,中心位置が $\vec{c}$ で位置が $\vec{p}$ とすると

(\vec{p}-\vec{c})\cdot(\vec{p}-\vec{c}) = r^2.

レイトレーシングでは光線と物体の衝突した位置が知りたいので,光線の方程式

\vec{p}(t) = \vec{o} + t\vec{d},

を代入すると

(\vec{p}(t)-\vec{c})\cdot (\vec{p}(t)-\vec{c}) = r^2.

つまり

(\vec{o} + t\vec{d} - \vec{c})\cdot (\vec{o} + t\vec{d} - \vec{c}) = r^2.

ここで $oc = \vec{o} - \vec{c}$ とすると

(oc + t\vec{d})\cdot (oc + t\vec{d}) - r^2 = 0.

展開すると

(\vec{d}\cdot\vec{d})t^2 + 2(\vec{d}\cdot oc)t + (oc\cdot oc) - r^2 = 0.

ここで「2次方程式の解の公式」を使います.
2次方程式 $ax^2 + bx + c = 0$ の解は

x = \frac{-b\pm\sqrt{b^2-4ac}}{2a}.

$oc$ を展開して,$ax^2 + bx + c = 0$ に当てはめてみると

\begin{align}
a &= (\vec{d}\cdot\vec{d}) \\[2ex]
b &= 2(\vec{d}\cdot (\vec{o} - \vec{c})) \\[2ex]
c &= ((\vec{o} - \vec{c})\cdot (\vec{o} - \vec{c})) - r^2,
\end{align}

となります.ここで,$D = b^2-4ac$ を判別式といい,解がいくつあるかがわかります.$D>0$ なら,2つの解が存在します.$D=0$ なら1つの解が存在し,$x = -\frac{b}{2a}$ で与えられます.残りの $D < 0$ では解は存在しません.

これは図で表すと

tuto-raytracing-sphere.png

光線と球が衝突したかどうかは判別式のみで十分です.それでは球をレンダリングしてみましょう.hit_sphere 関数を追加し,球との衝突処理を実装しています.そして color 関数に球との衝突コードを追加しています.

bool hit_sphere(const vec3& center, float radius, const rayt::Ray& r) {
    vec3 oc = r.origin() - center;
    float a = dot(r.direction(), r.direction());
    float b = 2.0f * dot(r.direction(), oc);
    float c = dot(oc,oc) - pow2(radius);
    float D = b*b-4*a*c;
    return (D > 0);
}

vec3 color(const rayt::Ray& r) {
    if (hit_sphere(vec3(0,0,-1), 0.5f, r)) {
        return vec3(1.0f, 0.0f, 0.0f);
    }
    vec3 d = normalize(r.direction());
    float t = 0.5f*(r.direction().getY() + 1.0f);
    return lerp(t, vec3(0.5f, 0.7f, 1.0f), vec3(1));
}

これは次のような画像になります.(rayt103.cpp)

tuto-raytracing-sphere-output.png

球体の法線

光を反射させるには衝突した位置の法線が必要です.法線は衝突した位置と球体の中心位置の差で求まります.

tuto-raytracing-sphere-normal.png

衝突したら,位置を算出する必要があるので解を求めます.正規化した法線は各要素が[-1,1]の範囲になるため,これをRGBの範囲[0-1]に変換するには 1 を足した後に 0.5 を掛けます.

hit_spherecolor を次のように書き換えます.

float hit_sphere(const vec3& center, float radius, const rayt::Ray& r) {
    vec3 oc = r.origin() - center;
    float a = dot(r.direction(), r.direction());
    float b = 2.0f * dot(r.direction(), oc);
    float c = dot(oc, oc) - pow2(radius);
    float D = b*b - 4 * a*c;
    if (D < 0) {
        return -1.0f;
    }
    else {
        return (-b - sqrtf(D)) / (2.0f*a);
    }
}

vec3 color(const rayt::Ray& r) {
    vec3 c(0,0,-1);
    float t = hit_sphere(c, 0.5f, r);
    if (t > 0.0f) {
        vec3 N = normalize(r.at(t) - c);
        return 0.5f*(N+vec3(1.0f));
    }
    vec3 d = normalize(r.direction());
    t = 0.5f*(r.direction().getY() + 1.0f);
    return lerp(t, vec3(1), vec3(0.5f, 0.7f, 1.0f));
}

これは次のような画像になります.(rayt104.cpp)

tuto-raytracing-sphere-normal-output.png

シーンクラス

ここで,少しコードを整理するためにシーンクラスを追加します.シーンはカメラや画像出力,球などの物体を管理します.これまで作成した hit_spherecolor 関数,カメラや画像クラスのインスタンスをシーンクラスにまとめます.また,背景色としてグラデーションのものは空っぽいので backgroundSky 関数,単色用に background 関数を追加し,color 関数から呼び出しています.好きな方を使ってください.render 関数はレンダリングを行う関数で,build 関数はレンダリングするときに一度呼ばれる関数でカメラや物体の生成などの初期化を行います.

#include "rayt.h"

namespace rayt {

    class Scene {
    public:
        Scene(int width, int height)
            : m_image(std::make_unique<Image>(width, height))
            , m_backColor(0.2f) {
        }

        void build() {

            // Camera

            vec3 w(-2.0f, -1.0f, -1.0f);
            vec3 u(4.0f, 0.0f, 0.0f);
            vec3 v(0.0f, 2.0f, 0.0f);
            m_camera = std::make_unique<Camera>(u, v, w);
        }

        float hit_sphere(const vec3& center, float radius, const rayt::Ray& r) const {
            vec3 oc = r.origin() - center;
            float a = dot(r.direction(), r.direction());
            float b = 2.0f * dot(r.direction(), oc);
            float c = dot(oc, oc) - pow2(radius);
            float D = b*b - 4 * a*c;
            if (D < 0) {
                return -1.0f;
            }
            else {
                return (-b - sqrtf(D)) / (2.0f*a);
            }
        }

        vec3 color(const rayt::Ray& r) {
            vec3 c(0, 0, -1);
            float t = hit_sphere(c, 0.5f, r);
            if (t > 0.0f) {
                vec3 N = normalize(r.at(t) - c);
                return 0.5f*(N + vec3(1.0f));
            }
            return backgroundSky(r.direction());
        }

        vec3 background(const vec3& d) const {
            return m_backColor;
        }

        vec3 backgroundSky(const vec3& d) const {
            vec3 v = normalize(d);
            float t = 0.5f * (v.getY() + 1.0f);
            return lerp(t, vec3(1), vec3(0.5f, 0.7f, 1.0f));
        }

        void render() {

            build();

            int nx = m_image->width();
            int ny = m_image->height();
            for (int j = 0; j<ny; ++j) {
                std::cerr << "Rendering (y = " << j << ") " << (100.0 * j / (ny - 1)) << "%" << std::endl;
                for (int i = 0; i<nx; ++i) {

                    float u = float(i + drand48()) / float(nx);
                    float v = float(j + drand48()) / float(ny);
                    Ray r = m_camera->getRay(u, v);
                    vec3 c = color(r);
                    m_image->write(i, (ny - j - 1), c.getX(), c.getY(), c.getZ());
                }
            }

            stbi_write_bmp("render.bmp", nx, ny, sizeof(Image::rgb), m_image->pixels());
        }

    private:
        std::unique_ptr<Camera> m_camera;
        std::unique_ptr<Image> m_image;
        vec3 m_backColor;
    };

} // namespace rayt

main 関数はシーンクラスのインスタンスを作成し,render 関数を呼びます.(rayt104.cpp)

int main()
{
    int nx = 200;
    int ny = 100;
    std::unique_ptr<rayt::Scene> scene(std::make_unique<rayt::Scene>(nx, ny));
    scene->render();

    return 0;
}

並列処理 - OpenMP

レイトレーシングは並列処理に向いています.今回の実装では各光線ごとに計算を並列で処理させることが可能です.Visual Studio には並列処理のフレームワークである OpenMP が使えるようになっています.OpenMP はコード中にディレクティブを挿入することで簡単に並列処理を実現することができます.render 関数に以下のディレクティブを入れます.

#pragma omp parallel for schedule(dynamic, 1) num_threads(NUM_THREAD)

実際のコードは以下のようになります.

// rayt.cpp に記載
class Scene {
    ...

    void render() {

        build();

        int nx = m_image->width();
        int ny = m_image->height();
#pragma omp parallel for schedule(dynamic, 1) num_threads(NUM_THREAD)
        for (int j = 0; j<ny; ++j) {
            std::cerr << "Rendering (y = " << j << ") " << (100.0 * j / (ny - 1)) << "%" << std::endl;
            for (int i = 0; i<nx; ++i) {

                float u = float(i + drand48()) / float(nx);
                float v = float(j + drand48()) / float(ny);
                Ray r = m_camera->getRay(u, v);
                vec3 c = color(r);
                m_image->write(i, (ny - j - 1), c.getX(), c.getY(), c.getZ());
            }
        }

        stbi_write_bmp("render.bmp", nx, ny, sizeof(Image::rgb), m_image->pixels());
    }
}

ディレクティブについての詳しい内容はネットで検索してみてください.NUM_THREAD は使用するスレッドの数です.お使いのPCに合わせたスレッド数を指定します.

// rayt.h に記載
#define NUM_THREAD 8 // 例

このままでは,OpenMP のディレクティブは有効になっていません.Visual Studio のプロジェクト設定で OpenMP を有効にする必要があります.

tuto-raytracing-openmp.png

あとはビルドして実行するだけで並列処理されます.また,OpenMP は無効になっていたらディレクティブは無視されるのでコードを変更する必要がありません.

もしレンダリングした画像がおかしくなったとき,並列処理が原因となっている場合があります.そのときは一度並列処理を疑ってみてください.例えば schedule(dynamic,1) を消してみたり,並列処理そのものを無効にするなどで確認します.

(rayt105.cpp, rayt105.h)

複数の物体への対応

今までは球が1つでしたが,複数の物体にも対応できるように拡張します.

衝突情報

まず,物体に衝突したときの情報を格納するクラス HitRec を用意します.

class HitRec {
public:
    float t;
    vec3 p;
    vec3 n;
};

t は光線のパラメータ,p は衝突した位置,n は衝突した位置の法線です.

物体の抽象クラス

次に,物体の抽象クラス Shape を追加します.これは衝突関数が仮想関数で定義されています.

class Shape {
public:
    virtual bool hit(const Ray& r, float t0, float t1, HitRec& hrec) const = 0;
};

t0t1 は光線の衝突範囲です.

球クラス

そして,球クラス Sphere を追加します.

class Sphere : public Shape {
public:
    Sphere() {}
    Sphere(const vec3& c, float r)
        : m_center(c)
        , m_radius(r) {
    }

    virtual bool hit(const Ray& r, float t0, float t1, HitRec& hrec) const override {
        vec3 oc = r.origin() - m_center;
        float a = dot(r.direction(), r.direction());
        float b = 2.0f*dot(oc, r.direction());
        float c = dot(oc, oc) - pow2(m_radius);
        float D = b*b - 4*a*c;
        if (D > 0) {
            float root = sqrtf(D);
            float temp = (-b - root) / (2.0f*a);
            if (temp < t1 && temp > t0) {
                hrec.t = temp;
                hrec.p = r.at(hrec.t);
                hrec.n = (hrec.p - m_center) / m_radius;
                return true;
            }
            temp = (-b + root) / (2.0f*a);
            if (temp < t1 && temp > t0) {
                hrec.t = temp;
                hrec.p = r.at(hrec.t);
                hrec.n = (hrec.p - m_center) / m_radius;
                return true;
            }
        }

        return false;
    }

private:
    vec3 m_center;
    float m_radius;
};

override は仮想関数をオーバーライドすることをコンパイラに知らせるキーワードです.引数の間違いや誤字などでオーバーライドが出来ていないことを防げます.

hit 関数での衝突判定は少し変更しています.解が2つあるとき,$t = \frac{-b-\sqrt{D}}{2a}$ の方が始点に近いため,近いほうを先に判定しています.

物体リスト

複数の物体を管理するクラスを追加します.これは Shape を継承することで,物体リストも単体の物体として動作するようになります.

class ShapeList : public Shape {
public:
    ShapeList() {}

    void add(const ShapePtr& shape) {
        m_list.push_back(shape);
    }

    virtual bool hit(const Ray& r, float t0, float t1, HitRec& hrec) const override {
        HitRec temp_rec;
        bool hit_anything = false;
        float closest_so_far = t1;
        for (auto& p : m_list) {
            if (p->hit(r, t0, closest_so_far, temp_rec)) {
                hit_anything = true;
                closest_so_far = temp_rec.t;
                hrec = temp_rec;
            }
        }
        return hit_anything;
    }

private:
    std::vector<ShapePtr> m_list;
};

add 関数で物体を追加できます.物体は std::shared_ptr<Shape> で渡します.これを打つのは面倒なので typedef で定義しておきます.

class Shape;
typedef std::shared_ptr<Shape> ShapePtr;

また,STL のコンテナを使っているので必要なヘッダをインクルードする必要があります.

// rayt.h
#include <vector>

それでは球を2つ表示してみます.そのためには Scene クラスをいくつか修正する必要があります.まず,物体リストのインスタンスを追加します.

std::unique_ptr<Shape> m_world;

次に build 関数で物体を生成します.

ShapeList* world = new ShapeList();
world->add(std::make_shared<Sphere>(vec3(0, 0, -1), 0.5f));
world->add(std::make_shared<Sphere>(vec3(0, -100.5, -1), 100));
m_world.reset(world);

そして,color 関数に物体を渡し,Shape::hit 関数を呼ぶようにします.

vec3 color(const rayt::Ray& r, const Shape* world) const {
    HitRec hrec;
    if (world->hit(r, 0, FLT_MAX, hrec)) {
        return 0.5f*(hrec.n + vec3(1.0f));
    }
    return backgroundSky(r.direction());
}

render 関数では color 関数に物体を渡します.

vec3 c = color(r, m_world.get());

最終的に Scene クラスは次のようになります.

namespace rayt {
    ...

    class Scene {
    public:
        Scene(int width, int height)
            : m_image(std::make_unique<Image>(width, height))
            , m_backColor(0.2f) {
        }

        void build() {

            // Camera

            vec3 w(-2.0f, -1.0f, -1.0f);
            vec3 u(4.0f, 0.0f, 0.0f);
            vec3 v(0.0f, 2.0f, 0.0f);
            m_camera = std::make_unique<Camera>(u, v, w);

            // Shapes

            ShapeList* world = new ShapeList();
            world->add(std::make_shared<Sphere>(
                vec3(0, 0, -1), 0.5f));
            world->add(std::make_shared<Sphere>(
                vec3(0, -100.5, -1), 100));
            m_world.reset(world);
        }

        vec3 color(const rayt::Ray& r, const Shape* world) const {
            HitRec hrec;
            if (world->hit(r, 0, FLT_MAX, hrec)) {
                return 0.5f*(hrec.n + vec3(1.0f));
            }
            return backgroundSky(r.direction());
        }

        vec3 background(const vec3& d) const {
            return m_backColor;
        }

        vec3 backgroundSky(const vec3& d) const {
            vec3 v = normalize(d);
            float t = 0.5f * (v.getY() + 1.0f);
            return lerp(t, vec3(1), vec3(0.5f, 0.7f, 1.0f));
        }

        void render() {

            build();

            int nx = m_image->width();
            int ny = m_image->height();
#pragma omp parallel for schedule(dynamic, 1) num_threads(NUM_THREAD)
            for (int j = 0; j<ny; ++j) {
                std::cerr << "Rendering (y = " << j << ") " << (100.0 * j / (ny - 1)) << "%" << std::endl;
                for (int i = 0; i<nx; ++i) {

                    float u = float(i) / float(nx);
                    float v = float(j) / float(ny);
                    Ray r = m_camera->getRay(u, v);
                    vec3 c = color(r, m_world.get());
                    m_image->write(i, (ny - j - 1), c.getX(), c.getY(), c.getZ());
                }
            }

            stbi_write_bmp("render.bmp", nx, ny, sizeof(Image::rgb), m_image->pixels());
        }

    private:
        std::unique_ptr<Camera> m_camera;
        std::unique_ptr<Image> m_image;
        std::unique_ptr<Shape> m_world;
        vec3 m_backColor;
    };

} // namespace rayt

実行すると次のような画像になります.(rayt106.cpp, rayt106.h)

tuto-raytracing-shape.png

アンチエイリアシング

これまでレンダリングした画像を見てみると球体の周りでジャギーが目立っています.これを取り除くために,各ピクセルごとにいくつも光線を飛ばして(このとき乱数でずらします),その平均を使用します.図で表すと次のようになります.

tuto-raytracing-antialias.png

この手法を「スーパーサンプリング」といいます.まずサンプリング数を指定できるようにするために Scene のコンストラクタを変更します.

Scene(int width, int height, int samples)
    : m_image(std::make_unique<Image>(width, height))
    , m_backColor(0.2f)
    , m_samples(samples) {
}

メンバに m_samples を追加します.

int m_samples;

そして,光線を飛ばしているところを次のように書き換えます.

void render() {
    ...

    int nx = m_image->width();
    int ny = m_image->height();
#pragma omp parallel for schedule(dynamic, 1) num_threads(NUM_THREAD)
    for (int j = 0; j<ny; ++j) {
        std::cerr << "Rendering (y = " << j << ") " << (100.0 * j / (ny - 1)) << "%" << std::endl;
        for (int i = 0; i<nx; ++i) {
            vec3 c(0);
            for (int s = 0; s<m_samples; ++s) {
                float u = (float(i) + drand48()) / float(nx);
                float v = (float(j) + drand48()) / float(ny);
                Ray r = m_camera->getRay(u, v);
                c += color(r, m_world.get());
            }
            c /= m_samples;
            m_image->write(i, (ny - j - 1), c.getX(), c.getY(), c.getZ());
        }
    }

    ...
}

main 関数でサンプリング数を指定します.

int main()
{
    int nx = 200;
    int ny = 100;
    int ns = 100;
    std::unique_ptr<rayt::Scene> scene(std::make_unique<rayt::Scene>(nx, ny, ns));
    scene->render();
    return 0;
}

実行すると次のようになります.(rayt107.cpp)

tuto-raytracing-antialias-output.png

拡散反射

光は物体の表面に当たると,反射するか,透過するか,もしくは両方の現象が起きます.それは物体の表面の材質によって変わります.ほとんどの物体の表面は凹凸があり,眼で確認できるものもあれば,顕微鏡などで見ないとわからないとても小さい凹凸があります.このような表面に光が当たると様々な方向に反射されます.そのような現象を乱反射といい,コンピュータグラフィックスでは拡散反射としてシミュレーションします.

tuto-raytracing-diffuse.png

この現象をシミュレーションするために,光線が物体に当たったとき,ランダムな方向に反射させればよいことになります.それでは,ランダムに反射する方向を決めるにはどうすればよいでしょうか.一様乱数を生成し,単位球内にある点を取り出して使うとよさそうです.図にすると $\vec{s}$ の位置に向かって反射させます.

tuto-raytracing-diffuse-unit-random.png

区間[0,1] の乱数を [-1,1] の区間に変換して,半径 1 の球内であればそれを使用します.

// rayt.h
inline vec3 random_vector() {
    return vec3(drand48(), drand48(), drand48());
}

inline vec3 random_in_unit_sphere() {
    vec3 p;
    do {
        p = 2.f * random_vector() - vec3(1.f);
    } while (lengthSqr(p) >= 1.f);
    return p;
}

random_vector は一様乱数を使ってベクトルを作成します.random_in_unit_sphere は単位球の中の任意の点を生成します.この関数では,ランダムに生成したベクトルの半径が 1 以下になるまで繰り返すようになっています.このように条件を満たすまで一様乱数を何度も生成して繰り返すことを「棄却法」といいます.これを使って拡散反射を実装してみます.color 関数を次のように書き換えます.

vec3 color(const rayt::Ray& r, const Shape* world) const {
    HitRec hrec;
    if (world->hit(r, 0, FLT_MAX, hrec)) {
        vec3 target = hrec.p + hrec.n + random_in_unit_sphere();
        return 0.5f * color(Ray(hrec.p, target - hrec.p), world);
    }
    return backgroundSky(r.direction());
}

ここでは反射率を 50% にしました.次のような画像になります.(rayt108.cpp, rayt108.h)

tuto-raytracing-diffuse-output.png

この画像暗い感じがしますね.これはディスプレイのガンマ補正の影響です.プログラム上では色をリニア空間で扱っています.例えば区間[0,1] の値はリニア空間つまり線形なので,次の図のように直線になります.

tuto-raytracing-linear-space.png

ディスプレイは入力した信号を基本的にガンマ補正するようになっています.特に設定を変更していなければガンマ係数は 2.2 になっています.この場合は下に歪んでいます.

tuto-raytracing-gamma-space.png

つまり,全体的に暗くなってしまいます.ではどうすればいいかというと,ガンマ補正がかかることを想定して,ガンマ補正されたらリニア空間になるように画像データの方を調整します.具体的にはリニア空間のデータを以下のようなカラー空間に変換します.

tuto-raytracing-srgb-space.png

このようなカラー空間を「sRGB空間」といいます.リニア空間とsRGB空間との変換式は次のようになっています.

C_{srgb} = C_{linear}^{1/2.2}, \quad C_{linear} = C_{srgb}^{2.2}

これは近似式で他に厳密な式が存在します.ではこの処理を入れてみます.
最初にリニア空間とsRGB空間との変換関数を追加します.

// rayt.h
inline vec3 linear_to_gamma(const vec3& v, float gammaFactor) {
    float recipGammaFactor = recip(gammaFactor);
    return vec3(
        powf(v.getX(), recipGammaFactor),
        powf(v.getY(), recipGammaFactor),
        powf(v.getZ(), recipGammaFactor));
}

inline vec3 gamma_to_linear(const vec3& v, float gammaFactor) {
    return vec3(
        powf(v.getX(), gammaFactor),
        powf(v.getY(), gammaFactor),
        powf(v.getZ(), gammaFactor));
}

linear_to_gamma がリニア空間からsRGB空間へ,gamma_to_linear がsRGB空間からリニア空間に変換します.

Image クラスの write で色を書き込んていますので,そこでフィルター処理するようにします.また,新しいフィルターを簡単に追加できるような仕組みにします.まず,フィルターの抽象クラスを定義します.

// rayt.h
class ImageFilter {
public:
    virtual vec3 filter(const vec3& c) const = 0;
};

filter 関数は色を受け取ってフィルター処理した色を返します.このフィルタークラスを派生して,ガンマ補正をかけるフィルタークラスを作成します.

// rayt.h
class GammaFilter : public ImageFilter {
public:
    GammaFilter(float factor) : m_factor(factor) {}
    virtual vec3 filter(const vec3& c) const override {
        return linear_to_gamma(c, m_factor);
    }
private:
    float m_factor;
};

ガンマ係数を指定できるようにしていますが,通常は 2.2 になります.
次に Image クラスで,フィルター処理を追加します.ImageFilter クラスを複数持てるようにリストでメンバに追加します.

std::vector< std::unique_ptr<ImageFilter> > m_filters;

コンストラクタでフィルターを生成します.

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<GammaFilter>(GAMMA_FACTOR));
}

そして,write 関数内でフィルター処理を行います.

void write(int x, int y, float r, float g, float b) {
    vec3 c(r, g, b);
    for (auto& f : m_filters) {
        c = f->filter(c);
    }
    int index = m_width*y + x;
    m_pixels[index].r = static_cast<unsigned char>(c.getX() * 255.99f);
    m_pixels[index].g = static_cast<unsigned char>(c.getY() * 255.99f);
    m_pixels[index].b = static_cast<unsigned char>(c.getZ() * 255.99f);
}

この結果,次のようになります.(rayt109.h)

tuto-raytracing-gamma-output.png

この画像,まだおかしなところがあります.このシーンでは光線が最終的に空に向かうのですが,球体の上部や,床に暗いところが多くあります.何か明るいところと暗いところが混じって斑模様のように見えませんか? これは計算誤差によるもので,反射した位置から次に衝突する位置を判定するときに,t が 0 に近い位置で当たったと判定されていて,反射方向がおかしくなっているようです.そのため,当たり判定のときに,t の区間を 0 からではなく例えば 0.001 のようにずらします.

if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
  ...
}

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

tuto-raytracing-gamma-output-corrected.png

先ほどの画像より暗い部分が減っていることがわかると思います.斑模様みたいなものは「シャドウアクネ」と呼ばれることがあります.

材質

物体の反射の特性は表面の材質によって決まります.今は拡散反射だけですが,これから鏡面反射などを追加するためにも,材質をクラスにしましょう.
材質に必要な機能は,その材質の表面に光が当たったときにどの方向に光が反射するかということと,入射した光に対してどれくらい反射するかを表す反射率です.光が物体表面に当たったときに反射する,つまり光の向きを変える現象を「散乱」といいます.また,物体の表面は光を吸収することがあり,吸収されなかった光が散乱して放出されます.反射率は入射した光が吸収されない光の割合とも言えます.この反射率を「アルベド」または「リフレクタンス」といいます.

材質クラスは抽象クラスで次のように定義します.

class Material {
public:
    virtual bool scatter(const Ray& r, const HitRec& hrec, ScatterRec& srec) const = 0;
};

scatter 関数は散乱をシミュレーションします.散乱後の光の向きと,反射率を ScatterRec クラスに格納します.ScatterRec は次のようになっています.

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

ray には散乱後の新しい光線,albedo は反射率です.前に実装した拡散反射のような材質を「ランバート」と言うことがあります.材質クラスとして次のように書き直します.

class Lambertian : public Material {
public:
    Lambertian(const vec3& c)
        : m_albedo(c) {
    }

    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, target - hrec.p);
        srec.albedo = m_albedo;
        return true;
    };

private:
    vec3 m_albedo;
};

このランバート材質はアルベドを指定するようになっています.
次にこの材質を処理するようにいくつか変更します.HitRec に材質を追加します.

class HitRec {
public:
    float t;
    vec3 p;
    vec3 n;
    MaterialPtr mat;
};

MaterialPtr は typedef で定義しています.

class Material;
typedef std::shared_ptr<Material> MaterialPtr;

次に Sphere に材質を指定できるようにし,hit 関数で当たったとき, HitRec の材質に設定するようにします.

class Sphere : public Shape {
public:
    Sphere() {}
    Sphere(const vec3& c, float r, const MaterialPtr& mat)
        : m_center(c)
        , m_radius(r)
        , m_material(mat) {
    }

    virtual bool hit(const Ray& r, float t0, float t1, HitRec& hrec) const override {
        vec3 oc = r.origin() - m_center;
        float a = dot(r.direction(), r.direction());
        float b = 2.0f*dot(oc, r.direction());
        float c = dot(oc, oc) - pow2(m_radius);
        float D = b*b - 4*a*c;
        if (D > 0) {
            float root = sqrtf(D);
            float temp = (-b - root) / (2.0f*a);
            if (temp < t1 && temp > t0) {
                hrec.t = temp;
                hrec.p = r.at(hrec.t);
                hrec.n = (hrec.p - m_center) / m_radius;
                hrec.mat = m_material;
                return true;
            }
            temp = (-b + root) / (2.0f*a);
            if (temp < t1 && temp > t0) {
                hrec.t = temp;
                hrec.p = r.at(hrec.t);
                hrec.n = (hrec.p - m_center) / m_radius;
                hrec.mat = m_material;
                return true;
            }
        }

        return false;
    }

private:
    vec3 m_center;
    float m_radius;
    MaterialPtr m_material;
};

そして,color 関数で材質を処理します.

vec3 color(const rayt::Ray& r, const Shape* world) const {
    HitRec hrec;
    if (world->hit(r, 0.001f, FLT_MAX, hrec)) {
        ScatterRec srec;
        if (hrec.mat->scatter(r, hrec, srec)) {
            return mulPerElem(srec.albedo, color(srec.ray, world));
        }
        else {
            return vec3(0);
        }
    }
    return backgroundSky(r.direction());
}

散乱させて,次の方向から受け取った色と反射率を乗算しています.
あとは作成する球体を追加し,材質を設定します.

ShapeList* world = new ShapeList();
world->add(std::make_shared<Sphere>(
    vec3(0.6, 0, -1), 0.5f,
    std::make_shared<Lambertian>(vec3(0.1f, 0.2f, 0.5f))));
world->add(std::make_shared<Sphere>(
    vec3(-0.6, 0, -1), 0.5f,
    std::make_shared<Lambertian>(vec3(0.8f, 0.0f, 0.0f))));
world->add(std::make_shared<Sphere>(
    vec3(0, -100.5, -1), 100,
    std::make_shared<Lambertian>(vec3(0.8f, 0.8f, 0.0f))));
m_world.reset(world);

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

tuto-raytracing-lambert.png

鏡面反射

金属のような材質のものは,拡散反射はせずに鏡面反射します.物体の表面が完全に平坦なとき,入射角と反射角は同じになります.このような反射ベクトルを求めるには次のように考えます.

tuto-raytracing-reflect.png

$\vec{r}$ が求める反射ベクトルで,$\vec{r} = \vec{v} + 2\vec{b}$ です.$\vec{v}$ は入射ベクトルで,$\vec{b}$ は $\vec{v}$ を $\vec{n}$ に投影したベクトルを反転したものです.$-(\vec{v}\cdot\vec{n})\cdot\vec{n}$ となります.$\vec{n}$ は単位ベクトルで,$\vec{v}$ は任意の長さのベクトルです.これをコードにすると

inline vec3 reflect(const vec3& v, const vec3& n) {
    return v - 2.f * dot(v, n)*n;
}

これを使って鏡面反射する材質を作成します.この材質は「金属」と呼びます.

class Metal : public Material {
public:
    Metal(const vec3& c, float fuzz)
        : m_albedo(c) {
    }

    virtual bool scatter(const Ray& r, const HitRec& hrec, ScatterRec& srec) const override {
        vec3 reflected = reflect(normalize(r.direction()), hrec.n);
        srec.ray = Ray(hrec.p, reflected);
        srec.albedo = m_albedo;
        return dot(srec.ray.direction(), hrec.n) > 0;
    }

private:
    vec3 m_albedo;
};

scatter の戻り値は反射ベクトルと法線ベクトルとのなす角度が 0 より大きいときに真とします.

シーンの球体1つを金属材質に変えてみます.

ShapeList* world = new ShapeList();
world->add(std::make_shared<Sphere>(
    vec3(0.6, 0, -1), 0.5f,
    std::make_shared<Lambertian>(vec3(0.1f, 0.2f, 0.5f))));
world->add(std::make_shared<Sphere>(
    vec3(-0.6, 0, -1), 0.5f,
    std::make_shared<Metal>(vec3(0.8f, 0.8f, 0.8f))));
world->add(std::make_shared<Sphere>(
    vec3(0, -100.5, -1), 100,
    std::make_shared<Lambertian>(vec3(0.8f, 0.8f, 0.0f))));
m_world.reset(world);

次のような画像になります.(rayt113.cpp, rayt113.h)

tuto-raytracing-metal-output.png

今は完全に平坦な表面を考えていましたが,表面は多少凹凸しています.凹凸のある表面に光線が当たると反射ベクトルは理想よりずれます.この現象をシミュレーションするためには,高度なものだと微小面モデルとかあるのですが,ここでは簡易な方法で行います.理想的な反射ベクトルの先を中心とした単位球内の点を無作為に選んでそこまでのベクトルを反射ベクトルとします.図にすると以下のようになります.

tuto-raytracing-reflect-fuzz.png

赤いベクトルは理想的な反射ベクトルです.それに対して青いベクトルはずらした反射ベクトルです.
ずらし具合をパラメータで調整できるようにします.このパラメータを fuzz とします.これを実装します.

class Metal : public Material {
public:
    Metal(const vec3& c, float fuzz)
        : m_albedo(c)
        , m_fuzz(fuzz) {
    }

    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;
        return dot(srec.ray.direction(), hrec.n) > 0;
    }

private:
    vec3 m_albedo;
    float m_fuzz;
};

通常通り反射ベクトルを計算したあとに,単位球から無作為に点を生成して,それに fuzz を乗算して,計算した反射ベクトルに加算します.

reflected += m_fuzz*random_in_unit_sphere();

球体の金属材質を調整します.

std::make_shared<Metal>(vec3(0.8f, 0.8f, 0.8f), 1.f)));

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

tuto-raytracing-metal-fuzz-output.png

ぼやけているように見えます.

さて,拡散反射や鏡面反射において,ランダムな方向に反射させるようになったので,反射が複雑になってきています.color 関数は反射される度に呼び出されていて再帰的です.どれくらい再帰的に呼び出されるか調査してみると 1 つの光線で最高 180~250 回呼ばれていました.光は反射される度に反射率で減衰していきます.今は反射率 0.8 とか設定しているのですが,その反射率で例えば 100 回反射したら最終的に $0.8^{100}$ となりかなり小さい値になり結局黒くなります.反射率と反射回数との関係をグラフにしてみると次のようになります.

tuto-raytracing-reflectance-relationship.png

反射率 0.8 なら反射回数 40 回ぐらいでほぼ 0 になります.先ほどのシーンでは,ランバート材質の物体にも当たることがあり,反射率は 0.5 とかさらに小さいので,これよりもっと早く 0 に収束していきます.このことから,ある程度の反射回数以上は計算の無駄なので,どこかで打ち切ります.この打ち切る条件を決めるために,「ロシアンルーレット」と呼ばれる手法がよく使われるのですが,今回は特定の反射回数で打ち切ります.

color 関数を次のように書き換えます.

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

depth は再起呼び出しの深さで,MAX_DEPTH 以上になったら打ち切ります.MAX_DEPTH は適当に決めます.とりあえず 50 回にしておきます.

// rayt.h
#define MAX_DEPTH 50

color 関数を呼び出しているところで,depth に 0 を渡します.

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

これでレンダリングすると次のようになります.結果はほとんど変わりません.(rayt115.cpp,rayt115.h)

tuto-raytracing-depth-limit.png

屈折

水やガラス,ダイヤモンドなどは見てみると,その先が歪んで見えたり,反射した方向の景色が映りこんで見えます.このような材質に光が当たると,光は鏡面反射する光と物体内部に透過する光に分かれます.光は電磁波の一種で,媒質の中を移動しています.空気も媒質の一つです.例えば光が空気中から別の媒質に入ると光の速度が変化して屈折します.どれくらい屈折するかは入る前の媒質の屈折率,侵入する媒質の屈折率,侵入するときの入射角によって求めることができます.この関係を表したのが「スネルの法則」です.それによれば,媒質 $A$ の屈折率を $\eta_1$, 入射角を $\theta_1$,媒質 $B$ の屈折率を $\eta_2$ ,出射角を $\theta_2$ とすると以下の式が成り立ちます.

\eta_1\sin\theta_1 = \eta_2\sin\theta_2

空気の屈折率は 1, ガラスは 1.3-1.7,ダイヤモンドは 2.4 です.他の屈折率はネットで検索すれば多くの情報が出てくると思います.

屈折ベクトルはスネルの法則から導くことができます.図のように $\vec{r}$ は屈折ベクトル,$\vec{v}$ は方向ベクトル,$\vec{n}$ は表面の法線です.

tuto-raytracing-snell.png

屈折ベクトル $\vec{r}$ は

\vec{r} = -\frac{\eta_1}{\eta_2}(\vec{v}-(\vec{v}\cdot\vec{n})\vec{n})-
\vec{n}\sqrt{
    1-\left(\frac{\eta_1}{\eta_2}\right)^2(1-(\vec{v}\cdot\vec{n}^2)
}

この式の導出については「スネルの法則」を参照してください.これをコードにすると次のようになります.

inline bool refract(const vec3& v, const vec3& n, float ni_over_nt, vec3& refracted) {
    vec3 uv = normalize(v);
    float dt = dot(uv, n);
    float D = 1.f - pow2(ni_over_nt) * (1.f - pow2(dt));
    if (D > 0.f) {
        refracted = -ni_over_nt * (uv - n*dt) - n*sqrt(D);
        return true;
    }
    else {
        return false;
    }
}

D は判別式です.一体何を判別しているのでしょうか.ここでスネルの法則を次のように変形します.

\sin\theta_2 = \frac{\eta_1}{\eta_2}\sin\theta_1

$\theta_1$ は [0,90] の範囲なので,$\sin\theta_1$ の取る範囲は $0 \leq \sin\theta_1 \leq 1$ です.$\eta_1 < \eta_2$ ならば $\frac{\eta_1}{\eta_2} < 1$ となるので,上記の式から $0 \leq \sin\theta_2 \leq 1$ となることがわかります.しかし,$\eta_1 > \eta_2$ の場合は $\frac{\eta_1}{\eta_2} > 1$ なので,$\sin\theta_1$ が大きいと,$\sin\theta_2 > 1$ となって,解がありません.これは,屈折光がなくなり,反射光のみになります.この現象を「全反射」といいます.

屈折をする材質を「誘導体」といいます.これを実装すると

class Dielectric : public Material {
public:
    Dielectric(float ri)
        : m_ri(ri) {
    }

    virtual bool scatter(const Ray& r, const HitRec& hrec, ScatterRec& srec) const override {

        vec3 outward_normal;
        vec3 reflected = reflect(r.direction(), hrec.n);
        float ni_over_nt;
        if (dot(r.direction(), hrec.n) > 0) {
            outward_normal = -hrec.n;
            ni_over_nt = m_ri;
        }
        else {
            outward_normal = hrec.n;
            ni_over_nt = recip(m_ri);
        }

        srec.albedo = vec3(1);

        vec3 refracted;
        if (refract(-r.direction(), outward_normal, ni_over_nt, refracted)) {
            srec.ray = Ray(hrec.p, refracted);
        }
        else {
            srec.ray = Ray(hrec.p, reflected);
            return false;
        }

        return true;
    }

private:
    float m_ri;
};

m_ri は屈折率です.屈折ベクトルは物体の内部に入射するときと,内部から外部に出射するときとでは屈折率が反転するので,内積を計算して判定しています.

この誘導体材質の球を追加します.

ShapeList* world = new ShapeList();
world->add(std::make_shared<Sphere>(
    vec3(0.6, 0, -1), 0.5f,
    std::make_shared<Lambertian>(vec3(0.1f, 0.2f, 0.5f))));
world->add(std::make_shared<Sphere>(
    vec3(-0.6, 0, -1), 0.5f,
    std::make_shared<Dielectric>(1.5f)));
world->add(std::make_shared<Sphere>(
    vec3(0, -0.35, -0.8f), 0.15f,
    std::make_shared<Metal>(vec3(0.8f, 0.8f, 0.8f), 0.2f)));
world->add(std::make_shared<Sphere>(
    vec3(0, -100.5, -1), 100,
    std::make_shared<Lambertian>(vec3(0.8f, 0.8f, 0.0f))));
m_world.reset(world);

次のような画像になります.(rayt116.cpp, rayt116.h)

tuto-raytracing-dieletric-output.png

表面の法線に対して光線の向きが表面に近い角度をグレージング角度といいます.金属や誘導体はグレージング角度に近づくと屈折光が少なくなり,材質によっては全反射します.このような材質が光が入射したときに,どれくらいの比率で反射光と屈折光に分配されるかは「フレネルの方程式」で求めることができます.フレネルの方程式は偏光されていないのであれば,Schlick の近似式がよく使われます.

F_r(\theta) \approx F_{0} + (1-F_{0})(1-\cos\theta)^5

このとき,$F_r(\theta)$ はフレネル反射係数で,$F_{0}$ は,表面に対して垂直に光が入射したときのフレネル反射係数です.$F_r(\theta)$ が,入射した光に対する鏡面反射率を表しています.

$F_{0}$ は屈折率を使って求められます.

F_{0} = \frac{(\eta_1 - \eta_2)^2}{(\eta_1 + \eta_2)^2} = \left(\frac{\eta_1 - \eta_2}{\eta_1 + \eta_2}\right)^2

これをコードにすると次のようになります.

// rayt.h
inline float schlick(float cosine, float ri) {
    float r0 = pow2((1.f - ri) / (1.f + ri));
    return r0 + (1.f - r0) * pow5(1.f - cosine);
}

基本的に大気中からある媒質に入ることを想定するので,片方の媒質は空気で屈折率はほぼ 1 です.
これを誘導体クラスに組み込みます.次のように書き換えます.

class Dielectric : public Material {
public:
    Dielectric(float ri)
        : m_ri(ri) {
    }

    virtual bool scatter(const Ray& r, const HitRec& hrec, ScatterRec& srec) const override {
        vec3 outward_normal;
        vec3 reflected = reflect(r.direction(), hrec.n);
        float ni_over_nt;
        float reflect_prob;
        float cosine;
        if (dot(r.direction(), hrec.n) > 0) {
            outward_normal = -hrec.n;
            ni_over_nt = m_ri;
            cosine = m_ri * dot(r.direction(), hrec.n) / length(r.direction());
        }
        else {
            outward_normal = hrec.n;
            ni_over_nt = recip(m_ri);
            cosine = -dot(r.direction(), hrec.n) / length(r.direction());
        }

        srec.albedo = vec3(1);

        vec3 refracted;
        if (refract(-r.direction(), outward_normal, ni_over_nt, refracted)) {
            reflect_prob = schlick(cosine, m_ri);
        }
        else {
            reflect_prob = 1;
        }

        if (drand48() < reflect_prob) {
            srec.ray = Ray(hrec.p, reflected);
        }
        else {
            srec.ray = Ray(hrec.p, refracted);
        }

        return true;
    }

private:
    float m_ri;
};

cosine は入射角です.全反射しなければ,近似式を使ってフレネル反射率を求めます.それを使って,反射するか屈折するかを決定します.半径を負にして少し小さくした球体を追加し,すでにある誘導体材質の球体に重ねると,水泡のような表現ができます.

ShapeList* world = new ShapeList();
world->add(std::make_shared<Sphere>(
    vec3(0.6, 0, -1), 0.5f,
    std::make_shared<Lambertian>(vec3(0.1f, 0.2f, 0.5f))));
world->add(std::make_shared<Sphere>(
    vec3(-0.6, 0, -1), 0.5f,
    std::make_shared<Dielectric>(1.5f)));
world->add(std::make_shared<Sphere>(
    vec3(-0.6, 0, -1), -0.45f,
    std::make_shared<Dielectric>(1.5f)));
world->add(std::make_shared<Sphere>(
    vec3(0, -0.35, -0.8f), 0.15f,
    std::make_shared<Metal>(vec3(0.8f, 0.8f, 0.8f), 0.2f)));
world->add(std::make_shared<Sphere>(
    vec3(0, -100.5, -1), 100,
    std::make_shared<Lambertian>(vec3(0.8f, 0.8f, 0.0f))));
m_world.reset(world);

この結果は次のようになります.(rayt117.cpp, rayt117.h)

tuto-raytracing-dieletric-fresnel-output.png

沢山表示してみる

物体を多く置いてみます.また,カメラを LookAt 方式で指定します.

// Camera

vec3 lookfrom(13,2,3);
vec3 lookat(0,0,0);
vec3 vup(0,1,0);
float aspect = float(m_image->width()) / float(m_image->height());
m_camera = std::make_unique<Camera>(lookfrom, lookat, vup, 20, aspect);

// Shapes

ShapeList* world = new ShapeList();

int N = 11;
for (int i=-N; i<N; ++i) {
    for (int j=-N; j<N; ++j) {
        float choose_mat = drand48();
        vec3 center(i+0.9f*drand48(), 0.2f, j+0.9f*drand48());
        if (length(center-vec3(4,0.2,0)) > 0.9f) {
            if (choose_mat < 0.8f) {
                world->add(std::make_shared<Sphere>(
                    center, 0.2f, 
                    std::make_shared<Lambertian>(mulPerElem(random_vector(),random_vector()))));
            }
            else if (choose_mat < 0.95f) {
                world->add(std::make_shared<Sphere>(
                    center, 0.2f,
                    std::make_shared<Metal>(0.5f * (random_vector()+vec3(1)), 0.5f*drand48())));
            }
            else {
                world->add(std::make_shared<Sphere>(
                    center, 0.2f,
                    std::make_shared<Dielectric>(1.5f)));
            }
        }
    }
}

world->add(std::make_shared<Sphere>(
    vec3(0, -1000, -1), 1000,
    std::make_shared<Lambertian>(vec3(0.5f, 0.5f, 0.5f))));
m_world.reset(world);

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

tuto-raytracing-scene1.png

次は「テクスチャとコーネルボックス」 です.


Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away