14 Light Transport 1: Surface Reflection

用蒙特卡洛方法来估算LTE

\[\begin{align*} L_o(p,\omega_o)&=\int_{s^2}f(p,\omega_o,\omega_i)L_i(p,\omega_i)cos\theta_i d\omega_i \\ &= \frac{1}{N}\sum_{j=1}^N \frac{f(p, \omega_o, \omega_j)L_i(p,\omega_j)cos\theta_j}{p(\omega_j)} \end{align*}\]

14.1 Sampling Reflection Functions

Spectrum BxDF::Sample_f(const Vector3f &wo, Vector *wi, const Point2f &u, Float *pdf, BxDFType *sampledType) const

输入一个出射方向wo，二位随机变量u，根据BxDF的分布采样一个wi，计算采样wi的pdf，并返回wi和wo计算得到的brdf的值。默认的实现是半球上的cosine分布。

14.1.3 Specular Reflection and Transmission (Handle Dirac Delta)

Dirac Delta 函数的定义: \(\delta = 0 \quad for \quad all \ x \ne 0 \\\) 且 \(\int_{-\inf}^{\inf}\delta(x)dx = 1\)

这个函数没法积分。但在Monte Carlo Estimator中，通常会同时出现在分子和分母中，通常的处理方式是在公式中直接消去。

14.1.5 Application: Estimating Reflectance

这里用Hemisphere-Direction Reflectance的计算作为Monte Carlo方法的一个例子。

\[\begin{align*} \rho_{hd}(\omega_o)&=\int_{H^2(n)}f_r(\omega_o, \omega_i) |cos\theta|d\omega_i \\ &=\frac{1}{N} \sum_{j}^N \frac{f_r(\omega_o, \omega_j)|cos\theta|}{p(\omega_j)} \end{align*}\]

BxDF::rho() 实现了这个过程。

也用Hemisphere-Hemisphere Reflectance ($\rho_{hh}$)做了例子，但不太明白 $\rho_{hh}$ 是用来计算什么东西的。

14.2 Sampling Lighting Sources

通过采样BxDF来命中场景中的光源，效率会很低。因此需要应用 Multiple Importance Sampling 方法，直接对光源进行采样。PBRT中采样光源的函数是：

virtual Spectrum Light::Sample_Li(const Interaction &ref, const Point2f &u, Vector3f *wi, Float *pdf, VisibilityTester *vis) const = 0

14.2.1 Lights with Singularities

Dirac Delta形式的光源（如点光源，方向光），采样到光源位置或方向的概率为0，处理方式也是在公式的分母分子中直接把Delta项消掉。

14.3 Direct Lighting

对于一根Camera Ray，可以一次采所有的光源，或者一根Ray只采一个光源。

一根Ray采所有的光源：

def UnifromSampleAllLights():
    L = Spectrum(0.0f)
    for light in all lights:
        Ld = Spectrum(0.0f)
        for i = 0..nSamples:
            Ld += EstimateDirect(light) # Monte carlo estimator.        
        Ld /= nSamples  # Normalize
    L += Ld
    return L

一根Ray只采一个光源:

def UniformSampleOneLight():
    # 随机选择场景中的一个light
    nLight = scene.getLightNum()
    lightIdx = random() * nLight
    light = scene.getLight(lightIdx)
    # 采这个light的pdf
    lightPdf = 1.0 / nLight    
    # 采样这个light
    L = EstimateDirect(light)
    # Monte Carlo Weighting
    L /= lightPdf   # == L * nLight
    return L

EstimateDirect的计算采用了MIS，分别采BxDF和光源。它要计算的项为： \(\frac{f(p,\omega_o,\omega_i)L_d(p,\omega_i)|cos \theta_i|}{p(\omega_i)}\) 其中 $L_d$ 的下标 $d$ 表示 direct（直接光）。

