问题四十九：怎么用ray tracing画supertoroid（超级圆环）

49.1 数学推导

49.2 看C++代码实现

----------------------------------------------supertoroid.h ------------------------------------------

supertoroid.h

 

#ifndef SUPERTOROID_H
#define SUPERTOROID_H

#include <hitable.h>
#include "material.h"
#include "log.h"

class supertoroid : public hitable
{
    public:
        supertoroid() {}
        supertoroid(vec3 cen, float a1, float a2, float a3, float a4, float e1, float e2, int in, float tol,  material *m) :
            center(cen), intercept_x(a1), intercept_y(a2), intercept_z(a3), radius(a4), p_e1(e1), p_e2(e2),
            initial_number(in), tolerance(tol), ma(m) {}
/*
f(x,y,z)=(( (x/a1)^(2/e2) + (z/a3)^(2/e2) )^(e2/2) - a4)^(2/e1) + (y/a2)^(2/e1) -1 = 0
a4=R/sqrt(a1*a1+a2*a2), R is the torus radius.
in: initial number
tol: tolerance
*/
        virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;
        vec3 center;
        float intercept_x, intercept_y, intercept_z, radius;
        float p_e1, p_e2;
        int initial_number;
        float tolerance;
        material *ma;
};

#endif // SUPERTOROID_H

----------------------------------------------supertoroid.cpp ------------------------------------------

supertoroid.cpp

 

#include "supertoroid.h"

#include <iostream>
#include <limits>
#include "float.h"
#include "log.h"

using namespace std;

