问题五十一：怎么用ray tracing画tear drop

tear drop是长这个样子的：

51.1 数学推导

在网上找到tear drop的参数方程：

51.2 看C++代码实现

----------------------------------------------vec3.h ------------------------------------------

vec3.h

/*重写了解一元二次方程、一元三次方程、一元四次方程的函数，将其所有数据类型由float换成double，同时使用double数据类型的最小正数“DBL_EPSILON”，绝对值小于该最小正数的数即为double数据类型的0*/
#define EPS 1e-10double* roots_quadratic_equation_rain(double a, double b, double c);
double* roots_cubic_equation_rain(double a, double b, double c, double d);
bool roots_quartic_equation2_rain(double a, double b, double c, double d, double e, double (&roots)[5]);

----------------------------------------------vec3.cpp ------------------------------------------

vec3.cpp

double* roots_quadratic_equation_rain(double a, double b, double c) {//the first element is the number of the real roots, and other elements are the real roots.double *roots = new double[3];if (fabs(a) < DBL_EPSILON) {if (fabs(b) < DBL_EPSILON) {roots[0] = 0.0;}else {roots[1] = -c/b;roots[0] = 1.0;}}else {double d = b*b - 4*a*c;if (d < -DBL_EPSILON) {roots[0] = 0.0;}else if (fabs(d) <= DBL_EPSILON) {roots[1] = double((-b) / (2*a));roots[2] = double((-b) / (2*a));roots[0] = double(1.0);}else {roots[1] = double((-b + sqrt(d)) / (2*a));roots[2] = double((-b - sqrt(d)) / (2*a));roots[0] = double(2.0);}}return roots;
}double* roots_cubic_equation_rain(double a, double b, double c, double d) {//the first element is the number of the real roots, and other elements are the real roots.//Shengjin's formuladouble *roots = new double[4];if (fabs(a) < DBL_EPSILON) {double *roots2;roots2 = roots_quadratic_equation_rain(b, c, d);for (int i=0; i<int(roots2[0])+1; i++) {roots[i] = roots2[i];}delete [] roots2;}else {double A = double(b*b - 3.0*a*c);double B = double(b*c - 9.0*a*d);double C = double(c*c - 3.0*b*d);double deita = double(B*B - 4.0*A*C);if ((A == B) && (fabs(A) < DBL_EPSILON)) {//the three roots are the sameif (!(fabs(a) < DBL_EPSILON)) {roots[1] = -b/(3.0*a);}else {if (!(fabs(b) < DBL_EPSILON)) {roots[1] = -c/b;}else {if (!(fabs(c) < DBL_EPSILON)) {roots[1] = -3.0*d/c;}}}roots[2] = roots[1];roots[3] = roots[1];roots[0] = double(3.0);}else if (deita > DBL_EPSILON) {//only one real rootdouble y1 = double(A*b + (3*a)*(-B + sqrt(deita))/2);double y2 = double(A*b + (3*a)*(-B - sqrt(deita))/2);double pow_y1, pow_y2;if (y1 < -DBL_EPSILON) {//for pow(a,b), when b is not int, a should not be negative.pow_y1 = - pow(double(-y1), double(1.0/3.0));}else if (fabs(y1) <= DBL_EPSILON) {pow_y1 = 0;}else {pow_y1 = pow(double(y1), double(1.0/3.0));}if (y2 < -DBL_EPSILON) {pow_y2 = - pow(double(-y2), double(1.0/3.0));}else if (fabs(y2) <= DBL_EPSILON) {pow_y2 = 0;}else {pow_y2 = pow(double(y2), double(1.0/3.0));}roots[1] = double((-b - pow_y1 - pow_y2) / (3.0*a));roots[0] = double(1.0);}else if (fabs(deita) <= DBL_EPSILON) {//three real roots and two of them are the samedouble K = B/A;roots[1] = double(-b/a + K);roots[2] = double(-K/2.0);roots[3] = double(-K/2.0);roots[0] = double(3.0);}else if (deita < -DBL_EPSILON) {//three different real rootsdouble theta = acos((2.0*A*b-3.0*a*B) / (2.0*pow(A, 1.5)));roots[1] = double((-b - 2.0*sqrt(A)*cos(theta/3.0)) / (3.0*a));roots[2] = double((-b + sqrt(A) * (cos(theta/3.0) + sqrt(double(3.0))*sin(theta/3.0))) / (3.0*a));roots[3] = double((-b + sqrt(A) * (cos(theta/3.0) - sqrt(double(3.0))*sin(theta/3.0))) / (3.0*a));roots[0] = double(3.0);}}return roots;
}bool roots_quartic_equation2_rain(double a, double b, double c, double d, double e, double (&roots)[5]) {//the first element is the number of the real roots, and other elements are the real roots.//Descartes's Method.if (fabs(a) < DBL_EPSILON) {double *roots3;roots3 = roots_cubic_equation_rain(b, c, d, e);for (int i=0; i<int(roots3[0])+1; i++) {roots[i] = roots3[i];}delete [] roots3;}else {double a1 = b/a;double b1 = c/a;double c1 = d/a;double d1 = e/a;double p = (-3*a1*a1)/8 + b1;double q = (a1*a1*a1)/8 - (a1*b1)/2 + c1;double r = (-3*a1*a1*a1*a1)/256 + (a1*a1*b1)/16 - (a1*c1)/4 + d1;double sa = 2*p;double sb = p*p - 4*r;double sc = -q*q;double k, m, t;double *roots_s = roots_cubic_equation_rain(1.0, sa, sb, sc);
#if VEC3_LOG == 1for (int i=0; i<(roots_s[0]+1); i++) {std::cout << "roots_s3[" << i << "]=" << roots_s[i] << endl;}
//        std::cout << "DBL_EPSILON:" << DBL_EPSILON << endl;
#endif // VEC3_LOGdouble temp;for (int i=1; i<int(roots_s[0]); i++) {for (int j=i+1; j<int(roots_s[0])+1; j++) {if (roots_s[i] > roots_s[j]) {temp = roots_s[i];roots_s[i] = roots_s[j];roots_s[j] = temp;}}}if (roots_s[int(roots_s[0])] > DBL_EPSILON) {k = sqrt(double(roots_s[int(roots_s[0])]));delete [] roots_s;}else {if (fabs(q) < DBL_EPSILON) {k = 0.0;delete [] roots_s;}else {roots[0] = 0.0;delete [] roots_s;return true;}}if (fabs(k) == 0.0) {if (fabs(q) < DBL_EPSILON) {double *roots2;roots2 = roots_quadratic_equation_rain(1.0, -p, r);if (roots2[0] == 0.0) {roots[0] = 0.0;delete [] roots2;return true;}else {m = roots2[1];t = roots2[2];}delete [] roots2;}else {roots[0] = 0.0;return true;}}else {m = (k*k*k + k*p + q) / (2*k);t = (k*k*k + k*p - q) / (2*k);}double *roots_y1;double *roots_y2;roots_y1 = roots_quadratic_equation_rain(1.0, k, t);
#if VEC3_LOG == 1for (int i=0; i<(roots_y1[0]+1); i++) {std::cout << "roots_y1[" << i << "]=" << roots_y1[i] << endl;}
#endif // VEC3_LOGroots_y2 = roots_quadratic_equation_rain(1.0, -k, m);
#if VEC3_LOG == 1for (int i=0; i<(roots_y2[0]+1); i++) {std::cout << "roots_y2[" << i << "]=" << roots_y2[i] << endl;}
#endif // VEC3_LOGif (roots_y1[0] != 0.0) {for (int i=1; i<roots_y1[0]+1; i++) {roots[i] = roots_y1[i] - a1/4;}}if (roots_y2[0] != 0.0) {int roots_y1_number = int(roots_y1[0]);for (int j=1; j<roots_y2[0]+1; j++) {roots[roots_y1_number+j] =roots_y2[j] - a1/4;}}roots[0] = roots_y1[0] + roots_y2[0];delete [] roots_y1;delete [] roots_y2;}return true;
}

