古典的レイトレーサーの実装(第2回)
古典的レイトレーサーを実装し, 簡単なシェーディングを行う

今日やること

  • Ray Tracerについて説明
  • ヘッダファイルについて
  • PPM画像について
  • 画像クラスを作ってみよう
  • 画像クラスで遊んでみよう
  • 座標系について
  • カメラクラスを作ってみよう
  • はじめてのレイトレ

Ray Tracerについて説明

https://speakerdeck.com/yumcyawiz/classical-ray-tracing

ヘッダファイルについて

今まで作ってきたクラスを一つのファイルにまとめると他のプログラムからそれらのクラスを利用することができます。この仕組みを実現するのが ヘッダファイル と呼ばれるもので、拡張子が.hとなっているファイルのことを言います。

実は今まで書いてきたプログラムで共通して使われているヘッダファイルがあります。

#include <cstdio>


int main() {
  ...
}

この#include <cstdio>という命令は、実は/usr/includeにあるstdio.hというヘッダファイルを読み込むという命令になっています。stdio.hの中にprintfscanfが定義されているので、それらの関数を使うことができるわけです。

定義と実装

stdio.hを開いてみると分かるのですが、ヘッダファイルの中に書かれているのは関数のプロトタイプ宣言だけです。

/* Write formatted output to STREAM.

   This function is a possible cancellation point and therefore not
   marked with __THROW.  */
extern int fprintf (FILE *__restrict __stream,
		    const char *__restrict __format, ...);
/* Write formatted output to stdout.

   This function is a possible cancellation point and therefore not
   marked with __THROW.  */
extern int printf (const char *__restrict __format, ...);
/* Write formatted output to S.  */
extern int sprintf (char *__restrict __s,
		    const char *__restrict __format, ...) __THROWNL;

実装は別のファイルで行われていて、その拡張子が.cppになっています。

digraph {
  node [shape=box]
  h
  cpp
}

一般的にはこのように定義と実装を別のファイルでやるのですが、面倒くさいのでこれからレイトレーサーを作っていく時にはヘッダファイルに定義と実装を一緒に書いていきます。

ヘッダファイルの作り方

自分で作ったVec3クラスをvec3.hというヘッダファイルにまとめてみましょう。まずはクラスを書きます。

class Vec3 {
  ...
};

さらに以下のような呪文でクラスをはさみます。

#ifndef VEC3_H
#define VEC3_H

class Vec3 {
  ...
};

#endif

#ifndef VEC3_HVEC3_Hが定義されていなかったら#endifまでを実行せよという意味です。その下の#define VEC3_HVEC3_Hが定義されています。

こうすることで、複数のプログラムからこのヘッダファイルが読み込まれても、Vec3クラスの定義は必ず最初の一回だけ行われるようになります。クラスの多重定義はコンパイルエラーの原因になるので、エラーを防ぐためにこのような呪文を書きます。

自分で作ったヘッダファイルを読み込むには新たにmain.cppを作り、以下のようにします。

#include "vec3.h"


int main() {
  Vec3 v(1, 1, 1);
  ...
}

#include "~"は同じディレクトリにあるヘッダファイルを読み込むときに使われる書き方です。/usr/includeにあるヘッダファイルを読み込むには#include <~>を使用します。

練習問題

  • 前回の宿題で作ったクラスをヘッダファイルにまとめてみよう

解答例

#ifndef VEC3_H
#define VEC3_H
#include <iostream>
#include <cmath>
class Vec3 {
    public:
        double x;
        double y;
        double z;

        Vec3() {
            x = y = z = 0;
        };
        Vec3(double _x, double _y, double _z) {
            x = _x;
            y = _y;
            z = _z;
        };

        double length() const {
            return std::sqrt(x*x + y*y + z*z);
        };
        double length2() const {
            return x*x + y*y + z*z;
        };
};


