问题五十七：怎么用ray tracing画translational sweeping图形

57.1 概述

我们这里考虑的translational sweeping是XOZ平面的封闭曲线沿着y轴平移得到的图形。类似于如下图形：我们把这种图形称为“prism”。我们这里的封闭曲线是由多段三次b-spline曲线组成的。

光线和prism相交有如下图中所示的几种情况：

第一种情况：光线和顶面所在的平面相交的交点在底面的投影在封闭曲线内，如红色光线所示。

第二种情况：光线和底面所在的平面相交的交点在封闭曲线内，如绿色光线所示。

第三种情况：光线和顶面所在的平面相交的交点在底面的投影、光线和底面所在平面相交的交点都不在封闭曲线内，如紫色光线所示。

所以，我们“求根”的算法以此分为如下几个步骤：

第一步：求光线R(t)和base plane、cap plane的交点t0、t1；

第二步：将光线R(t)投影到base plane得到R'(t)；

第三步：求出投影光线R’(t)与曲线的所有交点t’；

第四步：处理t0、t1有在曲线内的情况；

第五步：处理t0、t1都不在曲线内的情况；

57.2 求根

57.2.1求光线R(t)和base plane、cap plane的交点t0、t1

57.2.2将光线R(t)投影到base plane得到R'(t)

57.2.3 怎么用三次b-spline曲线段形成封闭曲线？

参考PPT截图如下：

我们接下来要画的图形，有17个控制点，则有14段（17-3=14）曲线段。示意图如下：

形成封闭曲线，则需要控制点收尾相连。所以，这17个控制点中前三个和后三个是重复的。

57.2.4 求R’(t)和封闭曲线的所有交点

要求R’(t)和封闭曲线的所有交点，我们先求R’(t)和每一曲线段的交点，然后将这14段曲线段对应的交点汇总则得到R’(t)和封闭曲线的所有交点。

接下来，我们求R’(t)和每一曲线段的交点。

57.2.5 判断t0和t1在base plane的投影是否在封闭曲线内

我们已经求得R’(t)与曲线的所有交点。示意如下图。

t0、t1在base plane的投影必定都在R’(t)上，如上图中三角形示意。

怎么判断点是否在曲线内呢？

若大于t0/t1的交点数为偶数，则t0/t1在曲线外；

若大于t0/t1的交点数为奇数，则t0/t1在曲线内；

我们现在不说t0或者t1，改说两者中较小的t0_t1_smaller，两者中较大的t0_t1_bigger。

若t0_t1_smaller在曲线内：

如“57.1”图示的红色光线。

最终要求的t是：t0_t1_smaller。

若t0_t1_bigger在曲线内：

如“57.1”图示的绿色光线。

最终要求的t是：所有有效交点中最小的那个对应的t

57.2.6处理t0、t1都不在曲线内的情况

针对t0、t1都不在曲线内的情况，如“57.1”图示的蓝色光线和紫色光线，我们直到蓝色光线没有撞击到prism，紫色光线撞击到了prism。

但是，具体怎么判断光线是否撞击到prism呢？

参考前人的成果，判断方式如下：

所有有效交点存在一个t’满足：“式子九”

|t’-t0|*Rd·n<h

其中Rd为光线方向向量，n为base plane的单位法向量，h为prism的高度。

将t’、t0都代入光线方程：R(t)=R0+Rd*t

然后，我们可以得到R(t0)-R(t’)= (R0+Rd*t0)-( R0+Rd*t’)=(t0-t’)*Rd

然后，再和n点积（n为单位向量，所以相当于长度乘以cos(theta)），得到高度。

具体如下图所示：

注意到：t1对应的那个高度为h。

上图中存在t’3和t’4满足“式子九”。

最终要求的t是：满足“式子九”的最小的那个t’

注意：

这个时候还需要判断对应的y值是否在[base_y, cap_y]之间。若在，该t即为所求；若不在，则光线没有撞击到prism。

具体做法：将t代入光线方程“式子一”求出y值，然后判断一下即可。

上图中为t’3.

57.3 求法向量

若t0_t1_smaller在曲线内：

法向量即为cap plane的法向量：（0, 1, 0）