def EstimateDiect(light, isect):
    
    # Sample light source.

    wi = Vector3f()
    lightPdf = 0.0
    Li = light.sample_Li(wi, lightPdf)   # 获取light的radiance，wi 和 pdf
    # 确定了wi，计算bsdf，和scatteringPdf
    f = isect.bsdf.f(isect.wo, wi)    
    scatteringPdf = isect.bsdf.pdf(isect.wo, wi)
    # the |cos(theta)| term
    f *= AbsDot(wi, isect.n)

    shadowRay = Ray(isect.p, wi)
    if scene.Occlude(shadowRay):
        # 遮挡了了
        Li = Spectrum(0.0)

    # 这里有个假设：如果bsdf带delta项，那采样光源得到的wi，不可能刚好撞上bsdf的delta方向。

    if light.isDelta():
        # 如果是Delta光源，不需要进行MIS，因为从概率上来讲采样bsdf得到wi不可能采到这个光源。
        return Li * f / lightPdf        # 普通Monte Carlo Estimate
    else:
        Ld = Spectrum(0.0)

        # 不是Delta光源，应用MIS。        
        weight = ligthPdf / (lightPdf + scatteringPdf)  # Balance heuristic
        Ld += Li * f * weight / lightPdf

        # 采样bsdf
        f = isect.bsdf.sample_f(wi, scatteringPdf, bsdfType)
        specularBsdf = bsdfType & BSDFTypeSpecular      # bsdf也有可能带有delta项

        # Spawn a ray to intersect with scene
        lightSect = Intersection()
        shadowRay = Ray(isect.p, wi)
        if scene.Occlude(shadowRay, lightSect):
            # Intersect with some geometry.
            if lightSect.getEmitter() == light:
                # 交到了当前正在计算的light
                Li = lightSect.Le(-wi)
            else:
                # 没有交到当前在估算的这个light
                Li = Spectrum(0.0)
        else:
            # Not intersect，这种情况下Le计算地是啥？方向光？
            Li = light.Le(-wi)

        weight = 1.0
        if specularBsdf:
            # 如果是带delta项的bsdf，也不需要MIS。且上面采样光源时计算得到的f必然是0，
            # 因此之前计算得到的Ld为0，值完全由采样bsdf的结果决定。
            pass            
        else:
            # 不带delta项的bsdf，计算bsdf的MIS的权重。
            weight = scatteringPdf / (lightPdf + scatteringPdf) # Balance heuristic
        
        Ld += Li * f * weight / scatteringPdf

        return Ld

14.4 The Light Transport Equation

14.4.1 Basic Derivation

不考虑光的波动属性，认为场景中光能的分布处于一个稳态（Equilibrium），因此对应场景中的任意一个位置： \(\Phi_o - \Phi_i = \Phi_e - \Phi_a\) 即从Flux的角度来说，输出的Flux减去输入的Flux的差值等于发射出的Flux减去吸收的Flux的差值。

由此： \(\Phi_o = \Phi_e + (\Phi_i - \Phi_a)\) 从Radiance的角度来看： \(L_o(p, \omega_o)=L_e(p,\omega_o)+\int_{S^2}f(p,\omega_o,\omega_i)L_i(p,\omega_i) |cos \theta_i|d\omega_i\)

用 ray-casting function \(t(p, \omega)\) 表示从 $p$ 点出发，沿着 $\omega$ 方向和场景交到的点，则： \(L_i(p,\omega_i) = L_o(t(p,\omega_i), -\omega_i)\)

这样可以将式子中的 $L$ 项都用 $L_o$ 来表示。进一步去掉下标 $o$ 后： \(L(p, \omega_o)=L_e(p,\omega_o)+\int_{S^2}f(p, \omega_o, \omega_i)L(t(p,\omega_i),-\omega_i)|cos \theta|d\omega_i\)

14.4.2 Analytic Solutions to LTE

只有在一些极简情况下有解析解，大多数情况下是没有解析解的。解析解可在极简情况下作为标准答案用于检查计算结果是否正确。

14.4.3 The Surface Form of the LTE

