问题四十一：怎么用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对应“数学推导”中的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