1 前言
这次的光线追踪，主要是基于 Radiometry 的一种实现，也就是基于物理的一种实现。本文对辐射度量学做了解释，同时给出了程序中的关键代码以及参考资料，实现了微表面模形，模形的代码正确无误，但是运用到路径追踪上会出现很强的噪声，这个项目无法解决。
2 辐射度量学
关于辐射度量学，这里给出一种理解方式，
- Radiant energy (barely used in CG) \(Q [J = Joule]\) ：The energy of electromagnetic radiation
Raiant energy说白了就是能量。可以类比于位置，假设初始的位置是 0 m， 最终的位置是 100 m。我们可以说，兔子一开始在0，最后在100， 我们也可以说乌龟一开始在0，最终在 100。很显然，位置没有什么用，我们必须要把位置和时间结合起来，才可以获得更加直观和直觉的信息。
我们不会说自行车和汽车都能走100公里，而是会说走这100公里的速度。
一个10w的灯泡可以放出1000J的能量，可这个时间就可能是开10个小时。一个100w的灯泡也可以放出1000J的能量，这个时间可能是1个小时。我们直觉的感到，10w的灯泡使用的时候会比100w的暗淡，如何描述这种性质呢？就引出了
- Radiant flux(power) $ \Phi \equiv \frac {dQ}{dt} $ [W=Watt][lm=lumen] ： Energy per unit time
我们可以采用量子化的方法，把一份能量看作一个小球。而 Radiant flux(power) 就是描述 一秒内 能射出出多少个小球。
这可以使用向量来表示，小球的射出方向就是向量的方向，小球的大小就是向量的大小。需要注意的一点是，在物理中，向量只有大小和方向，没有所谓的固定的起始点。
Radiant flux(power)就是包蓝色的小球。一包 power 里有许多 能量(黄色的小球)，我们用蓝色的时间包裹，就成了power。这是辐射度量学里常用的东西，也是基本的东西。而 power 比 energy更好描述这个世界，我们可以将它比喻为力量感或者什么东西，它代表了某种强度。
-
Radiant intensity $I(\omega) \equiv \frac {d\Phi }{d\omega } $ :power per unit solid angle。
-
Solid Angle $ \Omega $ = $ \frac {A}{r^ {2}} $ ：ratio of subtended area on sphere to radius squared
一个物体某个点发出的黄色的小球会朝着空间中四散而去, 在一个锥形的固体角中，一秒中的黄色小球数量。
需要注意的是，上述公式中用到了\(d\), 我们知道这是一个极限的过程，立体角的定义是球面上某处面积与整个球面面积之比，当处于极限过程时，一块面积趋近于一个非常非常小的点，而这时立体角就可以看作是一个箭头。有了这个思想，就接下来就好做多了。
光线可以理解为一个糖葫芦。穿着power。
- Irradiance
Definition: The irradiance is the power per unit area incident on a surface point.
可以理解为，四面八方来的光线这个面积上，可以引申为这个点上，由于极限的关系。
- Radiance
Definition: The radiance is the power emitted, reflected,transmitted or received by a surface, per unit solid angle, per projected unit area.
Radiance是光线Ray传播时(即不与物体交互时)的一个性质。他在物体发出光线，或者光线打到物体上时会有折损，这种折损与物体表面的法线方向与传播光线的方向有关** Lambert’s Cosine Law**。所以我们用Radiance表示光线的传播。
一个光线带有100power，打到物体上由于角度关系我们只得到了50power, 但是我们不能说光线是50power的。
接收到的 Irradiance 由于在接收时，会因为 Lambert’s Cosine Law而进行缩小， 所以在还原时，需要除以 \(\cos \theta\)进行放大，还原。
3 Coding
3.1 uniform random point in triangle in 3D
// Triangle::Sample()
void Sample(Intersection &pos, float &pdf){
float x = std::sqrt(get_random_float()), y = get_random_float();
pos.coords = v0 * (1.0f - x) + v1 * (x * (1.0f - y)) + v2 * (x * y);
pos.normal = this->normal;
pdf = 1.0f / area;
}
在Triangle::Sample()里给出了三角形内随机采样某一点的算法，来历如下。
https://www.cs.princeton.edu/~funk/tog02.pdf
https://math.stackexchange.com/questions/18686/uniform-random-point-in-triangle-in-3d
3.2 Multithreading
int spp = 64;
std::vector<std::thread> threads;
std::vector<Vector3f> colors(spp);
/* Muti-threads implenmetation for each Ray
Ray ray(eye_pos, dir);
for(int i = 0; i < spp; i++){
threads.emplace_back( std::thread(&Renderer::RT_thread, this, std::ref(colors[i]), std::ref(scene), ray) );
}
for(int i = 0; i < threads.size(); i++){
threads[i].join();
}
for (int k = 0; k < spp; k++){
framebuffer[m] += colors[k] / spp;
}
threads.clear();
*/
https://www.educative.io/blog/modern-multithreading-and-concurrency-in-cpp
3.3 Eval()
使用了 Diffuse / Lambertian Material 的BRDF Lecture17 P8。总计两条假设
-
入射的radiance为常数。Suppose the incident lighting is uniform
-
各个方向的出射radiance为常数 Light is equally reflected in each output direction
关于 ρ=1 时，就是所谓的入射与反射的能量相同。设一条入射光带着 x 能量打到表面，这些能量会均匀的向四面散开，这时某一方向的出射光能量很少，而有无数条从任意方向来的入射光带着 x 能量打到表面，某一方向的出射光 能量就等于 x。这建立在上述两条假设成立的情况下。
Vector3f Material::eval(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
switch(m_type){
case DIFFUSE:
{
// calculate the contribution of diffuse model
float cosalpha = dotProduct(N, wo);
if (cosalpha > 0.0f) {
/*
Vector3f diffuse = Kd / M_PI;
return diffuse;
*/
// wi is the cam_to_pos, wo is pos_to_light, so we need to reverse the wi let them all outward for the
// compute of the cos
return brdfMicrofacet(-wi, wo, N, 0.5f, 0.2, Kd, 0.5);
}
else
return Vector3f(0.0f);
break;
}
}
}
3.4 Sample()
https://www.zhihu.com/question/522815454
对于水平和垂直方向的极坐标角度θ和φ分别采样（利用三角函数将其控制在合适的区间内），然后再转化为直角坐标。
由于这个半球面是以 Normal 为中心的，所以需要进行一个变换,代码中为 toWorld()。由于代码中并没有给出这个变换的来源与依据, 经过测试, wi是一个单位向量。比较好的方法是，既然已经得到了一个球面的采样，与一个normal，要以这个 normal 定义的半球面采样，只需要判断这个球面的向量与normal的夹角，如果为正那很好，若为负数，直接取逆即可。
Vector3f Material::sample(const Vector3f &wi, const Vector3f &N){
switch(m_type){
case DIFFUSE:
{
// uniform sample on the hemisphere
float x_1 = get_random_float(), x_2 = get_random_float();
float z = std::fabs(1.0f - 2.0f * x_1);
float r = std::sqrt(1.0f - z * z), phi = 2 * M_PI * x_2;
Vector3f localRay(r*std::cos(phi), r*std::sin(phi), z);
return toWorld(localRay, N);
break;
}
}
}
4 Microfacet Model
GGX分布函数，可以理解为这个表面和 h 一样方向的有多少 [0-1]
G_Smith阴影遮蔽函数主要是由两部分的乘积组成，一个是 cam2pos 和 light2pos 这两个都会产生遮挡。
我们常会遇到这两个GGX分布函数的表达式，实际上他们两个是一样的, \(\theta_m\) 为 half-vector 与 normal 的夹角。
// For Microfacet Model 2007 Walter
// Author:GSN Composer from youtube
// lihgt = in, view = out, N = Macro Normal
// metallic use to define the dielectric(glass air and plastic with 0.0f) metal(Iron copper and gold with 1.0f)
// roughness using at GGX-D and G_Smith function
// relectance use to modifine the F0(Fresnel relectance), eg. if2 F0 = 0.04 the material is Glass
// baseColor use to the BRDF diffuse part and also F0(Fresnel relectance) whth metal
inline Vector3f brdfMicrofacet(Vector3f Light, Vector3f View, Vector3f N, float metallic, float roughness, Vector3f baseColor, float reflectance);
inline Vector3f fresnelSchlick(float cos, Vector3f F0);
inline float D_GGX(float NOH, float roughness);
inline float G_Smith(float NoV, float NoL, float roughness);
inline float G1_GGX_Schlick(float NoV, float roughness);
Vector3f Material::eval(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
switch(m_type){
case DIFFUSE:
{
// calculate the contribution of diffuse model
float cosalpha = dotProduct(N, wo);
if (cosalpha > 0.0f) {
/*
Vector3f diffuse = Kd / M_PI;
return diffuse;
*/
// wi is the cam_to_pos, wo is pos_to_light, so we need to reverse the wi let them all outward for the
// compute of the cos
return brdfMicrofacet(-wi, wo, N, 0.5f, 0.2, Kd, 0.5);
}
else
return Vector3f(0.0f);
break;
}
}
}
Vector3f Material::brdfMicrofacet(Vector3f Light, Vector3f View, Vector3f N, float metallic, float roughness, Vector3f baseColor, float reflectance)
{
// half vector
Vector3f H = normalize(Light + View);
float NoL = clamp(0.001f, 1.f, dotProduct(N, Light));
float NoV = clamp(0.001f, 1.f, dotProduct(N, View));
float NoH = clamp(0.001f, 1.f, dotProduct(N, H));
float VoH = clamp(0.001f, 1.f, dotProduct(View, H));
Vector3f F0(0.16 * (reflectance * reflectance));
F0 = lerp(F0, baseColor, metallic);
Vector3f F = fresnelSchlick(VoH, F0);
float D = D_GGX(NoH, roughness);
float G = G_Smith(NoV, NoL, roughness);
Vector3f spec = (F * D * G) / (4.0f * NoV * NoL);
Vector3f rolD = baseColor;
// only the transmitted part contribute to diffuse for the energy conservation;
rolD = rolD *(Vector3f(1.0f) - F);
// the diffuse part for the metal matriels must be equal to 0.0f;
rolD = rolD * (1.0f - metallic);
Vector3f diff = rolD / M_PI;
return diff + spec;
}
inline Vector3f Material::fresnelSchlick(float cos, Vector3f F0)
{
return F0 + ( Vector3f(1.0f) - F0 ) * std::pow(1.0 - cos, 5.0);
}
inline float Material::D_GGX(float NoH, float roughness)
{
float alpha = roughness * roughness;
float alpha2 = alpha * alpha;
float NoH2 = NoH *NoH;
float b = NoH2 * (alpha2 -1.f) + 1.0f;
return alpha2 / (M_PI * b * b);
}
inline float Material::G_Smith(float NoV, float NoL, float roughness)
{
return G1_GGX_Schlick(NoV, roughness) * G1_GGX_Schlick(NoL, roughness);
}
inline float Material::G1_GGX_Schlick(float NoV, float roughness)
{
float alpha = roughness * roughness;
float k = alpha * 0.5f;
return NoV /(NoV * (1.0f - k ) + k);
}
https://zhuanlan.zhihu.com/p/20119162
https://zhuanlan.zhihu.com/p/25539396
https://mitsuba2.readthedocs.io/en/latest/generated/plugins.html#rough-conductor-material-roughconductor
Microfacet BRDF: Theory and Implementation of Basic PBR Materials [Shaders Monthly #9]: https://www.youtube.com/watch?v=gya7x9H3mV0
GGX 2007 :https://www.cs.cornell.edu/~srm/publications/EGSR07-btdf.pdf