📌 光学
ここからはモンテカルロ法を使ってレイトレーシングを行う手法を説明します.そのためにまずは少し光学について触れておきます.これまで行ってきた光線を飛ばして光伝達をシミュレーションするような方法をレイトレーシングといいました.これは幾何光学という分野になります.なので,レイトレーシングは幾何光学で説明できる物理現象(例えば反射や屈折)を実現することが出来ます.それとは別に,光は電磁波の一種なので,波動の特性を持っています.光の波動による物理現象(例えば回折や干渉)は波動光学の分野になります.もちろん電磁波なので電磁学にも関連しています.例えば偏光がこの分野になります.光の物理現象をシミュレーションするには,こういった分野がすべて必要になってくるのですが,実際は計算が複雑になったり,よくわからん(切実)ので,基本的に幾何光学だけを考えます.フレネルの方程式のときに偏光が出てきましたが,無視して考えていたのはそういうことです.そのため,写実的なものをレンダリングするためには,まず幾何光学について知っておく必要があります.
📌 光の物理量
コンピュータグラフィックスにおいて,光の最も小さい単位は「光子」です.これを「フォトン(photon)」といいます.光のエネルギーというのはこの光子が集まったもので,放射エネルギーといいます.細かいことをいうと波長が関係してくるのですが,ここでは割愛します.
この放射エネルギーの時間的な変化を表したものを「放射束」といいます.これは放射エネルギーの仕事量のことでもあるので単位は「ワット(W)」です.電球とかで 60W とかありますよね? あれです.
これまで計算してきた数値は光学においてどのような値(量)なのでしょうか.それは「放射束」です.光伝達というのは放射束を輸送することを指します.
余談ですが,人間の眼が知覚する光の量を測光値といい,この分野を測光学といいます.測光学における放射束を「光束」といいます.
📌 モンテカルロレイトレーシング
古典的なレイトレーシングは,カメラから光線を飛ばして鏡面反射させていました.その後,鏡面反射の方向の周りに向かって複数の光線を飛ばして計算するように改善されました.これを分散レイトレーシングといいます.これから実装していくモンテカルロレイトレーシングはこの分散レイトレーシングをさらに改善させたものです.
光源から放射された放射束は反射を繰り返して物体表面に入射します.物体の表面は拡散反射や鏡面反射などを行ってあらゆる方向に放射束を反射します.そして,物体表面のある点を考えると,半球上のあらゆるところから放射束が入ってくることになります.つまり,半球上から入ってくるすべての放射束を積分する必要があります.これをシミュレーションしようとするのは困難なので,統計的にサンプリングを行って,結果を推定します.そこでモンテカルロ法を使います.積分範囲内で「無作為」に光線を発生させてサンプリングし,平均値を求めて積分の推定値とします.このように表面のある点に入ってくる放射束をモンテカルロ法で求めて,反射や屈折などを計算する方法をモンテカルロレイトレーシングといいます.
📌 モンテカルロ法について
モンテカルロ法,またそれに関する確率と統計については「CGのための確率・統計入門」を参照してください.また,立体角といった新しい用語が出てきます.それらについては「基礎からはじめる物理ベースレンダリング」で説明していますので,分からない場合は参照してください.これからはモンテカルロ法に基づいて,各種の物理現象をシミュレーションしていきます.前回の最後でコーネルボックスをレンダリングしました.そのコーネルボックスを使ってモンテカルロレイトレーシングを実装していきます.
📌 光の散乱
材質のところで散乱について説明しました.ここでもう一度定義すると,散乱とは物体表面に光線が当たったときに,光線の向きが変わる現象のことです.また,物体表面で吸収されずに,反射される光の割合を反射率(アルベド)といいます.ここで,光とは放射束のことです.
散乱したときの方向分布は,立体角の確率密度関数(以降,pdfと記述します)によって表され,これを
ある表面上の点に入射してくる全放射束は
ここで
ランバート面では
となり,今までの計算と変わらないことになります.
余談ですが,物理ベースレンダリングとかでよくでてくる双方向反射率分布関数(Bidirectional reflectance distribution function: BRDF)との関係は
であり,ランバート面の
これはランバートの BRDF と一致します.
また,関与媒質において
📌 重点的サンプリング
Importance sampling.
散乱 pdf と重み
材質クラスに散乱pdfを返す関数 scattering_pdf
を追加します.デフォルトで 0 を返します.
virtual float scattering_pdf(const Ray& r, const HitRec& hrec) const { return 0; }
ランバート材質の scattering_pdf
は
virtual float scattering_pdf(const Ray& r, const HitRec& hrec) const override {
return std::max(dot(hrec.n, normalize(r.direction())), 0.0f) / PI;
}
また,モンテカルロ法で積分するために 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
)
ノイズが多くありますね.これからこのノイズを無くしていきます.
半球上の無作為な方向の生成
重点的サンプリングを行うには表面の法線を中心に半球上の無作為な方向を生成する必要があります.まずは
方位角
確率変数
極角
単位球の面積は
極座標(
平方根の項を整理すると
なので,
これは単位球上の無作為な方向なので,半球上に制限します.半球上で一様分布な pdf は
デカルト座標に変換すると
次に pdf が
で置換積分します.
よって
デカルト座標で表すと
これをコードに書き起こしたものが以下になります.
// 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()
結果は次の通りです.
一段目の一番左を見ると半球上に点が生成されているのがわかります.一段目の中央,XZ のプロットを見てみるとZ値が高くなると点の数も増えているように見えませんか? ちょっとわかりづらいかもしれません.そんなときは2段目のヒストグラムを見てみます.一番右が Z値ですが,段階的になっていることがわかります.つまり,Z値が高くなるほど多くなっています.Zと比べてXYはどうでしょうか.1段目の一番右を見てみると円の中で一様分布しているように見えます.先ほどの関数はこのような分布をするベクトルを無作為に生成します.
正規直交基底
Orhto-normal bases.
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
)
光源のサンプリング
レイトレーシングでは,カメラから光線を飛ばして,最終的に光源まで辿りつかないと入射する放射束の精度が大きく下がります.光線の反射は基本的にランダムな方向なのですが,意図的に光源方向に光線を飛ばしてサンプリングすれば光源にたどり着く可能性が上がります(光線の開始点と光源との間に物体があると,光源に辿り着かないことがあります)
コーネルボックスでは光源は四角形です.この領域上をサンプリングするために適切な確率密度関数を決定する必要があります.
表面上のある点
点
球上の点をサンプリングする確率は
上の式から
この式を使って,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
)
光源の周りにつぶつぶがありますね.これは光源と天井の間に隙間があるためです.今の光源は両面から当たっても発光するので,片面のみにします.そのため,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
)
改善されているのがわかると思います.
合成PDF
これまでに,光源をサンプリングするための pdf と,
合成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
)
次に物体そのものを 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
)
それでは合成 pdf を実装します.今回は2つの 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
)
これまでは 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());
}
あとはコーネルボックスの箱を1つをアルミニウム(金属)にしてみます.
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
)
左上のところにノイズが強く発生してしまっています.これはアルミニウム箱へのサンプリングが減ってしまったためです.そのため,その箱に重点的サンプリングすると改善できるはずです.ただし,今回は箱ではなく球にして重点的サンプリングします.その理由は箱より簡単だからです.
ある点から球に対して重点的サンプリングする場合,球の立体角を一様分布でサンプリングすればいいと考えられます.一様分布に従う球のベクトルについてはすでにやりました.それによると,
今,
この
よって,
また,
これを先ほどの式に代入すると
これを使ってベクトルを生成できます.
ここで,
この図から
次に方向分布を表す pdf を算出する必要があります.球の全立体角は
これらを使って 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);
あとは,コーネルボックスの箱を1つ球に変更して,材質も 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
)
ノイズがひどいですね.これは重点的サンプリング用の物体が球だけだからです.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
)
これが目標の画像になります.
さて,とうとう最後になりましたが,どこかの計算でゼロ除算があったり,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 著の3冊をお勧めします.
- amazon: Ray Tracing in One Weekend
- amazon: Ray Tracing: the Next Week
- amazon: Ray Tracing: The Rest Of Your Life
この入門では取り扱っていない関与媒質,BVH,ノイズ,モーションブラー,被写界深度についても書かれています.一部名前を変えていますが,ほとんどが同じ名前になっているので,わかりやすくなっているのではないかと思います.冒頭でもお伝えした通り,この本の内容は現在インターネットから無料で手に入ります.
smallpt
とても小さいコードのパストレーシングです.次にパストレーシングに挑戦するなら参考になると思います.
smallpt の解説もありますのでこちらも参考にしてみてください.
Presentation slides about smallpt by David Cline
edupt
上記の smallpt を c++ で書き直し,ソースコードには日本語のコメントが沢山付いている物理ベースレンダラです.とりあえず,次のことに迷ったらこちらのソースコードを見ることをお勧めします.edupt の解説もありますので参考にしてみてください.(上記リンク先から参照できます)
また,作者の方は「Computer Graphics Gems JP 2015 コンピュータグラフィックス技術の最前線」という本で,「パストレーシングで始めるオフライン大域照明レンダリング入門」を執筆されているので,そちらもお勧めします.
memoRANDOM
レイトレーシングに関する情報が豊富にある日本語のサイトです.とても重宝しています.
「フォトンマッピングー実写に迫るコンピュータグラフィックス」
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」
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著,狩野智英訳「ゲームプログラミングのための3Dグラフィックス数学」ボーンデジタル社,2006
- 鈴木健太郎,「パストレーシングで始めるオフライン大域照明レンダリング入門」『Computer Graphics Gems JP 2015 コンピュータグラフィックス技術の最前線』ボーンデジタル,2015, p. 3