问题四十一：怎么用ray tracing画任意圆柱面（generalized cylinder）

我们之前在“35.2”章节中画过椭圆柱面：

我们还在“36.4”章节中画过圆柱面的Inverse Mapping图：

但是，这些柱面都是：底面与ZOX平面平行，中心线与Y轴（方向向量（0，1，0））平行。

这一章节，我们将画空间中任意位置、中心线方向任意的圆柱面。

41.1 数学推导

思路：根据圆柱面的中心线建立新的坐标系vuw，在vuw坐标系中，圆柱底面与wov平面平行，中心线与u轴平行；然后，将圆柱面中心点坐标、光线起点和方向向量的坐标转换到vuw坐标系。

这样就可以用之前画“特殊”圆柱面的方法来画了。

关键是建立新的坐标系vuw。

设圆柱面中心线的方向向量u的坐标为（Xu, Yu, Zu）。

ZOX平面内会有两个向量与u垂直：（-Zu, 0,Xu）和（Zu, 0, -Xu）

我们取v=（-Zu, 0, Xu）

w= cross(v, u)

这样很容易建立了vuw坐标系。

另外，我们希望引入theta角表示：圆柱面围绕中心线旋转的角度

（如果圆柱面是纯色的，则旋转是看不出来的；如果圆柱面有多种颜色，则可以看到旋转的效果）

在已经建立的vuw坐标系中：

旋转theta角后，v轴的方向向量转到v1位置，设v1为vuw坐标系中的单位向量，则v1在vuw坐标系中的坐标为（cos(theta), 0, sin(theta)）

将v1的坐标转换到xyz坐标系vector_trans_back(v1, v, u, w)

然后在xyz坐标系中w1 = cross(v1, u)

所以，v1, u, w1，即为旋转theta角后的最终的坐标系的基。

另外，补充一点：u.x()=0（即在YOZ平面）时，由于v固定为（1，0，0）

41.2 看C++代码实现

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

vec3.cpp

