问题五十八：怎么用ray tracing画conic sweeping图形

58.1 概述

在上一章节中画的translational sweeping图形是将曲线沿着y轴平移后得到的，在平移的过程中曲线保持不变。接下来，我们要画的是conic sweeping图形。这种图形也是曲线沿着y轴平移后得到的，但是，在平移的过程中对曲线进行了一定比例的缩放，这个比例一般是和y坐标成正比。示意图如下：

58.2 求根

后面求根的步骤和translational sweeping是完全一样了。

对应“57.2.5”和“57.2.6”

 

58.3 求法向量

若t0_t1_smaller在曲线内：

法向量即为cap plane的法向量：（0, 1, 0）

 

 

若t0_t1_bigger在曲线内

或者

t0、t1都不在曲线内

 

这两种情况的根是通过投影光线和曲线相交求得，对应的法向量即为曲线在该处的切向量和“该处和conic sweeping顶点连线向量”的叉乘。

Conic sweeping的法向量示意如下图：

58.4 看C++代码实现

下面只贴出在“问题五十七”的基础上有改动的部分代码。

 

----------------------------------------------sweeping_translational.h ------------------------------------------

sweeping_translational.h

 

class sweeping_translational : public hitable
{
    public:
        sweeping_translational() {}
        sweeping_translational(vec3 *cp, float by, float cy, material *m, int ty){
/*
ty=1: translational sweeping
ty=2: conic sweeping
*/
            base_y = by;
            cap_y = cy;
            ma = m;
            type = ty;

……

        material *ma;

        inttype;

};

 

----------------------------------------------sweeping_translational.cpp ------------------------------------------

sweeping_translational.cpp

 

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)) {
……
……
            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++) {

                if (type == 1) {
                    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;
                }
                if (type == 2) {
/*对应“式子58.6”*/
                    a3 = ((y0-base_y)*zd-z0*yd)*matrix_c_x[3][i] - ((y0-base_y)*xd-x0*yd)*matrix_c_z[3][i];
                    a2 = ((y0-base_y)*zd-z0*yd)*matrix_c_x[2][i] - ((y0-base_y)*xd-x0*yd)*matrix_c_z[2][i];
                    a1 = ((y0-base_y)*zd-z0*yd)*matrix_c_x[1][i] - ((y0-base_y)*xd-x0*yd)*matrix_c_z[1][i];
                    a0 = ((y0-base_y)*zd-z0*yd)*matrix_c_x[0][i] - ((y0-base_y)*xd-x0*yd)*matrix_c_z[0][i] + (cap_y-base_y)*(xd*z0 - zd*x0);
                }

                roots3_base = roots_cubic_equation(a3, a2, a1, a0);
                if (roots3_base[0] != 0.0) {
                    for (int j=1; j<(int(roots3_base[0])+1); j++) {
                        u = roots3_base[j];
                        if ((u>=0.0) && (u<=1.0)) {
                            xu = 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;

                            if (type == 1) {
                                roots_base[num_roots_base+1][0] = (xu-x0)/xd;
                                roots_base[num_roots_base+1][1] = i;
                                roots_base[num_roots_base+1][2] = u;
                                num_roots_base ++;
                            }
                            if (type == 2) {
/*对应“式子58.3”*/
                                rbt = ((y0-base_y)*xu-(cap_y-base_y)*x0)/((cap_y-base_y)*xd-yd*xu);
                                if (yd > 0) {
                                    rbt_s = (base_y-y0)/yd;
                                    rbt_b = (cap_y-y0)/yd;
                                }
                                if (yd < 0) {
                                    rbt_s = (cap_y-y0)/yd;
                                    rbt_b = (base_y-y0)/yd;
                                }
                                if ((rbt>rbt_s) && (rbt<rbt_b)) {
//                                if (1) {
                                    roots_base[num_roots_base+1][0] = rbt;
                                    roots_base[num_roots_base+1][1] = i;
                                    roots_base[num_roots_base+1][2] = u;
                                    num_roots_base ++;
                                }
                            }

                        }
                    }
                    roots_base[0][0] = float(num_roots_base);
                }
                delete [] roots3_base;
……
……
                            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;

                            if (type == 1) {
                                rec.normal = unit_vector(vec3(-nx, 0, nz));
                            }
                            if (type == 2) {
/*对应“58.3 求法向量”*/
                                vec3 va = vec3(nx, 0, nz);
                                float bx = matrix_c_x[0][num]+matrix_c_x[1][num]*nu+matrix_c_x[2][num]*nu*nu+matrix_c_x[3][num]*nu*nu*nu;
                                float bz = matrix_c_z[0][num]+matrix_c_z[1][num]*nu+matrix_c_z[2][num]*nu*nu+matrix_c_z[3][num]*nu*nu*nu;
                                vec3 vb = vec3(bx, (cap_y-base_y), bz);
                                rec.normal = unit_vector(cross(va, vb));
                            }

                        }
                        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;
}

----------------------------------------------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, -2.0, 3.0, new lambertian(vec3(1.0, 0.0, 0.0)), 2);
        hitable *world = new hitable_list(list,1);

        vec3 lookfrom(0.0001, 10, 20);
        vec3 lookat(0, 1, 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);

 

输出图形如下：（前边为conicsweeping，后边为translational sweeping）

(0.0001,10, 20)

 

如果将“式子58.3”和“式子58.6”中的“（y0-yb）”换成“y0”，则图形变成这个样子：

58.5 问题说明

如果在“58.2 求根”中不做如下验证：

 

----------------------------------------------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)), 2);
        hitable *world = new hitable_list(list,1);

        vec3 lookfrom(-5, 20, 12);
        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);

 

输出图形会多出一部分：

 

(-5, 20, 12)

 

(-5, 20, 8)

 

(-5, 20, 4)

(-5, 20, 4)（添加“验证投影光线是否有效”之后）