[Nori]作业5

Assignment 5

Part 1

第一部分，我们需要将microfacet.cpp里的微表面模型BRDF补充完整，用来模拟类似塑料的材质。它由漫反射BRDF和粗糙电介质BRDF线性组合而成。完整BRDF公式是：

fr(ωi,ωo)=kdπ+ksD(ωh) F((ωh⋅ωi),ηe,ηi) G(ωi,ωo,ωh)4cos⁡θicos⁡θocos⁡θh,  ωh=(ωi+ωo)∣∣ωi+ωo∣∣2f_r({\omega_i},{\omega_o}) = \frac{k_d}{\pi} + {k_s} \frac{D({\omega_{h}})~ F\left({({\omega_h} \cdot {\omega_i})}, \eta_{e},\eta_{i}\right)~ G({\omega_i},{\omega_o},{\omega_{h}})}{4 \cos{\theta_i} \cos{\theta_o}\cos\theta_h}, ~~ {\omega_{h}} = \frac{\left({\omega_i} + {\omega_o}\right)}{\left|\left|{\omega_i} + {\omega_o}\right|\right|_2}

其实就是Lambert漫反射和Cook-Torrance高光结合在一起。kd∈[0,1]3k_d \in [0,1]^3是RGB漫反射系数，为了保证能量守恒，这里简单将ks定义为：ks=1−max⁡(kd)k_s = 1 - \max(k_d)

F就是菲涅尔项，作业4介绍过了，直接拿过来用就行

nori选择的法线分布项和几何遮蔽项是Beckmann，对应公式是：

G(ωi,ωo,ωh)=G1(ωi,ωh) G1(ωo,ωh)G({\omega_i},{\omega_o},{\omega_{h}}) = G_1({\omega_i},{\omega_{h}})~G_1({\omega_o},{\omega_{h}})