bool get_vector_vw(vec3& u, float angle, vec3& v, vec3& w) {
/*determine the v axis and w axis of u-v-w space*/vec3 v1, w1;u = unit_vector(u);float theta = angle*M_PI/180;
if (u.x() == 0.0) {
//u在ZOY平面时，v1为x轴上的单位向量（这里的v1对应“数学推导”中的v）v1 = vec3(1.0, 0.0, 0.0);theta = (angle)*M_PI/180;}else {v1 = unit_vector(vec3(-u.z(), 0.0, u.x()));theta = (angle)*M_PI/180;}w1 = cross(v1, u);//这里的w1对应“数学推导”中的wfloat x = cos(theta);float z = sin(theta);
if (fabs(x) < 0.0001) {
/*这里是考虑到，计算机中计算结果中的0有时候表示为一个10-8数量级的数，而不是0，这样当theta为90*N度时，cos(theta)并不是真正的0，导致后面计算出错。*/x = 0.0;}if (fabs(z) < 0.0001) {z = 0.0;}vec3 v_uv1w1 = vec3(x, 0, z);v = unit_vector(vector_trans_back(v_uv1w1, v1, u, w1));w = cross(v, u);return true;
}

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

vec3.h

bool get_vector_vw(vec3& u, float angle, vec3& v, vec3& w);

----------------------------------------------quadratic_cylinder_all.h ------------------------------------------

quadratic_cylinder_all.h

#ifndef QUADRATIC_CYLINDER_ALL_H
#define QUADRATIC_CYLINDER_ALL_H#include <hitable.h>
#include "material.h"
#include "log.h"class quadratic_cylinder_all : public hitable
{public:quadratic_cylinder_all() {}quadratic_cylinder_all(vec3 cen, float a, float b, float c, float hy, material *m, vec3 u, float an) {center = cen;intercept_x = a;intercept_y = b;intercept_z = c;height_half_y = hy;ma = m;vector_u = u;get_vector_vw(vector_u, an, vector_v, vector_w);}
/*
(x-xc)^2/a^2 + 0*(y-yc)^2/b^2 + (z-zc)^2/c^2 = 1
*/virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;vec3 center;float intercept_x;float intercept_y;float intercept_z;float height_half_y;material *ma;vec3 vector_u;vec3 vector_v;vec3 vector_w;
};#endif // QUADRATIC_CYLINDER_ALL_H

----------------------------------------------quadratic_cylinder_all.cpp ------------------------------------------

quadratic_cylinder_all.cpp

#include "quadratic_cylinder_all.h"#include <iostream>
using namespace std;bool quadratic_cylinder_all::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
#if QUADRATIC_CYLINDER_ALL_LOG == 1std::cout << "-------------quadratic_cylinder_all::hit----------------" << endl;
#endif // QUADRATIC_CYLINDER_ALL_LOGvec3 direction = vector_trans(r.direction(), vector_v, vector_u, vector_w);vec3 origin = vector_trans(r.origin(), vector_v, vector_u, vector_w);vec3 center_vuw = vector_trans(center, vector_v, vector_u, vector_w);
/*这里就是将光线的起点和方向向量、圆柱面的中心从xyz坐标系转换到vuw坐标系*/float ab_square = intercept_x*intercept_x*intercept_y*intercept_y;float bc_square = intercept_y*intercept_y*intercept_z*intercept_z;float ac_square = 0.0*intercept_x*intercept_x*intercept_z*intercept_z;float abc_square = 1.0*intercept_x*intercept_x*intercept_y*intercept_y*intercept_z*intercept_z;vec3 inter_square = vec3(bc_square, ac_square, ab_square);vec3 rd_square = vec3(direction.x()*direction.x(),direction.y()*direction.y(),direction.z()*direction.z());float A = dot(inter_square, rd_square);vec3 r0_c = origin - center_vuw;vec3 r0_c_rd = vec3(r0_c.x()*direction.x(),r0_c.y()*direction.y(),r0_c.z()*direction.z());float B = 2*dot(r0_c_rd, inter_square);vec3 r0_c_square = vec3(r0_c.x()*r0_c.x(),r0_c.y()*r0_c.y(),r0_c.z()*r0_c.z());float C = dot(r0_c_square, inter_square) - abc_square;float temp, temp1, temp2;vec3 pc;if(A == 0) {if (B == 0) {return false;}else {temp = -C/B;if (temp < t_max && temp > t_min) {rec.t = temp;rec.p = vector_trans(r.point_at_parameter(rec.t), vector_v, vector_u, vector_w); //将撞击点坐标转换到vuw坐标系if (((rec.p.y()-center_vuw.y()) > -height_half_y) && ((rec.p.y()-center_vuw.y()) < height_half_y)) {pc = rec.p - center_vuw;rec.normal = unit_vector(vec3(2*bc_square*pc.x(), 2*ac_square*pc.y(), 2*ab_square*pc.z()));if (dot(rec.normal, direction) > 0) {rec.normal = -rec.normal;}rec.normal = vector_trans_back(rec.normal, vector_v, vector_u, vector_w);//最终的撞击点的法向量需要从vuw坐标系转回到xyz坐标系。rec.mat_ptr = ma;if ((intercept_x == intercept_y) && (intercept_y == intercept_z)) {//0/1:cylinderrec.v = (rec.p.y() - (center_vuw.y() - height_half_y)) / (2*height_half_y);vec3 pc1 = vec3(rec.p.x()-center_vuw.x(), 0, rec.p.z()-center_vuw.z());vec3 vx = vec3(0, 0, 1);float u = acos(dot(pc1, vx) / (pc1.length()*vx.length())) / (2*M_PI);
//                            float u = acos((rec.p.x() - center_vuw.x()) / intercept_x) / (2*M_PI);if ((rec.p.x() - center_vuw.x()) < 0) {rec.u = 1-u;}else {rec.u = u;}rec.p = vector_trans_back(rec.p, vector_v, vector_u, vector_w);return true;}else {rec.u = -1.0;rec.v = -1.0;rec.p = vector_trans_back(rec.p, vector_v, vector_u, vector_w);return true;}}else {return false;}}}}else {float discriminant = B*B - 4*A*C;if (discriminant >= 0) {temp1 = (-B - sqrt(discriminant)) / (2.0*A);temp2 = (-B + sqrt(discriminant)) / (2.0*A);if (temp1 > temp2) {//make sure that temp1 is smaller than temp2temp = temp1;temp1 = temp2;temp2 = temp;}if (temp1 < t_max && temp1 > t_min) {rec.t = temp1;rec.p = vector_trans(r.point_at_parameter(rec.t), vector_v, vector_u, vector_w); //将撞击点坐标转换到vuw坐标系if (((rec.p.y()-center_vuw.y()) > -height_half_y) && ((rec.p.y()-center_vuw.y()) < height_half_y)) {pc = rec.p - center_vuw;rec.normal = unit_vector(vec3(2*bc_square*pc.x(), 2*ac_square*pc.y(), 2*ab_square*pc.z()));if (dot(rec.normal, direction) > 0) {rec.normal = -rec.normal;}rec.normal = vector_trans_back(rec.normal, vector_v, vector_u, vector_w); //最终的撞击点的法向量需要从vuw坐标系转回到xyz坐标系。rec.mat_ptr = ma;if ((intercept_x == intercept_y) && (intercept_y == intercept_z)) {//cylinderrec.v = (rec.p.y() - (center_vuw.y() - height_half_y)) / (2*height_half_y);vec3 pc1 = vec3(rec.p.x()-center_vuw.x(), 0, rec.p.z()-center_vuw.z());vec3 vx = vec3(0, 0, 1);float u = acos(dot(pc1, vx) / (pc1.length()*vx.length())) / (2*M_PI);
//                            float u = acos((rec.p.x() - center_vuw.x()) / intercept_x) / (2*M_PI);if ((rec.p.x() - center_vuw.x()) < 0) {rec.u = 1-u;}else {rec.u = u;}rec.p = vector_trans_back(rec.p, vector_v, vector_u, vector_w);return true;}else {rec.u = -1.0;rec.v = -1.0;rec.p = vector_trans_back(rec.p, vector_v, vector_u, vector_w);return true;}}else {
//                        return false;}}if (temp2 < t_max && temp2 > t_min) {rec.t = temp2;rec.p = vector_trans(r.point_at_parameter(rec.t), vector_v, vector_u, vector_w); //将撞击点坐标转换到vuw坐标系if (((rec.p.y()-center_vuw.y()) > -height_half_y) && ((rec.p.y()-center_vuw.y()) < height_half_y)) {pc = rec.p - center_vuw;rec.normal = unit_vector(vec3(2*bc_square*pc.x(), 2*ac_square*pc.y(), 2*ab_square*pc.z()));if (dot(rec.normal, direction) > 0) {rec.normal = -rec.normal;}rec.normal = vector_trans_back(rec.normal, vector_v, vector_u, vector_w); //最终的撞击点的法向量需要从vuw坐标系转回到xyz坐标系。rec.mat_ptr = ma;if ((intercept_x == intercept_y) && (intercept_y == intercept_z)) {//cylinderrec.v = (rec.p.y() - (center_vuw.y() - height_half_y)) / (2*height_half_y);vec3 pc1 = vec3(rec.p.x()-center_vuw.x(), 0, rec.p.z()-center_vuw.z());vec3 vx = vec3(0, 0, 1);float u = acos(dot(pc1, vx) / (pc1.length()*vx.length())) / (2*M_PI);
//                            float u = acos((rec.p.x() - center_vuw.x()) / intercept_x) / (2*M_PI);if ((rec.p.x() - center_vuw.x()) < 0) {rec.u = 1-u;}else {rec.u = u;}rec.p = vector_trans_back(rec.p, vector_v, vector_u, vector_w);return true;}else {rec.u = -1.0;rec.v = -1.0;rec.p = vector_trans_back(rec.p, vector_v, vector_u, vector_w);return true;}}else {
//                        return false;}}}return false;}
}

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

main.cpp

hitable *list[2];

list[0] = new quadratic(vec3(-3.75, 3.25, 0), 1.25, 1.25, 1.25, 0, 1,1.25, new lambertian(vec3(0.8, 0.8, 0.0)));

list[1] = newquadratic_cylinder_all(vec3(3.75, 3.25, 0), 1.25, 1.25, 1.25, 1.25,

new lambertian(vec3(0.8,0.8, 0.0)), vec3(1, 1, 0), 0);

hitable *world = new hitable_list(list,2);

vec3 lookfrom(0, 5, 20);

vec3 lookat(0, 3, 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);

输出图片如下：

咦~咦~咦

我们设置的旋转角度为0，为什么会发生旋转呢？？？

41.3 补偿旋转角度

u向量在ZOX平面的投影向量为u_zox的坐标应该为（Xu, 0,Zu），u_zox向量和-Z轴的单位方向向量z_n（0，0，-1）的夹角为α，

v向量和与XOY平面的夹角为β，这β=α

β即为vuw坐标系相对于xyz坐标系旋转的角度。

如果要使圆柱面的实际效果是旋转角度为0，则需要补偿一个β角度。

具体做法：

u.x()=0（即在YOZ平面）时，由于v固定为（1，0，0），所以，不管β是0度还是180度，都不需要补偿。

u.x()>0，需要减去一个β

u.x()<0，需要加上一个β

代码修改如下：

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

vec3.cpp

bool get_vector_vw(vec3& u, float angle, vec3& v, vec3& w) {
/*determine the v axis and w axis of u-v-w space*/vec3 v1, w1;u = unit_vector(u);float theta = angle*M_PI/180;if (u.x() == 0.0) {v1 = vec3(1.0, 0.0, 0.0);theta = (angle)*M_PI/180;}else {v1 = unit_vector(vec3(-u.z(), 0.0, u.x()));vec3 z_n = vec3(0, 0, -1);vec3 u_xoz = vec3(u.x(), 0, u.z());float phi = acos(dot(z_n, u_xoz) / (z_n.length()*u_xoz.length()));if (u.x() > 0) {theta = (angle)*M_PI/180 - phi;}else {theta = (angle)*M_PI/180 + phi;}}w1 = cross(v1, u);float x = cos(theta);float z = sin(theta);if (fabs(x) < 0.0001) {x = 0.0;}if (fabs(z) < 0.0001) {z = 0.0;}vec3 v_uv1w1 = vec3(x, 0, z);v = unit_vector(vector_trans_back(v_uv1w1, v1, u, w1));w = cross(v, u);return true;
}bool roots_quadratic_equation2(float a, float b, float c, float (&roots)[3]) {//the first element is the number of the real roots, and other elements are the real roots.if (a == 0.0) {if (b == 0.0) {roots[0] = 0.0;}else {roots[1] = -c/b;roots[0] = 1.0;}}else {float d = b*b - 4*a*c;if (d < 0.0) {roots[0] = 0.0;}else {roots[1] = (-b + sqrt(d)) / (2*a);roots[2] = (-b - sqrt(d)) / (2*a);roots[0] = 2.0;}}return true;
}

旋转角度补偿前后图片对比：

补偿前：（之前已经贴出图片，这里再贴一次）

补偿后：

接下来，我们测几组图片：

vec3(0, 0, 1), 0

vec3(0, 0, -1), 0

vec3(1, 0, 0), 0

vec3(-1, 0, 0), 0

vec3(1, 1, 1), 0

vec3(1, 1, -1), 0

vec3(1, -1, 1), 0

vec3(1,- 1, -1), 0

vec3(-1, 1, 1), 0

vec3(-1, 1, -1), 0

vec3(-1, -1, 1), 0

vec3(-1, -1, -1), 0

vec3(-1, 1, 0), 0

vec3(-1, 0,-1), 0

vec3(-1, -1, 0), 0

vec3(0, 1, 1), 0

vec3(1, 1, 1), 0

vec3(1, 1, 1), 45

vec3(1, 1, 1), 90

vec3(1, 1, 1), 135