tear drop是长这个样子的：
51.1 数学推导
在网上找到tear drop的参数方程：
51.2 看C++代码实现
----------------------------------------------vec3.h ------------------------------------------
vec3.h
/*重写了解一元二次方程、一元三次方程、一元四次方程的函数，将其所有数据类型由float换成double，同时使用double数据类型的最小正数“DBL_EPSILON”，绝对值小于该最小正数的数即为double数据类型的0*/
#define EPS 1e-10
double* 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 formula
double *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 same
if (!(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 root
double 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 same
double 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 roots
double 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 == 1
for (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_LOG
double 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 == 1
for (int i=0; i<(roots_y1[0]+1); i++) {
std::cout << "roots_y1[" << i << "]=" << roots_y1[i] << endl;
}
#endif // VEC3_LOG
roots_y2 = roots_quadratic_equation_rain(1.0, -k, m);
#if VEC3_LOG == 1
for (int i=0; i<(roots_y2[0]+1); i++) {
std::cout << "roots_y2[" << i << "]=" << roots_y2[i] << endl;
}
#endif // VEC3_LOG
if (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))/2
y=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 == 1
std::cout << "-------------rain::hit---1-------------" << endl;
#endif // RAIN_LOG
vec3 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数据类型的部分常量。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
