问题六十：怎么用ray tracing画回旋体（rotational sweeping / revolution）

60.1 概述

回旋体，大概是长这个样子：

 

回旋体是指曲线（称为“基本曲线”）围绕y轴转一圈得到的图形。

（基本曲线是由多段b-spline曲线段连接而成）

 

这里先强调一下：

上图中有两个坐标系：熟悉的空间三维xyz坐标系和曲线的局部二维uv坐标系。

 

对于回旋体，有如下特点：

上面这两个坐标系有这种关系：

对于回旋体上任意一点P(x,y,z)都对应一条曲线，则对应一个曲线的局部二维uv坐标系，则对应这个局部坐标系中一个坐标(u,v)。

60.2 求根

60.3 求法向量

60.4 看C++代码实现

----------------------------------------------sweeping_rotational.h ------------------------------------------

sweeping_rotational.h

/*这个文件和sweeping_translational.h差别不大，可以对比参照来理解*/

#ifndef SWEEPING_ROTATIONAL_H
#define SWEEPING_ROTATIONAL_H

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

class sweeping_rotational : public hitable
{
    public:
        sweeping_rotational() {}
        sweeping_rotational(vec3 cen, vec3 *cp, material *m){
            center = cen;
            ma = m;
            float matrix_t_6[4][4] = {{ 1,  4,  1, 0},
                                      {-3,  0,  3, 0},
                                      { 3, -6,  3, 0},
                                      {-1,  3, -3, 1}};
            float matrix_t[4][4];
            for (int i=0; i<4; i++) {
                for (int j=0; j<4; j++) {
                    matrix_t[i][j] = matrix_t_6[i][j] / 6.0;
                }
            }
            float points_u[6], points_v[6], b_points_u[4][1], b_points_v[4][1], b_result_u[4][1], b_result_v[4][1];
            int b_num = 0;
            float max_u, min_v, max_v;
            for (int i=0; i<6; i++) {
                points_u[i] = cp[i].x();
                points_v[i] = cp[i].y();
            }
            max_u = points_u[0];
            min_v = points_v[0];
            max_v = points_v[0];
            for (int i=0; i<6; i++) {
                if (max_u < points_u[i]) {
                    max_u = points_u[i];
                }
                if (min_v > points_v[i]) {
                    min_v = points_v[i];
                }
                if (max_v < points_v[i]) {
                    max_v = points_v[i];
                }
            }
            sweeping_bl = vec3(-max_u+center.x(), min_v+center.y(), -max_u+center.z());
            sweeping_bh = vec3( max_u+center.x(), max_v+center.y(),  max_u+center.z());

            for (int i=0; i<3; i=i+1) {
                b_points_u[0][0] = points_u[i];
                b_points_u[1][0] = points_u[i+1];
                b_points_u[2][0] = points_u[i+2];
                b_points_u[3][0] = points_u[i+3];
                b_points_v[0][0] = points_v[i];
                b_points_v[1][0] = points_v[i+1];
                b_points_v[2][0] = points_v[i+2];
                b_points_v[3][0] = points_v[i+3];
                matrix_4_4_multiply_4_1(matrix_t, b_points_u, b_result_u);
                matrix_4_4_multiply_4_1(matrix_t, b_points_v, b_result_v);
                for (int j=0; j<4; j++) {
                    matrix_c_u[j][b_num] = ((fabs(b_result_u[j][0])<1e-6)? (0.0):(b_result_u[j][0]));
                    matrix_c_v[j][b_num] = ((fabs(b_result_v[j][0])<1e-6)? (0.0):(b_result_v[j][0]));
                }
                b_num ++;
            }

        }
        virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;
        float matrix_c_u[4][3], matrix_c_v[4][3];
        vec3 sweeping_bl, sweeping_bh, center;
        material *ma;
};

#endif // SWEEPING_ROTATIONAL_H

----------------------------------------------sweeping_rotational.cpp ------------------------------------------

sweeping_rotational.cpp

#include "sweeping_rotational.h"

#include <iostream>
#include <limits>
#include "float.h"

