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 == 1
std::cout <<"-------------sweeping_translational::hit---1-------------" <<endl;
#endif // SWEEPING_TRANSLATIONAL_LOG
float 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)
扫一扫
308
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