Vec3 operator+(const Vec3& v1, const Vec3& v2) {
    return Vec3(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
}
Vec3 operator-(const Vec3& v1, const Vec3& v2) {
    return Vec3(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
}

Vec3 operator*(double k, const Vec3& v) {
    return Vec3(k * v.x, k * v.y, k * v.z);
}
Vec3 operator*(const Vec3& v, double k) {
    return k*v;
}

Vec3 operator/(double k, const Vec3& v) {
    return Vec3(k / v.x, k / v.y, k / v.z);
}
Vec3 operator/(const Vec3& v, double k) {
    return Vec3(v.x / k, v.y / k, v.z / k);
}

std::ostream& operator<<(std::ostream& stream, const Vec3& v) {
    stream << "(" << v.x << ", " << v.y << ", " << v.z << ")";
    return stream;
}

double dot(const Vec3& v1, const Vec3& v2) {
    return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
}
Vec3 cross(const Vec3& v1, const Vec3& v2) {
    return Vec3(v1.y*v2.z - v1.z*v2.y, v1.z*v2.x - v1.x*v2.z, v1.x*v2.y - v1.y*v2.x);
}

Vec3 normalize(const Vec3& v) {
    return v/v.length();
}
#endif
#ifndef RAY_H
#define RAY_H
#include "vec3.h"
class Ray {
    public:
        Vec3 origin;
        Vec3 direction;

        Ray(const Vec3& origin, const Vec3& direction) : origin(origin), direction(direction) {};
};
#endif
#ifndef SPHERE_H
#define SPHERE_H
#include "vec3.h"
#include "ray.h"
#include "hit.h"
class Sphere {
    public:
        Vec3 center;
        double radius;

        Sphere(const Vec3& center, double radius) : center(center), radius(radius) {};

        bool intersect(const Ray& ray, Hit& hit) const {
            double a = ray.direction.length2();
            double b = 2*dot(ray.direction, ray.origin - center);
            double c = (ray.origin - center).length2() - radius*radius;
            double D = b*b - 4*a*c;
            if(D < 0) return false;

            double t1 = (-b - std::sqrt(D))/(2*a);
            double t2 = (-b + std::sqrt(D))/(2*a);
            
            double t = t1;
            if(t < 0) {
                t = t2;
                if(t < 0) return false;
            }

            hit.t = t;
            hit.hitPos = ray.origin + t*ray.direction;
            hit.hitNormal = normalize(hit.hitPos - center);
            return true;
        };
};
#endif
  • Vec3をターミナルに出力するメンバ関数void Vec3::print() constを追加してみよう
class Vec3 {
...
  void print() const {
    printf("(%f, %f, %f)\n", x, y, z);
  };

PPM画像について

以下を参照

https://www.wikiwand.com/ja/PNM_(%E7%94%BB%E5%83%8F%E3%83%95%E3%82%A9%E3%83%BC%E3%83%9E%E3%83%83%E3%83%88)

PPM画像は次のようなフォーマットになっています。

P3
横幅 縦幅
255(階調数)
0 0 255
0 255 0
255 0 0
...

画像データは一行ずつRGBとして格納されています。画像データの最初は最も左上のピクセルに対応し、そこから横に並んでいく形になっています。Wikipediaの図を参照してください。

ImageMagickのインストール

ImageMagickは種々のフォーマットの画像を表示したり、変換したりすることのできるツールです。PPM画像を表示したり、PPM画像をPNG画像に変換するときに今後使っていくので、インストールしておきましょう。

Mac

brew install imagemagick

Ubuntu

sudo apt-get install imagemagick

PPM画像を表示するには、displayコマンドを使います。

display output.ppm

PPM画像をPNG画像に変換するには、convertコマンドを使います

convert output.ppm output.png

PNG画像に変換しておけば、自分のレンダリング画像をTwitterなどに上げることができますね!!

PPM画像を作ってみよう

上記のフォーマットでファイルを作成して、拡張子をppmとすればそのままppm画像になります。C++で上のようなフォーマットでファイル出力することをやってみましょう。

#include <iostream>
#include <fstream>

int main() {
  //ファイルを開く
  std::ofstream file("output.ppm");
  
  //このようにしてファイルに文字列を出力できる
  file << "P3" << std::endl;
  //横幅と縦幅は512
  file << "512 512" << std::endl;
  //階調数は255
  file << "255" << std::endl;
  
  //画像データ出力
  //全てのピクセルを(255, 255, 255)=白色にする
  for(int i = 0; i < 512; i++) {
    for(int j = 0; j < 512; j++) {
      file << "255 255 255" << std::endl;
    }
  }
  
  //使い終わったらファイルを閉じる
  file.close();
  return 0;
}

これを実行するとプログラムと同じディレクトリにoutput.ppmというppm画像が生成されます。

すべて白色だと面白くないので、グラデーションをつけてみましょう。

#include <iostream>
#include <fstream>

int main() {
  //ファイルを開く
  std::ofstream file("output.ppm");
  
  //このようにしてファイルに文字列を出力できる
  file << "P3" << std::endl;
  //横幅と縦幅は512
  file << "512 512" << std::endl;
  //階調数は255
  file << "255" << std::endl;
  
  //画像データ出力
  //全てのピクセルを(255, 255, 255)=白色にする
  for(int i = 0; i < 512; i++) {
    for(int j = 0; j < 512; j++) {
      int r = i/512.0 * 255; //赤色
      int g = j/512.0 * 255; //緑色
      int b = 255; //青色
      file << r << " " << g << " " << b << std::endl;
    }
  }
  
  //使い終わったらファイルを閉じる
  file.close();
  return 0;
}

実行するとこんな画像が生成されます。

実行結果
実行結果

画像クラスを作ってみよう

次のようなImageクラスがあるとこの先なにかと便利なので作ってみましょう。

メンバ変数

  • 横幅width
  • 縦幅height
  • Vec3(RGBを表す)の配列dataを持ち、それがそのまま画像データを表す。RGBの各成分は0から1までの実数で表すことにする。

メンバ関数

  • void setPixel(int i, int j, const Vec3& color)dataの位置(i,j)(i, j)colorを書き込む
  • Vec3 getPixel(int i, int j) constdataの位置(i,j)(i, j)のRGBを取得できる。
  • void ppm_output() constdataから"output.ppm"というPPM画像を生成する

dataは一次元配列で表すことにします。(i,j)(i, j)の画素にアクセスするにはdata[widthi+j]data[width*i + j]とすればよいです。RGBを0から255ではなく、0から1の範囲で表すことに注意してください。こうすると後々色の掛け算をするときに便利になります。

Imageクラスをimage.hというヘッダファイルにまとめてみましょう。

雛形

#ifndef IMAGE_H
#define IMAGE_H

#include "vec3.h"


class Image {
  public:
    int width; //横幅
    int height; //縦幅
    Vec3* data; //配列先頭へのポインタ
    
    //コンストラクタ 横幅と縦幅を受け取り、dataを初期化する
    Image(int _width, int _height) {
      width = _width;
      height = _height;
      data = new Vec3[width*height]; //配列を確保
    };
    
    //デストラクタ 使い終わったメモリを解放する
    ~Image() {
      delete[] data;
    };
    
    Vec3 getPixel(int i, int j) const {
      ...
    };
    void setPixel(int i, int j, const Vec3& color) {
      ...
    };
    
    void ppm_output() const {
      ...
    };
};
#endif

練習問題

  • 必要なところを自分で実装してみよう

解答例

#ifndef IMAGE_H
#define IMAGE_H
#include <iostream>
#include <fstream>
#include "vec3.h"
class Image {
    public:
        int width; //横幅
        int height; //縦幅
        Vec3* data;

        Image(int _width, int _height) {
            width = _width;
            height = _height;
            data = new Vec3[width*height];
        };
        ~Image() {
            delete[] data;
        };

        //(i, j)のRGBを返す
        Vec3 getPixel(int i, int j) const {
            return data[width*i + j];
        };
        //(i, j)にRGBを書き込む
        void setPixel(int i, int j, const Vec3& color) {
            data[width*i + j] = color;
        };

        //PPM画像を出力する
        void ppm_output() const {
            std::ofstream file("output.ppm");
            
            file << "P3" << std::endl;
            file << width << " " << height << std::endl;
            file << "255" << std::endl;
            for(int i = 0; i < height; i++) {
                for(int j = 0; j < width; j++) {
                    Vec3 color = 255*this->getPixel(j, i);
                    int r = (int)color.x;
                    int g = (int)color.y;
                    int b = (int)color.z;
                    file << r << " " << g << " " << b << std::endl;
                }
            }
        };
};
#endif

画像クラスで遊んでみよう

ピクセル座標(i,j)(i, j)に対して (icx)2+(jcy)2<r2(i - c_x)^2 + (j - c_y)^2 < r^2 を満たす点を白く塗ると、中心座標(cx,cy)(c_x, c_y)、半径rrの円を描画することができます。

#include "image.h"

int main() {
  Image img(512, 512); //512*512の画像
  
  int cx = 256;
  int cy = 256;
  int r = 100;
  
  for(int i = 0; i < 512; i++) {
    for(int j = 0; j < 512; j++) {
      int dx = i - cx;
      int dy = j - cy;
      if(dx*dx + dy*dy < r*r) {
        img.setPixel(i, j, Vec3(255, 255, 255)); //白色
      }
      else {
        img.setPixel(i, j, Vec3(0, 0, 0)); //黒色
      }
    }
  }
  img.ppm_output();
}
実行結果
実行結果

練習問題

  • 楕円を描画してみよう

座標系について

これからVec3をバリバリ使って3次元空間上での演算をたくさんしていくのですが、その前に3次元空間の座標系について説明しておきます。

これから扱う座標系は x方向が右、y方向が上、z方向が前 に対応することにします。この座標系の取り方はコンピューターグラフィクスでよく使われます。

座標系
座標系

カメラクラスを作ってみよう

カメラの役割は画素の座標(i,j)(i, j)に対応するレイを計算して返すことです。カメラクラスをcamera.hに実装して、画素に対応するレイを取得できるようにしましょう。

カメラクラスはメンバ変数として次のようなデータを持ちます。

  • Vec3 camPos カメラの位置
  • Vec3 camForward カメラの前方向

メンバ関数として次のような関数を持ちます。

  • Ray getRay(double u, double v) const (i,j)(i, j)に対応するレイを計算して返す

ピンホールカメラ

ピンホールカメラは最も単純なカメラモデルです。

https://www.wikiwand.com/ja/%E3%83%94%E3%83%B3%E3%83%9B%E3%83%BC%E3%83%AB%E3%82%AB%E3%83%A1%E3%83%A9

ピンホールカメラ
ピンホールカメラ

ピンホールという小さな穴があり、光はそこを通ってカメラセンサーに到達します。穴の大きさが無限小ならボケのない像を得ることができます。現実では無限小の穴のピンホールカメラを作ることは不可能ですが、コンピューターの中では再現することが可能です。

最初に作るカメラクラスとして、このピンホールカメラを作っていきましょう。

カメラの位置(センサ中心の位置)をcp\vec{c_p}, カメラの前方向をcf\vec{c_f}とします。ピンホールがセンサ中心から距離ttのところにあるなら、ピンホールの位置はcp+tcf\vec{c_p} + t\vec{c_f}と表せます。

レイの生成
レイの生成

ピンホールを通過するRayを作るためには画素(i,j)(i, j)のシーン空間上での位置p(i,j)\vec{p_{(i, j)}}が必要です。カメラの横方向をcr\vec{c_r}, カメラの上方向をcu\vec{c_u}とすると

p(i,j)=cp+ucr+vcu\vec{p_{(i, j)}} = \vec{c_p} + u\vec{c_r} + v\vec{c_u}

と表すことができます。ここでuuはセンサの左端で-1、右端で1を取る実数で、同様にvvは下端で-1、上端で1を取る実数です。

(i,j)(i, j)(u,v)(u, v)に変換するには、カメラセンサーの横幅をww、縦幅をhhとして

u=2.0iwwv=2.0jhh\begin{aligned} u &= \frac{2.0i - w}{w} \\ v &= \frac{2.0j - h}{h} \\ \end{aligned}

とすればよいです。

カメラの横方向cr\vec{c_r}と上方向cu\vec{c_u}は外積を使って前方向cf\vec{c_f}から計算することができます。

cr=(cf)×(0,1,0)cu=cf×cr\vec{c_r} = (-\vec{c_f}) \times (0, 1, 0) \\ \vec{c_u} = \vec{c_f} \times \vec{c_r}

カメラクラスの実装

上の式を使ってRayの始点と方向を計算していきます。方向を正規化 するのをお忘れずに!!!

#ifndef CAMERA_H
#define CAMERA_H
class Camera {
    public:
        Vec3 camPos;
        Vec3 camForward;
        Vec3 camRight;
        Vec3 camUp;

        Camera(const Vec3& camPos, const Vec3& camForward) : camPos(camPos), camForward(camForward) {
            camRight = -1 * normalize(cross(camForward, Vec3(0, 1, 0))); //カメラの横方向
            camUp = normalize(cross(camForward, camRight)); //カメラの上方向
        };

        Ray getRay(double u, double v) const {
            Vec3 pinhole = camPos + camForward; //ピンホールの位置
            Vec3 sensorPos = camPos + u*camRight + v*camUp; //(i, j)の位置
            return Ray(sensorPos, normalize(pinhole - sensorPos));
        };
};
#endif

実際に使うには次のようにします。

int main() {
    Img img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));
    
    for(int i = 0; i < 512; i++) {
        for(int j = 0; j < 512; 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);
        }
    }
}

