问题三十五： 怎么用ray tracing画二次曲面(quadratic surfaces)（2）——单页双曲面、双页双曲面、椭圆锥面、椭圆柱面

35.2.1 数学推导

单页双曲面、双页双曲面、椭圆锥面、椭圆柱面。

这四个二次曲面方程共同形式：

但是，注意到，这些曲面都是开放曲面。在画图时，需要限制曲面的范围（以免曲面覆盖整个画面）。

我们在这里是限制曲面在y轴方向距离中心点的长度为height_y（引入该参数）。

所以，我们可以在根据实根t求得交点坐标后，对交点坐标作如下判断：

((rec.p.y()-center.y()) >-height_half_y) && ((rec.p.y()-center.y()) < height_half_y)

 

这里，要特别注意：

之前，球面的处理方式是：

所以，其一：我们需要对两个实根进行排序（先处理小的）

 

另外，由于，是开放曲面，也就是，光线有可能撞击到曲面的正反两面，所以，对于撞击点处的标准化之后的法向量，我们需要做如下判断：

                        if (dot(rec.normal,r.direction()) > 0) {

                            rec.normal =-rec.normal;

                        }//（法向量决定着反射光线和折射光线）

 

还有，由于我们引入了height_y参数来限制曲面的高度，但是，我们要注意到：小实根对应的交点超出高度范围时（之前的一贯做法：小根不在t_min、t_max范围，就直接return false），大实根是有可能在高度范围的（而且，如果在范围的话，光线撞击的曲面的内表面，这时的法向量是需要反向的）。

35.2.2 看C++代码实现

----------------------------------------------quadratic.h ------------------------------------------

quadratic.h

 

#ifndef QUADRATIC_H
#define QUADRATIC_H

#include <hitable.h>
#include "material.h"
#include "log.h"


class quadratic : public hitable
{
    public:
        quadratic() {}
        quadratic(vec3 cen, float a, float b, float c, int s1, int s2, int hy, material *m) : center(cen), intercept_x(a), intercept_y(b), intercept_z(c), sign1(s1), sign2(s2), height_half_y(hy), ma(m) {}
/*
(x-xc)^2/a^2 + s1*(y-yc)^2/b^2 + (z-zc)^2/c^2 = s2
s1=  -1,s2=  1:  hyperboloid of one sheet
s1=  -1,s2= -1:  hyperboloid of two sheets
s1=  -1,s2=  0:  elliptic cone
s1=   0,s2=  1:  elliptic cylinder
*/
        virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;
        vec3 center;
        float intercept_x;
        float intercept_y;
        float intercept_z;
        int sign1, sign2, height_half_y;
        material *ma;
};

#endif // QUADRATIC_H

----------------------------------------------quadratic.cpp ------------------------------------------

quadratic.cpp

 

#include "quadratic.h"

#include <iostream>
using namespace std;

