问题四十八:怎么用ray tracing画superhyperboloid(超级双曲面)
48.1 数学推导
Superhyperboloid的方程和“问题四十六”中描述的superellipsoid的如下方程非常接近。
到获得迭代格式,超级双曲面和超级椭圆面的推导是相似的,此处不表。
超级双曲面和超级椭圆面的差异之处在于:
超级双曲面是不封闭的面,所以,我们在y轴方向对曲面的范围进行限制。
其一:由于这个“限制”,导致迭代得到的解也必须在这个“限制”内。
其二:由于方程的改变,求迭代初值的几个box也有所变化。
其三:对迭代初值的选取要求更为苛刻(一不小心就少根了)。
48.2 看C++代码实现
----------------------------------------------superhyperboloid.h ------------------------------------------
superhyperboloid.h
#ifndef SUPERHYPERBOLOID_H
#define SUPERHYPERBOLOID_H#include <hitable.h>
#include "material.h"
#include "log.h"class superhyperboloid : public hitable
{public:superhyperboloid() {}superhyperboloid(vec3 cen, float a1, float a2, float a3, float e1, float e2, float s1, float s2, float hy, int in, float tol, material *m) :center(cen), intercept_x(a1), intercept_y(a2), intercept_z(a3), p_e1(e1), p_e2(e2),sign1(s1), sign2(s2), half_y(hy), initial_number(in), tolerance(tol), ma(m) {}
/*
f(x,y,z)=( (x/a1)^(2/e2) + s2*(z/a3)^(2/e2) )^(e2/e1) + s1*(y/a2)^(2/e1) -1 = 0
in: initial number//迭代初值个数
tol: tolerance
s1,s2: 1, 1: superellipsoid
s1,s2:-1, 1: superhyperboloids of one sheet
s1,s2:-1,-1: superhyperboloids of two sheets
hy: half height of the surface in y-direction
*/virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;vec3 center;float intercept_x, intercept_y, intercept_z;float p_e1, p_e2;float sign1, sign2;float half_y;int initial_number;float tolerance;material *ma;
};#endif // SUPERHYPERBOLOID_H
----------------------------------------------superhyperboloid.cpp ------------------------------------------
superhyperboloid.cpp
#include "superhyperboloid.h"#include <iostream>
#include <limits>
#include "float.h"
#include "log.h"using namespace std;bool ray_hit_box_h(const ray& r, const vec3& vertex_l, const vec3& vertex_h, float& t_near, float& 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 SUPERHYPERBOLOID_LOG == 1std::cout << "the ray is parallel to the planes and the origin X0 is not between the slabs. return false" <<endl;
#endif // SUPERHYPERBOLOID_LOGreturn 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 SUPERHYPERBOLOID_LOG == 1std::cout << "the ray is parallel to the planes and the origin Y0 is not between the slabs. return false" <<endl;
#endif // SUPERHYPERBOLOID_LOGreturn 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 SUPERHYPERBOLOID_LOG == 1std::cout << "the ray is parallel to the planes and the origin Z0 is not between the slabs. return false" <<endl;
#endif // SUPERHYPERBOLOID_LOGreturn 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 SUPERHYPERBOLOID_LOG == 1std::cout << "array1[" << i << "]:" << array1[i] <<endl;std::cout << "array1[" << i+1 << "]:" << array1[i+1] <<endl;
#endif // SUPERHYPERBOLOID_LOGif(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 SUPERHYPERBOLOID_LOG == 1std::cout << "No.(0=X;2=Y;4=Z):" << i << " :t_near > t_far. return false" <<endl;
#endif // SUPERHYPERBOLOID_LOGreturn false;}if(t_far < 0) {
#if SUPERHYPERBOLOID_LOG == 1std::cout << "No.(0=X;2=Y;4=Z):" << i << " :t_far < 0. return false" <<endl;
#endif // SUPERHYPERBOLOID_LOGreturn false;}}if (t_near != t_near) {t_near = t_near * 1;}return true;
}bool get_superellipsoid_function_and_derivative_h(float a1, float a2, float a3, float e1, float e2, float s1, float s2, 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 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;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+s2*pow_z) == 0) {
f = pow_y - 1;
fd = (2/e1_r)*pow_y_d*yd_a2;
}
elseif ((pow_x+s2*pow_z) < 0) {
//对于双页双曲面,s2为负,所以这里需要判断一下
f = -pow(-(pow_x+s2*pow_z), (e2_r/e1_r)) + s1*pow_y - 1;
fd = (2/e1_r)*(-pow(-(pow_x+s2*pow_z), ((e2_r/e1_r)-1))
*(pow_x_d*xd_a1+s2*pow_z_d*zd_a3)+ s1*pow_y_d*yd_a2);
}
else {f = pow(pow_x+s2*pow_z, (e2_r/e1_r)) + s1*pow_y - 1;fd = (2/e1_r)*(pow(pow_x+s2*pow_z, ((e2_r/e1_r)-1))
*(pow_x_d*xd_a1+s2*pow_z_d*zd_a3) + s1*pow_y_d*yd_a2);}return true;
}bool get_roots_by_newton_iteration_h(const ray& r, vec3 c, float a1, float a2, float a3, float e1, float e2, float s1, float s2, int in, float tol, float *x0, float (&roots)[3]) {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_h(a1, a2, a3, e1, e2, s1, s2, 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])) {
/*传进该函数的x0序列是从t_near到t_far排序好的,由于y轴方向的高度限制,所求的解必须在(t_near, t_far)范围内。By the way, 我们在这里值求一个解,这里求得的解不确定是大根还是小根。*/roots[j+1] = float(t_k1);j++;break;}else {break;}}else {t_k = t_k1;}}else {break;}}else {break;}}if (j == 1) {break;}}roots[0] = float(j);return true;
}bool superhyperboloid::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
#if SUPERHYPERBOLOID_LOG == 1std::cout << "-------------superhyperboloid::hit----------------" << endl;
#endif // SUPERHYPERBOLOID_LOGvec3 vertex_l[1], vertex_h[1];float intercept_hyper_x, intercept_hyper_z;float pow_y = pow((half_y/intercept_y), (2/p_e1));intercept_hyper_x = intercept_x*pow((pow_y+1), (p_e1/2));intercept_hyper_z = intercept_z*pow((pow_y+1), (p_e1/2));
/*y轴方向的半长度为half_y。在方程中,设y=half_y, x=0,可以求得intercept_hyper_z;在方程中,设y=half_y, z=0,可以求得intercept_hyper_x 。这里的intercept_hyper_x, half_y, intercept_hyper_z就相当于“问题四十六”中的a1, a2, a3*/vertex_l[0] = vec3(center.x()-intercept_hyper_x, center.y()-half_y, center.z()-intercept_hyper_z);vertex_h[0] = vec3(center.x()+intercept_hyper_x, center.y()+half_y, center.z()+intercept_hyper_z);float roots[3] = {0.0, -1.0, -1.0};float x0[initial_number+1];float t_near = 0;float t_far = 0;if (ray_hit_box_h(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的值作为迭代初值(序列)*/}x0[0] = float(initial_number);get_roots_by_newton_iteration_h(r, center, intercept_x, intercept_y, intercept_z, p_e1, p_e2, sign1, sign2, 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_d;float d_a1 = intercept_x;float d_a2 = intercept_y;float d_a3 = intercept_z;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);if (((sign1 == -1)&&(sign2 == 1)) || ((sign1 == -1)&&(sign2 == -1))) {if (!(((rec.p.y()-center.y()) >= -half_y)
&& ((rec.p.y()-center.y()) <= half_y))) {
/*其实这句是多余的,因为在迭代函数中已经用t_near, t_far对解进行了限制*/continue;}}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+sign2*pow_z) == 0) {pow_x_z_d = 0;}if ((pow_x+sign2*pow_z) < 0) {//当sign2为负时,需要判断一下pow_x_z_d = -pow(-(pow_x+sign2*pow_z), double(p_e2/p_e1-1));}else {pow_x_z_d = pow(pow_x+sign2*pow_z, double(p_e2/p_e1-1));}nx = double(2/p_e1) * pow_x_z_d * pow_x_d / d_a1;ny = double(2/p_e1) *sign1* pow_y_d / d_a2;nz = double(2/p_e1) * 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;
}
48.2.1 Superhyperboloid one sheet
----------------------------------------------main.cpp ------------------------------------------
main.cpp
hitable *list[1];list[0] = new superhyperboloid(vec3(0, 3, 0), 3, 4.5, 3, 3, 0.1, -1, 1, 4.5, 16, 0.001, new lambertian(vec3(0.0, 1.0, 0.0)));
/*half_y为4.5;迭代初值个数为16,只有将t_near, t_far等分为这么多份,才能保证图片基本不缺块,这就是之前说的“苛刻”的地方;误差为0.001*/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=0.1
e1=1, e2=0.1
e1=2, e2=0.1
e1=3, e2=0.1
e1=0.1, e2=1
e1=1, e2=1
e1=2, e2=1
e1=3, e2=1
e1=0.1, e2=2
e1=1, e2=2
e1=2, e2=2
e1=3, e2=2
e1=0.1, e2=3
e1=1, e2=3
e1=2, e2=3
e1=3, e2=3
48.2.2 Superhyperboloid two sheet
----------------------------------------------main.cpp ------------------------------------------
main.cpp
hitable *list[1];
list[0] = newsuperhyperboloid(vec3(0, 3, 0), 3, 4.5, 3, 0.1, 0.1, -1, -1,4.5, 16, 0.001, new lambertian(vec3(0.0, 1.0, 0.0)));
/*half_y为4.5;迭代初值个数为16,只有将t_near, t_far等分为这么多份,才能保证图片基本不缺块,这就是之前说的“苛刻”的地方;误差为0.001*/
hitable *world = newhitable_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);
输出图片如下:
双页双曲面的实现还有问题,估计迭代初值和误差等等还需要调。对于学习进度的考虑,在此不调了。
问题四十八:怎么用ray tracing画superhyperboloid(超级双曲面)相关推荐
- 问题四十七:怎么用ray tracing画superellipsoid (2)
网络上有另外一种更为简单的superellipsoid的方程形式:(这种方程,更快些) 数学推导和"问题四十六"类似(应该是更为简单),此处不表. 47.1 看C++代码实现 -- ...
- 问题五十:怎么用ray tracing画blobs
这一节,画这个: 参考文献: Blinn, J.F., A generalization of algebraic surfacedrawing. ACM Trans. Graph. 1(3) , 2 ...
- 问题六十:怎么用ray tracing画回旋体(rotational sweeping / revolution)
60.1 概述 回旋体,大概是长这个样子: 回旋体是指曲线(称为"基本曲线")围绕y轴转一圈得到的图形. (基本曲线是由多段b-spline曲线段连接而成) 这里先强调一下: 上图 ...
- 问题四十一:怎么用ray tracing画任意圆柱面(generalized cylinder)
我们之前在"35.2"章节中画过椭圆柱面: 我们还在"36.4"章节中画过圆柱面的Inverse Mapping图: 但是,这些柱面都是:底面与ZOX平面平行, ...
- OpenCV学习笔记(四十六)——FAST特征点检测features2D OpenCV学习笔记(四十七)——VideoWriter生成视频流highgui OpenCV学习笔记(四十八)——PCA算
OpenCV学习笔记(四十六)--FAST特征点检测features2D 特征点检测和匹配是计算机视觉中一个很有用的技术.在物体检测,视觉跟踪,三维常年关键等领域都有很广泛的应用.这一次先介绍特征点检 ...
- 【零基础学Java】—自定义异常(四十八)
[零基础学Java]-自定义异常(四十八) 自定义异常类:java提供的异常类,不够我们使用,需要自己定义一个异常类 格式: public class XXXException extends Exc ...
- JavaScript学习(四十八)—原型对象的增删改查
JavaScript学习(四十八)-原型对象的增删改查 一.构造方法与原型对象的图解 二.型对象的增删改查 (一).什么是原型 每个函数都会有一个属性--prototype属性,这个属性都会有一个对象 ...
- 罗永浩:我今年四十八岁,还可以承受无数次的失败;iOS14 或将推出系统级「小程序」功能; PyCharm新版发布| 极客头条...
整理 | 屠敏 头图 | CSDN 下载自视觉中国 快来收听极客头条音频版吧,智能播报由标贝科技提供技术支持. 「极客头条」-- 技术人员的新闻圈! CSDN 的读者朋友们早上好哇,「极客头条」来啦, ...
- 【Visual C++】游戏开发四十八 浅墨DirectX教程十六 三维地形系统的实现
分享一下我老师大神的人工智能教程!零基础,通俗易懂!http://blog.csdn.net/jiangjunshow 也欢迎大家转载本篇文章.分享知识,造福人民,实现我们中华民族伟大复兴! 本系列文 ...
最新文章
- Redis实际应用:限流
- 用javascript模拟分子扩散并思考熵与序
- shell获取git最近一次提交信息_Git修改commit提交信息
- 深入了解asp.net框架。生命周期以及事件处理机制
- JavaScript 中最​​重要的保留字
- sqlserver 2005 数据库的差异备份与还原
- 基于JAVA+SpringMVC+Mybatis+MYSQL的高校勤工助学管理系统
- java解析未知key json_获取JsonObject某一未知key的值操作
- optparse命令解析模块
- UE4之UMG用户界面
- Hadoop源码分析28 JobTracker 处理JobClient请求
- Android 基于监听的事件处理机制
- 华为手机6130失效_华为手机的拨号键这5个功能,用过的人都拍手叫好,绝不虚吹...
- java搭建rtmp服务器,利用docker搭建RTMP直播流服务器实现直播
- 学习Spring,这篇就够了
- 电脑关机程序(源码)
- 想更快成长更应该关注的博客
- 类似宝塔linux面板,类似宝塔面板的软件有没有呢?
- 二叉树的顺序存储和三种遍历(二)
- AES解密报错,Input length must be multiple of 16 when decrypting with padded cipher