ちゃんとRayが飛んでいるか確認してみましょう

#include "vec3.h"
#include "image.h"
#include "sphere.h"
#include "camera.h"


int main() {
    Image img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));

    for(int i = 0; i < 512; i++) {
        for(int j = 0; j < 512; 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);
            img.setPixel(i, j, (ray.direction + Vec3(1, 1, 1))/2.0);
        }
    }
    img.ppm_output();
}
実行結果
実行結果

はじめてのレイトレ

Cameraクラスが実装できたので、いよいよRayを飛ばして球を描画してみましょう。

シーン構成
シーン構成

Sphereを一つ作成し、飛ばしたRayがそれに衝突したら白色、しなかったら黒色にする処理にしてみましょう。

#include "vec3.h"
#include "image.h"
#include "sphere.h"
#include "camera.h"


int main() {
    Image img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));

    Sphere sphere(Vec3(0, 0, 0), 1.0);

    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);

            if(sphere.intersect(ray)) {
                img.setPixel(i, j, Vec3(1, 1, 1));
            }
            else {
                img.setPixel(i, j, Vec3(0, 0, 0));
            }
        }
    }

    img.ppm_output();
}

実行結果
実行結果
球が出ました!!

法線で色付けをする

球の衝突計算において、衝突したかどうかだけではなく、衝突距離・衝突位置・法線といった情報を返すようにしましょう。

