Games101 基于蒙特卡洛积分的光线路径追踪 作业7 框架解读

目录

1 前言

这次的光线追踪，主要是基于 Radiometry 的一种实现，也就是基于物理的一种实现。本文对辐射度量学做了解释，同时给出了程序中的关键代码以及参考资料，实现了微表面模型，模型的代码正确无误，但是运用到路径追踪上会出现很强的噪声，以目前的实力无法解决。
image

image

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

image

我们可以采用量子化的方法，把一份能量看作一个小球。而 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
    image

一个物体某个点发出的黄色的小球会朝着空间中四散而去, 在一个锥形的固体角中，一秒中的黄色小球数量。

需要注意的是，上述公式中用到了\(d\), 我们知道这是一个极限的过程，立体角的定义是球面上某处面积与整个球面面积之比，当处于极限过程时，一块面积趋近于一个非常非常小的点，而这时立体角就可以看作是一个箭头。有了这个思想，就接下来就好做多了。

光线可以理解为一个糖葫芦。穿着power。

  • Irradiance
    Definition: The irradiance is the power per unit area incident on a surface point.

\[E(x)= \frac {d\Phi(x)}{dA} \]

\[[\frac {W}{m^ {2}}] [ \frac {lm}{m^ {2}} =lux] \]

可以理解为，四面八方来的光线这个面积上，可以引申为这个点上，由于极限的关系。

  • Radiance
    Definition: The radiance is the power emitted, reflected,transmitted or received by a surface, per unit solid angle, per projected unit area.

\[L(p, \omega)= \frac {d^ {2}\Phi (p,\omega )}{d\omega dA\cos \theta } \]

Radiance是光线Ray传播时(即不与物体交互时)的一个性质。他在物体发出光线，或者光线打到物体上时会有折损，这种折损与物体表面的法线方向与传播光线的方向有关** Lambert’s Cosine Law**。所以我们用Radiance表示光线的传播。

一个光线带有100power，打到物体上由于角度关系我们只得到了50power, 但是我们不能说光线是50power的。

\[L(p, \omega )= \frac {dE(p)}{d\omega \cos \theta } \]

\[L(p, \omega ) \cos \theta d\omega = dE(p) \]

接收到的 Irradiance 由于在接收时，会因为 Lambert’s Cosine Law而进行缩小， 所以在还原时，需要除以 \(\cos \theta\)进行放大，还原。

3 Coding

Pseudo Code

shade(p, wo)
	Uniformly sample the light at xx (pdf_light = 1 / A)
	Shoot a ray from p to x
	If the ray is not blocked in the middle
		L_dir = L_i * f_r * cos_theta * cos_theta_x / |x-p|^2 / pdf_light
	L_indir = 0.0
	Test Russian Roulette with probability P_RR
	Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
	Trace a ray r(p, wi)
	If ray r hit a non -emitting object at q
		L_indir = shade(q, wi) * f_r * cos_theta / pdf_hemi / P_RR
	Return L_dir + L_indir

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。总计两条假设

  1. 入射的radiance为常数。Suppose the incident lighting is uniform

  2. 各个方向的出射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 的夹角。

\[D_{GGX}(\mathbf{m}) = \frac{\alpha^2}{\pi((\mathbf{n} \cdot \mathbf{m})^2 (\alpha^2 - 1) + 1)^2} \\ \]

\[D_{GGX}(\mathbf{m}) = \frac{\alpha^2}{\pi\cos^4\theta_m (\alpha^2 +\tan ^2\theta_m)^2} \\ \]

\[\mathbf{n} \cdot \mathbf{m} = \cos \theta_m \]

\[\tan \theta_m = \frac{\sin \theta_m}{ \cos \theta_m} \]

    // 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

posted @ 2023-07-17 22:04  Dba_sys  阅读(82)  评论(0)  编辑  收藏  举报