G1(ωv,ωh)=χ+(ωv⋅ωhωv⋅n){3.535b+2.181b21+2.276b+2.577b2,b<1.6,1,b≥1.6G_1({\omega_v},{\omega_h}) =\chi^+\left(\frac{\omega_v\cdot\omega_h}{\omega_v\cdot n}\right) \begin{cases} \frac{3.535b+2.181b^2}{1+2.276b+2.577b^2}, & b \lt 1.6, \\ 1, & b \ge 1.6 \end{cases}

b=(αtan⁡θv)−1,  χ+(c)={1,c>0,0,c≤0,b = (\alpha \tan{\theta_v})^{-1}, ~~ \chi^+(c) = \begin{cases} 1, & c > 0, \\ 0, & c \le 0, \end{cases} \\

D(ωh)=e−tan2θhα2πα2cos4θhD(\omega_h) = \frac{e^\frac{-{tan}^2\theta_h}{\alpha^2}}{\pi\alpha^2{cos}^4\theta_h}

其中，G1G1的θv\theta_v是表面法线nn和ωv\omega_v之间的夹角。

知道了BRDF计算公式后，还需要推导出采样公式，这部分中我们需要使用以下概率密度函数生成样本（这个就是pdf公式）：

ks D(ωh) Jh+(1−ks)cos⁡θoπk_s ~ D(\omega_h) ~ J_h + (1-k_s) \frac{\cos{\theta_o}}{\pi}

其中，Jh=(4(ωh⋅ωo))−1J_h = (4 (\omega_h \cdot \omega_o))^{-1}是半程向量映射的雅可比矩阵（wtf这是啥），采样结果就是出射方向ωo\omega_o，具体采样步骤：

  • 通过比较均匀分布的随机变量ξ1\xi_1和ksk_s来决定是漫反射还是高光反射
  • 变换并偏移随机变量ξ1\xi_1以便它可以在以后的采样步骤中重复使用（参考DiscretePDF::sampleReuse）（意思就是直接CV233）
  • 如果是漫反射，按照 src/diffuse.cpp 中的方法在球体上生成余弦加权方向
  • 如果是高光反射
    • 使用作业3中实现的Warp::squareToBeckmann采样一个法线
    • 使用此法线反射入射方向以生成出射方向

现在，所有拼图都有了，我们可以将microfacet.cpp补充完整了

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
static float DistributeBeckmann(const Vector3f& wh, float alpha) {//Beckmann法线分布项
  float tanTheta = Frame::tanTheta(wh);
  float cosTheta = Frame::cosTheta(wh);
  float a = std::exp(-(tanTheta * tanTheta) / (alpha * alpha));
  float b = M_PI * alpha * alpha * std::pow(cosTheta, 4.0f);
  return a / b;
}

static float G1(const Vector3f& wv, const Vector3f& wh, float alpha) {//Beckmann几何遮蔽项
  float c = wv.dot(wh) / Frame::cosTheta(wv);
  if (c <= 0) {
    return 0;
  }
  float b = 1.0f / (alpha * Frame::tanTheta(wv));
  return b < 1.6f ? (3.535f * b + 2.181f * b * b) / (1.f + 2.276f * b + 2.577f * b * b) : 1;
}

Color3f eval(const BSDFQueryRecord& bRec) const {
  Vector3f wh = (bRec.wi + bRec.wo).normalized();
  auto cosThetaI = Frame::cosTheta(bRec.wi);
  auto cosThetaO = Frame::cosTheta(bRec.wo);
  auto d = DistributeBeckmann(wh, m_alpha);
  auto g = G1(bRec.wi, wh, m_alpha) * G1(bRec.wo, wh, m_alpha);
  auto f = fresnel(wh.dot(bRec.wi), m_extIOR, m_intIOR);
  //kd + ks
  return m_kd / M_PI + m_ks * ((d * f * g) / (4.0f * cosThetaI * cosThetaO));
}

float pdf(const BSDFQueryRecord& bRec) const {
  if (bRec.wo.z() <= 0) {
    return 0;
  }
  Vector3f wh = (bRec.wi + bRec.wo).normalized();
  float d = DistributeBeckmann(wh, m_alpha);
  float jacobian = 1 / (4.0f * abs(wh.dot(bRec.wo)));
  return m_ks * d * Frame::cosTheta(wh) * jacobian + (1 - m_ks) * Frame::cosTheta(bRec.wo) * INV_PI;
}

Color3f sample(BSDFQueryRecord& bRec, const Point2f& _sample) const {
  if (Frame::cosTheta(bRec.wi) <= 0) {
    return Color3f(0.0f);
  }
  if (_sample.x() > m_ks) {  //漫反射
    Point2f sample((_sample.x() - m_ks) / (1.f - m_ks), _sample.y());
    bRec.wo = Warp::squareToCosineHemisphere(sample);
  } else {  //高光
    Point2f sample(_sample.x() / m_ks, _sample.y());
    Vector3f wh = Warp::squareToBeckmann(sample, m_alpha);
    bRec.wo = ((2.0f * wh.dot(bRec.wi) * wh) - bRec.wi).normalized();
  }
  if (bRec.wo.z() < 0.f) {
    return Color3f(0.0f);
  }
  // Note: Once you have implemented the part that computes the scattered
  // direction, the last part of this function should simply return the
  // BRDF value divided by the solid angle density and multiplied by the
  // cosine factor from the reflection equation, i.e.
  // 应该返回 BRDF 计算结果除以立体角密度并乘以反射系数（注释提示的，我也不知道为啥）
  return eval(bRec) * Frame::cosTheta(bRec.wo) / pdf(bRec);
}

首先来验证下代码正确性


BRDF结果可视化
看起来都是正确的，再来跑一下示例场景


还是很像塑料的

Part 2

第二部分，我们需要实现最原始的暴力路径追踪，说是暴力，其实是因为对材质表面采样，光线能不能碰到光源全靠运气，采样结果会有很多噪点

在whitted-style光追中我们只实现了直接光照，即从摄像机发射一条射线，反射一次到光源。现在我们要加上环境光，也就是从摄像机发射一条射线，在场景中多次碰撞后最终射入光源，这种只生成一条光线，连结视点和光源的路径的算法就是路径追踪。

路径追踪是一个递归算法，nori建议我们使用循环代替递归。我们可以用一个临时变量保存光线和表面交互时，该表面对最终结果的贡献

递归算法需要终止条件，不然会无限递归下去，但是不可以直接停止递归而什么都不做，会造成能量丢失使结果不准确。这里引入俄罗斯轮盘赌（Russian Roulette）来决定光线什么时候停止递归。如果还能继续递归，要将路径贡献除以概率保证结果正确。nori建议至少递归3次再进行轮盘赌，以减少噪音

在src中创建path_mats.cpp

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
#include <nori/integrator.h>
#include <nori/scene.h>
#include <nori/emitter.h>
#include <nori/sampler.h>
#include <nori/bsdf.h>

NORI_NAMESPACE_BEGIN

class PathMats : public Integrator {
 public:
  PathMats(const PropertyList& props) {}

  Color3f Li(const Scene* scene, Sampler* sampler, const Ray3f& _ray) const override {
    Color3f color = 0;//最终结果
    Color3f t = 1;//本次弹射对最终结果的贡献
    Ray3f rayRecursive = _ray;
    float probability;//本次继续的概率
    int depth = 1;//递归深度
    while (true) {
      Intersection its;
      if (!scene->rayIntersect(rayRecursive, its))
        break;
      //光源贡献
      if (its.mesh->isEmitter()) {
        EmitterQueryRecord lRecE(rayRecursive.o, its.p, its.shFrame.n);
        color += t * its.mesh->getEmitter()->eval(lRecE);
      }
      //俄罗斯轮盘赌
      if (depth >= 3) {
        probability = std::min(t.maxCoeff(), 0.99f);//选择路径贡献中最大项
        if (sampler->next1D() > probability)
          break;
        t /= probability;
      }
      BSDFQueryRecord bRec(its.shFrame.toLocal(-rayRecursive.d));
      Color3f f = its.mesh->getBSDF()->sample(bRec, sampler->next2D());
      t *= f;//将贡献叠加到路径上
      //继续递归
      rayRecursive = Ray3f(its.p, its.toWorld(bRec.wo));
      depth++;
    }
    return color;
  }

  std::string toString() const {
    return "PathMats[]";
  }
};

NORI_REGISTER_CLASS(PathMats, "path_mats");
NORI_NAMESPACE_END

看看结果
512次采样
256次采样

Part 3

暴力追踪的结果噪音很大，一个解决办法就是whitted-style里已经实现过的光源重要性采样，这同样适用于路径追踪，基本原理在作业4分析过了，直接上代码吧

在path_ems.cpp

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
#include <nori/integrator.h>
#include <nori/scene.h>
#include <nori/emitter.h>
#include <nori/sampler.h>
#include <nori/bsdf.h>

NORI_NAMESPACE_BEGIN

class PathEms : public Integrator {
 public:
  PathEms(const PropertyList& props) {}

  Color3f Li(const Scene* scene, Sampler* sampler, const Ray3f& _ray) const override {
    Color3f color = 0;//最终结果
    Color3f t = 1;//本次弹射对最终结果的贡献
    Ray3f rayRecursive = _ray;
    float probability;
    int depth = 1;
    int isDelta = 1;//判断是不是diffuse材质
    while (true) {
      Intersection its;
      if (!scene->rayIntersect(rayRecursive, its))
        break;
      if (its.mesh->isEmitter()) {//光源贡献
        EmitterQueryRecord lRecE(rayRecursive.o, its.p, its.shFrame.n);
        color += t * its.mesh->getEmitter()->eval(lRecE) * isDelta;
      }
      if (its.mesh->getBSDF()->isDiffuse()) {
        auto light = scene->getRandomEmitter(sampler->next1D());//均匀采样一个光源
        EmitterQueryRecord lRec(its.p);
        Color3f Li = light->getEmitter()->sample(light, lRec, sampler);//采样出射方向
        if (scene->rayIntersect(lRec.shadowRay)) {//出射方向有没有碰到物体
          Li = 0;
        }
        float cosTheta = Frame::cosTheta(its.shFrame.toLocal(lRec.wi));
        BSDFQueryRecord bRec(its.toLocal(-rayRecursive.d), its.toLocal(lRec.wi), ESolidAngle);
        Color3f f = its.mesh->getBSDF()->eval(bRec);//计算bsdf
        color += Li * f * cosTheta * t;//叠加到最终结果上
        isDelta = 0;
      } else {
        isDelta = 1;//不是diffuse就下次递归再说
      }
      //俄罗斯轮盘赌
      if (depth >= 3) {
        probability = std::min(t.maxCoeff(), 0.99f);
        if (sampler->next1D() > probability)
          break;
        t /= probability;
      }
      //采样下一个方向
      BSDFQueryRecord bRec(its.shFrame.toLocal(-rayRecursive.d));
      Color3f f = its.mesh->getBSDF()->sample(bRec, sampler->next2D());
      t *= f;//叠加路径贡献
      rayRecursive = Ray3f(its.p, its.toWorld(bRec.wo));
      depth++;
    }
    return color;
  }

  std::string toString() const {
    return "PathEms[]";
  }
};

NORI_REGISTER_CLASS(PathEms, "path_ems");
NORI_NAMESPACE_END

512次采样
256次采样
可以看到，同样的采样次数，整个画面噪点比暴力少了很多，但是，第二张图，明显看到越光滑的板，高光的噪音越明显，造成这个现象的原因是长板越光滑，它就越类似镜面反射的brdf，也就是说只有一个小范围才能对结果有贡献，而光源很大，均匀采样的时候很难找到有贡献的范围。

Part 4

既然对BRDF采样和对光源采样都有缺陷，我们可以将它们结合起来。

https://zhuanlan.zhihu.com/p/360420413

个人理解的多重重要性采样，其实是将两种采样方法都试一次，然后按照启发公式来叠加两种pdf

nori建议使用的启发公式：

wLight(pLight,pBRDF)=pLightpLight+pBRDFw_\mathrm{Light}(p_\mathrm{Light}, p_\mathrm{BRDF}) = \frac{p_\mathrm{Light}}{p_\mathrm{Light} + p_\mathrm{BRDF}}

wBRDF(pLight,pBRDF)=pBRDFpLight+pBRDFw_\mathrm{BRDF}(p_\mathrm{Light}, p_\mathrm{BRDF}) = \frac{p_\mathrm{BRDF}}{p_\mathrm{Light} + p_\mathrm{BRDF}}

翻了很多文章和资料，关于为什么它是正确的（看不懂），不过不影响我们使用它（

还有一个重要的点，两种pdf必须以同样的度量表示，不能一个是立体角，一个是单位面积

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
75
76
77
78
79
80
81
82
83
84
#include <nori/integrator.h>
#include <nori/scene.h>
#include <nori/emitter.h>
#include <nori/sampler.h>
#include <nori/bsdf.h>

NORI_NAMESPACE_BEGIN

class PathMisIntegrator : public Integrator {
 public:
  PathMisIntegrator(const PropertyList& props) {}

  Color3f Li(const Scene* scene, Sampler* sampler, const Ray3f& ray) const override {
    Color3f color = 0;
    Color3f t = 1;
    Ray3f rayRecursive = ray;
    float probability;
    float w_mats = 1.0f;//BRDF权重，放在外面是给下一轮路径使用的
    int depth = 1;
    Intersection its;
    if (!scene->rayIntersect(rayRecursive, its)) {//如果没有碰撞直接返回
      return color;
    }
    while (true) {
      if (its.mesh->isEmitter()) {
        EmitterQueryRecord lRec(rayRecursive.o, its.p, its.shFrame.n);
        lRec.uv = its.uv;
        //计算光源贡献时叠加BRDF权重
        color += t * w_mats * its.mesh->getEmitter()->eval(lRec);
      }
      //光源重要性采样
      const Mesh* mesh = scene->getRandomEmitter(sampler->next1D());
      const Emitter* light = mesh->getEmitter();
      EmitterQueryRecord lRec(its.p);
      lRec.uv = its.uv;
      Color3f Li = light->sample(mesh, lRec, sampler) * scene->getEmitters().size();
      float pdf_em = light->pdf(mesh, lRec);//光源pdf
      if (!scene->rayIntersect(lRec.shadowRay)) {//没有遮挡物
        float cosTheta = std::max(0.f, Frame::cosTheta(its.shFrame.toLocal(lRec.wi)));
        BSDFQueryRecord bRec(its.toLocal(-rayRecursive.d), its.toLocal(lRec.wi), ESolidAngle);
        Color3f f = its.mesh->getBSDF()->eval(bRec);
        float pdf_mat = its.mesh->getBSDF()->pdf(bRec);//BRDF pdf
        //使用启发公式计算光源权重
        float w_ems = pdf_mat + pdf_em > 0.0f ? pdf_em / (pdf_mat + pdf_em) : pdf_em;
        color += Li * f * cosTheta * w_ems * t;//叠加到结果
      }
      if (depth >= 3) {//俄罗斯轮盘赌
        probability = std::min(t.maxCoeff(), 0.99f);
        if (sampler->next1D() > probability) {
          return color;
        }
        t /= probability;
      }
      //BRDF采样
      BSDFQueryRecord bRec(its.shFrame.toLocal(-rayRecursive.d));
      Color3f f = its.mesh->getBSDF()->sample(bRec, sampler->next2D());
      t *= f;
      rayRecursive = Ray3f(its.p, its.toWorld(bRec.wo));
      float pdf_mat = its.mesh->getBSDF()->pdf(bRec);//BRDF pdf
      Point3f origin = its.p;
      if (!scene->rayIntersect(rayRecursive, its)) {
        return color;
      }
      if (its.mesh->isEmitter()) {//如果是光源就更新BRDF权重
        EmitterQueryRecord newLRec = EmitterQueryRecord(origin, its.p, its.shFrame.n);
        lRec.uv = its.uv;
        float new_pdf_em = its.mesh->getEmitter()->pdf(its.mesh, newLRec);
        w_mats = pdf_mat + new_pdf_em > 0.f ? pdf_mat / (pdf_mat + new_pdf_em) : pdf_mat;
      }
      if (bRec.measure == EDiscrete) {//如果PDF不是立体角上计算
        w_mats = 1.0f;
      }
      depth++;
    }
    return color;
  }

  std::string toString() const {
    return "PathMisIntegrator[]";
  }
};

NORI_REGISTER_CLASS(PathMisIntegrator, "path_mis");
NORI_NAMESPACE_END

256次采样
128次采样
看着好像噪点还是有点多？实际采样数只有前面的一半，但效果几乎差不多，而且第二张图中，光滑平板的高光问题也没有了

再来张水杯，效果太棒了！


再跑一下测试

作业 > Nori
#图形学
[Nori]作业5
https://ksgfk.github.io/2021/08/06/Nori-作业5/
作者
ksgfk
发布于
2021年8月6日
许可协议