レイをトレースしたい　第三回　球の描画

前回まででレイの定義ができたのでいよいよ球の描画をしてみたいと思います。

点p(a,b,c)を中心とする球の方程式

$$(x-a)^2 + (y-b)^2 + (z-c)^2 = C^2$$

に対して、位置ベクトル$r = (x,y,z)$の内積を用いて表し、定数Cを移項すると

$$(r-p)\cdot (r-p) - C^2 = 0$$

と表せます。これとカメラの位置$a$及びレイのベクトル$b$によって表される直線の方程式

$$ r = a + tb$$

が交点を持てば良いので、これを連立して整理することで、

$$b\cdot b t^2 + 2*b \cdot (a-p) t + (a-p) \cdot (a-p) -C^2  = 0$$

の方程式がtについて実数解を持つとき、レイと球が交差することがわかります。また、このとき２次方程式解の公式よりtの値を求めることができ、代入することで交点の座標を求める事もできます。

これを利用して球の衝突判定をするプログラムがこのようになります

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

衝突したときにはtの値を返すようにしています。これを用いて衝突判定をして球を描くプログラムが次のようになります。

```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){
    double a = ray.dot(ray);
    double b = ray.dot(camera_pos-center)*2.0;
    double c = (camera_pos-center).dot(camera_pos-center)-size*size;
    double D = b*b-4*a*c;
    if (D < 0){
        return -1.0;
    }else{
        return (-b-sqrt(D))/(2.0*a);
    }
}


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

}

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*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,150)));
    camera cam ;
    Vector3d sp_pos;
    sp_pos << 0,0,0;
    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 t = hit_sphere(sp_pos,ray,cam.position,1);
            Vector3d hit_point = cam.position + ray*t;
            if(t>0){
                pic[i][j][0] = 0;
                pic[i][j][1] = 255;
                pic[i][j][2] = 0;
            }else{
                pic[i][j][0] = 255;
                pic[i][j][1] = 255;
                pic[i][j][2] = 255;
            }

        }
    }
    write_ppm(pic);
    return 0;
}
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244678/b2615110-0d9b-fb22-94de-b2b338531eac.png)

ついでに法線の導出もしてしまいます。
球の場合には法線はかなり単純で球の中心からレイの交点までのベクトルを取れば良いです。それを実際に表示すると次のようになります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/244678/3c390a91-bdc7-58cf-039c-3ffb4d786e47.png)