----------------------------------------------rain.h ------------------------------------------

rain.h

#ifndef RAIN_H
#define RAIN_H#include "hitable.h"
#include "material.h"
#include "log.h"class rain : public hitable
{public:rain() {}rain(vec3 cen, float a1, float a2, float a3, material *m) : center(cen), x_a1(a1), y_a2(a2), z_a3(a3), ma(m) {}
/*
parametric equation:x=a1*cos(u)*sin(v)*(1-cos(v))/2y=a2*cos(v)z=a3*sin(u)*sin(v)*(1-cos(v))/2
imlipcit equation:4*(x/a1)^2+4*(z/a3)^2+(y/a2)^4-2*(y/a2)^3+2*(y/a2)-1=0
*/virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;vec3 center;float x_a1, y_a2, z_a3;material *ma;
};
#endif // RAIN_H

----------------------------------------------rain.cpp ------------------------------------------

rain.cpp

#include "rain.h"
#include <iostream>
using namespace std;bool rain::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
#if RAIN_LOG == 1std::cout << "-------------rain::hit---1-------------" << endl;
#endif // RAIN_LOGvec3 oc1 = r.origin() - center;vec3 dir1 = r.direction();vec3 oc = vec3(double(oc1.x()), double(oc1.y()), double(oc1.z()));vec3 dir = vec3(double(dir1.x()), double(dir1.y()), double(dir1.z()));double A = 4*z_a3*z_a3*y_a2*y_a2*y_a2*y_a2;double B = 4*x_a1*x_a1*y_a2*y_a2*y_a2*y_a2;double C = x_a1*x_a1*z_a3*z_a3;double D = -2*x_a1*x_a1*z_a3*z_a3*y_a2;double E = 2*x_a1*x_a1*z_a3*z_a3*y_a2*y_a2*y_a2;double F = -x_a1*x_a1*z_a3*z_a3*y_a2*y_a2*y_a2*y_a2;double a4 = dir.y()*dir.y()*dir.y()*dir.y()*C;double a3 = 4*oc.y()*dir.y()*dir.y()*dir.y()*C + dir.y()*dir.y()*dir.y()*D;double a2 = 6*oc.y()*oc.y()*dir.y()*dir.y()*C + 3*oc.y()*dir.y()*dir.y()*D + dir.z()*dir.z()*B + dir.x()*dir.x()*A;double a1 = 4*oc.y()*oc.y()*oc.y()*dir.y()*C + 3*oc.y()*oc.y()*dir.y()*D + 2*oc.z()*dir.z()*B + 2*oc.x()*dir.x()*A + dir.y()*E;double a0 = oc.y()*oc.y()*oc.y()*oc.y()*C + oc.y()*oc.y()*oc.y()*D + oc.y()*E + F + oc.z()*oc.z()*B + oc.x()*oc.x()*A;double roots[5];roots_quartic_equation2_rain(a4, a3, a2, a1, a0, roots);double 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;}}}vec3 pc;double nx, ny, nz;for (int k=1; k<int(roots[0])+1; k++) {if ((float(roots[k]) < t_max) && (float(roots[k]) > t_min)) {rec.t = float(roots[k]);rec.p = r.point_at_parameter(rec.t);pc = rec.p - center;if (rec.p.y() < 1) {
//this is a bad trick, there are something wrong in the progam when r.direction().y() superimposes the drop picture.continue;}nx = 2.0*A*double(pc.x());ny = 4.0*C*double(pc.y())*double(pc.y())*double(pc.y()) + 3.0*D*double(pc.y())*double(pc.y()) + E;nz = 2.0*B*double(pc.z());double length = sqrt(nx*nx+ny*ny+nz*nz);if (fabs(length) >= EPS) {/* if the length is very small,when we trasform the type from double to float,we will get: nx=ny=nz=0.*/nx = nx/length;ny = ny/length;nz = nz/length;}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;
}

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