若t0_t1_bigger在曲线内

或者

t0、t1都不在曲线内

这两种情况的根是通过投影光线和曲线相交求得，对应的法向量即为曲线在该处的切向量的垂直向量。

我们可以通过对“式子六”求对u的微分来求得切向量。“式子十”

57.4 看C++代码实现

----------------------------------------------vec3.h ------------------------------------------

vec3.h

bool matrix_4_4_multiply_4_1(const float matrix1[4][4], const float matrix2[4][1], float (&result)[4][1]);

----------------------------------------------vec3.cpp ------------------------------------------

vec3.cpp

bool matrix_4_4_multiply_4_1(const float matrix1[4][4], const float matrix2[4][1], float (&result)[4][1]) {for (int k=0; k<1; k++) {for (int i=0; i<4; i++) {result[i][k] = 0.0;for (int j=0; j<4; j++) {result[i][k] = result[i][k] + matrix1[i][j]*matrix2[j][k];}}}return true;
}

----------------------------------------------sweeping_translational.h ------------------------------------------

sweeping_translational.h

#ifndef SWEEPING_TRANSLATIONAL_H
#define SWEEPING_TRANSLATIONAL_H#include <hitable.h>
#include <log.h>class sweeping_translational : public hitable
{public:sweeping_translational() {}sweeping_translational(vec3 *cp, float by, float cy, material *m){base_y = by;cap_y = cy;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_x[17], points_z[17], b_points_x[4][1], b_points_z[4][1], b_result_x[4][1], b_result_z[4][1];int b_num = 0;float min_x, max_x, min_z, max_z;for (int i=0; i<17; i++) {points_x[i] = cp[i].x();points_z[i] = cp[i].z();}min_x = points_x[0];max_x = points_x[0];min_z = points_z[0];max_z = points_z[0];for (int i=0; i<17; i++) {if (min_x > points_x[i]) {min_x = points_x[i];}if (max_x < points_x[i]) {max_x = points_x[i];}if (min_z > points_z[i]) {min_z = points_z[i];}if (max_z < points_z[i]) {max_z = points_z[i];}}sweeping_bl = vec3(min_x, base_y, min_z);sweeping_bh = vec3(max_x, cap_y, max_z);
/*获取包围prism的最小的长方体的参数*/for (int i=0; i<14; i=i+1) {b_points_x[0][0] = points_x[i];b_points_x[1][0] = points_x[i+1];b_points_x[2][0] = points_x[i+2];b_points_x[3][0] = points_x[i+3];b_points_z[0][0] = points_z[i];b_points_z[1][0] = points_z[i+1];b_points_z[2][0] = points_z[i+2];b_points_z[3][0] = points_z[i+3];matrix_4_4_multiply_4_1(matrix_t, b_points_x, b_result_x);matrix_4_4_multiply_4_1(matrix_t, b_points_z, b_result_z);for (int j=0; j<4; j++) {matrix_c_x[j][b_num] = ((fabs(b_result_x[j][0])<1e-6)? (0.0):(b_result_x[j][0]));matrix_c_z[j][b_num] = ((fabs(b_result_z[j][0])<1e-6)? (0.0):(b_result_z[j][0]));
/*这里是求“式子六”中提到的[Cx][Cz]*/}b_num ++;}}virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;float matrix_c_x[4][14], matrix_c_z[4][14];vec3 sweeping_bl, sweeping_bh;float base_y, cap_y;material *ma;
};#endif // SWEEPING_TRANSLATIONAL_H

----------------------------------------------sweeping_translational.cpp ------------------------------------------

sweeping_translational.cpp

