今日やること
- 複数の物体を扱えるようにする
- 影の計算
- 乱数の生成
- スーパーサンプリング(SSAA)
- 並列化計算
- アンビエントオクルージョン(AO)
複数の物体を扱えるようにする
前回までで、一つの球に影をつけて表示することができるようになりました。今回は複数の球を扱えるように拡張してみましょう。
Accelクラス
球の集合を保持するAccelクラスを作りましょう。
メンバ変数
std::vector<std::shared_ptr<Sphere>> spheres
球の集合
メンバ関数
void add(std::shared_ptr<Sphere> p)
球を追加するbool intersect(const Ray& ray, Hit& hit) const
rayと衝突計算を行い、結果をhitに格納する
#ifndef ACCEL_H
#define ACCEL_H
#include <vector>
#include <memory>
#include "ray.h"
#include "hit.h"
#include "sphere.h"
class Accel {
public:
std::vector<std::shared_ptr<Sphere>> spheres;
Accel() {};
void add(std::shared_ptr<Sphere> p) {
...
};
bool intersect(const Ray& ray, Hit& hit) const {
...
};
};
#endif
intersectの実装
まず全物体に対して衝突計算します。
Hit hit_each;
for(auto p : spheres) {
if(p->intersect(ray, hit_each)) {
//当たった
}
}
各物体の衝突情報はhit_each
に格納されます。
複数の物体との衝突計算を行う場合、返す必要のある衝突情報は 最も衝突距離が近い物体の衝突情報 だけです。
これを踏まえるとAccel::intersect()
は次のようになるでしょう。
bool intersect(const Ray& ray, Hit& hit) const {
bool isHit = false;
//hitの距離を大きな値で初期化しておく
hit.t = 10000000;
//全ての物体に対し衝突計算をする
Hit hit_each;
for(auto p : spheres) {
if(p->intersect(ray, hit_each)) {
isHit = true;
//衝突距離が短いものが見つかったらhitをそれで置き換える
if(hit_each.t < hit.t) {
hit = hit_each;
}
}
}
return isHit;
};
最初にhit.t
に大きな値をセットしておき、もしhit_each.t
がhit.t
より小さかったらhit
をhit_each
で更新するという処理です。
以上のことを踏まえてAccelクラスを自分で実装してみましょう。
練習問題
- Accelクラスを実装してみよう
複数の球で法線レンダリング
これで複数の球を扱えるようになったので、さっそく球を2つ追加して法線表示を行ってみましょう。
#include <iostream>
#include <memory>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"
int main() {
Image img(512, 512);
Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
Accel accel;
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0))); //球1
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000))); //球2
for(int i = 0; i < img.width; i++) {
for(int j = 0; j < img.height; j++) {
double u = (2.0*i - img.width)/img.width;
double v = (2.0*j - img.height)/img.height;
Ray ray = cam.getRay(u, v);
Hit hit;
if(accel.intersect(ray, hit)) {
img.setPixel(i, j, (hit.hitNormal + Vec3(1, 1, 1))/2.0); //法線表示
}
else {
img.setPixel(i, j, Vec3(0, 0, 0));
}
}
}
img.ppm_output();
}
画像が逆さまになってしまっていますね。でもこれはバグではありません。ピンホールカメラ を使っているので、逆さまの画像が出力されるのです。
しかし逆さまに表示されるのは不便ですね。ピンホールカメラの実装を変更して直しましょう。
Ray getRay(double u, double v) const {
v = -v; //変更
Vec3 pinhole = camPos + camForward;
Vec3 sensorPos = camPos + u*camRight + v*camUp;
return Ray(camPos, normalize(pinhole - sensorPos));
}
影の計算
複数の物体を扱えるようになったので、いよいよ影の計算をやってみましょう。
まず前回と同じように、光源の方向を指定します。
Vec3 lightDir = normalize(Vec3(1, 1, -1))
前回と同じようにレイトレーシングします。
for(int i = 0; i < img.width; i++) {
for(int j = 0; j < img.height; j++) {
double u = (2.0*i - img.width)/img.width;
double v = (2.0*j - img.height)/img.height;
Ray ray = cam.getRay(u, v);
Hit hit;
if(accel.intersect(ray, hit)) {
//物体に色をつける処理
}
else {
img.setPixel(i, j, Vec3(0, 0, 0));
}
}
}
物体上の点に影がつくということは、点から光源の方向を見たときに何らかの物体がそこに存在するので、光が遮られるということです。
つまり、光源の方向に物体があるかどうか判定し、物体があれば影をつけるという処理にすればよいのです。
Ray shadowRay = Ray(hit.hitPos, lightDir);
Hit hit2;
//光源との間に物体がある
if(accel.intersect(shadowRay, hit2)) {
//影をつける
}
else {
//影をつけない
}
これを実装すると次のようになります。
#include <iostream>
#include <memory>
#include <algorithm>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"
int main() {
Image img(512, 512);
Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
Accel accel;
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0)));
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000)));
Vec3 lightDir = normalize(Vec3(1, 1, -1));
for(int i = 0; i < img.width; i++) {
for(int j = 0; j < img.height; j++) {
double u = (2.0*i - img.width)/img.width;
double v = (2.0*j - img.height)/img.height;
Ray ray = cam.getRay(u, v);
Hit hit;
//最初に飛ばすRayが当たった場合
if(accel.intersect(ray, hit)) {
Ray shadowRay = Ray(hit.hitPos, lightDir); //光源の方向に飛ばすRay
Hit hit_shadow;
//もし当たらなかったら色をつける
if(!accel.intersect(shadowRay, hit_shadow)) {
double I = std::max(dot(hit.hitNormal, lightDir), 0.0); //光の強さ
img.setPixel(i, j, I*Vec3(1, 1, 1));
}
//当たったら影にする
else {
img.setPixel(i, j, Vec3(0, 0, 0));
}
}
else {
img.setPixel(i, j, Vec3(0, 0, 0));
}
}
}
img.ppm_output();
}
ノイズだらけの不気味な画像が出てきました.
なぜこんな結果になってしまったのでしょう。
コンピューターによる浮動小数点数の演算は誤差がつきものという話をしましたね。この結果は実はそこからきています。最初のRayの衝突点を求めるという計算において、わずかな誤差が発生するので、衝突点が床の下側にめり込んでしまっているのです。そのために、そこから光源の方向に伸びるRayが必ず衝突してしまうことになり、影がついて黒くなってしまうのです。
これを防ぐには、衝突点を法線方向にわずかにずらすという処理を加えます。
Ray shadowRay = Ray(hit.hitPos + 0.001*hit.hitNormal, lightDir);
Hit hit2;
//光源との間に物体がある
if(accel.intersect(shadowRay, hit2)) {
//影をつける
}
else {
//影をつけない
}
プログラム全体は次のようになります。
#include <iostream>
#include <memory>
#include <algorithm>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"
int main() {
Image img(512, 512);
Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
Accel accel;
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0)));
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000)));
Vec3 lightDir = normalize(Vec3(1, 1, -1));
for(int i = 0; i < img.width; i++) {
for(int j = 0; j < img.height; j++) {
double u = (2.0*i - img.width)/img.width;
double v = (2.0*j - img.height)/img.height;
Ray ray = cam.getRay(u, v);
Hit hit;
if(accel.intersect(ray, hit)) {
Ray shadowRay = Ray(hit.hitPos + hit.hitNormal*0.001, lightDir); //法線方向にわずかにずらす
Hit hit_shadow;
if(!accel.intersect(shadowRay, hit_shadow)) {
double I = std::max(dot(hit.hitNormal, lightDir), 0.0);
img.setPixel(i, j, I*Vec3(1, 1, 1));
}
else {
img.setPixel(i, j, Vec3(0, 0, 0));
}
}
else {
img.setPixel(i, j, Vec3(0, 0, 0));
}
}
}
img.ppm_output();
}
やった!! ちゃんとした影がつきました!!
乱数の生成
ここで急に話の流れが変わりますが、乱数の生成方法について説明しましょう。
C言語の方法
古いC言語の方法では次のようにして乱数を生成していました。
#include <cstdio>
#include <cstdlib>
int main() {
//シード値
srand(1);
for(int i = 0; i < 100; i++) {
printf("%f\n", (double)rand()/RAND_MAX);
}
return 0;
}
実行結果
0.840188
0.394383
0.783099
0.798440
0.911647
0.197551
0.335223
0.768230
...
srand()
には シード値(seed) という値をセットします。コンピューターで生成できる乱数はあくまで疑似乱数であり、何らかの漸化式などに従って数列を次々と生成していきます。このとき初期値が必要になるのですが、これがシード値です。
rand()
は0からRAND_MAX
までの整数の疑似乱数を返します。プログラムでは、それをRAND_MAX
で割ることで、0から1までの疑似乱数にしています。
シード値が同じなら生成される疑似乱数列は全く同じものになります。なので、通常はシード値に時刻などの値をセットします。
#include <cstdio>
#include <cstdlib>
#include <ctime>
int main() {
//シード値
srand((unsigned long)time(NULL));
for(int i = 0; i < 100; i++) {
printf("%f\n", (double)rand()/RAND_MAX);
}
return 0;
}
こうすると実行するたびに表示される疑似乱数列が変わるはずです。
この疑似乱数は 線形合同法 という方法で生成されているのですが、科学技術計算に用いようとする場合には次のような問題点があります。
- 周期が短い
- 分布が一様でない
メルセンヌ・ツイスタ
STLにはrandom
という乱数を扱うライブラリーが存在します。その中にあるメルセンヌ・ツイスタと呼ばれる疑似乱数生成方法は
- 周期が長い
- 分布が一様
という優れた特徴を持つために、科学技術計算において使用される疑似乱数生成方法です。これからレイトレにおいて乱数が必要になった場合にはこのメルセンヌ・ツイスタを使用していきます。
メルセンヌ・ツイスタを用いて乱数を生成するには次のようにします。
#include <cstdio>
#include <random>
std::random_device dev_rnd;
std::mt19937 mt(dev_rnd());
std::uniform_real_distribution<> dist(0, 1);
double rnd() {
return dist(mt);
}
int main() {
for(int i = 0; i < 100; i++) {
printf("%f\n", rnd());
}
return 0;
}
少し複雑なので一つ一つ見ていきます。
std::random_device dev_rnd;
std::random_device
はデバイスの物理的な情報(熱など)から乱数を生成してくれるものです。動作が遅いので、多くの乱数生成に使用することはできません。その代わりに擬似乱数生成器の初期値として使用します。
std::mt19937 mt(dev_rnd())
std::mt19937
はメルセンヌ・ツイスタによる擬似乱数生成器です。初期値をデバイス乱数で初期化しています。
std::uniform_real_distribution<> dist(0, 1);
std::uniform_real_distribution<>
は疑似乱数生成器による乱数を指定された範囲の一様分布に従うように変換してくれるものです。この場合は0から1までの一様分布になるように変換されます。
double rnd() {
return dist(mt);
}
std::uniform_real_distribution<>
の引数に疑似乱数生成器を入れてあげれば、0から1までの一様分布に従う乱数が得られます。
スーパーサンプリング(SSAA)
先ほどのレンダリング画像をもう一度見てみましょう。 物体の境界、影の境界がギザギザしていますね。これは本来連続的な曲線である物体の境界などを、離散的な画素で表現しているために現れるもので、エイリアシング(Aliasing) と呼ばれています。
これを防ぐ方法が スーパーサンプリング(SSAA) と呼ばれるものです。
一つの画素内で始点をわずかにずらして複数のRayを飛ばし、平均した色をその画素に書き込みます。こうすることで、物体の境界線などがやんわりとした感じになり、美しく見えるようになります。
メルセンヌ・ツイスタを使ってこれを実装すると次のようになります。
#include <iostream>
#include <memory>
#include <algorithm>
#include <random>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"
std::random_device dev_rnd;
std::mt19937 mt(dev_rnd());
std::uniform_real_distribution<> dist(0, 1);
double rnd() {
return dist(mt);
}
int main() {
Image img(512, 512);
Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
Accel accel;
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0)));
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000)));
Vec3 lightDir = normalize(Vec3(1, 1, -1));
for(int k = 0; k < 100; k++) {
for(int i = 0; i < img.width; i++) {
for(int j = 0; j < img.height; j++) {
//ピクセル内の位置をランダムにずらす
double u = (2.0*(i + rnd()) - img.width)/img.width;
double v = (2.0*(j + rnd()) - img.height)/img.height;
Ray ray = cam.getRay(u, v);
Vec3 color;
Hit hit;
if(accel.intersect(ray, hit)) {
Ray shadowRay = Ray(hit.hitPos + hit.hitNormal*0.001, lightDir);
Hit hit_shadow;
if(!accel.intersect(shadowRay, hit_shadow)) {
double I = std::max(dot(hit.hitNormal, lightDir), 0.0);
color = I*Vec3(1, 1, 1);
}
else {
color = Vec3(0, 0, 0);
}
}
else {
color = Vec3(0, 0, 0);
}
img.setPixel(i, j, img.getPixel(i, j) + 1/100.0*color); //平均をとる
}
}
}
img.ppm_output();
}
先ほどよりだいぶ綺麗になりましたね!!
並列化計算
SSAAもやるようになってくると、だいぶ計算に時間がかかるようになってきます。もう少し早く計算できるようになると嬉しいですね。
現在のCPUにはたいてい複数のコアが搭載されているのですが、今まで書いてきたプログラムは一つのコアのみを用いて計算を行っています。
早く計算できるようにするために、CPUのコアを全て使って計算を行うことを考えてみましょう。
このようにfor文を幾つかの領域に分けて、それぞれの領域を別々のCPUコアで計算することができるようになります。
openMP
openMPは並列化計算を簡単に実現できるようにしてくれるプラットフォームです。
openMPのインストール
https://qiita.com/kaityo256/items/ae9329dae24ea8828ae0
openMPの使い方
int main() {
for(int i = 0; i < 100; i++) {
c[i] = a[i] + b[i];
}
}
並列化
#include <omp.h> //これが必要
int main() {
#pragma omp parallel for
for(int i = 0; i < 100; i++) {
c[i] = a[i] + b[i];
}
}
#pragma omp parallel for
を並列化したいfor文の直前につけるだけです。
さらにコンパイルするときに次のようなオプションを追加してあげます。
g++ -fopenmp main.cpp
これでCPUのコアを全て使って並列化計算できるようになります。
並列化計算の問題点
並列化計算が正しく行われるためには、それぞれの計算タスクが独立 している必要があります。独立していない場合、計算が行われる順序によって計算結果が変わってきてしまいます。
しかしそれを防ごうとして計算順序を指定すると、その計算が終わるまで他の計算を進めることができなくなるので、結果的に並列化の恩恵を受けられなくなります。
つまり、何でもかんでも並列化すればいいというわけではなく、計算タスクが独立している部分だけに使用するべきです。
幸い、レイトレーシングは画素ごとに色を求める計算は互いに独立しているので、この部分に対して並列化を適用することができます。
並列化計算によるレイトレ
#include <iostream>
#include <memory>
#include <algorithm>
#include <random>
#include <omp.h>
#include "vec3.h"
#include "ray.h"
#include "hit.h"
#include "sphere.h"
#include "image.h"
#include "camera.h"
#include "accel.h"
std::random_device dev_rnd;
std::mt19937 mt(dev_rnd());
std::uniform_real_distribution<> dist(0, 1);
double rnd() {
return dist(mt);
}
int main() {
Image img(512, 512);
Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
Accel accel;
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, 0, 0), 1.0)));
accel.add(std::make_shared<Sphere>(Sphere(Vec3(0, -10001, 0), 10000)));
Vec3 lightDir = normalize(Vec3(1, 1, -1));
#pragma omp parallel for
for(int k = 0; k < 100; k++) {
for(int i = 0; i < img.width; i++) {
for(int j = 0; j < img.height; j++) {
//ピクセル内の位置をランダムにずらす
double u = (2.0*(i + rnd()) - img.width)/img.width;
double v = (2.0*(j + rnd()) - img.height)/img.height;
Ray ray = cam.getRay(u, v);
Vec3 color;
Hit hit;
if(accel.intersect(ray, hit)) {
Ray shadowRay = Ray(hit.hitPos + hit.hitNormal*0.001, lightDir);
Hit hit_shadow;
if(!accel.intersect(shadowRay, hit_shadow)) {
double I = std::max(dot(hit.hitNormal, lightDir), 0.0);
color = I*Vec3(1, 1, 1);
}
else {
color = Vec3(0, 0, 0);
}
}
else {
color = Vec3(0, 0, 0);
}
img.setPixel(i, j, img.getPixel(i, j) + 1/100.0*color); //平均をとる
}
}
}
img.ppm_output();
}
並列計算が正しく動いていれば、計算時間が先ほどより短くなっているはずです。
アンビエントオクルージョン(AO)
ある点から空がどれくらいの割合で見えるかを表すのが アンビエントオクルージョン(AO) です。物体が隣接している場所は空が見える範囲が狭くなるので割合が小さくなり、逆に何もないところでは0.5に近づきます。この値を可視化してみましょう。
ランダムな方向の生成
ランダムに単位球上の点への方向を生成する関数Vec3 randomDir()
を作ってみましょう。
天球全体へのランダムな方向を作るには、単位球上の点をランダムに生成して、原点からその点への方向を代わりに使えば良いです。
そこで、単位球上の点をランダムに生成することを考えます。球面座標系のとをランダムに生成すると、単位球上の点は
となります。これをそのまま方向ベクトルとして使えばよいです。(実は一様サンプリングにはなっていない)
Vec3 randomDir() {
double theta = M_PI*rnd();
double phi = 2*M_PI*rnd();
double x = std::cos(phi)*std::sin(theta);
double y = std::cos(theta);
double z = std::sin(phi)*std::sin(theta);
return Vec3(x, y, z);
}
AOの計算
あとはshadowRayの方向をrandomDir()
で生成し、結果の平均を取るようにすればよいです。
#pragma omp parallel for
for(int k = 0; k < 100; k++) {
for(int i = 0; i < img.width; i++) {
for(int j = 0; j < img.height; j++) {
//ピクセル内の位置をランダムにずらす
double u = (2.0*(i + rnd()) - img.width)/img.width;
double v = (2.0*(j + rnd()) - img.height)/img.height;
Ray ray = cam.getRay(u, v);
Vec3 color;
Hit hit;
if(accel.intersect(ray, hit)) {
Vec3 lightDir = randomDir(); //ランダムな方向を生成
Ray shadowRay = Ray(hit.hitPos + hit.hitNormal*0.001, lightDir);
Hit hit_shadow;
if(!accel.intersect(shadowRay, hit_shadow)) {
color = Vec3(1, 1, 1);
}
else {
color = Vec3(0, 0, 0);
}
}
else {
color = Vec3(1, 1, 1);
}
img.setPixel(i, j, img.getPixel(i, j) + 1/100.0*color); //平均をとる
}
}
}
影がとてもリアルに見えますね!!
AOの計算でやっていることを数式で表す
天球全体を、点から方向が見えるということをで表すと、点におけるAOの計算は次のように表すことができます。
はからの方向に何もないときに、それ以外のときはになる関数です。方向のことを 立体角(Solid Angle) といいます。全球の立体角はになるので、積分をで割っています。
上の計算では、ランダムな方向列を生成し、
とすることで、上の積分を近似していました。一般的には上の積分と、この近似式の値はを大きくしても一致することはありません。
上の積分を正しく近似して計算した結果は以下のようになります。 先ほどの近似計算による結果とはかなり異なります。
では、上の積分に一致するように近似式を作るにはどうすればよいのでしょうか?この続きは物理ベースレンダリング理論のところでやります。
練習問題
問1
絶対屈折率から絶対屈折率の媒質に入射するレイがある。媒質境界面の法線をとしたときに、レイが屈折した方向を求めよ。なおとする。
補足
スネルの法則により、入射角と屈折角には の関係がある。
問2
講義でランダムな方向を生成する関数randomDir()
を作ったが、実はこの関数は単位球面上の点の一様サンプリングにはなっていない。例えば、北極、南極点に近いところでは点が密集する傾向にある。単位球面上での一様サンプリングを行うにはどうすればよいだろうか?計算方法を考えてみよ。
問3
円周率を乱数を用いて評価してみよう。 の正方形内の点を個ランダムに生成する。これらの点が中心座標で半径の円に入った回数をとする。すると
で推定することができる。この方法で円周率を求めるプログラムを作成せよ。