问题六十三：怎么用ray tracing画sphere sweeping图形

63.1 概述

Translational sweeping、conic sweeping、rotational sweeping都是任意曲线以某种固定方式移动后形成的图形。

接下来，我们要画的sphere sweeping则是球沿着任意曲线形成的图形：球心沿着三次b-spline曲线移动，在移动过程中球的半径也随之改变。

示意图如下：

球心的变化轨迹：c(u)；半径的变化轨迹：r(u)

Sphere sweeping上任一点P的坐标为（x,y,z），则sphere sweeping可以表示为：

63.2 求光线和sphere sweeping的交点

如上图，两个红球与光线相切，则光线与sphere sweeping的交点为：

光线与两个红球之间的无数个绿球的交点的总和。但是，真正有效的交点只有最前和最后两个交点（如上如绿色点t(u2),t(u3)表示）。

数学上，不难理解，u2、u3为函数t(u)的极值点（即导数值为零的点）：

63.3 求法向量

在求得交点坐标之后，同时我们也已经求得交点对应的球的球心坐标：

63.4 看C++代码实现

----------------------------------------------sweeping_ sphere.h ------------------------------------------

sweeping_ sphere.h

/*这个初始化文件中，需要注意的是：半径r的变化规律也是三次b样条曲线，所以其对应的系数求法和球心对应的系数求法是一样的。*/#ifndef SWEEPING_SPHERE_H
#define SWEEPING_SPHERE_H#include <hitable.h>
#include <log.h>class sweeping_sphere : public hitable
{public:sweeping_sphere() {}sweeping_sphere(vec3 *ccp, vec3 *rcp, material *m){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_cx[8], points_cy[8], points_cz[8], points_r[8];float b_points_cx[4][1], b_points_cy[4][1], b_points_cz[4][1], b_points_r[4][1];float b_result_cx[4][1], b_result_cy[4][1], b_result_cz[4][1], b_result_r[4][1];int b_num = 0;float min_cx, max_cx, min_cy, max_cy, min_cz, max_cz, max_r;for (int i=0; i<8; i++) {//8points_cx[i] = ccp[i].x();points_cy[i] = ccp[i].y();points_cz[i] = ccp[i].z();points_r[i]  = rcp[i].x();}min_cx = points_cx[0];max_cx = points_cx[0];min_cy = points_cy[0];max_cy = points_cy[0];min_cz = points_cz[0];max_cz = points_cz[0];max_r  = points_r[0];for (int i=0; i<8; i++) {//8if (min_cx > points_cx[i]) {min_cx = points_cx[i];}if (max_cx < points_cx[i]) {max_cx = points_cx[i];}if (min_cy > points_cy[i]) {min_cy = points_cy[i];}if (max_cy < points_cy[i]) {max_cy = points_cy[i];}if (min_cz > points_cz[i]) {min_cz = points_cz[i];}if (max_cz < points_cz[i]) {max_cz = points_cz[i];}if (max_r < points_r[i]) {max_r = points_r[i];}}sweeping_bl = vec3(min_cx-max_r, min_cy-max_r, min_cz-max_r);sweeping_bh = vec3(max_cx+max_r, max_cy+max_r, max_cz+max_r);for (int i=0; i<5; i=i+1) {//5b_points_cx[0][0] = points_cx[i];b_points_cx[1][0] = points_cx[i+1];b_points_cx[2][0] = points_cx[i+2];b_points_cx[3][0] = points_cx[i+3];b_points_cy[0][0] = points_cy[i];b_points_cy[1][0] = points_cy[i+1];b_points_cy[2][0] = points_cy[i+2];b_points_cy[3][0] = points_cy[i+3];b_points_cz[0][0] = points_cz[i];b_points_cz[1][0] = points_cz[i+1];b_points_cz[2][0] = points_cz[i+2];b_points_cz[3][0] = points_cz[i+3];b_points_r[0][0]  = points_r[i];b_points_r[1][0]  = points_r[i+1];b_points_r[2][0]  = points_r[i+2];b_points_r[3][0]  = points_r[i+3];matrix_4_4_multiply_4_1(matrix_t, b_points_cx, b_result_cx);matrix_4_4_multiply_4_1(matrix_t, b_points_cy, b_result_cy);matrix_4_4_multiply_4_1(matrix_t, b_points_cz, b_result_cz);matrix_4_4_multiply_4_1(matrix_t, b_points_r, b_result_r);for (int j=0; j<4; j++) {matrix_c_cx[j][b_num] = ((fabs(b_result_cx[j][0])<1e-6)? (0.0):(b_result_cx[j][0]));matrix_c_cy[j][b_num] = ((fabs(b_result_cy[j][0])<1e-6)? (0.0):(b_result_cy[j][0]));matrix_c_cz[j][b_num] = ((fabs(b_result_cz[j][0])<1e-6)? (0.0):(b_result_cz[j][0]));matrix_c_r[j][b_num]  = ((fabs(b_result_r[j][0])<1e-6)? (0.0):(b_result_r[j][0]));}b_num ++;}}virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;float matrix_c_cx[4][5], matrix_c_cy[4][5], matrix_c_cz[4][5], matrix_c_r[4][5];vec3 sweeping_bl, sweeping_bh;material *ma;
};#endif // SWEEPING_SPHERE_H

----------------------------------------------sweeping_ sphere.cpp ------------------------------------------

sweeping_ sphere.cpp

#include "sweeping_sphere.h"bool sweeping_sphere::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
#if SWEEPING_SPHERE_LOG == 1std::cout << "-------------sweeping_sphere::hit---1-------------" << endl;
#endif // SWEEPING_SPHERE_LOGfloat t_near, t_far;if (ray_hit_box_general(r, sweeping_bl, sweeping_bh, t_near, t_far)) {float xx0 = r.origin().x();float yy0 = r.origin().y();float zz0 = r.origin().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;}float f20 = xxd*xxd + yyd*yyd + zzd*zzd;if (fabs(f20) < 1e-6) {f20 = 1e-6;}float f13, f12, f11, f10, aa2, aa1, aa0, cc5, cc4, cc3, cc2, cc1, cc0, dd4, dd3, dd2, dd1, dd0;float f06, f05, f04, f03, f02, f01, f00, bb5, bb4, bb3, bb2, bb1, bb0;float cx3, cx2, cx1, cx0;float cy3, cy2, cy1, cy0;float cz3, cz2, cz1, cz0;float cr3, cr2, cr1, cr0;float gg10[11], hh10[11], ff10[11], ee10[11], roots[11], roots_t[31][3], uuu, f0d, f1d, ttt, temp0, temp1, temp2;float tol = 1e-6;int num_roots_t = 0;float px, py, pz, cx, cy, cz;int num;for (int i=0; i<5; i++) {//5cx3 = matrix_c_cx[3][i];cx2 = matrix_c_cx[2][i];cx1 = matrix_c_cx[1][i];cx0 = matrix_c_cx[0][i];cy3 = matrix_c_cy[3][i];cy2 = matrix_c_cy[2][i];cy1 = matrix_c_cy[1][i];cy0 = matrix_c_cy[0][i];cz3 = matrix_c_cz[3][i];cz2 = matrix_c_cz[2][i];cz1 = matrix_c_cz[1][i];cz0 = matrix_c_cz[0][i];cr3 = matrix_c_r[3][i];cr2 = matrix_c_r[2][i];cr1 = matrix_c_r[1][i];cr0 = matrix_c_r[0][i];f13 = -2*(xxd*cx3 + yyd*cy3 + zzd*cz3);f12 = -2*(xxd*cx2 + yyd*cy2 + zzd*cz2);f11 = -2*(xxd*cx1 + yyd*cy1 + zzd*cz1);f10 = -2*(xxd*cx0 + yyd*cy0 + zzd*cz0) + 2*(xx0*xxd+yy0*yyd+zz0*zzd);aa2 = 3*f13;aa1 = 2*f12;aa0 = f11;
/*上方系数对应“式子63.9.2”*/f06 = cx3*cx3 + cy3*cy3 + cz3*cz3 - cr3*cr3;f05 = 2*(cx2*cx3 + cy2*cy3 + cz2*cz3 - cr2*cr3);f04 = 2*(cx1*cx3 + cy1*cy3 + cz1*cz3 - cr1*cr3) + (cx2*cx2 + cy2*cy2 + cz2*cz2 - cr2*cr2);f03 = 2*(cx1*cx2 + cy1*cy2 + cz1*cz2 - cr1*cr2) + 2*((cx0-xx0)*cx3 + (cy0-yy0)*cy3 + (cz0-zz0)*cz3 - cr0*cr3);f02 = 2*((cx0-xx0)*cx2 + (cy0-yy0)*cy2 + (cz0-zz0)*cz2 - cr0*cr2) + (cx1*cx1 + cy1*cy1 + cz1*cz1 - cr1*cr1);f01 = 2*((cx0-xx0)*cx1 + (cy0-yy0)*cy1 + (cz0-zz0)*cz1 - cr0*cr1);f00 = (cx0-xx0)*(cx0-xx0) + (cy0-yy0)*(cy0-yy0) + (cz0-zz0)*(cz0-zz0) - cr0*cr0;bb5 = 6*f06;bb4 = 5*f05;bb3 = 4*f04;bb2 = 3*f03;bb1 = 2*f02;bb0 = f01;
/*上方系数对应“式子63.9.3”*/hh10[0]  = f20*(bb5*bb5);hh10[1]  = f20*(2*bb4*bb5);hh10[2]  = f20*(2*bb3*bb5 + bb4*bb4);hh10[3]  = f20*(2*bb2*bb5 + 2*bb3*bb4);hh10[4]  = f20*(2*bb1*bb5 + 2*bb2*bb4 + bb3*bb3);hh10[5]  = f20*(2*bb0*bb5 + 2*bb1*bb4 + 2*bb2*bb3);hh10[6]  = f20*(2*bb0*bb4 + 2*bb1*bb3 + bb2*bb2);hh10[7]  = f20*(2*bb0*bb3 + 2*bb1*bb2);hh10[8]  = f20*(2*bb0*bb2 + bb1*bb1);hh10[9]  = f20*(2*bb0*bb1);hh10[10] = f20*(bb0*bb0);
/*上方系数对应“式子63.9.4”*/cc5 = f13*aa2;cc4 = f13*aa1 + f12*aa2;cc3 = f13*aa0 + f12*aa1 + f11*aa2;cc2 = f12*aa0 + f11*aa1 + f10*aa2;cc1 = f11*aa0 + f10*aa1;cc0 = f10*aa0;ff10[0]  = bb5*cc5;ff10[1]  = bb5*cc4 + bb4*cc5;ff10[2]  = bb5*cc3 + bb4*cc4 + bb3*cc5;ff10[3]  = bb5*cc2 + bb4*cc3 + bb3*cc4 + bb2*cc5;ff10[4]  = bb5*cc1 + bb4*cc2 + bb3*cc3 + bb2*cc4 + bb1*cc5;ff10[5]  = bb5*cc0 + bb4*cc1 + bb3*cc2 + bb2*cc3 + bb1*cc4+ bb0*cc5;ff10[6]  = bb4*cc0 + bb3*cc1 + bb2*cc2 + bb1*cc3 + bb0*cc4;ff10[7]  = bb3*cc0 + bb2*cc1 + bb1*cc2 + bb0*cc3;ff10[8]  = bb2*cc0 + bb1*cc1+ bb0*cc2;ff10[9]  = bb1*cc0 + bb0*cc1;ff10[10] = bb0*cc0;
/*上方系数对应“式子63.9.5”*/dd4 = aa2*aa2;dd3 = 2*aa1*aa2;dd2 = 2*aa0*aa2 + aa1*aa1;dd1 = 2*aa0*aa1;dd0 = aa0*aa0;ee10[0]  = dd4*f06;ee10[1]  = dd4*f05 + dd3*f06;ee10[2]  = dd4*f04 + dd3*f05 + dd2*f06;ee10[3]  = dd4*f03 + dd3*f04 + dd2*f05 + dd1*f06;ee10[4]  = dd4*f02 + dd3*f03 + dd2*f04 + dd1*f05 + dd0*f06;ee10[5]  = dd4*f01 + dd3*f02 + dd2*f03 + dd1*f04 + dd0*f05;ee10[6]  = dd4*f00 + dd3*f01 + dd2*f02 + dd1*f03 + dd0*f04;ee10[7]  = dd3*f00 + dd2*f01 + dd1*f02 + dd0*f03;ee10[8]  = dd2*f00 + dd1*f01 + dd0*f02;ee10[9]  = dd1*f00 + dd0*f01;ee10[10] = dd0*f00;
/*上方系数对应“式子63.9.6”*/for (int k=0; k<11; k++) {gg10[k] = hh10[k] - ff10[k] + ee10[k];
/*上方系数对应“式子63.9.7”*/if (fabs(gg10[k])<1e-6) {gg10[k] = 1e-6;}}roots_equation_10th(gg10, 0, 1, tol, roots);if (int(roots[0]) >= 1) {for (int j=1; j<(int(roots[0])+1); j++) {uuu = roots[j];//the real roots u which is the parameter of cubic b-spline.if (fabs(uuu)>1e-6) {f0d = bb5*pow(uuu,5) + bb4*pow(uuu,4) + bb3*pow(uuu,3) + bb2*pow(uuu,2) + bb1*uuu + bb0;f1d = aa2*pow(uuu,2) + aa1*uuu + aa0;if (fabs(f1d)>1e-6) {ttt = -f0d/f1d;
//note that b-spline is the trajectory of sphere center and the distance from hit-point to origin of ray is determined by this fomula
/*上方系数对应“式子63.6”*/roots_t[num_roots_t+1][0] = ttt;
//store the distance from hit-point to origin.if (roots_t[num_roots_t+1][0] > 0) {roots_t[num_roots_t+1][1] = i;
//store the b-spline curve number that the corresponding sphere center lies onroots_t[num_roots_t+1][2] = roots[j];
//store the b-spline parameter of the center of the hitted sphere.num_roots_t ++;}}}}}}roots_t[0][0] = float(num_roots_t);if (roots_t[0][0] != 0.0) {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);px = rec.p.x();py = rec.p.y();pz = rec.p.z();num = int(roots_t[1][1]);uuu = roots_t[1][2];cx = matrix_c_cx[0][num] + matrix_c_cx[1][num]*uuu + matrix_c_cx[2][num]*uuu*uuu + matrix_c_cx[3][num]*uuu*uuu*uuu;cy = matrix_c_cy[0][num] + matrix_c_cy[1][num]*uuu + matrix_c_cy[2][num]*uuu*uuu + matrix_c_cy[3][num]*uuu*uuu*uuu;cz = matrix_c_cz[0][num] + matrix_c_cz[1][num]*uuu + matrix_c_cz[2][num]*uuu*uuu + matrix_c_cz[3][num]*uuu*uuu*uuu;rec.normal = unit_vector(vec3(px-cx, py-cy, pz-cz));
/*上方系数对应“式子63.10”*/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_sc[8] = {vec3(-4.0, -4.0,  8.0), vec3(-4.0,  0.0, -4.0),vec3(-4.0,  8.0, -4.0), vec3(-4.0,  8.0,  4.0),vec3( 4.0,  8.0,  4.0), vec3( 4.0,  8.0, -4.0),vec3( 4.0,  0.0, -4.0), vec3( 4.0, -4.0,  8.0)};vec3 ctrl_points_sr[8] = {vec3( 2.0,  0.0,  0.0), vec3( 2.0,  0.0,  0.0),vec3( 1.0,  0.0,  0.0), vec3( 0.5,  0.0,  0.0),vec3( 0.5,  0.0,  0.0), vec3( 1.0,  0.0,  0.0),vec3( 2.0,  0.0,  0.0), vec3( 2.0,  0.0,  0.0)};hitable *list[1];list[0] = new sweeping_sphere(ctrl_points_sc, ctrl_points_sr, new lambertian(vec3(1.0, 0.0, 0.0)));hitable *world = new hitable_list(list,1);vec3 lookfrom(0, 4, 20);vec3 lookat(0, 4, 0);float dist_to_focus = (lookfrom - lookat).length();float aperture = 0.0;camera cam(lookfrom, lookat, vec3(0,1,0), 40, float(nx)/float(ny), aperture, 0.7*dist_to_focus);

输出图形如下：

vec3 lookfrom(0, 4, 20);

vec3 lookfrom(20, 4, 20);

vec3 lookfrom(10, 20, 10);