[图形学]smallpt代码详解(2)
一、简介
本文紧接在[图形学]smallpt代码详解(1)之后,继续详细讲解smallpt中的代码,包括自定义函数(第41到47行)
和递归路径跟踪函数(第48到74行)
部分。
二、smallpt代码详解
1.自定义函数(第41到47行)
smallpt源代码在第41到47行声明定义了clamp()
、toInt()
和intersect()
三个函数:
double clamp(double x)
:inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
clamp()
函数将输入参数x
限制在[0,1]
范围内,避免浮点数计算出错,出现小于0.0或者大于1.0的计算结果。int toInt(double x)
:inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }
toInt()
函数将输入的在[0,1]
范围内的输入参数x
先进行伽马矫正(Gamma Correction),伽马校正线性颜色空间转换为 sRGB 颜色空间。伽马校正有助于提高图像在显示器上的视觉效果。,然后映射到范围[0,255]
。bool intersect(const Ray &r, double &t, int &id)
:inline bool intersect(const Ray &r, double &t, int &id){ // n 为场景中的球面物体个数 // d 用来记录场景中与光线相交的第一个物体id // t 用来记录场景中与光线相交的第一个交点到光线原点的距离 double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20; // 遍历场景中的所有球面对象,然后使用 sphere[i] 对象中的函数interscet() 函数计算 // 光线是否与 球面对象 sphere[i] 相交,并找到距离光线原点最近的交点距离 t 和球面对象 i for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;} return t<inf; }
intersect()
函数用于计算光线r
与整个渲染场景中所有球面对象的相交计算,如果光线r与场景存在交点,则返回true,并且将交点到光线原点的距离t,写入输入的变量double &t中,并且返回相交球面物体的id,写入输入的变量int &id中。
代码首先确定场景中的球面物体个数n:
然后遍历场景中的所有球面对象,然后使用double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20;
sphere[i]
对象中的interscet()
函数计算光线是否与球面对象sphere[i]
相交,对于sphere[i]
对象中的intersect()
函数的讲解请查看[图形学]smallpt代码详解(1)中的三、smallpt代码详解-1.自定义数据结构部分(第4到第22行)部分,并找到距离光线原点最近的交点距离t
和球面对象i
。for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;}
2.递归路径跟踪函数(第48到74行)
radiance()
函数是光线跟踪(路径跟踪)中最主要的代码,该函数使用递归的方式实现路径跟踪,对场景进行渲染。
函数Vec radiance(const Ray &r, int depth, unsigned short *Xi)
的返回值和参数介绍如下:
- 返回值
Vec
:radiance() 函数计算的光线在此次传播中的辐射值; - 参数
r
:要处理的光线 - 参数
depth
:此光线的传播深度 - 参数
Xi
:随机数,作为随机采样种子值
radiance()函数
整体可以分为四个部分是否结束递归(与场景无交点或者俄罗斯轮盘赌失败)(红色)
、处理漫反射(绿色)
、处理镜面反射(蓝色)
和处理折射(橙色)
部分。流程如下:
radiance()函数
代码流程:
下面介绍各部分的详细代码实现:
2.1) 是否结束递归(与场景无交点或者俄罗斯轮盘赌失败)
:
double t; // distance to intersection
int id=0; // id of intersected object
// 计算当前光线 r 和场景的交点,如果没有交点,则返回 Vec(),即返回 (0,0,0)
if (!intersect(r, t, id)) return Vec(); // if miss, return black
// 如果 r 与场景有交点,那么记录交点对应的球面对象 obj
const Sphere &obj = spheres[id]; // the hit object
// x 为光线 r 与 obj的交点位置
// n 为交点 x 处的球面法向
// nl 为定向法向,用来描述光线r是从球外射向球内,还是从球外射向球内
// 假如 r 为从球外射向球内,则 nl 为朝向球外的法向
// 假如 r 为从球内射向球外,则 nl 为朝向球内的法向
// f 为 obj 的漫辐射颜色
Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
/* 使用俄罗斯轮盘赌策略确定是否结束递归 */
double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
if (++depth>5) if (erand48(Xi)<p) f=f*(1/p); else return obj.e; //R.R.
代码中首先判断光线r是否与场景相交:
- 1). 如果不相交则说明光线射出到场景之外,因此直接返回
(0,0,0)
。 - 2). 如果光线r与场景相交,则计算得到交点位置
x
,交点所处的球面对象obj
,和交点处的法向n
,与r在球面同边的法向nl
。- 2.1) 如果此时的递归深度大于5,那么使用俄罗斯轮盘赌确定是否结束递归。在俄罗斯轮盘赌策略中,smallpt使用交点处颜色f的(r,g,b)三个颜色分量中的最大值作为阈值
p
,然后使用erand48()
函数生成一个在[0,1]范围内的随机数。此处之所以使用颜色(r,g,b)
中最大值作为阈值p
是考虑假如交点处材质颜色接近黑色,可以尽可能地结束递归,因为此时(r,g,b)
三个分量值都很小,因此该点最终反射的radiance很小,对渲染地结果影响也很小。反之,该交点颜色接近白色,说明该交点最终反射的radiance可能很大,因此要尽可能地让光线继续传播下去。- 2.1.1) 假如该随机数小于
p
,说明赢得此次俄罗斯轮盘赌,将交点处的颜色值f变为原来的(1/p)
倍数,然后继续往下运行; - 2.1.2) 假如该随机数大于
p
,说明败于此次俄罗斯轮盘赌,直接返回交点处的自发光颜色obj.e
。
- 2.1.1) 假如该随机数小于
- 2.2) 如果此时的递归深度小于等于5,则继续往下运行。
- 2.1) 如果此时的递归深度大于5,那么使用俄罗斯轮盘赌确定是否结束递归。在俄罗斯轮盘赌策略中,smallpt使用交点处颜色f的(r,g,b)三个颜色分量中的最大值作为阈值
2.2) 处理漫反射
:
if (obj.refl == DIFF){ // Ideal DIFFUSE reflection
// r1 是在范围 [0,2Pi] 范围内的随机数
// r2 是在范围 [0,1] 范围内的随机数
double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
return obj.e + f.mult(radiance(Ray(x,d),depth,Xi));
在处理漫反射材质表面时,在交点处的上半球面上的cos加权采样
,采样得到光线r对应的入射方向d
,然后递归地计算光线下一次传播radiance(Ray(x,d),...)
的结果。
半球面上的cos加权采样
,即各方向的采样概率密度与
c
o
s
(
θ
)
cos(\theta)
cos(θ)成正比,
θ
\theta
θ为采样方向与交点处法向nl
的夹角。
采样步骤如下:
- 首先采样方位角
r1
和半球对应的圆平面半径长度r2s
。 - 然后计算以交点处法向
nl
为(0,0,1)轴的坐标系(u,v,w)
。 - 再之后根据采样的方位角
r1
和圆平面半径长度r2s
在u-v平面的圆平面上采样得到采样点p,该方法即在平面u-v上的圆平面内均匀采样得到采样p。采样点p对应到半球面上的映射点P的向量即为目标方向向量d
。
该方法与半球面上的cos加权采样
等效的原因请查看博客[图形学]在半球面上均匀采样和cos加权采样中的3.在半球面上的cos加权采样:方法二
部分。
采样点p在(u,v,w)坐标系下的坐标为:
(
c
o
s
(
r
1
)
∗
r
2
s
,
s
i
n
(
r
1
)
∗
r
2
s
,
0
)
(cos(r1)*r2s, sin(r1)*r2s, 0)
(cos(r1)∗r2s,sin(r1)∗r2s,0),
映射点P在(u,v,w)坐标系下的坐标即为:
(
c
o
s
(
r
1
)
∗
r
2
s
,
s
i
n
(
r
1
)
∗
r
2
s
,
s
q
r
t
(
1
−
r
2
s
2
)
)
=
(
c
o
s
(
r
1
)
∗
r
2
s
,
s
i
n
(
r
1
)
∗
r
2
s
,
s
q
r
t
(
1
−
r
2
)
)
(cos(r1)*r2s, sin(r1)*r2s, sqrt(1-r2s^2))=(cos(r1)*r2s, sin(r1)*r2s, sqrt(1-r2))
(cos(r1)∗r2s,sin(r1)∗r2s,sqrt(1−r2s2))=(cos(r1)∗r2s,sin(r1)∗r2s,sqrt(1−r2))
那么在世界坐标系(x,y,z)下,半球面上的映射点P坐标即为:
(
u
∗
c
o
s
(
r
1
)
∗
r
2
s
,
v
∗
s
i
n
(
r
1
)
∗
r
2
s
,
w
∗
s
q
r
t
(
1
−
r
2
)
)
(u*cos(r1)*r2s,v*sin(r1)*r2s,w*sqrt(1-r2))
(u∗cos(r1)∗r2s,v∗sin(r1)∗r2s,w∗sqrt(1−r2))
由于我们是使用交点处的法向nl
作为坐标系w
轴,暗含了以该交点为坐标系原点,因此映射点P的坐标即采样向量,即代码中的变量d
。
最后返回:该交点处的自发光obj.e
+该交点处的颜色值f
乘以光线下一次传播计算得到的radiance值。
2.3) 处理镜面反射
:
} else if (obj.refl == SPEC) // Ideal SPECULAR reflection
return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
处理镜面反射相对比较简单,根据光线r
方向和交点处法向n
可以方便地计算得到入射光线方向。
根据反射光方向和法向计算入射光方向(或者根据入射光方向和法向计算反射光方向)的算法示意图如下所示:
最后返回:该交点处的自发光obj.e
+该交点处的颜色值f
乘以光线下一次传播计算得到的radiance值。
2.4) 处理折射
:
Ray reflRay(x, r.d-n*2*n.dot(r.d)); // Ideal dielectric REFRACTION
bool into = n.dot(nl)>0; // Ray from outside going in?
double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0) // Total internal reflection
return obj.e + f.mult(radiance(reflRay,depth,Xi));
Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ? // Russian roulette
radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) :
radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
处理折射时,首先判断光线是否发生全内反射
,如果发生全内反射,那么就直接根据理想的反射方向继续在玻璃内传播。如果没有发生全内反射
则基于斯涅尔定律计算折射方向
,然后使用 菲涅尔方程和Schlick近似 计算折射和反射的比例
。
-
全内反射判断
:
根据全内反射定律,当入射角 θ a > θ c \theta_{a}>\theta_{c} θa>θc时发生全内反射。而 θ c = a r c s i n ( n b n a ) \theta_{c}=arcsin(\frac{n_b}{n_a}) θc=arcsin(nanb);
即当:
θ a > θ c = a r c s i n ( n b n a ) θ a > a r c s i n ( n b n a ) s i n ( θ a ) > n b n a s i n 2 ( θ a ) > ( n b n a ) 2 1 − c o s 2 ( θ a ) > ( n b n a ) 2 ( n a n b ) 2 ( 1 − c o s 2 ( θ a ) ) > 1 1 − ( n a n b ) 2 ( 1 − c o s 2 ( θ a ) ) < 0 (1) \theta_{a} > \theta_{c}=arcsin(\frac{n_b}{n_a}) \\\ \theta_{a} > arcsin(\frac{n_b}{n_a}) \\ sin(\theta_{a}) > \frac{n_b}{n_a} \\ sin^2(\theta_{a}) > (\frac{n_b}{n_a})^2 \\ 1-cos^{2}(\theta_{a})>(\frac{n_b}{n_a})^2 \\ (\frac{n_a}{n_b})^2(1-cos^{2}(\theta_{a}))>1 \\ 1 - (\frac{n_a}{n_b})^2(1-cos^{2}(\theta_{a})) < 0 \tag{1} θa>θc=arcsin(nanb) θa>arcsin(nanb)sin(θa)>nanbsin2(θa)>(nanb)21−cos2(θa)>(nanb)2(nbna)2(1−cos2(θa))>11−(nbna)2(1−cos2(θa))<0(1)
时,发生全内反射。
而代码中nnt
即 n a n b \frac{n_a}{n_b} nbna,ddn
等于 − c o s ( θ a ) -cos(\theta_{a}) −cos(θa),ddn*ddn
即 c o s 2 ( θ a ) cos^{2}(\theta_{a}) cos2(θa)。
因此,当cos2t=1 - nnt * nnt * (1 - ddn * ddn)<0
时,发生全内反射。 -
计算折射方向
:
根据斯涅尔定律,假如入射光方向为 D D D,折射光方向为 T T T,那么应该满足如下公式:
n a ∗ s i n ( θ a ) = n b ∗ s i n ( θ b ) (2) n_{a}*sin(\theta_{a}) = n_{b}*sin(\theta_{b}) \tag{2} na∗sin(θa)=nb∗sin(θb)(2)
其中 θ a \theta_{a} θa为入射角, t h e t a b theta_{b} thetab为折射角, n a n_a na为入射介质折射率, n b n_b nb为出射介质折射率。
如上图所示,我们可以得到:
s i n ( θ a ) = 1 − c o s 2 ( θ a ) = 1 − ( D ∗ N ) 2 s i n ( θ b ) = n a n b s i n ( θ a ) = n a n b 1 − ( D ∗ N ) 2 c o s ( θ b ) = 1 − s i n 2 ( θ b ) = 1 − n a 2 n b 2 ( 1 − ( D ∗ N ) 2 ) (3) sin(\theta_a)=\sqrt{1-cos^2(\theta_a)}=\sqrt{1-(D*N)^2}\\ sin(\theta_b)=\frac{n_a}{n_b}sin(\theta_a) = \frac{n_a}{n_b}\sqrt{1-(D*N)^2} \\ cos(\theta_b) = \sqrt{1-sin^2(\theta_b)}=\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D*N)^2)} \tag{3} sin(θa)=1−cos2(θa)=1−(D∗N)2sin(θb)=nbnasin(θa)=nbna1−(D∗N)2cos(θb)=1−sin2(θb)=1−nb2na2(1−(D∗N)2)(3)
单位向量 B B B等于:
B = n o r m ( D − ∣ c o s ( θ a ) ∣ N ) = D − ∣ c o s ( θ a ) ∣ N s i n ( θ a ) = D + ( D ⋅ N ) N 1 − ( D ⋅ N ) 2 (4) B = norm(D - |cos(\theta_a)| N) = \frac{D - |cos(\theta_a)| N}{sin(\theta_a)}=\frac{D +(D\cdot N)N}{\sqrt{1-(D\cdot N)^2}} \tag{4} B=norm(D−∣cos(θa)∣N)=sin(θa)D−∣cos(θa)∣N=1−(D⋅N)2D+(D⋅N)N(4)
目标折射光向量 T T T等于:
T = B s i n ( θ b ) − N c o s ( θ b ) = D + ( D ⋅ N ) N 1 − ( D ⋅ N ) 2 ∗ n a n b 1 − ( D ⋅ N ) 2 − N 1 − n a 2 n b 2 ( 1 − ( D ⋅ N ) 2 ) = n a n b ∗ ( D + D ⋅ N ) N − N 1 − n a 2 n b 2 ( 1 − ( D ⋅ N ) 2 ) = n a n b D + n a n b N ( D ⋅ N ) − N 1 − n a 2 n b 2 ( 1 − ( D ⋅ N ) 2 ) = n a n b D + N ( n a n b ( D ⋅ N ) − 1 − n a 2 n b 2 ( 1 − ( D ⋅ N ) 2 ) ) (5) T = Bsin(\theta_b)-Ncos(\theta_b) = \frac{D +(D\cdot N)N}{\sqrt{1-(D\cdot N)^2}} * \frac{n_a}{n_b}\sqrt{1-(D \cdot N)^2} - N\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2)} \\ = \frac{n_a}{n_b}*(D+D\cdot N)N-N\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2)} \\ = \frac{n_a}{n_b}D+\frac{n_a}{n_b}N(D\cdot N)-N\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2)} \\ = \frac{n_a}{n_b}D+N \left( \frac{n_a}{n_b}(D\cdot N)-\sqrt{1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2) } \right) \tag{5} T=Bsin(θb)−Ncos(θb)=1−(D⋅N)2D+(D⋅N)N∗nbna1−(D⋅N)2−N1−nb2na2(1−(D⋅N)2)=nbna∗(D+D⋅N)N−N1−nb2na2(1−(D⋅N)2)=nbnaD+nbnaN(D⋅N)−N1−nb2na2(1−(D⋅N)2)=nbnaD+N(nbna(D⋅N)−1−nb2na2(1−(D⋅N)2))(5)
因为代码中nnt
= n a n b \frac{n_a}{n_b} nbna,ddn
= − c o s ( θ a ) = − ( D ⋅ N ) -cos(\theta_{a})=-(D\cdot N) −cos(θa)=−(D⋅N),cos2t=1 - nnt * nnt * (1 - ddn * ddn)
= 1 − n a 2 n b 2 ( 1 − ( D ⋅ N ) 2 ) 1-\frac{n_{a}^2}{n_b^2}(1-(D\cdot N)^2) 1−nb2na2(1−(D⋅N)2)。
因此,结果折射光线 T T T就等于代码中的:Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
-
计算折射和反射的比例
:
根据菲涅尔方程和Schlick近似,可以计算得到光线发生折射和反射的比例,Schlick近似公式如下:
n = n a n b R 0 = ( n − 1 ) 2 ( n + 1 ) 2 = ( n a − n b ) / n b ( n a + n b ) / n b = n a − n b n a + n b R r ( θ a ) = R o + ( 1 − R o ) ( 1 − c o s ( θ a ) ) 5 (6) n=\frac{n_a}{n_b} \\ R_{0}=\frac{(n-1)^2}{(n+1)^2}=\frac{(n_a-n_b)/n_b}{(n_a+n_b)/n_b}=\frac{n_a-n_b}{n_a+n_b} \\ R_{r}(\theta_{a}) = R_{o} + (1-R_{o})(1-cos(\theta_{a}))^{5} \tag{6} n=nbnaR0=(n+1)2(n−1)2=(na+nb)/nb(na−nb)/nb=na+nbna−nbRr(θa)=Ro+(1−Ro)(1−cos(θa))5(6)
其中 R 0 R_{0} R0为光线垂直入射(即入射角 θ a = 0 \theta_{a}=0 θa=0)时的反射率,即代码中的变量R0
, R r ( θ a ) R_{r}(\theta_{a}) Rr(θa)为入射角为 θ a \theta_{a} θa时的反射率,即代码中的变量Re
。double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
-
其他
:return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ? // Russian roulette radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) : radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
在结束折射部分代码时,由于smallpt中规定折射材质的颜色
c=(1.0,1.0,1.0)
,在radiance()
函数开始,使用p=max(f.r, f.g,f.b)
进行俄罗斯轮盘赌时肯定无法失败、结束递归。因此在此处增加了一次俄罗斯轮盘赌,使用P=0.25+0.5*Re
作为俄罗斯轮盘赌的阈值。
当递归深度大于2时使用P
作为阈值进行俄罗斯轮盘赌,决定是否结束递归。如果不结束,就将Re
比例的光线用于接下来的镜面反射,Tr=1-Re
比例的光线用于接下来的折射。
在接下来的 [图形学]smallpt代码详解(3) 中,将继续讲解 smallpt
中的main函数部分
,包括渲染场景的定义
、相机的设置
、光线的生成
、滤波处理减少噪点
和保存最终渲染结果
几部分。
三、参考
[1].smallpt: Global Illumination in 99 lines of C++
[2].smallpt: Global Illumination in 99 lines of C+±Presentation slides
[3].光线跟踪smallpt详解 (一)
[4].斯涅尔定律
[5].菲涅尔方程
[6].Schlick近似
原文地址:https://blog.csdn.net/Strengthennn/article/details/142711240
免责声明:本站文章内容转载自网络资源,如本站内容侵犯了原著者的合法权益,可联系本站删除。更多内容请关注自学内容网(zxcms.com)!