问题三十五：怎么用ray tracing画二次曲面(quadratic surfaces)（3）—

35.3 椭球抛物面

35.3.1 数学推导

椭球抛物面的方程如下：

所以，其一：我们需要对两个实根进行排序（先处理小的）

另外，由于，是开放曲面，也就是，光线有可能撞击到曲面的正反两面，所以，对于撞击点处的标准化之后的法向量，我们需要做如下判断：

if (dot(rec.normal,r.direction()) > 0) {

rec.normal =-rec.normal;

}//（法向量决定着反射光线和折射光线）

还有，由于我们引入了height_y参数来限制曲面的高度，但是，我们要注意到：小实根对应的交点超出高度范围时（之前的一贯做法：小根不在t_min、t_max范围，就直接return false），大实根是有可能在高度范围的（而且，如果在范围的话，光线撞击的曲面的内表面，这时的法向量是需要反向的）。

35.3.2 看C++代码实现

----------------------------------------------quadratic_paraboloid.h ------------------------------------------

quadratic_paraboloid.h

#ifndef QUADRATIC_PARABOLOID_H
#define QUADRATIC_PARABOLOID_H#include <hitable.h>class quadratic_paraboloid : public hitable
{public:quadratic_paraboloid() {}quadratic_paraboloid(vec3 cen, float p, float q, int hy, material *m) : center(cen), focus_directrix_p(p), focus_directrix_q(q), height_half_y(hy), 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;material *ma;
};#endif // QUADRATIC_PARABOLOID_H

----------------------------------------------quadratic_paraboloid.cpp ------------------------------------------

quadratic_paraboloid.cpp

#include "quadratic_paraboloid.h"#include <iostream>
using namespace std;bool quadratic_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)) {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 temp2temp = 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)) {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)) {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_paraboloid(vec3(0.5, 1.5, 0), 1, 2, 1, new lambertian(vec3(0.9, 0.1, 0.5)));hitable *world = new hitable_list(list,2);vec3 lookfrom(0,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);

输出图片：