bool ray_hit_box_t(const ray& r, const vec3& vertex_l, const vec3& vertex_h, float& t_near, float& t_far) {
/*光线撞击包围超级圆环的box，得到t_near, 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())) {
#if SUPERTOROID_LOG == 1
                std::cout << "the ray is parallel to the planes and the origin X0 is not between the slabs. return false" <<endl;
#endif // SUPERTOROID_LOG
                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())) {
#if SUPERTOROID_LOG == 1
                std::cout << "the ray is parallel to the planes and the origin Y0 is not between the slabs. return false" <<endl;
#endif // SUPERTOROID_LOG
                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())) {
#if SUPERTOROID_LOG == 1
                std::cout << "the ray is parallel to the planes and the origin Z0 is not between the slabs. return false" <<endl;
#endif // SUPERTOROID_LOG
                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 SUPERTOROID_LOG == 1
            std::cout << "array1[" << i << "]:" << array1[i] <<endl;
            std::cout << "array1[" << i+1 << "]:" << array1[i+1] <<endl;
#endif // SUPERTOROID_LOG
            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) {
#if SUPERTOROID_LOG == 1
                std::cout << "No.(0=X;2=Y;4=Z):" << i << "  :t_near > t_far. return false" <<endl;
#endif // SUPERTOROID_LOG
                return false;
            }
            if(t_far < 0) {
#if SUPERTOROID_LOG == 1
                std::cout << "No.(0=X;2=Y;4=Z):" << i << "  :t_far < 0. return false" <<endl;
#endif // SUPERTOROID_LOG
                return false;
            }
        }
        if (t_near != t_near) {
            t_near = t_near * 1;
        }
        return true;
}


bool get_superellipsoid_function_and_derivative_t(float a1, float a2, float a3, float a4, float e1, float e2, float xo, float yo, float zo, float xd, float yd, float zd, double t, double& f, double& fd) {
/*求函数值和函数导数值*/
        double a1_r = double(a1);
        double a2_r = double(a2);
        double a3_r = double(a3);
        double a4_r = double(a4);
        double e1_r = double(e1);
        double e2_r = double(e2);
        double xo_r = double(xo);
        double yo_r = double(yo);
        double zo_r = double(zo);
        double xd_r = double(xd);
        double yd_r = double(yd);
        double zd_r = double(zd);
        double pow_x, pow_y, pow_z, pow_x_d, pow_y_d, pow_z_d, pow_x_z, pow_x_z_d, pow_x_z_a4, pow_x_z_a4_d;
        double xd_a1, yd_a2, zd_a3;

        if ((xo_r+xd_r*t) < 0) {
            xd_a1 = -xd_r/a1_r;
        }
        else {
            xd_a1 = xd_r/a1_r;
        }
        if ((yo_r+yd_r*t) < 0) {
            yd_a2 = -yd_r/a2_r;
        }
        else {
            yd_a2 = yd_r/a2_r;
        }
        if ((zo_r+zd_r*t) < 0) {
            zd_a3 = -zd_r/a3_r;
        }
        else {
            zd_a3 = zd_r/a3_r;
        }

        if ((xo_r+xd_r*t) == 0) {
            pow_x = 0;
            pow_x_d = 0;
        }
        else {
            pow_x = pow(fabs(xo_r/a1_r + xd_r*t/a1_r), (2/e2_r));
            pow_x_d = pow(fabs(xo_r/a1_r + xd_r*t/a1_r), ((2/e2_r)-1));
        }
        if ((yo_r+yd_r*t) == 0) {
            pow_y = 0;
            pow_y_d = 0;
        }
        else {
            pow_y = pow(fabs(yo_r/a2_r + yd_r*t/a2_r), (2/e1_r));
            pow_y_d = pow(fabs(yo_r/a2_r + yd_r*t/a2_r), ((2/e1_r)-1));
        }
        if ((zo_r+zd_r*t) == 0) {
            pow_z = 0;
            pow_z_d = 0;
        }
        else {
            pow_z = pow(fabs(zo_r/a3_r + zd_r*t/a3_r), (2/e2_r));
            pow_z_d = pow(fabs(zo_r/a3_r + zd_r*t/a3_r), ((2/e2_r)-1));
        }

        if ((pow_x+pow_z) == 0) {
            pow_x_z = 0;
            pow_x_z_d = 0;
        }
        else {
            pow_x_z = pow(pow_x+pow_z, (e2_r/2));
            pow_x_z_d = pow(pow_x+pow_z, (e2_r/2)-1);
        }
        if ((pow_x_z-a4_r) == 0) {
            pow_x_z_a4 = 0;
            pow_x_z_a4_d = 0;
        }
        else {
            pow_x_z_a4 = pow(fabs(pow_x_z-a4_r), (2/e1_r));
            pow_x_z_a4_d = pow(fabs(pow_x_z-a4_r), (2/e1_r)-1);
        }
/*如上对幂函数的基数取了绝对值。实际情况，基数可正可负，对于不同的指数，幂函数的奇偶性不同，所以，为了计算方便，在此取其绝对值。这种做法的结果和真实情况是有出入的。*/
        f = pow_x_z_a4+pow_y - 1;
        fd = (2/e1_r)*(pow_x_z_a4_d*pow_x_z_d
*(pow_x_d*xd_a1+pow_z_d*zd_a3)+pow_y_d*yd_a2);

        return true;
}

bool get_roots_by_newton_iteration_t(const ray& r, vec3 c, float a1, float a2, float a3, float a4, float e1, float e2, int in, float tol, float *x0, float (&roots)[2]) {
/*牛顿迭代*/
        float xo = r.origin().x() - c.x();
        float yo = r.origin().y() - c.y();
        float zo = r.origin().z() - c.z();
        float xd = r.direction().x();
        float yd = r.direction().y();
        float zd = r.direction().z();
        double t_k, t_k1, ft_k, ft_d_k;
        int j=0, in_r;
        if (in > int(x0[0])) {
            in_r = int(x0[0]);
        }
        else {
            in_r = in;
        }

        for (int i=1; i<in_r; i++) {
            t_k = double(x0[i]);
            for (int k=0; k<50; k++) {
                if (!(isnan(t_k))) {
                    get_superellipsoid_function_and_derivative_t(a1, a2, a3, a4, e1, e2, xo, yo, zo, xd, yd, zd, t_k, ft_k, ft_d_k);
                    if ((ft_d_k != 0) && !(isnan(ft_k)) && !(isnan(ft_d_k))) {
                        t_k1 = t_k - ft_k/ft_d_k;
//                        if (fabs(t_k1) >= 1) {
                            if (fabs((t_k1 - t_k)/t_k1) < tol) {
                                if ((t_k1 >= x0[1]) && (t_k1 <= x0[in_r])) {
                                    roots[j+1] = float(t_k1);
                                    j++;
                                    break;
                                }
                                else {
                                    break;
                                }
                            }
                            else {
                                t_k = t_k1;
                            }
/*
                        }
                        else {
                            if (fabs(t_k1 - t_k) < tol) {
                                roots[j+1] = float(t_k1);
                                j++;
                                break;
                            }
                            else {
                                t_k = t_k1;
                            }
                        }
*/
                    }
                    else {
                        break;
                    }
                }
                else {
                    break;
                }
            }

            if (j == 1) {
                break;
            }

        }
        roots[0] = float(j);
        if (j == 0) {

        }
        return true;
}