main.cpp

    vec3 color(const ray& r, hitable *world, int depth) {hit_record rec;if (world->hit(r, 0.001, (numeric_limits<float>::max)(), rec)) {ray scattered;vec3 attenuation;if (depth < 50 && rec.mat_ptr->scatter(r, rec, attenuation, scattered)) {
//                std::cout << "-------------color()----------------depth:" << depth << endl;return attenuation*color(scattered, world, depth+1);}else {return vec3(0,0,0);}}else {vec3 unit_direction = unit_vector(r.direction());float t = 0.5*(unit_direction.y() + 1.0);

if (isnan(t)) {

t = 0;

}

/*经过这个判断之后，再也不用担心像素点值溢出了*/

//            std::cout << "-------------color()--5----------------" << endl;return (1.0-t)*vec3(1.0, 1.0, 1.0) + t*vec3(0.5, 0.7, 1.0);//white, light blue}}int main(){int nx = 200;int ny = 100;int ns = 100;//        ofstream outfile( ".\\results\\superhyperboloid-two.txt", ios_base::out);ofstream outfile( ".\\results\\rain.txt", ios_base::out);outfile << "P3\n" << nx << " " << ny << "\n255\n";std::cout << "P3\n" << nx << " " << ny << "\n255\n";hitable *list[5];list[0] = new rain(vec3(0, 3, 0), 2, 2, 2, new lambertian(vec3(1.0, 0.0, 0.0)));list[1] = new rain(vec3(-3, 4, -4), 2, 2, 2, new lambertian(vec3(0.0, 1.0, 0.0)));list[2] = new rain(vec3(3, 4, -4), 2, 2, 2, new lambertian(vec3(0.0, 1.0, 0.0)));list[3] = new rain(vec3(-6.5, 4.5, -6), 2, 2, 2, new lambertian(vec3(0.0, 0.0, 1.0)));list[4] = new rain(vec3(6.5, 4.5, -6), 2, 2, 2, new lambertian(vec3(0.0, 0.0, 1.0)));hitable *world = new hitable_list(list,5);vec3 lookfrom(0, 0, 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);

输出图片如下：

51.3 问题说明

在rain.cpp的hit()函数中若没有这段code:

if (rec.p.y() < 1) {

//this is abad trick, there are something wrong in the progam when r.direction().y()superimposes the drop picture.

continue;

}

输出图片会出现多余的像素点（如下方红色框标出）：

原因是：一元四次方程的实根多了。设置断点拦截了一个，如下：

但是通过网页上一元四次方程计算器，这个方程是没有实根的：

方程为什么会多出实根呢？解方程的函数看了又看，查了又查，没有找到彻底解决方法。

问题原因是：

出现多根的位置是，y=r.direction().y()附近。

vec3 lookfrom(0, 0, 20);

即光线起点附近，对应一元四次方程的最高次系数非常小，如上拦截的方程对应最高次系数为10^-13数量级（double的最小正数DBL_EPSILON为10^-16数量级）。

取巧解决问题：

将图形在y轴方向与y=r.direction().y()，然后在程序中过滤掉的y=r.direction().y()实根。

临时解决。

顺便，下方贴出C++中double、float数据类型的部分常量。