bool ray_hit_box_sr(const ray& r, const vec3& vertex_l, const vec3& vertex_h, float& t_near, float& t_far) {
/*罪恶般地又抄了一遍这个函数，光线撞击包围回旋体外面的最小盒子*/
        t_near = (numeric_limits<float>::min)();
        t_far = (numeric_limits<float>::max)();
        vec3 direction = r.direction();
        vec3 origin = r.origin();
        vec3 bl = vertex_l;
        vec3 bh = vertex_h;
        float array1[6];

        if(direction.x() == 0) {
            if((origin.x() < bl.x()) || (origin.x() > bh.x())) {
                return false;
            }
            array1[0] = (numeric_limits<float>::min)();
            array1[1] = (numeric_limits<float>::max)();
        }
        if(direction.y() == 0) {
            if((origin.y() < bl.y()) || (origin.y() > bh.y())) {
                return false;
            }
            array1[2] = (numeric_limits<float>::min)();
            array1[3] = (numeric_limits<float>::max)();
        }
        if(direction.z() == 0) {
            if((origin.z() < bl.z()) || (origin.z() > bh.z())) {
                return false;
            }
            array1[4] = (numeric_limits<float>::min)();
            array1[5] = (numeric_limits<float>::max)();
        }

        if((direction.x() != 0) && (direction.y() != 0) && (direction.z() != 0)) {
            array1[0] = (bl.x()-origin.x())/direction.x();
            array1[1] = (bh.x()-origin.x())/direction.x();
            array1[2] = (bl.y()-origin.y())/direction.y();
            array1[3] = (bh.y()-origin.y())/direction.y();
            array1[4] = (bl.z()-origin.z())/direction.z();
            array1[5] = (bh.z()-origin.z())/direction.z();
        }

        for (int i=0; i<6; i=i+2){
            if(array1[i] > array1[i+1]) {
                float t = array1[i];
                array1[i] = array1[i+1];
                array1[i+1] = t;
            }
            if(array1[i] >= t_near) {t_near = array1[i];}
            if(array1[i+1] <= t_far) {t_far = array1[i+1];}
            if(t_near > t_far) {
                return false;
            }
            if(t_far < 0) {
                return false;
            }
        }
        if (t_near != t_near) {
            t_near = t_near * 1;
        }
        return true;
}

bool sweeping_rotational::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
#if SWEEPING_ROTATIONAL_LOG == 1
        std::cout << "-------------sweeping_rotational::hit---1-------------" << endl;
#endif // SWEEPING_ROTATIONAL_LOG
        float t_near, t_far;
        if (ray_hit_box_sr(r, sweeping_bl, sweeping_bh, t_near, t_far)) {
            float xx0 = r.origin().x() - center.x();
            float yy0 = r.origin().y() - center.y();
            float zz0 = r.origin().z() - center.z();
            float xxd = r.direction().x();
            float yyd = r.direction().y();
            float zzd = r.direction().z();
            if (fabs(xxd) < 1e-6) {
                xxd = 1e-6;
            }
            if (fabs(yyd) < 1e-6) {
                yyd = 1e-6;
            }
            if (fabs(zzd) < 1e-6) {
                zzd = 1e-6;
            }
/*对应“式子60.3”*/
            float aaa = xxd*xxd + zzd*zzd;
            float bbb = 2 * ((xx0*xxd+zz0*zzd)*yyd - aaa*yy0);
            float ccc = (xx0*yyd-xxd*yy0)*(xx0*yyd-xxd*yy0) + (zz0*yyd-zzd*yy0)*(zz0*yyd-zzd*yy0);
            float ddd = yyd*yyd;
            float ss6[7], roots[7], roots_t[19][3], yyv, temp0, temp1, temp2;
            float tol = 1e-6;
            int num_roots_t = 0;

            for (int i=0; i<3; i++) {
/*“基本曲线”是有三段b-spline曲线段组成，分别对每一段求交。对应“式子60.5”*/
                ss6[0] = aaa*matrix_c_v[3][i]*matrix_c_v[3][i] - ddd*matrix_c_u[3][i]*matrix_c_u[3][i];
                ss6[1] = aaa*2*matrix_c_v[2][i]*matrix_c_v[3][i] - ddd*2*matrix_c_u[2][i]*matrix_c_u[3][i];
                ss6[2] = aaa*(matrix_c_v[2][i]*matrix_c_v[2][i]+2*matrix_c_v[1][i]*matrix_c_v[3][i]) -
                      ddd*(matrix_c_u[2][i]*matrix_c_u[2][i]+2*matrix_c_u[1][i]*matrix_c_u[3][i]);
                ss6[3] = aaa*2*(matrix_c_v[0][i]*matrix_c_v[3][i]+matrix_c_v[1][i]*matrix_c_v[2][i]) + bbb*matrix_c_v[3][i] -
                      ddd*2*(matrix_c_u[0][i]*matrix_c_u[3][i]+matrix_c_u[1][i]*matrix_c_u[2][i]);
                ss6[4] = aaa*(matrix_c_v[1][i]*matrix_c_v[1][i]+2*matrix_c_v[0][i]*matrix_c_v[2][i]) + bbb*matrix_c_v[2][i] -
                      ddd*(matrix_c_u[1][i]*matrix_c_u[1][i]+2*matrix_c_u[0][i]*matrix_c_u[2][i]);
                ss6[5] = aaa*2*matrix_c_v[0][i]*matrix_c_v[1][i] + bbb*matrix_c_v[1][i]  - ddd*2*matrix_c_u[0][i]*matrix_c_u[1][i];
                ss6[6] = aaa*matrix_c_v[0][i]*matrix_c_v[0][i] + bbb*matrix_c_v[0][i]  - ddd*matrix_c_u[0][i]*matrix_c_u[0][i] + ccc;
                roots_equation_6th(ss6, 0, 1, tol, roots);
                for (int j=1; j<(int(roots[0])+1); j++) {
                    yyv = matrix_c_v[0][i]+matrix_c_v[1][i]*roots[j]+matrix_c_v[2][i]*roots[j]*roots[j]+matrix_c_v[3][i]*roots[j]*roots[j]*roots[j];
                    roots_t[num_roots_t+1][0] = (yyv-yy0)/yyd;
                    roots_t[num_roots_t+1][1] = i;
                    roots_t[num_roots_t+1][2] = roots[j];
                    num_roots_t ++;
                }
            }
            roots_t[0][0] = float(num_roots_t);

            if (roots_t[0][0] != 0.0) {
/*对所有的交点t从小到大排序，所以最小的值跑到了最前面*/
                for (int i=1; i<int(roots_t[0][0]); i++) {
                    for (int j=i+1; j<int(roots_t[0][0])+1; j++) {
                        if (roots_t[i][0] > roots_t[j][0]) {
                            temp0 = roots_t[i][0];
                            roots_t[i][0] = roots_t[j][0];
                            roots_t[j][0] = temp0;
                            temp1 = roots_t[i][1];
                            roots_t[i][1] = roots_t[j][1];
                            roots_t[j][1] = temp1;
                            temp2 = roots_t[i][2];
                            roots_t[i][2] = roots_t[j][2];
                            roots_t[j][2] = temp2;
                        }
                    }
                }
                if ((roots_t[1][0] < t_max) && (roots_t[1][0] > t_min)) {
                    rec.t = roots_t[1][0];
                    rec.p = r.point_at_parameter(rec.t);
                    float nx = rec.p.x() - center.x();
                    float nz = rec.p.z() - center.z();
                    int num = int(roots_t[1][1]);
                    float sss = roots_t[1][2];
                    float nu = -(matrix_c_v[1][num]+2.0*matrix_c_v[2][num]*sss+3.0*matrix_c_v[3][num]*sss*sss);
                    float nv =  (matrix_c_u[1][num]+2.0*matrix_c_u[2][num]*sss+3.0*matrix_c_u[3][num]*sss*sss);
                    float ny = sqrt(nx*nx+nz*nz)*nv/nu - center.y();
/*求法向量，对应“式子60.6”*/
                    rec.normal = unit_vector(vec3(nx, ny, nz));
                    if(dot(r.direction(), rec.normal) > 0) {
                        rec.normal = - rec.normal;
                    }
                    rec.mat_ptr = ma;
                    rec.u = -1.0;
                    rec.v = -1.0;
                    return true;
                }
            }
        }
        return false;
}

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

