• 【GAMES101】作业 7: 路径追踪

    文章目录

    作业 7: 路径追踪

    作业要求

    在之前的练习中，我们实现了 Whitted-Style Ray Tracing 算法，并且用 BVH 等加速结构对于求交过程进行了加速。在本次实验中，我们将在上一次实验的基 础上实现完整的 Path Tracing 算法。至此，我们已经来到了光线追踪版块的最后一节内容。

    你需要从上一次编程练习中直接拷贝以下函数到对应位置： -

    • Triangle::getIntersection in Triangle.hpp: 将你的光线-三角形相交函数 粘贴到此处，请直接将上次实验中实现的内容粘贴在此。
    • IntersectP(const Ray& ray, const Vector3f& invDir, const std::array& dirIsNeg) in the Bounds3.hpp: 这个函数的作用是判断包围盒 BoundingBox 与光线是否相交，请直接将上次实验中实现的内容粘贴在此处，并且注意检查 t_enter = t_exit 的时候的判断是否正确。
    • getIntersection(BVHBuildNode* node, const Ray ray) in BVH.cpp: BVH查找过程，请直接将上次实验中实现的内容粘贴在此处

    在本次实验中，你只需要修改这一个函数:

    • castRay(const Ray ray, int depth) in Scene.cpp: 在其中实现 Path Tracing 算法

    可能用到的函数有：

    • intersect(const Ray ray)in Scene.cpp: 求一条光线与场景的交点
    • sampleLight(Intersection pos, float pdf) in Scene.cpp: 在场景的所有 光源上按面积 uniform 地 sample 一个点，并计算该 sample 的概率密度
    • sample(const Vector3f wi, const Vector3f N) in Material.cpp: 按照该 材质的性质，给定入射方向与法向量，用某种分布采样一个出射方向
    • pdf(const Vector3f wi, const Vector3f wo, const Vector3f N) in Material.cpp: 给定一对入射、出射方向与法向量，计算 sample 方法得到该出射 方向的概率密度
    • eval(const Vector3f wi, const Vector3f wo, const Vector3f N) in Material.cpp: 给定一对入射、出射方向与法向量，计算这种情况下的 f_r 值

    可能用到的变量有：

    • RussianRoulette in Scene.cpp: P_RR, Russian Roulette 的概率

    代码迁移

    Triangle::getIntersection()

    inline Intersection Triangle::getIntersection(Ray ray)
{
    Intersection inter;

    if (dotProduct(ray.direction, normal) > 0)
        return inter;
    double u, v, t_tmp = 0;
    Vector3f pvec = crossProduct(ray.direction, e2);
    double det = dotProduct(e1, pvec);
    if (fabs(det) < EPSILON)
        return inter;

    double det_inv = 1. / det;
    Vector3f tvec = ray.origin - v0;
    u = dotProduct(tvec, pvec) * det_inv;
    if (u < 0 || u > 1)
        return inter;
    Vector3f qvec = crossProduct(tvec, e1);
    v = dotProduct(ray.direction, qvec) * det_inv;
    if (v < 0 || u + v > 1)
        return inter;
    t_tmp = dotProduct(e2, qvec) * det_inv;

    if (t_tmp < 0) return inter;

    // TODO find ray triangle intersection
    inter.happened = true;
    inter.distance = t_tmp;
    inter.normal = normal;
    inter.obj = this;
    inter.m = m;
    inter.coords = ray(t_tmp);

    return inter;
}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35

    IntersectP()

    inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
    const std::array<int, 3>& dirIsNeg) const
{
    // invDir: ray direction(x,y,z), invDir=(1.0/x,1.0/y,1.0/z), use this because Multiply is faster that Division
    // dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x>0),int(y>0),int(z>0)], use this to simplify your logic
    // TODO test if ray bound intersects

    float tmin_x = ((dirIsNeg[0] ? pMin.x : pMax.x) - ray.origin.x) * invDir.x;
    float tmax_x = ((dirIsNeg[0] ? pMax.x : pMin.x) - ray.origin.x) * invDir.x;
    float tmin_y = ((dirIsNeg[1] ? pMin.y : pMax.y) - ray.origin.y) * invDir.y;
    float tmax_y = ((dirIsNeg[1] ? pMax.y : pMin.y) - ray.origin.y) * invDir.y;
    float tmin_z = ((dirIsNeg[2] ? pMin.z : pMax.z) - ray.origin.z) * invDir.z;
    float tmax_z = ((dirIsNeg[2] ? pMax.z : pMin.z) - ray.origin.z) * invDir.z;

    float t_enter = std::max(tmin_x, std::max(tmin_y, tmin_z));
    float t_exit = std::min(tmax_x, std::min(tmax_y, tmax_z));

    if (t_enter <= t_exit && t_exit >= 0) return true;
    return false;
}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20

    逐一这里最后的条件判断 t_enter <= t_exit 课上说不考虑边界情况，但是 cornell box 模型的墙壁是一个平面，所以包围盒没有宽度，如果不加等号，将判断为不相交，导致结果一片黑。所以应当加上等号。

    getIntersection()

    Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    // TODO Traverse the BVH to find intersection
    Intersection inter, hit1, hit2;
    if (!node->bounds.IntersectP(
        ray, Vector3f(1.0 / ray.direction.x, 1.0 / ray.direction.y, 1.0 / ray.direction.z),
        std::array<int, 3>{int(ray.direction.x > 0), int(ray.direction.y > 0), int(ray.direction.z > 0)}))
        return inter;
    if (!node->left && !node->right)
        return node->object->getIntersection(ray);
    hit1 = getIntersection(node->left, ray);
    hit2 = getIntersection(node->right, ray);
    return hit1.distance < hit2.distance ? hit1 : hit2;
}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14

    路径追踪 castRay()

    根据下面算法实现：

    shade(p, wo)
    # Contribution from the light source.
    L_dir = 0.0
	Uniformly sample the light at x’ (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 θ * cos θ’ / |x’ - p|^2 / pdf_light	
    
    # Contribution from other reflectors.
    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 θ / pdf_hemi / P_RR
    Return L_dir + L_indir

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16

    首先判断光线是否投射到物体上，如果不是直接返回，如果投射到光源上也直接返回光源的亮度

    	Intersection inter = intersect(ray);
	if (!inter.happened) return Vector3f(0.0f, 0.0f, 0.0f);
	if (inter.m->hasEmission()) return inter.m->getEmission();
...

    • 1
    • 2
    • 3
    • 4

    随后计算光直接照射的辐射亮度，注意这里的定义： w i w_i wi​ 方向从 p \mathbf p p 到 x ′ \mathbf x\prime x′ ， w o w_o wo​ 方向为光线方向，与课上定义相反。

    ...	
	// 计算直接光照
    Vector3f L_dir;
    Intersection lightInter;
    float pdf_light;
    sampleLight(lightInter, pdf_light);
	if (pdf_light < EPSILON) pdf_light = EPSILON;	// 防止白点
    Vector3f p = inter.coords;
    Vector3f n = inter.normal;
    Vector3f xp = lightInter.coords;
    Vector3f np = lightInter.normal;
    Vector3f p2xp = xp - p;
    float distance = p2xp.norm();
    float distance_square = dotProduct(p2xp, p2xp);
    Vector3f wi = normalize(p2xp);   // 与课程上方向定义相反
    Vector3f wo = ray.direction;   // 与课程上方向定义相反
    float cos_theta = dotProduct(n, wi);
    Intersection p2light = intersect(Ray(p, wi));
    //
    // 下面如果用EPSILON = 0.00001会出现黑线
    //
    if (p2light.happened && distance - p2light.distance < 0.0001) {
        float cos_thetap = dotProduct(np, -wi);
        Vector3f f_r = inter.m->eval(wo, wi, n);
        Vector3f L_i = lightInter.emit;
        L_dir = L_i * f_r * cos_theta * cos_thetap / distance_square / pdf_light;
    }
...

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28

    这里当 p d f \mathrm {pdf} pdf 过小时，辐射亮度会非常大，所以会出现白点，所以我们要限制 p d f \mathrm {pdf} pdf 大于 ε \varepsilon ε ，下面也同理。

    注意这里如果直接判断 distance 与 p2light.distance 相等再计算 L_dir 会导致全暗，因为浮点数经度问题相等时很难的，所以应当设置一定的 ε \varepsilon ε。这里设置 ε = 0.00001 \varepsilon = 0.00001 ε=0.00001 会出现黑线，增大才可消除。

    最后计算其他的反射光：

    ...
	// 计算其他反射光线
    if (get_random_float() > RussianRoulette) return L_dir;
    Vector3f L_indir;
    wi = normalize(inter.m->sample(wo, n));
    Ray newRay = Ray(p, wi);
    Intersection newInter = intersect(newRay);
    if (newInter.happened && !newInter.m->hasEmission()) {
        Vector3f f_r = inter.m->eval(wo, wi, n);
        cos_theta = dotProduct(n, wi);
        float pdf_hemi = inter.m->pdf(wo, wi, n);
        if (pdf_hemi < EPSILON) pdf_hemi = EPSILON;	// 防止白点
        L_indir = castRay(newRay, depth + 1) * f_r * cos_theta / pdf_hemi / RussianRoulette;
    }
    
    return L_dir + L_indir;

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16

    完整代码：

    Vector3f Scene::castRay(const Ray& ray, int depth) const
{
    Intersection inter = intersect(ray);
    if (!inter.happened) return Vector3f(0.0f, 0.0f, 0.0f);
    if (inter.m->hasEmission()) return inter.m->getEmission();
    

    // 计算直接光照
    Vector3f L_dir;
    Intersection lightInter;
    float pdf_light;
    sampleLight(lightInter, pdf_light);
    if (pdf_light < EPSILON) pdf_light = EPSILON;
    Vector3f p = inter.coords;
    Vector3f n = inter.normal;
    Vector3f xp = lightInter.coords;
    Vector3f np = lightInter.normal;
    Vector3f p2xp = xp - p;
    float distance = p2xp.norm();
    float distance_square = dotProduct(p2xp, p2xp);
    Vector3f wi = normalize(p2xp);   // 与课程上方向定义相反
    Vector3f wo = ray.direction;   // 与课程上方向定义相反
    float cos_theta = dotProduct(n, wi);
    Intersection p2light = intersect(Ray(p, wi));
    //
    // 下面如果用EPSILON = 0.00001会出现黑线
    //
    if (p2light.happened && distance - p2light.distance < 0.0001) {
        float cos_thetap = dotProduct(np, -wi);
        Vector3f f_r = inter.m->eval(wo, wi, n);
        Vector3f L_i = lightInter.emit;
        L_dir = L_i * f_r * cos_theta * cos_thetap / distance_square / pdf_light;
    }

    // 计算其他反射光线
    if (get_random_float() > RussianRoulette) return L_dir;
    Vector3f L_indir;
    wi = normalize(inter.m->sample(wo, n));
    Ray newRay = Ray(p, wi);
    Intersection newInter = intersect(newRay);
    if (newInter.happened && !newInter.m->hasEmission()) {
        Vector3f f_r = inter.m->eval(wo, wi, n);
        cos_theta = dotProduct(n, wi);
        float pdf_hemi = inter.m->pdf(wo, wi, n);
        if (pdf_hemi < EPSILON) pdf_hemi = EPSILON;
        L_indir = castRay(newRay, depth + 1) * f_r * cos_theta / pdf_hemi / RussianRoulette;
    }
    
    return L_dir + L_indir;

}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48
    • 49
    • 50
    • 51

    这里的 get_random_float() 框架代码有些问题，更改为

    static std::random_device dev;
static std::mt19937 rng(dev());
static std::uniform_real_distribution<float> dist(0.f, 1.f); // distribution in range [1, 6]

inline float get_random_float()
{
 return dist(rng);
}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8

    多线程

    我们知道渲染的过程是逐一遍历屏幕像素，然后从视点向像素中心投射 spp 条光线，大致流程如下：

    for (uint32_t j = 0; j < scene.height; ++j)
    for (uint32_t i = 0; i < scene.width; ++i) 
        for (int k = 0; k < spp; k++){
            framebuffer[m] += castRay(eye_pos) / spp;
            m++;
        }

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6

    我们可以用多线程的方法，同时计算一些像素。这里我们按行来划分线程，每一线程处理一定的行数。

    首先定义处理一定行的渲染函数，这里用 lambda 表达式来写：

    auto renderInRows = [&](uint32_t beginRow, uint32_t endRow){
    for (uint32_t j = beginRow; j < endRow; ++j){
        for (uint32_t i = 0; i < scene.width; ++i){
			float x = (2 * (i + 0.5) / (float)scene.width - 1) *
                      imageAspectRatio * scale;
            float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;

            Vector3f dir = normalize(Vector3f(-x, y, 1));
            for (int k = 0; k < spp; k++){
                framebuffer[m] += scene.castRay(Ray(eye_pos, dir), 0) / spp;  
            }
            m++;
        }
        mtx.lock();
        process++; // 用来更新进度条，表示已经处理的行数
        UpdateProgress((float)process / scene.height);
        mtx.unlock();
	}
}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19

    定义线程

    const int numOfThreads = 16;
std::thread th[numOfThreads];
std::mutex mtx;						// 防止冲突
std::atomic_int process = 0;        // 用来更新进度条，表示已经处理的行数

    • 1
    • 2
    • 3
    • 4

    多线程运行

    for (int i = 0; i < numOfThreads; ++i)
    th[i] = std::thread(renderInRows, i * scene.height / numOfThreads, (i + 1) * scene.height / numOfThreads);
for (int i = 0; i < numOfThreads; ++i)
    th[i].join();

    • 1
    • 2
    • 3
    • 4

    完整代码：

    void Renderer::Render(const Scene& scene)
{
    std::vector<Vector3f> framebuffer(scene.width * scene.height);

    float scale = tan(deg2rad(scene.fov * 0.5));
    float imageAspectRatio = scene.width / (float)scene.height;
    Vector3f eye_pos(278, 273, -800);
    int m = 0;

    // change the spp value to change sample ammount
    int spp = 16;
    std::cout << "SPP: " << spp << "\n";
    bool isMultithreading = true;
    if (isMultithreading) {
        // 定义线程
        const int numOfThreads = 16;
        std::thread th[numOfThreads];
        std::atomic_int process = 0;        // 用来更新进度条，表示已经处理的行数
        //
        // lambda表达式定义函数，渲染一定行
        //
        auto renderInRows = [&](uint32_t beginRow, uint32_t endRow) {
            for (uint32_t j = beginRow; j < endRow; ++j) {
                for (uint32_t i = 0; i < scene.width; ++i) {
                    float x = (2 * (i + 0.5) / (float)scene.width - 1) *
                        imageAspectRatio * scale;
                    float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;

                    Vector3f dir = normalize(Vector3f(-x, y, 1));
                    for (int k = 0; k < spp; k++) {
                        framebuffer[j * scene.width + i] += scene.castRay(Ray(eye_pos, dir), 0) / spp;
                    }
                }
                process++;
                UpdateProgress((float)process / scene.height);
            }
        };

        for (int i = 0; i < numOfThreads; ++i)
            th[i] = std::thread(renderInRows, i * scene.height / numOfThreads, (i + 1) * scene.height / numOfThreads);
        for (int i = 0; i < numOfThreads; ++i)
            th[i].join();
        UpdateProgress(1.f);
    }
    else {
        for (uint32_t j = 0; j < scene.height; ++j) {
            for (uint32_t i = 0; i < scene.width; ++i) {
                float x = (2 * (i + 0.5) / (float)scene.width - 1) *
                    imageAspectRatio * scale;
                float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;

                Vector3f dir = normalize(Vector3f(-x, y, 1));
                for (int k = 0; k < spp; k++) {
                    framebuffer[m] += scene.castRay(Ray(eye_pos, dir), 0) / spp;
                }
                m++;
            }
            UpdateProgress(j / (float)scene.height);
        }
        UpdateProgress(1.f);
    }

    // save framebuffer to file
    FILE* fp = fopen("binary.ppm", "wb");
    (void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);
    for (auto i = 0; i < scene.height * scene.width; ++i) {
        static unsigned char color[3];
        color[0] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].x), 0.6f));
        color[1] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].y), 0.6f));
        color[2] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].z), 0.6f));
        fwrite(color, 1, 3, fp);
    }
    fclose(fp);    
}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48
    • 49
    • 50
    • 51
    • 52
    • 53
    • 54
    • 55
    • 56
    • 57
    • 58
    • 59
    • 60
    • 61
    • 62
    • 63
    • 64
    • 65
    • 66
    • 67
    • 68
    • 69
    • 70
    • 71
    • 72
    • 73
    • 74

    实现 M i c r o f a c e t \mathrm{Microfacet} Microfacet 材质

    对于粗糙表面：

    • 宏观尺度：平坦又粗糙
    • 微观尺度：颠簸又光滑

    表面的每个元素像一个镜子：

    • 称为微表面
    • 每个微表面都有自己的法线

    微表面的法线分布决定粗糙程度：

    分布集中就是又光泽的，分布分散就是漫反射的。

    微表面模型下，镜面反射的 BRDF 通过下式计算：
    f ( i , o ) = F ( i , h ) G ( i , o , h ) D ( h ) 4 ( n , i ) ( n , o ) f(\mathbf{i}, \mathbf{o})=\frac{\mathbf{F}(\mathbf{i}, \mathbf{h}) \mathbf{G}(\mathbf{i}, \mathbf{o}, \mathbf{h}) \mathbf{D}(\mathbf{h})}{4(\mathbf{n}, \mathbf{i})(\mathbf{n}, \mathbf{o})} f(i,o)=4(n,i)(n,o)F(i,h)G(i,o,h)D(h)​
    其中 F ( i , h ) \mathbf{F}(\mathbf{i}, \mathbf{h}) F(i,h) 为菲涅尔项； G ( i , o , h ) \mathbf{G}(\mathbf{i}, \mathbf{o}, \mathbf{h}) G(i,o,h) 为阴影遮蔽的项，解决微表面间互相遮挡的问题；$ \mathbf{D}(\mathbf{h})$ 为法线分布。

    菲涅尔项 Fresnel Reflection / Term

    菲涅尔项可以解释有多少能量发生反射，有多少能量发生折射
    R s = ∣ n 1 cos ⁡ θ i − n 2 cos ⁡ θ t n 1 cos ⁡ θ i + n 2 cos ⁡ θ t ∣ 2 = ∣ n 1 cos ⁡ θ i − n 2 1 − ( n 1 n 2 sin ⁡ θ i ) 2 n 1 cos ⁡ θ i + n 2 1 − ( n 1 n 2 sin ⁡ θ i ) 2 ∣ 2 , R p = ∣ n 1 cos ⁡ θ t − n 2 cos ⁡ θ i n 1 cos ⁡ θ t + n 2 cos ⁡ θ i ∣ 2 = ∣ n 1 1 − ( n 1 n 2 sin ⁡ θ i ) 2 − n 2 cos ⁡ θ i n 1 1 − ( n 1 n 2 sin ⁡ θ i ) 2 + n 2 cos ⁡ θ i ∣ . R e f f = 1 2 ( R s + R p )

    Rs=|n1cos⁡θi−n2cos⁡θtn1cos⁡θi+n2cos⁡θt|2=|n1cos⁡θi−n21−(n1n2sin⁡θi)2n1cos⁡θi+n21−(n1n2sin⁡θi)2|2,Rp=|n1cos⁡θt−n2cos⁡θin1cos⁡θt+n2cos⁡θi|2=|n11−(n1n2sin⁡θi)2−n2cos⁡θin11−(n1n2sin⁡θi)2+n2cos⁡θi|." role="presentation" style="position: relative;">Rs=|n1cos⁡θi−n2cos⁡θtn1cos⁡θi+n2cos⁡θt|2=|n1cos⁡θi−n21−(n1n2sin⁡θi)2n1cos⁡θi+n21−(n1n2sin⁡θi)2|2,Rp=|n1cos⁡θt−n2cos⁡θin1cos⁡θt+n2cos⁡θi|2=|n11−(n1n2sin⁡θi)2−n2cos⁡θin11−(n1n2sin⁡θi)2+n2cos⁡θi|.
    \\R_{\mathrm{eff}} =\frac{1}{2} (R_s+R_p) Rs​=∣ ∣​n1​cosθi​+n2​cosθt​n1​cosθi​−n2​cosθt​​∣ ∣​2=∣ ∣​n1​cosθi​+n2​1−(n2​n1​​sinθi​)2 ​n1​cosθi​−n2​1−(n2​n1​​sinθi​)2 ​​∣ ∣​2,Rp​=∣ ∣​n1​cosθt​+n2​cosθi​n1​cosθt​−n2​cosθi​​∣ ∣​2=∣ ∣​n1​1−(n2​n1​​sinθi​)2 ​+n2​cosθi​n1​1−(n2​n1​​sinθi​)2 ​−n2​cosθi​​∣ ∣​.​Reff​=21​(Rs​+Rp​)
    上面的算法过于复杂，可以用 Schlick’s approximation 近似：
    R ( θ ) = R 0 + ( 1 − R 0 ) ( 1 − cos ⁡ θ ) 5 R 0 = ( n 1 − n 2 n 1 + n 2 ) 2
    R(θ)=R0+(1−R0)(1−cos⁡θ)5R0=(n1−n2n1+n2)2" role="presentation" style="position: relative;">R(θ)=R0+(1−R0)(1−cos⁡θ)5R0=(n1−n2n1+n2)2
    R(θ)R0​​=R0​+(1−R0​)(1−cosθ)5=(n1​+n2​n1​−n2​​)2​
    代码框架没有进行近似，也直接给出了求解菲涅尔项的函数

    void fresnel(const Vector3f &I, const Vector3f &N, const float &ior, float &kr) const
    {
        float cosi = clamp(-1, 1, dotProduct(I, N));
        float etai = 1, etat = ior;
        if (cosi > 0) {  std::swap(etai, etat); }
        // Compute sini using Snell's law
        float sint = etai / etat * sqrtf(std::max(0.f, 1 - cosi * cosi));
        // Total internal reflection
        if (sint >= 1) {
            kr = 1;
        }
        else {
            float cost = sqrtf(std::max(0.f, 1 - sint * sint));
            cosi = fabsf(cosi);
            float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost));
            float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost));
            kr = (Rs * Rs + Rp * Rp) / 2;
        }
        // As a consequence of the conservation of energy, transmittance is given by:
        // kt = 1 - kr;
    }

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21

    这里也可以用 Schlick’s approximation 近似

        float Fresnel_Reflection(const Vector3f& wi, const Vector3f& N, float n1, float n2) {
        float R0 = (n1 - n2) / (n1 + n2);
        R0 *= R0;
        float cos_theta = dotProduct(-wi, N);
        float R = R0 + (1 - R0) * std::pow(1 - cos_theta, 5);
        return R;
    }

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7

    实际上差别不大。

    法线分布 Normal Distribution Function (NDF)

    微表面法线分布定义为
    ∫ Ω D ( h ) cos ⁡ θ h d h = 1 \int_\Omega D(h)\cos \theta_h\mathrm d h = 1 ∫Ω​D(h)cosθh​dh=1

    Beckmann NDF

    与高斯函数相似
    D ( h ) = e − tan ⁡ 2 θ h α 2 π α 2 cos ⁡ 4 θ h D(h)=\frac{e^{-\frac{\tan ^{2} \theta_{h}}{\alpha^{2}}}}{\pi \alpha^{2} \cos ^{4} \theta_{h}} D(h)=πα2cos4θh​e−α2tan2θh​​​
    α \alpha α ：表面的粗糙程度（越小越像镜面）关于表面粗糙度的取值，参考文章中是这么说的：“在迪士尼原理着色模型Disney principled shading model中，推荐将粗糙度控制以 α = r 2 \alpha = r^2 α=r2 暴露给用户，其中 r r r 是 0 到 1 的粗糙度参数值，表示 r o u g h n e s s roughness roughness。

    θ h \theta_h θh​：法线 n n n 与向量 $h $ 的夹角， h h h 为半程向量，即入射出射角平分线

    Beckmann NDF 与 θ \theta θ 有关，与 ϕ \phi ϕ 无关，所以它只能表示各项同性的表面材质。

    Beckmann NDF 是定义在坡度空间 slope space 上的，这就解释了分子为什么为 tan ⁡ 2 θ \tan^2\theta tan2θ

    在这里插入图片描述

    这样定义的好处是在无限大的平面上定义高斯函数，并且可以保证不出现面朝下的微表面。

        float Beckmann_NDF(const Vector3f& wi, const Vector3f& wo, const Vector3f& N, float roughness) {
        Vector3f h = (-wi + wo).normalized();
        float cos_theta_h = dotProduct(h, N);
        float cos_theta_h_square = cos_theta_h * cos_theta_h;
        float alpha = roughness * roughness;
        float alpha_square = alpha * alpha;
        float tan2cos_theta_h = (1 - cos_theta_h_square) / cos_theta_h_square;
        float denom = M_PI * alpha_square * cos_theta_h_square * cos_theta_h_square;
        if (denom < EPSILON) denom = EPSILON;
        float D = std::exp(-tan2cos_theta_h / alpha_square) / denom;
        return D;
    }


    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13

    GGX NDF（TR）

    D ( h ) = α 2 π ( ( n ⋅ h ) 2 ( α 2 − 1 ) + 1 ) 2 D(h) = \frac{\alpha^2}{\pi((\mathbf{n\cdot h})^2(\alpha^2 -1)+1)^2} D(h)=π((n⋅h)2(α2−1)+1)2α2​

    一个特点：长尾 long tail

    光晕现象

       float GGX_NDF(const Vector3f& wi, const Vector3f& wo, const Vector3f& N, float roughness) {
        Vector3f h = (-wi + wo).normalized();
        float alpha = roughness * roughness;
        float alpha_square = alpha * alpha;
        float nDoth = dotProduct(N, h);
        float denom = nDoth * nDoth * (alpha_square - 1) + 1;
        denom = denom * denom * M_PI;
        if (denom < EPSILON) denom = EPSILON;
        float D = alpha_square / denom;
        return D;
    }


    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12

    Extending GGX——GTR(Generalized Trowbridge-Reitz)

    有更长的尾巴

    γ = 2 \gamma =2 γ=2 为 GGX 模型， γ \gamma γ 越大，就越接近 Beckmann 模型

    Shadowing-Masking Term

    阴影遮蔽的项，解决微表面间互相遮挡的问题

    • 从光线出发，出现阴影（左图）叫做 Shadowing
    • 从眼睛出发，看不到一些微表面（右图）叫做 Masking

    从上往下看应该不会有遮挡现象，而当以掠视角 grazing angles 观察，遮挡现象将会很明显

    Smith Shadowing-Masking Term

    将 shadowing 与 masking 分离
    G ( i , o , m ) ≈ G 1 ( i , m ) G 1 ( o , m ) G ( i , m ) = i ⋅ m ( i ⋅ m ) ( 1 − k ) + k G(\mathbf{i}, \mathbf{o}, \mathbf{m}) \approx G_{1}(\mathbf{i}, \mathbf{m}) G_{1}(\mathbf{o}, \mathbf{m})\\G(\mathbf{i}, \mathbf{m})=\frac{\mathbf{i}\cdot \mathbf{m}}{(\mathbf{i}\cdot \mathbf{m})(1-k)+k} G(i,o,m)≈G1​(i,m)G1​(o,m)G(i,m)=(i⋅m)(1−k)+ki⋅m​
    其中 k = ( r o u g h n e s s + 1 ) 2 / 8 k=(roughness + 1)^2/8 k=(roughness+1)2/8
    在这里插入图片描述

    其中 m \mathbf m m 是微表面法线。

        float Smith_SMT(const Vector3f& wi, const Vector3f& wo, const Vector3f& N, float roughness) {
        float k = (roughness + 1) * (roughness + 1) / 8;
        float woDotN = dotProduct(wo, N);
        float wiDotN = dotProduct(-wi, N);
        float denomm = wiDotN * (1 - k) + k;
        if (denomm < EPSILON) denomm = EPSILON;
        float G1 = wiDotN / denomm;
        denomm = woDotN * (1 - k) + k;
        if (denomm < EPSILON) denomm = EPSILON;
        float G2 = woDotN / denomm;
        return G1 * G2; 
    }


    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13

    Cook-Torrance BRDF

    从微观的尺度去理解了物体表面镜面反射与漫反射的现象，分开考虑镜面反射与漫反射现象，得到
    f r = k d f lambert  + k s f cook-torrance  f_{r}=k_{d} f_{\text {lambert }}+k_{s} f_{\text {cook-torrance }} fr​=kd​flambert ​+ks​fcook-torrance ​
    其中 k d , k s k_d,k_s kd​,ks​分别指入射光线中被折射部分的能量所占的比率与被反射部分的比率(一般由菲涅尔项决定， k s k_s ks​ 就是菲涅尔项的值， k d k_d kd​ 简单设为 1 − k s 1-k_s 1−ks​)，而 k f l a m b e r t k_{flambert} kflambert​ 指漫反射的BRDF， f c o o k − t o r r a n c e f_{cook−torrance} fcook−torrance​ 指镜面反射的BRDF。

    其中
    f l a m b e r t = ρ π ,    f c o o k − t o r r a n c e = F ( i , h ) G ( i , o , h ) D ( h ) 4 ( n , i ) ( n , o ) f_{lambert}=\frac{\rho}{\pi},\ \ f_{cook-torrance}=\frac{\mathbf{F}(\mathbf{i}, \mathbf{h}) \mathbf{G}(\mathbf{i}, \mathbf{o}, \mathbf{h}) \mathbf{D}(\mathbf{h})}{4(\mathbf{n}, \mathbf{i})(\mathbf{n}, \mathbf{o})} flambert​=πρ​,  fcook−torrance​=4(n,i)(n,o)F(i,h)G(i,o,h)D(h)​

    这个 cook-torrance 相当于一种补偿，补偿多次弹射损失的能量，把 diffuse 补进来了。但是按照闫老师说这么补是不对的（不严谨的），应该使用 Kulla-Conty 考虑光线多次弹射近似这种结果是正确的做法。

    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;
            }
            else
                return Vector3f(0.0f);
            break;
        }
        case MICROFACET:
        {
            // 微表面模型
            float cosalpha = dotProduct(N, wo);
            if (cosalpha <= 0.0f)
                return Vector3f(0.0f);
            //
            // 只处理表面朝上的微表面
            // 

            float roughness = 0.0f;

            // 法线分布
            float D = Beckmann_NDF(wi, wo, N, 0.8);
            //float D = GGX_NDF(wi, wo, N, roughness);

            // Shadowing-Masking Term
            float G = Smith_SMT(wi, wo, N, roughness);

            // 菲涅尔项
            //float F = Fresnel_Reflection(wi, N, 1.0f, 1.85f);

            float F;
            float etat = 1.85;
            fresnel(wi, N, etat, F);

            // 计算 BRDF
            float numer = D * G * F;
            float denom = 4 * dotProduct(wo, N) * dotProduct(-wi, N);
            if (denom < EPSILON) denom = EPSILON;

            Vector3f fl = numer / denom;
            Vector3f fc = Kd / M_PI;

            float ks = F;
            float kd = 1.0f - ks;

            return ks * fl + kd * fc;

        }
    }
}


    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48
    • 49
    • 50
    • 51
    • 52
    • 53
    • 54
    • 55
    • 56
    • 57

    同时，Material::getColorAt(),Material::pdf()只需复制 DIFFUSE 处。

    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;
        }
        case MICROFACET:
        {
            // 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;
        }
    }
}


    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    float Material::pdf(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
    switch(m_type){
        case DIFFUSE:
        {
            // uniform sample probability 1 / (2 * PI)
            if (dotProduct(wo, N) > 0.0f)
                return 0.5f / M_PI;
            else
                return 0.0f;
            break;
        }
        case MICROFACET:
        {
            // uniform sample probability 1 / (2 * PI)
            if (dotProduct(wo, N) > 0.0f)
                return 0.5f / M_PI;
            else
                return 0.0f;
            break;
        }
    }
}


    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23

    效果

    2048 SPP Cornell Box (DIFFUSE)

    多线程一小时，分辨率 784x784

    1024 SPP Cornell Box + Bunny (Microfacet)

    多线程，分辨率：1568x1568

    在这里插入图片描述

    请添加图片描述

  • 相关阅读:
    【优雅】Java List转二维Map
    6.eureka服务发现实例（springcloud）
    java不同类加载器重复加载一个类
    攻防世界数据逆向 2023
    重读百度移动生态：“第一曲线”的创新“延长线”
    Python函数装饰器的深入解析
    密码学 | RC4算法Native层分析
    Git实战技巧-如何同时撤回远程和本地分支合并操作
    一个例子形象的理解协程和线程的区别
    【OS】进程通信
  • 原文地址：https://blog.csdn.net/m0_57714486/article/details/126291616