#include "sweeping_translational.h"#include <iostream>
#include <limits>
#include "float.h"bool ray_hit_box_st(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_translational::hit(const ray& r, float t_min, float t_max,hit_record& rec) const {
#if SWEEPING_TRANSLATIONAL_LOG == 1std::cout <<"-------------sweeping_translational::hit---1-------------" <<endl;
#endif // SWEEPING_TRANSLATIONAL_LOGfloat t_near, t_far;if (ray_hit_box_st(r,sweeping_bl, sweeping_bh, t_near, t_far)) {float x0 = r.origin().x();float y0 =r.origin().y();float z0 =r.origin().z();float xd =r.direction().x();float yd =r.direction().y();float zd =r.direction().z();if (fabs(xd) < 1e-6) {xd = 1e-6;}if (fabs(yd) < 1e-6) {yd = 1e-6;}if (fabs(zd) < 1e-6) {zd = 1e-6;}
/*非常接近0的值都用1e-6替代*/float cos_theta =dot(r.direction(), vec3(0, -1, 0))/(r.direction().length());float t0, t1,t0_t1_smaller, t0_t1_bigger, a3, a2, a1, a0, *roots3_base, roots_base[21][3],temp0, temp1, temp2, root[5], u, check_y;
/*
root[0]:    0表示没有根；1表示1个根
root[1]:   存放根t
root[2]:    1表示t0_t1_smaller在曲线内；2表示t0_t1_bigger在曲线内；3表示t0、t1都不在曲线内。（不同情况，法向量的求解方式不一样）
root[3]:    存放曲线段编号（法向量求解时需要）
root[4]:    存放根t对应的u的值（法向量求解时需要）
*/roots_base[0][0] =0.0;root[0] = 0.0;root[2] = 0.0;root[3] = -1.0;root[4] = -1.0;int num_roots_base= 0;intnum_t0_t1_smaller = 0;intnum_t0_t1_bigger = 0;intnum_roots_base_valid = 0;t0 = (base_y - y0) / yd;t1 = (cap_y - y0) / yd;
/*对应“式子三”*/t0_t1_smaller = (t0>t1)? t1:t0;t0_t1_bigger =(t0>t1)? t0:t1;for (int i=0;i<14; i++) {
/*14段曲线段。求解14个一元三次方程。对应“式子七”*/a3 =zd*matrix_c_x[3][i] - xd*matrix_c_z[3][i];a2 =zd*matrix_c_x[2][i] - xd*matrix_c_z[2][i];a1 =zd*matrix_c_x[1][i] - xd*matrix_c_z[1][i];a0 =zd*matrix_c_x[0][i] - xd*matrix_c_z[0][i] + xd*z0 - zd*x0;roots3_base =roots_cubic_equation(a3, a2, a1, a0);if(roots3_base[0] != 0.0) {for (intj=1; j<(int(roots3_base[0])+1); j++) {u =roots3_base[j];if((u>=0.0) && (u<=1.0)) {roots_base[num_roots_base+1][0] =(matrix_c_x[0][i]+matrix_c_x[1][i]*u+matrix_c_x[2][i]*u*u+matrix_c_x[3][i]*u*u*u-x0)/xd;
/*对应“式子八”，将u转化为t*/roots_base[num_roots_base+1][1] = i;roots_base[num_roots_base+1][2] = u;
/*同时保存曲线段编号和u，为了后面求法向量*/num_roots_base ++;}}roots_base[0][0] = float(num_roots_base);}delete []roots3_base;}if (roots_base[0][0] != 0.0) {
/*对所有根进行从小到大排序*/for (int i=1;i<int(roots_base[0][0]); i++) {for (intj=i+1; j<int(roots_base[0][0])+1; j++) {if(roots_base[i][0] > roots_base[j][0]) {temp0 =roots_base[i][0];roots_base[i][0] = roots_base[j][0];roots_base[j][0] = temp0;temp1 = roots_base[i][1];roots_base[i][1] = roots_base[j][1];roots_base[j][1] = temp1;temp2 = roots_base[i][2];roots_base[i][2] = roots_base[j][2];roots_base[j][2] = temp2;}}}for (int i=1;i<(int(roots_base[0][0])+1); i++) {
/*判断大于t0_t1_smaller的根的个数，为了判断其是否在曲线内，对应57.2.5*/if(roots_base[i][0] > (t0_t1_smaller)) {num_t0_t1_smaller = int(roots_base[0][0])+1-i;break;}}if((num_t0_t1_smaller%2) != 0) {root[1] =t0_t1_smaller;root[0] =1.0;root[2] =1.0;}else {for (inti=1; i<(int(roots_base[0][0])+1); i++) {
/*判断大于t0_t1_bigger的根的个数，为了判断其是否在曲线内，对应57.2.5*/if(roots_base[i][0] > (t0_t1_bigger)) {num_t0_t1_bigger = int(roots_base[0][0])+1-i;break;}}if((num_t0_t1_bigger%2) != 0) {root[1] = roots_base[1][0];//fabs(cos_theta);root[0] = 1.0;root[2] = 2.0;root[3] = roots_base[1][1];root[4] = roots_base[1][2];}else {
/*如下处理t0,t1都不在曲线内的情况，对应“57.2.6”*/for(int i=1; i<(int(roots_base[0][0])+1); i++) {if(fabs(roots_base[i][0]-t0_t1_bigger)*fabs(dot(r.direction(), vec3(0, -1,0)))<fabs(cap_y-base_y)) {
/*对应“式子九”*/
//                           if (roots_base[i][0] > t0_t1_smaller) {num_roots_base_valid = i;
/*保存大于t0_t1_smaller的最小的那个根的编号*/break;}}if(num_roots_base_valid != 0) {root[1] = roots_base[num_roots_base_valid][0];//fabs(cos_theta);check_y = y0 + yd * root[1];

if ((check_y > base_y) &&(check_y < cap_y)) {

/*这里就是“57.2.6”中提到的需要注意的地方，必须加这个限制条件。后续我们对比有无这个限制条件时的输出图形*/

                                root[0] = 1.0;root[2] = 3.0;root[3] = roots_base[num_roots_base_valid][1];root[4] = roots_base[num_roots_base_valid][2];}}}}if (root[0] != 0.0) {if ((root[1] < t_max) && (root[1] > t_min)) {rec.t = root[1];rec.p = r.point_at_parameter(rec.t);if (root[2] == 1.0) {
/*若t0_t1_smaller在曲线内：法向量即为cap plane的法向量：（0, 1, 0）*/rec.normal = vec3(0, 1, 0);}else {float nx, nz, nu;int num = int(root[3]);nu = root[4];nx = matrix_c_x[1][num]+2.0*matrix_c_x[2][num]*nu+3.0*matrix_c_x[3][num]*nu*nu;nz = matrix_c_z[1][num]+2.0*matrix_c_z[2][num]*nu+3.0*matrix_c_z[3][num]*nu*nu;
/*对应“式子十”*/rec.normal = unit_vector(vec3(nx, 0, 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;}else {return false;}}else {return false;}}else {return false;}}else {return false;}
}

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

main.cpp

        vec3 ctrl_points[17] = {vec3(-4.0,  0.0, -1.0), vec3(-3.0,  0.0, -4.0), vec3(-2.0,  0.0, -4.0), vec3(-1.0,  0.0,  -1.0),vec3( 1.0,  0.0, -1.0), vec3( 1.0,  0.0, -5.0), vec3( 3.0,  0.0, -5.0), vec3( 4.0,  0.0,   1.0),vec3( 2.0,  0.0,  2.0), vec3( 0.5,  0.0,  2.0), vec3(-2.0,  0.0,  1.0), vec3(-2.0,  0.0,  -2.0),vec3(-2.5,  0.0, -3.0), vec3(-3.0,  0.0, -1.0), vec3(-4.0,  0.0, -1.0), vec3(-3.0,  0.0,  -4.0), vec3(-2.0,  0.0, -4.0)};hitable *list[1];list[0] = new sweeping_translational(ctrl_points, 0.0, 3.0, new lambertian(vec3(1.0, 0.0, 0.0)));hitable *world = new hitable_list(list,1);vec3 lookfrom(0.0001, 20, 20);vec3 lookat(0, 1.5, 0);float dist_to_focus = (lookfrom - lookat).length();float aperture = 0.0;camera cam(lookfrom, lookat, vec3(0,1,0), 30, float(nx)/float(ny), aperture, 0.7*dist_to_focus);

输出图片如下：（前边是没有加“57.2.6”中提到的限制条件时；后边是加了该限制条件时）

(0.0001,20, 20)

(0.0001,20, 10)

(0.0001,20, 5)

(-5, 20, 10)