将基于方位角的公式变为基于场景中点和点之间radiance传输的公式： \(L(p' \rightarrow p) = L(p',\omega)\) \(f(p''\rightarrow p' \rightarrow p) = f(p', \omega_o, \omega_i)\) \(G(p\leftrightarrow p')=V(p\leftrightarrow p')\frac{|cos\theta|\ |cos\theta'|}{||p'-p||^2}\) G 这个项中的： \(\frac{|cos\theta'|}{||p'-p||^2}\) 是将积分由对方位角积分转换为对场景中的面积积分所需的Jacobian项。

这样原来的LTE变成了： \(L(p'\rightarrow p)=L_e(p'\rightarrow p)+\int_Af(p'' \rightarrow p' \rightarrow p)L(p''\rightarrow p')G(p''\leftrightarrow p')dA(p'')\)

14.4.4 Integral over Paths

将 L 递归的代入：

\[\begin{align*} L(p_1\rightarrow p_0)&=L_e(p_1\rightarrow p_0)+\int_A L(p_2\rightarrow p_1)f(p_2 \rightarrow p_1 \rightarrow p_0)G(p_2\leftrightarrow p_1)dA(p_2) \\ &=L_e(p_1\rightarrow p_0)+\int_A [ L_e(p_2\rightarrow p_1)+ \\ & \quad \int_Af(p_3 \rightarrow p_2 \rightarrow p_1)L(p_3\rightarrow p_2)G(p_3\leftrightarrow p_2)dA(p_3) ] f(p_2 \rightarrow p_1 \rightarrow p_0 ) G(p_2\leftrightarrow p_1)dA(p_2)\\ &=L_e(p_1\rightarrow p_0)+\int_A L_e(p_2\rightarrow p_1) f(p_2 \rightarrow p_1 \rightarrow p_0 )G(p_2\leftrightarrow p_1)dA(p_2) + \\ & \quad \int_A \int_A L(p_3\rightarrow p_2) f(p_3 \rightarrow p_2 \rightarrow p_1) G(p_3\leftrightarrow p_2) f(p_2 \rightarrow p_1 \rightarrow p_0 ) G(p_2\leftrightarrow p_1) dA(p_3) dA(p_2) \end{align*}\]

可以一直无穷的展开下去。

若定义

\[\bar{p}_n=p_0, p_1, p_2, ... \ p_n\]

为 n+1 个顶点组成的路径，且

\[P(\bar{p}_n) = \begin{matrix} \underbrace{ \int_A \int_A ... \int_A } & L_e(p_n \rightarrow p_{n-1}) \left[ \prod_{i=1}^{n-1}f(p_{i+1},p_i,p_{i-1})G(p_{i+1}\leftrightarrow p_i) \right] dA(p_2)dA(p_3)...dA(p_n)\\ n-1 \\ \end{matrix}\]

定义：

\[T(\bar{p}_n) = \prod_{i=1}^{n-1}f(p_{i+1},p_i,p_{i-1})G(p_{i+1}\leftrightarrow p_i)\]

则上式可以简化为：

\[\begin{matrix} P(\bar{p}_n) & = & \underbrace{ \int_A \int_A ... \int_A } & L_e(p_n \rightarrow p_{n-1}) T(\bar{p}_n) dA(p_2)dA(p_3)...dA(p_n)\\ & & n-1 \\ \end{matrix}\]

上面的 LTE 可改写为：

\[L(p_1 \rightarrow p_0) = \sum_{n=1}^{\infty}P(\bar{p}_n)\]

这样做有啥好处？将LTE的积分转换到了path-space中。

14.4.5 Delta Distributions in Integrand

同样是通过在分子分母中消去Delta项的方式来处理。

14.4.6 Partitioning the Integrand

按照Bounce分： \(\begin{matrix} L(p_1\rightarrow p_0) & = & \underbrace{P(\bar{p}_1)} & + & \underbrace{P(\bar{p}_2)} & + & \underbrace{\sum_{n=3}^{\infty}P(\bar{p}_n)} \\ & & Emission & & Direct & & Indirect \end{matrix}\)