レイをトレースしたい　第四回　球を増やす

球が一個だけだと流石に寂しいので、複数個の球を描画する方法を考えてみます。

複数の物体があるとき、まずその全てに対して衝突判定をする必要が出てきます。そこまでは特に問題はないと思いますが、複数の物体を書くときにはその奥行関係が必要となります。

この奥行きの判定をする方法は様々あると思いますが、私は交差判定方程式を解いた結果出てくるtの値の大小で判定しようと思います。tの値が小さいほど、レイを伸ばさずに物体と衝突するためです。

また、加えてレイが何とも衝突しなかった場合の処理をするために、衝突しなかった場合はtに負の大きな値を入れることにします。こうすることで衝突したかしてないかの判断はtの符号で判定できるのに加えて、tの値の大小で判定先のオブジェクトを判断するときも絶対値をとれば衝突していないオブジェクトに関しても同じように大小関係で判断できるからです。

すなわち、オブジェクトA,B,Cに対してレイを飛ばして、Aのみ交差せず、B,Cにはある値$t_B,t_C$でレイが交差するとしましょう。交差しない$t_A$については十分な負の値(たとえば-1e-30)が代入されるとします

1. $t_A > 0 \ || \ t_B > 0 \  || \  t_C > 0$を判定する。これによって交差するオブジェクトがあるか調べる。
2. abs(t)の中で最も小さいtをとる物体が衝突対象のオブジェクト

とこのように判定することにします。この衝突判定は当然min()のような配列に対しての最小値を求めるような関数を用いたり、オブジェクト自身に色を付けるメンバをもたせることによってスッキリかけますが、今回はプリミティブなif文による実装でひとまず済ませたいと思います。

今回のコードは以下のとおりです。

```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){
    int height = 512;
    int width = 512;
    double xstep = 2.0/width;
    double ystep = 2.0/height;
    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;
    for (int i = 0; i<height; ++i){
        for(int j = 0; j < width; ++ j){
            double u = xstep*j-1.0;
            double v = 1.0-ystep*i;
            Vector3d ray = cam.ray(u,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);
                    pic[i][j][0] = saturate(norm(0))*255;
                    pic[i][j][1] = saturate(norm(1))*255;
                    pic[i][j][2] = saturate(-norm(2))*255;
                }else{
                    pic[i][j][0] = 0;
                    pic[i][j][1] = 255;
                    pic[i][j][2] = 0;
                }
            }

        }
    }
    write_ppm(pic);
    return 0;
}

```

特に以下の部分が複数物体の衝突判定をしている部分になっています。

```cpp
            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);
                    pic[i][j][0] = saturate(norm(0))*255;
                    pic[i][j][1] = saturate(norm(1))*255;
                    pic[i][j][2] = saturate(-norm(2))*255;
                }else{
                    pic[i][j][0] = 0;
                    pic[i][j][1] = 255;
                    pic[i][j][2] = 0;
                }
            }
```

出力画像はこんな感じです。ちゃんと前後関係が反映されているのがわかると思います。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244678/c2b24c1b-31b8-278c-971c-ef001bffbd30.png)