法線
法線

C++では関数に複数の返り値をもたせることはできないので、参照渡しされた引数に直接返り値をセットするようにします。そのためにまずは衝突情報を格納するHitクラスを作りましょう。

#ifndef HIT_H
#define HIT_H
#include "vec3.h"
class Hit {
    public:
        double t; //衝突点までの距離
        Vec3 hitPos; //衝突位置
        Vec3 hitNormal; //衝突位置の法線

        Hit() {};
        Hit(double t, const Vec3& hitPos, const Vec3& hitNormal) : t(t), hitPos(hitPos), hitNormal(hitNormal) {};
};
#endif

これをhit.hとして保存します。

Sphereのintersectを次のような形に変えます。

bool intersect(const Ray& ray, Hit& hit)

練習問題

  • Sphere::intersectが正しく衝突情報をhitに格納するように書き直してみよう

あとは法線を使って色付けするだけです。

#include "vec3.h"
#include "image.h"
#include "sphere.h"
#include "camera.h"


int main() {
    Image img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));

    Sphere sphere(Vec3(0, 0, 0), 1.0);

    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(sphere.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();
}

実行結果
実行結果
先ほどより面白い見た目になりましたね!!

影をつける

簡単な陰影づけを行ってみましょう。

光源の方向をl\vec{l}、物体の法線をn\vec{n}とします。光の強さIII=max{ln,0}I = \max \{\vec{l}\cdot\vec{n}, 0\} で計算してみましょう。これは ランバートの拡散反射モデル と呼ばれます。