main.cpp

        vec3 ctrl_points[6] = {vec3(-1.0,  5.0,  0.0), vec3( 2.0,  4.0,  0.0),
                               vec3( 2.0,  1.0,  0.0), vec3(-0.5,  1.0,  0.0),
                               vec3( 1.5, -3.0,  0.0), vec3( 3.0,  0.0,  0.0)};
/*“基本曲线的控制点，6个点对应3条曲线段”*/

        hitable *list[1];
        list[0] = new sweeping_rotational(vec3( 0.0, 0.0,  0.0), ctrl_points, new lambertian(vec3(1.0, 0.0, 0.0)));
        hitable *world = new hitable_list(list,1);

        vec3 lookfrom(0, 10, 20);
        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);

 

输出图形如下：

 

        vec3 lookfrom(0, 10, 20);

 

        vec3 lookfrom(10, 10, 10);

来张组合图：

        vec3 ctrl_points[6] = {vec3(-1.0,  5.0,  0.0), vec3( 2.0,  4.0,  0.0),
                               vec3( 2.0,  1.0,  0.0), vec3(-0.5,  1.0,  0.0),
                               vec3( 1.5, -3.0,  0.0), vec3( 3.0,  0.0,  0.0)};
        hitable *list[4];
        list[0] = new sweeping_rotational(vec3( 0.0, 0,  0.0), ctrl_points, new dielectric(4.5));
        list[1] = new sweeping_rotational(vec3(-5.0, 0,  0.0), ctrl_points, new lambertian(vec3(1.0, 0.0, 0.0)));
        list[2] = new sweeping_rotational(vec3( 5.0, 0,  0.0), ctrl_points, new lambertian(vec3(1.0, 0.0, 1.0)));
        list[3] = new sphere(vec3(0,-103,-1), 100, new metal(vec3(0.0, 1.0, 0.5), 0.0), 0);
        hitable *world = new hitable_list(list,4);

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

输出图形如下：

 

一个玻璃材质，两个漫射材质，最下面的大球是镜面材质的

 

三个全是红色漫射材质时的图形如下：