レイをトレースしたい　第六回　diffuse material(Lambert)

いよいよマテリアルの実装をしていきたいと思います。
マテリアルの実装には様々な方法があると思いますが、大きく分けると計算負荷の軽い方法と重たい方法があります。

diffuse　materialは入射してきた光が当確率で表面半球上の任意の点方向へ反射するという仮定によって構成されます。

この等確率の反射を表現する方法として、一つは等確率の反射を、入射光が全方向に対して均等に分配されるという考え方に基づいて計算するモデルがあります。これはいわゆるLambertのモデルに基づく計算方法です。この方法では入射光の強さを求める計算以外が必要なく軽量に計算でき、また実装もかなり簡単であることがメリットです。ただしこの場合は物体表面で反射した光が更に周囲の物体に届く、いわば二次反射以降の影響を計算することは難しくなります。(私が知らないだけかもしれませんが)

それに対して、アンチエイリアスで使った複数のレイを投げる手法に近い方法があります。モンテカルロレイトレーシングと呼ばれる手法です。
この場合は1pxあたりに複数のレイを投げ、衝突のたびに実際にランダムに方向を判定し再びレイをトレースするという手法になります。この手法ではレイが間で通る全ての物体を考慮できるため環境光等を表現することに長けていますが、反射回数及びレイの本数によって負荷がどんどんと増大していくためこれを考える必要が出てきます。また、かなり十分なレイの本数を投げないとノイズが乗るという問題があります。

モンテカルロレイトレーシング実装もかなり大変なのでとりあえずLambertモデルの実装をしたいと思います。

Lambertモデルに関してはこちらを参考に実装します。
https://wgld.org/d/glsl/g010.html

球の法線と定義したレイベクトルの内積をとってそれを光の強さとして採用するだけになっています。

```cpp
#include <bits/stdc++.h>
#include "Eigen/Core"
using namespace std;
using namespace Eigen;
/*
Vector3d
Matrix3d
(i,j)の形で要素アクセス
A.array().sin()で要素ごと演算
A.norm()
A.normalize()
a.dot(b)
*/

const double PI = 3.141592653589793238;
const double PI2 = 6.28318530718;
const double eps = 1e-6;




double hit_sphere(Vector3d center, Vector3d ray, Vector3d camera_pos, double size){
    Vector3d oc = camera_pos-center;
    double a = ray.dot(ray);
    double b = ray.dot(oc)*2.0;
    double c = (oc).dot(oc)-size*size;
    double D = b*b-4*a*c;
    if (D < 0){
        return -1e30;
    }else{
        return (-b-sqrt(D))/(2.0*a);
    }
}

Vector3d norm_sphere(Vector3d center, Vector3d hit_point){
    Vector3d temp = (hit_point-center).normalized();
    return temp;
}

double saturate(double d){
    double ans = 0;

    if(d < 0){
        ans = 0.0;
    }else if(d > 1){
        ans = 1.0;
    }else{
        ans = d;
    }
    return ans;
}

void write_ppm(vector<vector<vector<int>>> pic){
    ofstream out;
    out.open("test.pnm");
    out << "P3" << endl;
    out << pic.size() << " " << pic[0].size() << endl;
    out << 255 << endl;
    for (int i = 0 ; i < pic.size(); ++i){
        for (int j = 0; j < pic[0].size(); ++j){
            out << pic[i][j][0] << " " << pic[i][j][1] << " " << pic[i][j][2] << endl;
        }
    }

}

Matrix3d Rod_rot(Vector3d axis, double angle){
    double theta = angle*PI/180;
    double c = cos(theta);
    double s = sin(theta);
    double x = axis(0);
    double y = axis(1);
    double z = axis(2);
    
    Matrix3d ret;

    ret << c+(1-c)*x*x    , x*y*(1-c)-z*s , x*z*(1-c)+y*s,
           y*x*(1-c)+z*s  , c+(1-c)*y*y   , y*z*(1-c)-x*s,
           z*x*(1-c)-y*s  , z*y*(1-c)+x*s , c+z*z*(1-c);
    return ret;
}

class camera{
    public:
        camera();
        Vector3d direction;
        Vector3d position;
        double angle;
        double rad_ang_half;
    Vector3d ray(double x, double y){
        Vector3d temp;
        temp << sin(rad_ang_half)*x, sin(rad_ang_half)*y, cos(rad_ang_half);
        temp = temp.normalized();
        return temp;
    }
    private:
        
};

camera::camera(){
    angle = 60;
    rad_ang_half = angle*0.5*PI/180;
    direction << 0,0,1;
    position << 0,0,-3;
}

int main(void){
    random_device rd{};
    int seed = rd();
    mt19937 mt(rd());
    uniform_real_distribution<> dist(0.01,0.99);
    int height = 512;
    int width = 512;
    double xstep = 2.0/width;
    double ystep = 2.0/height;
    const int ray_per_px = 10;
    vector<vector<vector<int>>> pic;
    pic = vector<vector<vector<int>>>(height,vector<vector<int>>(width,vector<int>(3,255)));
    camera cam ;
    Vector3d sp_pos;
    sp_pos << 0,0,2;
    Vector3d sp_pos2;
    sp_pos2 << 0,-20,10;
    Vector3d light_Dir;
    light_Dir << -1,1,-1;
    light_Dir.normalize();
    Vector3d albedo;
    albedo << 1.0,0,0;
    for (int i = 0; i<height; ++i){
        for(int j = 0; j < width; ++ j){
            vector<double> buf = vector<double>(3,0);
            double u = xstep*j-1.0;
            double v = 1.0-ystep*i;
            for (int k = 0; k < ray_per_px; ++ k){
                double dx = dist(mt)-0.5;
                double dy = dist(mt)-0.5;
                double temp_u = u+xstep*dx;
                double temp_v = v+ystep*dy;
                Vector3d ray = cam.ray(temp_u,temp_v);
                double t1 = hit_sphere(sp_pos,ray,cam.position,1);
                double t2 = hit_sphere(sp_pos2,ray,cam.position,20);
                Vector3d hit_point = cam.position + ray*t1;
                Vector3d hit_point2 = cam.position + ray*t2;
                if(t1 > 0 || t2 > 0){
                if (abs(t1) < abs(t2)){
                    Vector3d norm = norm_sphere(sp_pos,hit_point);
                    double ref_pow = light_Dir.dot(norm);
                    buf[0] += saturate(ref_pow)*albedo(0)*255;
                    buf[1] += saturate(ref_pow)*albedo(1)*255;
                    buf[2] += saturate(ref_pow)*albedo(2)*255;
                }else{
                    buf[0] += 0;
                    buf[1] += 255;
                    buf[2] += 0;
                }
                }else{
                    buf[0] += 255;
                    buf[1] += 255;
                    buf[2] += 255;
                } 

            }
            pic[i][j][0] = buf[0]/ray_per_px;
            pic[i][j][1] = buf[1]/ray_per_px;
            pic[i][j][2] = buf[2]/ray_per_px;
        }
    }
    write_ppm(pic);
    return 0;
}

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244678/2c7d3e78-899c-b08b-1030-7f8af438fc9a.png)

なんとなくのライティングはできていると思います(小さい方の球だけ)