问题三十五： 怎么用ray tracing画二次曲面(quadratic surfaces)（4）——双曲抛物面（马鞍面）

35.4 双曲抛物面（马鞍面）

35.4.1 数学推导

双曲抛物面的方程如下：

35.4.2 看C++代码实现

---------------------------------------------quadratic_hyperbolic_paraboloid.h ------------------------------------

quadratic_hyperbolic_paraboloid.h

 

#ifndef QUADRATIC_HYPERBOLIC_PARABOLOID_H#define QUADRATIC_HYPERBOLIC_PARABOLOID_H#include <hitable.h>class quadratic_hyperbolic_paraboloid : public hitable{    public:        quadratic_hyperbolic_paraboloid() {}        quadratic_hyperbolic_paraboloid(vec3 cen, float p, float q, int hy, int wx, material *m) : center(cen), focus_directrix_p(p), focus_directrix_q(-q), height_half_y(hy), width_x(wx), ma(m) {}/*(x-xc)^2/2*p + (z-zc)^2/2*q = (y-yc)(p>0, q>0)*/        virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;        vec3 center;        float focus_directrix_p;        float focus_directrix_q;        int height_half_y, width_x;        material *ma;};#endif // QUADRATIC_HYPERBOLIC_PARABOLOID_H

-------------------------------------------quadratic_hyperbolic_paraboloid.cpp ----------------------------------

quadratic_hyperbolic_paraboloid.cpp

 

#include "quadratic_hyperbolic_paraboloid.h"#include <iostream>using namespace std;bool quadratic_hyperbolic_paraboloid::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {        float A = focus_directrix_q*r.direction().x()*r.direction().x()                + focus_directrix_p*r.direction().z()*r.direction().z();        float B = 2*focus_directrix_q*r.direction().x()*(r.origin().x()-center.x())                + 2*focus_directrix_p*r.direction().z()*(r.origin().z()-center.z())                - 2*focus_directrix_p*focus_directrix_q*r.direction().y();        float C = focus_directrix_q*(r.origin().x()-center.x())*(r.origin().x()-center.x())                + focus_directrix_p*(r.origin().z()-center.z())*(r.origin().z()-center.z())                - 2*focus_directrix_p*focus_directrix_q*(r.origin().y()-center.y());        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 = r.point_at_parameter(rec.t);                    if ((((rec.p.y()-center.y()) > -height_half_y) && ((rec.p.y()-center.y()) < height_half_y)) && (((rec.p.x()-center.x()) > -width_x) && ((rec.p.x()-center.x()) < width_x)) ) {                        pc = rec.p - center;                        rec.normal = unit_vector(vec3(2*focus_directrix_q*pc.x(), -2*focus_directrix_p*focus_directrix_q, 2*focus_directrix_p*pc.z()));                        if (dot(rec.normal, r.direction()) > 0) {                            rec.normal = -rec.normal;                        }                        rec.mat_ptr = ma;                        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 = r.point_at_parameter(rec.t);                    if ((((rec.p.y()-center.y()) > -height_half_y) && ((rec.p.y()-center.y()) < height_half_y)) && (((rec.p.x()-center.x()) > -width_x) && ((rec.p.x()-center.x()) < width_x)) ) {                        pc = rec.p - center;                        rec.normal = unit_vector(vec3(2*focus_directrix_q*pc.x(), -2*focus_directrix_p*focus_directrix_q, 2*focus_directrix_p*pc.z()));                        if (dot(rec.normal, r.direction()) > 0) {                            rec.normal = -rec.normal;                        }                        rec.mat_ptr = ma;                        return true;                    }                    else {//                        return false;                    }                }                if (temp2 < t_max && temp2 > t_min) {                    rec.t = temp2;                    rec.p = r.point_at_parameter(rec.t);                    if ((((rec.p.y()-center.y()) > -height_half_y) && ((rec.p.y()-center.y()) < height_half_y)) && (((rec.p.x()-center.x()) > -width_x) && ((rec.p.x()-center.x()) < width_x)) ) {                        pc = rec.p - center;                        rec.normal = unit_vector(vec3(2*focus_directrix_q*pc.x(), -2*focus_directrix_p*focus_directrix_q, 2*focus_directrix_p*pc.z()));                        if (dot(rec.normal, r.direction()) > 0) {                            rec.normal = -rec.normal;                        }                        rec.mat_ptr = ma;                        return true;                    }                    else {//                        return false;                    }                }            }            return false;        }}

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

main.cpp

 

        hitable *list[2];        list[0] = new sphere(vec3(0.0,-100.5,-1), 100, new lambertian(vec3(0.8, 0.8, 0.0)));        list[1] = new quadratic_hyperbolic_paraboloid(vec3(0, 1.5, 0), 0.5, 2, 1.5, 2, new lambertian(vec3(0.0, 0.1, 0.5)));        hitable *world = new hitable_list(list,2);        vec3 lookfrom(2,5,10);        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);

输出图片：

如果改成这样vec3 lookfrom(5,5,10);，图是这样的：

如果改成这样vec3 lookfrom(10,5,10);，图是这样的：