这一节，画这个：
参考文献：
Blinn, J.F., A generalization of algebraic surfacedrawing. ACM Trans. Graph.
1(3) , 235-256, July 1982.
50.1 数学推导
引入密度函数（正态分布）：
函数图像如下：
该函数是关于x=u对称的，即任意到u距离相等的两个x值，其函数值是相等的。
现在将该函数拓展到三维空间。得到：
其实，量子力学中，原子中的电子表示为空间位置的密度函数，就是这个公式。球心为原子核，电子出现在某个位置的概率由该位置到原子核的距离决定。反过来，对于某个给定的概率值T，则其对应的到球心的距离（半径R）就确定，即，该概率下，电子只会出现在以该距离为半径的球面上。离球心越近，函数值越大；离球心越远，函数值越小。所以，该球面里面的点对应的函数值大于T，该球面外面的点对应的函数值小于T，该球面上的点对应的函数值等于T。
所以，方程
50.2 看C++代码实现
----------------------------------------------blobs.h ------------------------------------------
blobs.h
#ifndef BLOBS_H
#define BLOBS_H
#include <hitable.h>
#include "material.h"
#include "log.h"
class blobs : public hitable
{
public:
blobs() {}
blobs(vec3 cen1, vec3 cen2, float b1, float r1, float b2, float r2, float s, int in, float tol, material *m) :
center1(cen1), center2(cen2), blob_p1(b1), radius1(r1), blob_p2(b2), radius2(r2), sum(s),
initial_number(in), tolerance(tol), ma(m) {}
/*
f(x,y,z)= exp((B1/(R1^2) * ((x-x1)^2+(y-y1)^2+(z-z1)^2) - B1)
+ exp((B2/(R2^2) * ((x-x2)^2+(y-y2)^2+(z-z2)^2) - B2) - s = 0
//对应“式子十”
NOTE: in our program, x1=x2, z1=z2
s: should be bigger than 1
in: initial number
tol: tolerance
*/
virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;
vec3 center1, center2;
float blob_p1, radius1, blob_p2, radius2, sum;
int initial_number;
float tolerance;
material *ma;
};
#endif // BLOBS_H
----------------------------------------------blobs.cpp ------------------------------------------
blobs.cpp
#include "blobs.h"
#include <iostream>
#include <limits>
#include "float.h"
#include "log.h"
using namespace std;
bool ray_hit_box_b(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 BLOBS_LOG == 1
std::cout << "the ray is parallel to the planes and the origin X0 is not between the slabs. return false" <<endl;
#endif // BLOBS_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 BLOBS_LOG == 1
std::cout << "the ray is parallel to the planes and the origin Y0 is not between the slabs. return false" <<endl;
#endif // BLOBS_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 BLOBS_LOG == 1
std::cout << "the ray is parallel to the planes and the origin Z0 is not between the slabs. return false" <<endl;
#endif // BLOBS_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 BLOBS_LOG == 1
std::cout << "array1[" << i << "]:" << array1[i] <<endl;
std::cout << "array1[" << i+1 << "]:" << array1[i+1] <<endl;
#endif // BLOBS_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 BLOBS_LOG == 1
std::cout << "No.(0=X;2=Y;4=Z):" << i << " :t_near > t_far. return false" <<endl;
#endif // BLOBS_LOG
return false;
}
if(t_far < 0) {
#if BLOBS_LOG == 1
std::cout << "No.(0=X;2=Y;4=Z):" << i << " :t_far < 0. return false" <<endl;
#endif // BLOBS_LOG
return false;
}
}
if (t_near != t_near) {
t_near = t_near * 1;
}
return true;
}
bool get_blobs_function_and_derivative_b(const ray& r, vec3 c1, vec3 c2, float b1, float b2, float r1, float r2, float s, double t, double& f, double& fd) {
/*求函数值和函数导数值*/
double xo1 = double(r.origin().x() - c1.x());
double yo1 = double(r.origin().y() - c1.y());
double zo1 = double(r.origin().z() - c1.z());
double xo2 = double(r.origin().x() - c2.x());
double yo2 = double(r.origin().y() - c2.y());
double zo2 = double(r.origin().z() - c2.z());
double xd = double(r.direction().x());
double yd = double(r.direction().y());
double zd = double(r.direction().z());
double b1_r1 = double(b1/(r1*r1));
double b2_r2 = double(b2/(r2*r2));
double xo1_t = xo1+xd*t;
double yo1_t = yo1+yd*t;
double zo1_t = zo1+zd*t;
double xo2_t = xo2+xd*t;
double yo2_t = yo2+yd*t;
double zo2_t = zo2+zd*t;
double r1_2 = (xo1_t*xo1_t+yo1_t*yo1_t+zo1_t*zo1_t);
double r2_2 = (xo2_t*xo2_t+yo2_t*yo2_t+zo2_t*zo2_t);
double e1 = exp(b1_r1*(xo1_t*xo1_t+yo1_t*yo1_t+zo1_t*zo1_t)-double(b1));
double e2 = exp(b2_r2*(xo2_t*xo2_t+yo2_t*yo2_t+zo2_t*zo2_t)-double(b2));
f = e1 + e2 - double(s);
fd = e1*b1_r1*2*(xo1_t*xd+yo1_t*yd+zo1_t*zd)
+ e2*b2_r2*2*(xo2_t*xd+yo2_t*yd+zo2_t*zd);
/*分别对应“式子十二”，“式子十三”*/
if (e1 == 1.0) {
f = f*1;
}
return true;
}
bool get_roots_by_newton_iteration_b(const ray& r, vec3 c1, vec3 c2, float b1, float b2, float r1, float r2, float s, int in, float tol, float *x0, float (&roots)[2]) {
/*牛顿迭代*/
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_blobs_function_and_derivative_b(r, c1, c2, b1, b2, r1, r2, s, 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])) {
roots[j+1] = float(t_k1);
j++;
break;
}
else {
break;
}
}
else {
t_k = t_k1;
}
/*
}
else {
if (fabs(t_k1 - t_k) < tol) {
roots[j+1] = float(t_k1);
j++;
break;
}
else {
t_k = t_k1;
}
}
*/
}
else {
break;
}
}
else {
break;
}
}
if (j == 1) {
break;
}
}
roots[0] = float(j);
if (j == 0) {
}
return true;
}
bool blobs::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
#if BLOBS_LOG == 1
std::cout << "-------------blobs::hit----------------" << endl;
#endif // BLOBS_LOG
float box_blobs_x, box_blobs_y, box_blobs_z, center_y;
box_blobs_x = ((radius1 > radius2)? radius1:radius2);
box_blobs_y = (fabs(center1.y()-center2.y())+radius1+radius2)/2;
box_blobs_z = ((radius1 > radius2)? radius1:radius2);
if (center1.y() > center2.y()) {
center_y = ((center1.y()+radius1)+(center2.y()-radius2))/2;
}
else {
center_y = ((center2.y()+radius2)+(center1.y()-radius1))/2;
}
/*确定包围两个球的长方体，这些参数在“数学推导”部分有说明*/
vec3 vertex_l[1], vertex_h[1];
vertex_l[0] = vec3(center1.x()-box_blobs_x, center_y-box_blobs_y, center1.z()-box_blobs_z);
vertex_h[0] = vec3(center1.x()+box_blobs_x, center_y+box_blobs_y, center1.z()+box_blobs_z);
float roots[2] = {0.0, -1.0};
float x0[initial_number+1];
float t_near = 0;
float t_far = 0;
if (ray_hit_box_b(r, vertex_l[0], vertex_h[0], t_near, t_far)) {
if (initial_number == 1) {
x0[1] = t_near;
}
else {
for (int i=0; i<initial_number; i++) {
x0[i+1] = t_near + i*(t_far - t_near)/(initial_number-1);
}
}
x0[0] = float(initial_number);
get_roots_by_newton_iteration_b(r, center1, center2, blob_p1, blob_p2, radius1, radius2, sum, 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;
}
}
}
vec3 pc1, pc2;
double b1_r1 = double(blob_p1/(radius1*radius1));
double b2_r2 = double(blob_p2/(radius2*radius2));
double a1, a2, nx, ny, nz;
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);
pc1 = rec.p - center1;
pc2 = rec.p - center2;
a1 = b1_r1*2*exp(b1_r1*double(dot(pc1, pc1))-double(blob_p1));
a2 = b2_r2*2*exp(b2_r2*double(dot(pc2, pc2))-double(blob_p2));
nx = a1*double(pc1.x())+a2*double(pc2.x());
ny = a1*double(pc1.y())+a2*double(pc2.y());
nz = a1*double(pc1.z())+a2*double(pc2.z());
/*对应“式子十五”*/
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;
}
----------------------------------------------main.cpp ------------------------------------------
main.cpp
hitable *list[1];
list[0] = new blobs(vec3(5.2, 1.75, 0), vec3(5.2, 4.25, 0), -0.25, 1.2, -0.25, 1.2, 1.55, 5, 0.0001, new lambertian(vec3(0.0, 1.0, 0.0)));
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);
先解释一下这些参数：
list[0] = new blobs(vec3(5.2, 1.75, 0), vec3(5.2, 4.25, 0), -0.25, 1.2, -0.25, 1.2, 1.55, 5, 0.0001, new lambertian(vec3(0.0, 1.0, 0.0)));
球1：
中心：vec3(5.2, 1.75, 0)
半径：1.2
blob参数：-0.25
球2：
中心：vec3(5.2, 4.25, 0),
半径：1.2
blob参数：-0.25
分布叠加和（s）:1.55
迭代初值个数：5
容忍误差：0.0001
这里特别强调s值的调节，根据分布叠加图：
s值直接决定最终显示的图像（红色框内的部分）。
看一组图片：
vec3(5.2, 1.75, 0), vec3(5.2, 4.25, 0),-0.25, 1.2, -0.25, 1.2, 1.4,5, 0.0001,
vec3(5.2, 1.75, 0), vec3(5.2, 4.25, 0),-0.25, 1.2, -0.25, 1.2, 1.5,5, 0.0001,
vec3(5.2, 1.75, 0), vec3(5.2, 4.25, 0),-0.25, 1.2, -0.25, 1.2, 1.55,5, 0.0001,
如果S值太小，则图像会出现空白。
所以，对于图像中的空白处理，这里可以调节两个参数：
其一，增大S 值；
其二，增大迭代初值个数in。
接下来，看几组效果图：
第一组：半径为1，对应S值为1.35。
从左至右，blob参数一次为-4，-2，-1，-0.5，-0.25
vec3(5.2, 1.75, 0), vec3(5.2, 4.25, 0), -4, 1, -4, 1, 1.35, 5, 0.0001,
第二组：半径为1.2，对应S值为1.55。
从左至右，blob参数一次为-4，-2，-1，-0.5，-0.25
vec3(5.2, 1.75, 0), vec3(5.2, 4.25, 0), -4, 1.2, -4, 1.2, 1.55, 5, 0.0001,
第三组：半径为1.4，对应S值为1.75。
从左至右，blob参数一次为-4，-2，-1，-0.5，-0.25
vec3(5.2, 1.75, 0), vec3(5.2, 4.25, 0), -4, 1.4, -4, 1.4, 1.75, 5, 0.0001,
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
