我们之前在“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对应“数学推导”中的w
float 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 == 1
std::cout << "-------------quadratic_cylinder_all::hit----------------" << endl;
#endif // QUADRATIC_CYLINDER_ALL_LOG
vec3 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:cylinder
rec.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 temp2
temp = 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)) {//cylinder
rec.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)) {//cylinder
rec.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