bool supertoroid::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
#if SUPERELLIPSOID_LOG == 1
        std::cout << "-------------supertoroid::hit----------------" << endl;
#endif // SUPERELLIPSOID_LOG
        float intercept_tori_x, intercept_tori_y, intercept_tori_z;
        float pow_y = pow(1+pow(radius, (2/p_e1)), (p_e1/2));
        intercept_tori_x = (radius+1)*intercept_x;
        intercept_tori_y = pow_y*intercept_y;
        intercept_tori_z = (radius+1)*intercept_z;
/*这三个值，是分别令y, z为0，求得intercept_tori_x；令x, z为0，求得intercept_tori_y；令x, y为0，求得intercept_tori_z*/
        vec3 vertex_l[4], vertex_h[4];
        vertex_l[0] = vec3(center.x()-intercept_tori_x, center.y()-intercept_tori_y, center.z()-intercept_tori_z);
        vertex_h[0] = vec3(center.x()+intercept_tori_x, center.y()+intercept_tori_y, center.z()+intercept_tori_z);
        vertex_l[1] = vec3(center.x()-intercept_tori_x/2, center.y()-intercept_tori_y, center.z()-intercept_tori_z);
        vertex_h[1] = vec3(center.x()+intercept_tori_x/2, center.y()+intercept_tori_y, center.z()+intercept_tori_z);
        vertex_l[2] = vec3(center.x()-intercept_tori_x, center.y()-intercept_tori_y/2, center.z()-intercept_tori_z);
        vertex_h[2] = vec3(center.x()+intercept_tori_x, center.y()+intercept_tori_y/2, center.z()+intercept_tori_z);
        vertex_l[3] = vec3(center.x()-intercept_tori_x, center.y()-intercept_tori_y, center.z()-intercept_tori_z/2);
        vertex_h[3] = vec3(center.x()+intercept_tori_x, center.y()+intercept_tori_y, center.z()+intercept_tori_z/2);

        float roots[2] = {0.0, -1.0};
        float x0[initial_number+1];
        float t_near = 0;
        float t_far = 0;
        if (ray_hit_box_t(r, vertex_l[0], vertex_h[0], t_near, t_far)) {
            for (int i=0; i<initial_number; i++) {
                x0[i+1] = t_near + i*(t_far - t_near)/(initial_number-1);
            }
/*将t_near, t_far等分若干份，得到若干个迭代初值*/
            x0[0] = float(initial_number);
            get_roots_by_newton_iteration_t(r, center, intercept_x, intercept_y, intercept_z, radius, p_e1, p_e2, initial_number, tolerance, x0, roots);
        }
        else {
            return false;
        }

        float temp;
        if (roots[0] > 0.0001) {
            for (int i=1; i<int(roots[0]); i++) {
                for (int j=i+1; j<int(roots[0])+1; j++) {
                    if (roots[i] > roots[j]) {
                        temp = roots[i];
                        roots[i] = roots[j];
                        roots[j] = temp;
                    }
                }
            }
            double nx, ny, nz, pow_x, pow_z, pow_x_d, pow_z_d, pow_y_d, pow_x_z, pow_x_z_d, pow_x_z_a4_d;
            float d_a1 = intercept_x;
            float d_a2 = intercept_y;
            float d_a3 = intercept_z;
            float d_a4 = radius;
            vec3 pc;
            for (int k=1; k<int(roots[0])+1; k++) {
                if (roots[k] < t_max && roots[k] > t_min) {
                    rec.t = roots[k];
                    rec.p = r.point_at_parameter(rec.t);
                    pc = rec.p - center;
                    if (pc.x() < 0) {d_a1 = -intercept_x;}
                    if (pc.y() < 0) {d_a2 = -intercept_y;}
                    if (pc.z() < 0) {d_a3 = -intercept_z;}
                    if (pc.x() == 0) {
                        pow_x = 0;
                        pow_x_d = 0;
                    }
                    else {
                        pow_x = pow(double(fabs(pc.x()/d_a1)), double(2/p_e2));
                        pow_x_d = pow(double(fabs(pc.x()/d_a1)), double(2/p_e2-1));
                    }
                    if (pc.y() == 0) {
                        pow_y_d = 0;
                    }
                    else {
                        pow_y_d = pow(double(fabs(pc.y()/d_a2)), double(2/p_e1-1));
                    }
                    if (pc.z() == 0) {
                        pow_z = 0;
                        pow_z_d = 0;
                    }
                    else {
                        pow_z = pow(double(fabs(pc.z()/d_a3)), double(2/p_e2));
                        pow_z_d = pow(double(fabs(pc.z()/d_a3)), double(2/p_e2-1));
                    }
                    if ((pow_x+pow_z) == 0) {
                        pow_x_z = 0;
                        pow_x_z_d = 0;
                    }
                    else {
                        pow_x_z = pow(pow_x+pow_z, (p_e2/2));
                        pow_x_z_d = pow(pow_x+pow_z, (p_e2/2)-1);
                    }
                    if ((pow_x_z-d_a4) == 0) {
                        pow_x_z_a4_d = 0;
                    }
                    else {
                        pow_x_z_a4_d = pow(fabs(pow_x_z-d_a4), (2/p_e1)-1);
                    }

                    nx = double(2/p_e1) * pow_x_z_a4_d * pow_x_z_d  * pow_x_d / d_a1;
                    ny = double(2/p_e1) * pow_y_d / d_a2;
                    nz = double(2/p_e1) * pow_x_z_a4_d * pow_x_z_d  * pow_z_d / d_a3;

                    if (isnan(nx)) {
                        nx = nx * 1;
                    }

                    nx = nx/sqrt(nx*nx+ny*ny+nz*nz);
                    ny = ny/sqrt(nx*nx+ny*ny+nz*nz);
                    nz = nz/sqrt(nx*nx+ny*ny+nz*nz);

                    rec.normal = unit_vector(vec3(float(nx), float(ny), float(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;
                }
            }
            return false;
        }
        else {
            return false;
        }
        return false;
}

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

main.cpp

 

        hitable *list[1];
        list[0] = new supertoroid(vec3(0, 3, 0), 1, 1, 1, 6, 2, 0.1, 16, 0.0001, new lambertian(vec3(0.0, 1.0, 0.0)));

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

        vec3 lookfrom(10, 20, 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), 30, float(nx)/float(ny), aperture, 0.7*dist_to_focus);

 

输出结果如下：

 

e1=0.1, e2=1, initial_number=10, tol=0.0001

e1=0.1, e2=1, initial_number=16, tol=0.0001

将迭代初值个数调到16，还是有未能迭代到的实根，只有再增加迭代初值的个数。此处不再调试，太费时间了（调一张图要花差不多半个小时，太慢了）。

 

e1=1, e2=1, initial_number=10, tol=0.0001

e1=2, e2=1, initial_number=16, tol=0.0001

 

e1=4, e2=1, initial_number=16, tol=0.0001

 

e1=0.1, e2=0.1, initial_number=16, tol=0.0001

将采样次数设为1，调了比较接近的图:

e1=0.1, e2=1, initial_number=56, tol=0.002

e1=0.1, e2=1, initial_number=56, tol=0.004

 

e1=1, e2=0.1, initial_number=16, tol=0.0001

 

e1=2, e2=0.1, initial_number=16, tol=0.0001

 

e1=0.1, e2=2, initial_number=16, tol=0.0001

 

e1=1, e2=2, initial_number=16, tol=0.0001

 

e1=0.1, e2=3, initial_number=16, tol=0.0001

 

e1=1, e2=3, initial_number=16, tol=0.0001

 

最后贴一张截图，主要是看看生成图片是多花时间。