问题六十：怎么用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 == 1std::cout << "-------------sweeping_rotational::hit---1-------------" << endl;
#endif // SWEEPING_ROTATIONAL_LOGfloat 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);

输出图形如下：

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

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