首先，以球面为例学习怎么用ray tracing画参数方程表示的曲面；然后，再画一个牛角面。
特别说明：这一章节所画的曲面只是示意性的，所以先不care图片上的瑕疵。
53.1 数学推导
球面的参数方程f(φ,θ)如下：
画参数方程表示的球面的步骤如下：
第一步：
将φ分成n份，每一份为Δφ；将θ分成m份，每一份为Δθ。
每一对Δφ、Δθ对应球面上一小块曲面（如上图蓝色表示）。
所以，球面一共分成了n*m个这样的蓝色小块。
第二步：
根据蓝色小块的四个“顶点”的位置确定包围蓝色小块的最小长方体（如上图紫色长方体表示）。
所以，对应的小长方体个数也为n*m个。
“最小长方体”bl、bh两个顶点坐标的确定：
根据球面参数方程求得蓝色小块的四个“顶点”的坐标，选出其中最小的x、最小的y、最小的z作为bl的坐标，选出其中最大的x、最大的y、最大的z作为bh的坐标。
第三步：
光线撞击小长方体。只有撞击到小长方体的光线才有可能撞击到对应的蓝色小块。
所以，判断光线是否撞击球面，需要经过n*m次判断（光线是否撞击到这n*m个小长方体）。
光线撞击到长方体时，可以求得对应的t_near。
所以，每个长方体对应着一组：
第四步：
光线撞击到小长方体，不一定能够撞击到对应的蓝色小块。这里我们用三维牛顿迭代法来判断光线是否撞击到蓝色小块。若撞击到，求得其近似坐标值。（二维牛顿迭代对应的是切线，三维牛顿迭代对应的则是切平面）。
在牛顿迭代过程中会涉及到迭代终止条件。
下面重点分析其中的“第四步”
53.2 球面参数方程的各种微分的计算
53.3 牛角面参数方程的各种微分的计算
53.4 看C++代码实现
----------------------------------------------vec3.h ------------------------------------------
vec3.h
bool get_matrix_3_1(const vec3 a, float (&m)[3][1]);
bool get_matrix_3_3(const vec3 a, const vec3 b, const vec3 c, float (&m)[3][3]);
bool get_matrix_inverse_3_3(const float m[3][3], float (&inverse)[3][3]);
bool matrix_3_3_multiply_3_1(const float m1[3][3], const float m2[3][1], float (&m)[3][1]);
bool matrix_3_1_minus_3_1(const float m1[3][1], const float m2[3][1], float (&m)[3][1]);
----------------------------------------------vec3.cpp ------------------------------------------
vec3.cpp
bool get_matrix_3_1(const vec3 a, float (&m)[3][1]) {
/*将向量转换为矩阵*/
m[0][0] = a.x();
m[1][0] = a.y();
m[2][0] = a.z();
return true;
}
bool get_matrix_3_3(const vec3 a, const vec3 b, const vec3 c, float (&m)[3][3]) {
/*将向量转换为矩阵*/
m[0][0] = a.x();
m[1][0] = a.y();
m[2][0] = a.z();
m[0][1] = b.x();
m[1][1] = b.y();
m[2][1] = b.z();
m[0][2] = c.x();
m[1][2] = c.y();
m[2][2] = c.z();
return true;
}
bool get_matrix_inverse_3_3(const float m[3][3], float (&inverse)[3][3]) {
/*求3*3矩阵的逆矩阵*/
float det_m = m[0][0]*m[1][1]*m[2][2] + m[0][1]*m[1][2]*m[2][0] +
m[0][2]*m[1][0]*m[2][1] - m[0][2]*m[1][1]*m[2][0] -
m[0][0]*m[1][2]*m[2][1] - m[0][1]*m[1][0]*m[2][2];
/*求对应行列式的值*/
if (fabs(det_m) < 1e-6) {/*对应行列式的值不为零，才存在逆矩阵*/
return false;
}
else {
/*下面是求各种代数余子式，同时除以行列式的值*/
vec3 a = (1/det_m)*vec3(m[1][1]*m[2][2] - m[1][2]*m[2][1],
m[1][2]*m[2][0] - m[1][0]*m[2][2],
m[1][0]*m[2][1] - m[1][1]*m[2][0]);
vec3 b = (1/det_m)*vec3(m[0][2]*m[2][1] - m[0][1]*m[2][2],
m[0][0]*m[2][2] - m[0][2]*m[2][0],
m[0][1]*m[2][0] - m[0][0]*m[2][1]);
vec3 c = (1/det_m)*vec3(m[0][1]*m[1][2] - m[0][2]*m[1][1],
m[0][2]*m[1][0] - m[0][0]*m[1][2],
m[0][0]*m[1][1] - m[0][1]*m[1][0]);
get_matrix_3_3(a, b, c, inverse);
return true;
}
}
bool matrix_3_3_multiply_3_1(const float m1[3][3], const float m2[3][1], float (&m)[3][1]) {
/*3*3矩阵与3*1矩阵相乘*/
m[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0] + m1[0][2]*m2[2][0];
m[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0] + m1[1][2]*m2[2][0];
m[2][0] = m1[2][0]*m2[0][0] + m1[2][1]*m2[1][0] + m1[2][2]*m2[2][0];
return true;
}
bool matrix_3_1_minus_3_1(const float m1[3][1], const float m2[3][1], float (&m)[3][1]) {
/*矩阵相减*/
m[0][0] = m1[0][0] - m2[0][0];
m[1][0] = m1[1][0] - m2[1][0];
m[2][0] = m1[2][0] - m2[2][0];
return true;
}
----------------------------------------------parametric_surface.h ------------------------------------------
parametric_surface.h
#ifndef PARAMETRIC_SURFACE_H
#define PARAMETRIC_SURFACE_H
#include <hitable.h>
#include "material.h"
#include "log.h"
class parametric_surface : public hitable
{
public:
parametric_surface() {}
parametric_surface(vec3 cen, float a, float b, float c, float hy, material *m, float r, float t, float n) : center(cen), intercept_x(a), intercept_y(b), intercept_z(c), height_half_y(hy), ma(m), rho(r), theta(t), num(n) {}
/*rho即为ρ，一般为0.05；theta为相机的张角；num为长宽像素数目值中较大的那个值*/
virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;
vec3 center;
float intercept_x;
float intercept_y;
float intercept_z;
float height_half_y;
material *ma;
float rho, theta, num;
};
#endif // PARAMETRIC_SURFACE_H
----------------------------------------------parametric_surface.cpp ------------------------------------------
parametric_surface.cpp
#include "parametric_surface.h"
#include <iostream>
#include <limits>
#include "float.h"
#include "log.h"
using namespace std;
#define TYPE 2
/*
TYPE=1: sphere;
TYPE=2: horn;
*/
bool ray_hit_box_para(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 PARAMETRIC_SURFACE_LOG == 1
std::cout << "the ray is parallel to the planes and the origin X0 is not between the slabs. return false" <<endl;
#endif // PARAMETRIC_SURFACE_LOG
return 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 PARAMETRIC_SURFACE_LOG == 1
std::cout << "the ray is parallel to the planes and the origin Y0 is not between the slabs. return false" <<endl;
#endif // PARAMETRIC_SURFACE_LOG
return 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 PARAMETRIC_SURFACE_LOG == 1
std::cout << "the ray is parallel to the planes and the origin Z0 is not between the slabs. return false" <<endl;
#endif // PARAMETRIC_SURFACE_LOG
return 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 PARAMETRIC_SURFACE_LOG == 1
std::cout << "array1[" << i << "]:" << array1[i] <<endl;
std::cout << "array1[" << i+1 << "]:" << array1[i+1] <<endl;
#endif // PARAMETRIC_SURFACE_LOG
if(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 PARAMETRIC_SURFACE_LOG == 1
std::cout << "No.(0=X;2=Y;4=Z):" << i << " :t_near > t_far. return false" <<endl;
#endif // PARAMETRIC_SURFACE_LOG
return false;
}
if(t_far < 0) {
#if PARAMETRIC_SURFACE_LOG == 1
std::cout << "No.(0=X;2=Y;4=Z):" << i << " :t_far < 0. return false" <<endl;
#endif // VEC3_LOG
return false;
}
}
if (t_near != t_near) {
t_near = t_near * 1;
}
return true;
}
#if TYPE == 2
bool horn_parametric(const vec3 cen, const float a, const float b, const float c, const float u, const float v, vec3& parametric) {
/*求牛角曲面的参数方程*/
float sin_u = sin(2*M_PI*u);
float sin_v = sin(v);
float cos_u = cos(2*M_PI*u);
float cos_v = cos(v);
sin_u = (fabs(sin_u)<1e-6)? 0.0:sin_u;
sin_v = (fabs(sin_v)<1e-6)? 0.0:sin_v;
cos_u = (fabs(cos_u)<1e-6)? 0.0:cos_u;
cos_v = (fabs(cos_v)<1e-6)? 0.0:cos_v;
parametric = vec3(a*(2.0+u*cos_v)*sin_u+cen.x(),
c*((2+u*cos_v)*cos_u+2*u)+cen.y(),
b*u*sin_v+cen.z());
return true;
}
bool horn_parametric_f1(const ray& r, const float a, const float b, const float c, const float u, const float v, vec3& fu, vec3& fv, vec3& am) {
/*求牛角曲面的微分，为求Jacobi矩阵*/
float sin_u = sin(2*M_PI*u);
float sin_v = sin(v);
float cos_u = cos(2*M_PI*u);
float cos_v = cos(v);
sin_u = (fabs(sin_u)<1e-6)? 0.0:sin_u;
sin_v = (fabs(sin_v)<1e-6)? 0.0:sin_v;
cos_u = (fabs(cos_u)<1e-6)? 0.0:cos_u;
cos_v = (fabs(cos_v)<1e-6)? 0.0:cos_v;
fu = vec3(a*(2.0*M_PI*(2.0+u*cos_v)*cos_u+cos_v*sin_u),
c*(-2.0*M_PI*(2.0+u*cos_v)*sin_u+cos_v*cos_u+2.0),
b*(sin_v));
fv = vec3(a*(-u*sin_u*sin_v),
c*(-u*cos_u*sin_v),
b*(u*cos_v));
am = vec3(-r.direction().x(), -r.direction().y(), -r.direction().z());
return true;
}
bool horn_parametric_f1_f2(const float a, const float b, const float c, const float u, const float v, vec3& fu, vec3& fv, vec3& fuu, vec3& fuv, vec3& fvv) {
/*求牛角曲面的各种微分*/
float sin_u = sin(2*M_PI*u);
float sin_v = sin(v);
float cos_u = cos(2*M_PI*u);
float cos_v = cos(v);
sin_u = (fabs(sin_u)<1e-6)? 0.0:sin_u;
sin_v = (fabs(sin_v)<1e-6)? 0.0:sin_v;
cos_u = (fabs(cos_u)<1e-6)? 0.0:cos_u;
cos_v = (fabs(cos_v)<1e-6)? 0.0:cos_v;
fu = vec3(a*(2.0*M_PI*(2.0+u*cos_v)*cos_u+cos_v*sin_u),
c*(-2.0*M_PI*(2.0+u*cos_v)*sin_u+cos_v*cos_u+2.0),
b*(sin_v));
fv = vec3(a*(-u*sin_u*sin_v),
c*(-u*cos_u*sin_v),
b*(u*cos_v));
fuu = vec3(a*(-4.0*M_PI*M_PI*(2.0+u*cos_v)*sin_u+4.0*M_PI*cos_v*cos_u),
c*(-4.0*M_PI*(cos_v*sin_u+M_PI*(2.0+u*cos_v)*cos_u)),
0.0);
fuv = vec3(-a*sin_v*(sin_u+2*M_PI*u*cos_u),
-c*sin_v*(cos_u-2*M_PI*u*sin_u),
b*cos_v);
fvv = vec3(-a*u*sin_u*cos_v,
-c*u*cos_u*cos_v,
-b*u*sin_v);
return true;
}
bool horn_parametric_bl_bh(const vec3 cen, const float a, const float b, const float c, const float u1, const float u2, const float v1, const float v2, vec3& bl, vec3& bh) {
/*求牛角曲面的小长方体的bl、bh*/
float x[4], z[4], y[2];
float sin_u1 = sin(2*M_PI*u1);
float sin_u2 = sin(2*M_PI*u2);
float sin_v1 = sin(v1);
float sin_v2 = sin(v2);
float cos_u1 = cos(2*M_PI*u1);
float cos_u2 = cos(2*M_PI*u2);
float cos_v1 = cos(v1);
float cos_v2 = cos(v2);
sin_u1 = (fabs(sin_u1)<1e-6)? 0.0:sin_u1;
sin_u2 = (fabs(sin_u2)<1e-6)? 0.0:sin_u2;
sin_v1 = (fabs(sin_v1)<1e-6)? 0.0:sin_v1;
sin_v2 = (fabs(sin_v2)<1e-6)? 0.0:sin_v2;
cos_u1 = (fabs(cos_u1)<1e-6)? 0.0:cos_u1;
cos_u2 = (fabs(cos_u2)<1e-6)? 0.0:cos_u2;
cos_v1 = (fabs(cos_v1)<1e-6)? 0.0:cos_v1;
cos_v2 = (fabs(cos_v2)<1e-6)? 0.0:cos_v2;
float temp;
x[0] = a*((2+u1*cos_v1)*sin_u1)+cen.x();
x[1] = a*((2+u1*cos_v2)*sin_u1)+cen.x();
x[2] = a*((2+u2*cos_v1)*sin_u2)+cen.x();
x[3] = a*((2+u2*cos_v2)*sin_u2)+cen.x();
for (int i=0; i<3; i++) {
for (int j=i+1; j<4; j++) {
if (x[i] > x[j]) {
temp = x[i];
x[i] = x[j];
x[j] = temp;
}
}
}
y[0] = c*((2+u1*cos_v1)*cos_u1+2*u1)+cen.y();
y[1] = c*((2+u1*cos_v2)*cos_u1+2*u1)+cen.y();
y[2] = c*((2+u2*cos_v1)*cos_u2+2*u2)+cen.y();
y[3] = c*((2+u2*cos_v2)*cos_u2+2*u2)+cen.y();
for (int i=0; i<3; i++) {
for (int j=i+1; j<4; j++) {
if (y[i] > y[j]) {
temp = y[i];
y[i] = y[j];
y[j] = temp;
}
}
}
z[0] = b*u1*sin_v1+cen.z();
z[1] = b*u1*sin_v2+cen.z();
z[2] = b*u2*sin_v1+cen.z();
z[3] = b*u2*sin_v2+cen.z();
for (int i=0; i<3; i++) {
for (int j=i+1; j<4; j++) {
if (z[i] > z[j]) {
temp = z[i];
z[i] = z[j];
z[j] = temp;
}
}
}
bl = vec3(x[0], y[0], z[0]);
bh = vec3(x[3], y[3], z[3]);
return true;
}
#endif // TYPE
#if TYPE == 1
bool sphere_parametric(const vec3 cen, const float a, const float b, const float c, const float u, const float v, vec3& parametric) {
/*求球面的参数方程*/
float sin_u = sin(u);
float sin_v = sin(v);
float cos_u = cos(u);
float cos_v = cos(v);
sin_u = (fabs(sin_u)<1e-6)? 0.0:sin_u;
sin_v = (fabs(sin_v)<1e-6)? 0.0:sin_v;
cos_u = (fabs(cos_u)<1e-6)? 0.0:cos_u;
cos_v = (fabs(cos_v)<1e-6)? 0.0:cos_v;
parametric = vec3(a*sin_u*sin_v+cen.x(), c*cos_u+cen.y(), b*sin_u*cos_v+cen.z());
return true;
}
bool sphere_parametric_f1(const ray& r, const float a, const float b, const float c, const float u, const float v, vec3& fu, vec3& fv, vec3& am) {
/*求球面的微分，为求Jacobi矩阵*/
float sin_u = sin(u);
float sin_v = sin(v);
float cos_u = cos(u);
float cos_v = cos(v);
sin_u = (fabs(sin_u)<1e-6)? 0.0:sin_u;
sin_v = (fabs(sin_v)<1e-6)? 0.0:sin_v;
cos_u = (fabs(cos_u)<1e-6)? 0.0:cos_u;
cos_v = (fabs(cos_v)<1e-6)? 0.0:cos_v;
fu = vec3(a*cos_u*sin_v, -c*sin_u, b*cos_u*cos_v);
fv = vec3(a*sin_u*cos_v, 0, -b*sin_u*sin_v);
am = vec3(-r.direction().x(), -r.direction().y(), -r.direction().z());
return true;
}
bool sphere_parametric_f1_f2(const float a, const float b, const float c, const float u, const float v, vec3& fu, vec3& fv, vec3& fuu, vec3& fuv, vec3& fvv) {
/*求球面的各种微分*/
float sin_u = sin(u);
float sin_v = sin(v);
float cos_u = cos(u);
float cos_v = cos(v);
sin_u = (fabs(sin_u)<1e-6)? 0.0:sin_u;
sin_v = (fabs(sin_v)<1e-6)? 0.0:sin_v;
cos_u = (fabs(cos_u)<1e-6)? 0.0:cos_u;
cos_v = (fabs(cos_v)<1e-6)? 0.0:cos_v;
fu = vec3(a*cos_u*sin_v, -c*sin_u, b*cos_u*cos_v);
fv = vec3(a*sin_u*cos_v, 0, -b*sin_u*sin_v);
fuu = vec3(-a*sin_u*sin_v, -c*cos_u, -b*sin_u*cos_v);
fuv = vec3(a*cos_u*cos_v, 0, -b*cos_u*sin_v);
fvv = vec3(-a*sin_u*sin_v, 0, -b*sin_u*cos_v);
return true;
}
bool sphere_parametric_bl_bh(const vec3 cen, const float a, const float b, const float c, const float u1, const float u2, const float v1, const float v2, vec3& bl, vec3& bh) {
/*求球面的小长方体的bl、bh*/
float x[4], z[4], y[2];
float sin_u1 = sin(u1);
float sin_u2 = sin(u2);
float sin_v1 = sin(v1);
float sin_v2 = sin(v2);
float cos_u1 = cos(u1);
float cos_u2 = cos(u2);
float cos_v1 = cos(v1);
float cos_v2 = cos(v2);
sin_u1 = (fabs(sin_u1)<1e-6)? 0.0:sin_u1;
sin_u2 = (fabs(sin_u2)<1e-6)? 0.0:sin_u2;
sin_v1 = (fabs(sin_v1)<1e-6)? 0.0:sin_v1;
sin_v2 = (fabs(sin_v2)<1e-6)? 0.0:sin_v2;
cos_u1 = (fabs(cos_u1)<1e-6)? 0.0:cos_u1;
cos_u2 = (fabs(cos_u2)<1e-6)? 0.0:cos_u2;
cos_v1 = (fabs(cos_v1)<1e-6)? 0.0:cos_v1;
cos_v2 = (fabs(cos_v2)<1e-6)? 0.0:cos_v2;
float temp;
x[0] = a*sin_u1*sin_v1+cen.x();
x[1] = a*sin_u1*sin_v2+cen.x();
x[2] = a*sin_u2*sin_v1+cen.x();
x[3] = a*sin_u2*sin_v2+cen.x();
for (int i=0; i<3; i++) {
for (int j=i+1; j<4; j++) {
if (x[i] > x[j]) {
temp = x[i];
x[i] = x[j];
x[j] = temp;
}
}
}
z[0] = b*sin_u1*cos_v1+cen.z();
z[1] = b*sin_u1*cos_v2+cen.z();
z[2] = b*sin_u2*cos_v1+cen.z();
z[3] = b*sin_u2*cos_v2+cen.z();
for (int i=0; i<3; i++) {
for (int j=i+1; j<4; j++) {
if (z[i] > z[j]) {
temp = z[i];
z[i] = z[j];
z[j] = temp;
}
}
}
y[0] = ((c*cos_u1+cen.y())<(c*cos_u2+cen.y()))? (c*cos_u1+cen.y()):(c*cos_u2+cen.y());
y[1] = ((c*cos_u1+cen.y())<(c*cos_u2+cen.y()))? (c*cos_u2+cen.y()):(c*cos_u1+cen.y());
bl = vec3(x[0], y[0], z[0]);
bh = vec3(x[3], y[1], z[3]);
return true;
}
#endif // TYPE
int get_root_by_newton_iteration(const ray& r, vec3 cen, float a, float b, float c, vec3 x0, float tao, float uid, float vid, vec3& root) {
/*三维牛顿迭代*/
vec3 fun0, dis0, fu, fv, am, fun1, dis1, fuu, fuv, fvv, normal;
float m_d0[3][1], m_j[3][3], m_j_i[3][3], m_j_i_d0[3][1], m_x_k0[3][1], m_x_k1[3][1];
float d0, d1, d2, h;
get_matrix_3_1(x0, m_x_k0);
#if TYPE == 1
sphere_parametric(cen, a, b, c, m_x_k0[0][0], m_x_k0[1][0], fun0);
#endif // TYPE
#if TYPE == 2
horn_parametric(cen, a, b, c, m_x_k0[0][0], m_x_k0[1][0], fun0);
#endif // TYPE
dis0 = fun0 - r.point_at_parameter(m_x_k0[2][0]);
for (int i=0; i<20; i++) {
#if TYPE == 1
sphere_parametric_f1(r, a, b, c, m_x_k0[0][0], m_x_k0[1][0], fu, fv, am);
#endif // TYPE
#if TYPE == 2
horn_parametric_f1(r, a, b, c, m_x_k0[0][0], m_x_k0[1][0], fu, fv, am);
#endif // TYPE
get_matrix_3_1(dis0, m_d0);
get_matrix_3_3(fu, fv, am, m_j);
if (get_matrix_inverse_3_3(m_j, m_j_i)) {
/*迭代格式。对应上面推导过程中的“式子七”*/
matrix_3_3_multiply_3_1(m_j_i, m_d0, m_j_i_d0);
matrix_3_1_minus_3_1(m_x_k0, m_j_i_d0, m_x_k1);
// if (m_x_k1[2][0] > 1e-16)
{
#if TYPE == 1
sphere_parametric(cen, a, b, c, m_x_k1[0][0], m_x_k1[1][0], fun1);
#endif // TYPE
#if TYPE == 2
horn_parametric(cen, a, b, c, m_x_k1[0][0], m_x_k1[1][0], fun1);
#endif // TYPE
dis1 = fun1 - r.point_at_parameter(m_x_k1[2][0]);
if (dis1.length() < tao) {
#if TYPE == 1
sphere_parametric_f1_f2(a, b, c, m_x_k1[0][0], m_x_k1[1][0], fu, fv, fuu, fuv, fvv);
#endif // TYPE
#if TYPE == 2
horn_parametric_f1_f2(a, b, c, m_x_k1[0][0], m_x_k1[1][0], fu, fv, fuu, fuv, fvv);
#endif // TYPE
normal = unit_vector(cross(fu, fv));
d0 = dot(normal, fuu);
d1 = dot(normal, fuv);
d2 = dot(normal, fvv);
h = 0.5*(d0*uid*uid+2*d1*uid*vid+d2*vid*vid);
/*对应上面推导过程中的“式子八”*/
if (fabs(h) < tao) {
if (root.z() == 0.0) {
root = vec3(m_x_k1[0][0], m_x_k1[1][0], m_x_k1[2][0]);
}
else {
if (m_x_k1[2][0] < root.z()) {
root = vec3(m_x_k1[0][0], m_x_k1[1][0], m_x_k1[2][0]);
}
}
return 1;
}
else {
root = vec3(m_x_k1[0][0], m_x_k1[1][0], m_x_k1[2][0]);
return 2;//0;
}
}
else {
m_x_k0[0][0] = m_x_k1[0][0];
m_x_k0[1][0] = m_x_k1[1][0];
m_x_k0[2][0] = m_x_k1[2][0];
dis0 = dis1;
}
}
}
else {
return 0;
}
}
return 0;
}
bool parametric_surface::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
#if PARAMETRIC_SURFACE_LOG == 1
std::cout << "-------------parametric_surface::hit----------------" << endl;
#endif // PARAMETRIC_SURFACE_LOG
#if TYPE == 1
/*球面只需要分成4*8份，即可画出可接受的图片*/
#define UN 4
#define VN 8
#endif // TYPE
#if TYPE == 2
/*在当前程序中，牛角面分成80*80份，输出的图片还是有较多瑕疵*/
#define UN 80
#define VN 80
#endif // TYPE
float ui1, ui2, vi1, vi2, t_near, t_far, tao, u1, v1, uid, vid;
vec3 bl, bh, fu, fv, am, x0;
vec3 root = vec3(0.0, 0.0, 0.0);
int rn = 2;
int root_num = 0;
bool hitted1 = 0;
float t_near_a, t_far_a;
t_near_a = -1.0;
t_far_a = t_near_a;
bool hitted2 = 0;
#if TYPE == 1
#endif // TYPE
#if TYPE == 3
float delta_t[4] = = {-1.0, -1.0, -1.0, -1.0};
int min_num = 0;
#endif // TYPE
float ui1_a[4], ui2_a[4], vi1_a[4], vi2_a[4];
#if TYPE == 1
/*包围整个球面的长方体*/
bl = vec3(center.x()-intercept_x, center.y()-intercept_y, center.z()-intercept_z);
bh = vec3(center.x()+intercept_x, center.y()+intercept_y, center.z()+intercept_z);
#endif // TYPE
#if TYPE == 2
/*包围整个牛角面的长方体*/
bl = vec3(center.x()-intercept_x*2, center.y()-intercept_y*4, center.z()-intercept_z*2);
bh = vec3(center.x()+intercept_x*2, center.y()+intercept_y*4, center.z()+intercept_z*2);
#endif // TYPE
if(ray_hit_box_para(r, bl, bh, t_near, t_far)) {
for (int i=0; i<UN; i++) {
/*这两个for循环即是将曲面分成n*m个蓝色小块*/
#if TYPE == 1
ui1 = float(i)*M_PI/float(UN);
ui2 = float(i+1)*M_PI/float(UN);
#endif // TYPE
#if TYPE == 2
ui1 = float(i)/float(UN);
ui2 = float(i+1)/float(UN);
#endif // TYPE
for (int j=0; j<VN; j++) {
vi1 = float(j)*2*M_PI/float(VN);
vi2 = float(j+1)*2*M_PI/float(VN);
u1 = (ui1+ui2)/2;
v1 = (vi1+vi2)/2;
#if TYPE == 1
sphere_parametric_bl_bh(center, intercept_x, intercept_y, intercept_z, ui1, ui2, vi1, vi2, bl, bh);
#endif // TYPE
#if TYPE == 2
horn_parametric_bl_bh(center, intercept_x, intercept_y, intercept_z, ui1, ui2, vi1, vi2, bl, bh);
#endif // TYPE
hitted1 = ray_hit_box_para(r, bl, bh, t_near, t_far);
if(hitted1 && (t_near > 1e-7)) {
while (rn == 2) {
tao = fabs(2*rho*t_near*sin(theta*(M_PI/180)/2)/num);
x0 = vec3(u1, v1, t_near);
uid = ui2-ui1;
vid = vi2-vi1;
rn = get_root_by_newton_iteration(r, center, intercept_x, intercept_y, intercept_z, x0, tao, uid, vid, root);
if (rn == 1) {
root_num ++ ;
}
else if (rn == 2) {
ui1_a[0] = ui1;
ui2_a[0] = ui1+uid/2;
vi1_a[0] = vi1;
vi2_a[0] = vi1+vid/2;
ui1_a[1] = ui1;
ui2_a[1] = ui1+uid/2;
vi1_a[1] = vi1+vid/2;
vi2_a[1] = vi2;
ui1_a[2] = ui1+uid/2;
ui2_a[2] = ui2;
vi1_a[2] = vi1;
vi2_a[2] = vi1+vid/2;
ui1_a[3] = ui1+uid/2;
ui2_a[3] = ui2;
vi1_a[3] = vi1+vid/2;
vi2_a[3] = vi2;
for (int j=0; j<4; j++) {
ui1 = ui1_a[j];
ui2 = ui2_a[j];
vi1 = vi1_a[j];
vi2 = vi2_a[j];
#if TYPE == 1
sphere_parametric_bl_bh(center, intercept_x, intercept_y, intercept_z, ui1, ui2, vi1, vi2, bl, bh);
#endif // TYPE
#if TYPE == 2
horn_parametric_bl_bh(center, intercept_x, intercept_y, intercept_z, ui1, ui2, vi1, vi2, bl, bh);
#endif // TYPE
hitted2 = ray_hit_box_para(r, bl, bh, t_near_a, t_far_a);
if(hitted2 && (t_near_a > 1e-7)) {
// if (t_near_a == t_near) {
break;
// }
}
}
if (hitted2) {
t_near = root.z(); //t_near_a;
u1 = root.x();
v1 = root.y();
}
else {
break;
}
}
else if (rn == 0) {
return false;
}
}
if (root_num == 2) {
break;
}
}
}
if (root_num == 2) {
break;
}
}
}
else {
return false;
}
if (root.z() < t_max && root.z() > t_min) {
rec.t = root.z();
rec.p = r.point_at_parameter(rec.t);
vec3 pc = rec.p - center;
#if TYPE == 1
sphere_parametric_f1(r, intercept_x, intercept_y, intercept_z, root.x(), root.y(), fu, fv, am);
#endif // TYPE
#if TYPE == 2
horn_parametric_f1(r, intercept_x, intercept_y, intercept_z, root.x(), root.y(), fu, fv, am);
#endif // TYPE
/*在上面数学推导过程中的“式子八”求h时，有提到法向量的求法*/
rec.normal = unit_vector(cross(fu, fv));
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
/*先画个球面*/
int nx = 200;
int ny = 100;
int ns = 1;
// ofstream outfile( ".\\results\\superhyperboloid-two.txt", ios_base::out);
ofstream outfile( ".\\results\\parametric_surface.txt", ios_base::out);
outfile << "P3\n" << nx << " " << ny << "\n255\n";
std::cout << "P3\n" << nx << " " << ny << "\n255\n";
hitable *list[1]; list[0] = new parametric_surface(vec3(0.0, 3.0, 0.0), 3.0, 3.0, 3.0, 8.0, new lambertian(vec3(1.0, 0.0, 0.0)), 0.05, 20, 200);
hitable *world = new hitable_list(list,1);
vec3 lookfrom(0, 3, 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);
输出图片如下：
将采样次数ns设为100后的图片：
有明显的瑕疵，先不管了，重要的是方法，细节方面先不care。
----------------------------------------------main.cpp ------------------------------------------
main.cpp
/*再画个牛角面*/
int nx = 200;
int ny = 100;
int ns = 1;
// ofstream outfile( ".\\results\\superhyperboloid-two.txt", ios_base::out);
ofstream outfile( ".\\results\\parametric_surface.txt", ios_base::out);
outfile << "P3\n" << nx << " " << ny << "\n255\n";
std::cout << "P3\n" << nx << " " << ny << "\n255\n";
hitable *list[1];
list[0] = new parametric_surface(vec3(0.0, 3.0, 0.0), 3.0, 3.0, 3.0, 8.0, new lambertian(vec3(1.0, 0.0, 0.0)), 0.05, 70, 200);
hitable *world = new hitable_list(list,1);
vec3 lookfrom(0, 6, 20);
vec3 lookat(0, 6, 0);
float dist_to_focus = (lookfrom - lookat).length();
float aperture = 0.0;
camera cam(lookfrom, lookat, vec3(0,1,0), 70, float(nx)/float(ny), aperture, 0.7*dist_to_focus);
输出图片如下：
分成40*40份
分成80*40份
分成80*80份
分成100*100份
分成200*100份
分成300*100份
每张图片都是有的地方多了像素点，有的地方少了像素点。Anyway，“方法比细节重要”at this moment.
最后贴一张球面和牛角面的大小对比图（是分两次输出的两张图片，然后用画图工具合到了一块）
扫一扫
3179
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