Lambert Shading
Lambert Shading
#include <algorithm>
#include "vec3.h"
#include "image.h"
#include "sphere.h"
#include "camera.h"


int main() {
    Image img(512, 512);
    Camera cam(Vec3(0, 0, -3), Vec3(0, 0, 1));

    Sphere sphere(Vec3(0, 0, 0), 1.0);

    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(sphere.intersect(ray, hit)) {
                double I = std::max(dot(lightDir, hit.hitNormal), 0.0);
                img.setPixel(i, j, I*Vec3(1, 1, 1));
            }
            else {
                img.setPixel(i, j, Vec3(0, 0, 0));
            }
        }
    }

    img.ppm_output();
}

実行結果
実行結果
影がついて3DCG感が出てきました!!

練習問題

問1

Rayの方向d\vec{d}を法線n\vec{n}で反射させた方向を返す関数Vec3 reflect(const Vec3& d, const Vec3& n)を作れ。

反射ベクトル
反射ベクトル

問2

Phongの反射モデルを使って陰影計算を行ってみましょう。

Phongの反射モデルでは、光源の方向をl\vec{l}、Rayの方向をd\vec{d}、衝突点の法線をn\vec{n}とし、l-\vec{l}n\vec{n}で反射させたものをr\vec{r}とすると、光の強さII

I=kdmax{ln,0}+ksmax{(d)r,0}αI = k_d\max\{ \vec{l}\cdot\vec{n}, 0 \} + k_s\max \{ (-\vec{d})\cdot\vec{r} , 0\}^\alpha

で与えられます。

Phong Shading
Phong Shading

kd,ksk_d, k_skd+ks=1k_d + k_s = 1を満たす[0,1][0, 1]の区間に含まれる実数です。α\alphaはハイライトの強さを制御する0以上の実数です。これらのパラメーターを色々と変化させて画像を生成してみてください。