bool quadratic::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
        float ab_square = intercept_x*intercept_x*intercept_y*intercept_y;
        float bc_square = intercept_y*intercept_y*intercept_z*intercept_z;
        float ac_square = sign1*intercept_x*intercept_x*intercept_z*intercept_z;
        float abc_square = sign2*intercept_x*intercept_x*intercept_y*intercept_y*intercept_z*intercept_z;

        vec3 inter_square = vec3(bc_square, ac_square, ab_square);
        vec3 rd_square = vec3(r.direction().x()*r.direction().x(),
                              r.direction().y()*r.direction().y(),
                              r.direction().z()*r.direction().z());
        float A = dot(inter_square, rd_square);
        vec3 r0_c = r.origin() - center;
        vec3 r0_c_rd = vec3(r0_c.x()*r.direction().x(),
                            r0_c.y()*r.direction().y(),
                            r0_c.z()*r.direction().z());
        float B = 2*dot(r0_c_rd, inter_square);
        vec3 r0_c_square = vec3(r0_c.x()*r0_c.x(),
                                r0_c.y()*r0_c.y(),
                                r0_c.z()*r0_c.z());
        float C = dot(r0_c_square, inter_square) - abc_square;
        float temp, temp1, temp2;
        vec3 pc;

        if(A == 0) {
            if (B == 0) {
                return false;
            }
            else {
                temp = -C/B;
                if (temp < t_max && temp > t_min) {
                    rec.t = temp;
                    rec.p = r.point_at_parameter(rec.t);
                    if (((rec.p.y()-center.y()) > -height_half_y) && ((rec.p.y()-center.y()) < height_half_y)) {
                        pc = rec.p - center;
                        rec.normal = unit_vector(vec3(2*bc_square*pc.x(), 2*ac_square*pc.y(), 2*ab_square*pc.z()));
                        if (dot(rec.normal, r.direction()) > 0) {
                            rec.normal = -rec.normal;
                        }
                        rec.mat_ptr = ma;
                        return true;
                    }
                    else {
                        return false;
                    }
                }
            }
        }
        else {
            float discriminant = B*B - 4*A*C;
            if (discriminant >= 0) {
                temp1 = (-B - sqrt(discriminant)) / (2.0*A);
                temp2 = (-B + sqrt(discriminant)) / (2.0*A);
                if (temp1 > temp2) {//make sure that temp1 is smaller than temp2
                    temp = temp1;
                    temp1 = temp2;
                    temp2 = temp;
                }
/*对两个实根进行排序，t1<t2，先处理t1*/
                if (temp1 < t_max && temp1 > t_min) {
                    rec.t = temp1;
                    rec.p = r.point_at_parameter(rec.t);
                    if (((rec.p.y()-center.y()) > -height_half_y) && ((rec.p.y()-center.y()) < height_half_y)) {//限制曲面的高度
                        pc = rec.p - center;
                        rec.normal = unit_vector(vec3(2*bc_square*pc.x(), 2*ac_square*pc.y(), 2*ab_square*pc.z()));
						/*任何时候，都最好将法向量进行标准化*/
                        if (dot(rec.normal, r.direction()) > 0) {
                            rec.normal = -rec.normal;
                        }
/*如果光线和前面内表面相交，需要将法向量反向*/
                        rec.mat_ptr = ma;
                        return true;
                    }
                    else {
//                        return false;
/*这里要特别注意，一个实根不在高度范围内，不能直接return false，还需在判断另一个实根*/
                    }
                }
                if (temp2 < t_max && temp2 > t_min) {
                    rec.t = temp2;
                    rec.p = r.point_at_parameter(rec.t);
                    if (((rec.p.y()-center.y()) > -height_half_y) && ((rec.p.y()-center.y()) < height_half_y)) {
                        pc = rec.p - center;
                        rec.normal = unit_vector(vec3(2*bc_square*pc.x(), 2*ac_square*pc.y(), 2*ab_square*pc.z()));
                        if (dot(rec.normal, r.direction()) > 0) {
                            rec.normal = -rec.normal;
                        }
                        rec.mat_ptr = ma;
                        return true;
                    }
                    else {
//                        return false;
                    }
                }
            }
            return false;
        }
}

----------------------------------------------main.cpp ------------------------------------------

main.cpp

 

        hitable *list[5];
        list[0] = new sphere(vec3(0.0,-100.5,-1), 100, new lambertian(vec3(0.8, 0.8, 0.0)));
        list[1] = new quadratic(vec3(-2.5, 1.5, 0), 0.5, 0.7, 1, -1, 1, 1, new lambertian(vec3(0.0, 0.1, 0.5)));
        list[2] = new quadratic(vec3(-0.4, 1.5, 0), 0.15, 0.2, 0.2, -1, -1, 1, new lambertian(vec3(0.3, 0.1, 0.5)));
        list[3] = new quadratic(vec3(1.2, 1.5, 0), 0.5, 0.7, 1, -1, 0, 1, new lambertian(vec3(0.6, 0.1, 0.5)));
        list[4] = new quadratic(vec3(2.8, 1.5, 0), 0.5, 0.7, 1, 0, 1, 1, new lambertian(vec3(0.9, 0.1, 0.5)));

        hitable *world = new hitable_list(list,5);

        vec3 lookfrom(0,5,10);
        vec3 lookat(0,1,0);
        float dist_to_focus = (lookfrom - lookat).length();
        float aperture = 0.0;
        camera cam(lookfrom, lookat, vec3(0,1,0), 20, float(nx)/float(ny), aperture, 0.7*dist_to_focus);

输出图片：

从左到右依次是：单页双曲面、双页双曲面、椭圆锥面、椭圆柱面