PBRT筆記(10)——體積散射

體散射處理過程

3個影響參與介質在環境中的輻射度分佈的主要因素：

1. 吸取：減小光能，並將其轉化爲別的能量，例如熱量。
2. 發光：由光子發射光能至環境中。
3. 散射：因爲粒子碰撞，使得一個方向的輻射度散射至其餘方向。

吸取

吸取被描述爲一段介質截面區域。每單位距離媒介密度與吸取光能的比被定義爲\(\sigma_a\)。截面區域或隨着位置與方向的變化而變化，儘管如此，一般仍是會定義爲位置函數。它一般也被定義爲一個光譜變化量，\(\sigma_a\)與距離成反比關係\(m^{-1}\)。這意味這個值無需被判斷是否處於(0,1)。

沿長度爲dt的微分光線的輻射度變化可用微分方程來描述：
\(L_o(p,\omega)-L_i(p,-\omega)=dL_o(p,\omega)=-\sigma_a(p,\omega)L_i(p,-\omega)dt\)

求解分爲方程後得積分方程。若是假設光線於介質內p點處，沿方向ω行進距離d，則所有吸取量爲：
\(e^{-\int_0^d\sigma_a(p+t\omega,\omega)dt}\)

自發光

與吸取相對應，自發光過程是將能量轉換成可見光的化學、熱或核過程。自發光的輻射度變化微分方程爲：
\(dL_o(p,\omega)=L_e(p,\omega)dt\)

但這是基於入射光不會對發光產生影響的前提下創建的假設。

外散射和衰減

當光線穿過某一介質時，因爲粒子碰撞，致使其會往別的方向散射。它減小了出射位置處出射的輻射度，這種現象被稱爲外散射，與之對應的被稱爲內散射。

每單位距離外散射現象的發生機率由散射係數\(\sigma_s\)決定。微分長度dt上的外散射輻射度減小量爲：
\(dL_o(p,\omega)=-\sigma_s(p,\omega)L_i(p,-\omega)dt\)

總共減小的輻射度由吸取量與外散射量求和獲得。這兩種效果綜合被稱爲衰減或者消光。
\(\sigma_t(p,\omega)=\sigma_a(p,\omega)+\sigma_s(p,\omega)\)二者相集合能夠獲得如下兩個變量：

1. albedo：
   \(\rho=\frac{\sigma_s}{\sigma_t}\)
   他描述了發生散射現象的機率
2. 平均自由路徑:1/σt即發生散射現象前，光線在介質中行進的平均路徑。

對於給定的衰減係數σt，描述整個衰減的微分方程爲：\(\frac{dL_o(p,\omega)}{dt}=-\sigma_t(p,\omega)L_i(p,-\omega)\)

求解以後可得：
\(T_r(p\rightarrow p')=e^{-\int^d_0 \sigma_t(p+t\omega,\omega)dt}\)

d表示p到p'的距離，Tr表示透射率。
如圖11.8所示：\(T_r(p\rightarrow p'')=T_r(p\rightarrow p')+T_r(p'\rightarrow p'')\)

Tr中的逆置屬性τ被記做光學厚度：
\(\tau(p\rightarrow p')=\int^d_0 \sigma_t(p+t\omega,-\omega)dt\)

在均一介質中，σt是一個常量，根據Beer定律，\(T_r(p\rightarrow p')=e^{-\sigma_td}\)

內散射

內散射是由於其餘方向上發生的散射現象，從而增長了出射方向的輻射度的現象。

在忽略介質粒子間互相做用的狀況下，phase函數p(ω, ω')描述了在這一點輻射度的角度分佈狀況。它用於體積模擬有點相似BSDF。phase函數有一個歸一化條件，對於全部ω必須遵照：
\(\int_{g2}p(\omega,\omega')d\omega '=1\)這就意味着phase函數定義特定方向的散射的分佈機率。

根據內散射，內散射每單位距離增長的的輻射度爲：
\(dL_o(p,\omega)=L_s(p,\omega)dt\)

他包含了自發光與內散射現象：
\(L_s(p,\omega)=L_e(p,\omega)+\sigma_s(p,\omega)\int_{g2}L_i(p,\omega_i)d\omega_i\)

相位函數

目前有許多相位函數，通常分爲參數化模型(擁有少許變量的近似擬合函數)和散射輻射度解析模型。

大多數自然介質，其相位函數爲角度在ωo與ωi之間的一維函數，這些被常常寫做p(cos θ)。
由於在旋轉時，不會對入射光產生影響，因此被稱爲各向同性。

除了歸一化以外，另外一個重要特徵是:入射方向與出射方向是能夠互換的，且phase函數的值不變。

在由排列成凝聚結構的粒子組成的各向異性介質中，它的相位函數是帶兩個方向參數的4維函數。它具備較爲複雜的互相做用關係。（主要是晶體之類的）

相位函數自身既能夠呈現爲各項同性也呈現爲各項異形的。由於與方向無關，因此能夠定義爲：
\(p(\omega_o,\omega_i)=\frac{1}{4\pi}\)

PhaseFunction類位於medium.h中，
p()返回相位函數對於指定方向的值。

Henyey與Greenstein在1941年開發了一種相位函數，現今依然被普遍使用。它的設計目的是爲了便於與實測散射數據進行擬合。參數g（被稱爲非對稱變量）用於控制散射分佈。PhaseHG()對其進行了實現，其公式爲：
\(P_{HG}(cos\theta)=\frac{1}{4}\frac{1-g^2}{(1+g^2+2g(cos\theta))^{3/2}}\)

這個模型的g值畢竟在(-1,1)範圍內，g值爲負時表明後向散射，光線主要散射回入射方向，g值爲正時表明前向散射。g值越大散射現象就越進階入射或者出射方向。

HenyeyGreenstein類實現了這個模型。

非對稱變量g在這個模型中擁有明確意義，即爲相位函數估算值與cos(ω',ω)乘積的均值。

對於任意相位函數p,g的值能夠計算爲：

\(g=\int_{g^2}p(-\omega,\omega')(\omega\cdot \omega')d\omega '=2\pi\int^\pi_0p(-cos\theta)cos\theta sin\theta d\theta\)

因此各向同性介質的相位函數的g值爲0。

任意數量的相位函數均可以知足這個方程;光靠一個g值不足以描述散射分佈。然而方便地將一個複雜的散射分佈轉換成一個簡單的參數化模型每每更加劇要。

複雜的相位函數使用一個非對稱變量是不可以很好地描述分佈狀況的，因此通常會採用相位函數加權和的方式來建模：

\(p(\omega,\omega')=\sum_{i=1}^{n}w_i p_i (\omega \rightarrow \omega')\)

其中wi的和須要被保持標準化爲1。

介質

這段的徹底機翻

介質基類實現了多種特性的體積散射屬性，在複雜場景中可能有多個介質實例，用於表現不一樣的散射效果。

一個關鍵操做是介質類必須計算光束透射率。對應的函數爲Tr()，這個方法返回光線起點到終點的透射率的估算值。

Medium-aware積分器負責計算與光線相交的介質幾何圖元的透射率，具體會使用Tr()與蒙特卡洛方法解積分方程的方式。所以咱們假設傳遞給Tr()的射線是沒有被阻擋，且光線徹底處於當前介質中。

pbrt會將Medium實例與攝像機、燈光或者幾何圖元結合使用，用於模擬燈光霧等空間分佈效果。

在pbrt中，爲了表現場景中一塊區域的介質，會使用幾何圖元來表現，具體的方式是使用MediumInterface類，它存儲了外部與內部介質的指針。而不會只存儲一種介質。（由於不存在什麼介質都沒有的狀況）固然PBRT容許將測試設爲nullptr，即不對光線產生任何影響。

bool GeometricPrimitive::Intersect(const Ray &r,
                                   SurfaceInteraction *isect) const {
    Float tHit;
    if (!shape->Intersect(r, &tHit, isect)) return false;
    r.tMax = tHit;
    isect->primitive = this;
    CHECK_GE(Dot(isect->n, isect->shading.n), 0.);
    // Initialize _SurfaceInteraction::mediumInterface_ after _Shape_
    // intersection
    if (mediumInterface.IsMediumTransition())
        isect->mediumInterface = mediumInterface;
    else
        isect->mediumInterface = MediumInterface(r.medium);
    return true;
}

可是這樣有可能會出現不一樣物體所指定的外部介質不同的狀況。在這種狀況下，離開圖元朝向相機的光線與離開相機朝向圖元的光線將被視爲處於不一樣的介質中。所以光傳輸算法將沒法計算一致的結果。pbrt認爲用戶可以在場景中指定一致的介質配置，而且不值得爲此增長代碼的複雜性。（畢竟是渲染核心，指定參數這種操做仍是編輯器來操做比較好）

對於場景空間中的區域,咱們將實現Scene::IntersectTr()方法,若是相交成功成功，它返回第一個相交的會發生光散射現象的SurfaceInteraction的指針。否則就返回false。

bool Scene::IntersectTr(Ray ray, Sampler &sampler, SurfaceInteraction *isect,
                        Spectrum *Tr) const {
    *Tr = Spectrum(1.f);
    while (true) {
        bool hitSurface = Intersect(ray, isect);
        // Accumulate beam transmittance for ray segment
        if (ray.medium) *Tr *= ray.medium->Tr(ray, sampler);

        // Initialize next ray segment or terminate transmittance computation
        if (!hitSurface) return false;
        if (isect->primitive->GetMaterial() != nullptr) return true;
        ray = isect->SpawnRay(ray.d);
    }
}

均勻介質

均勻介質是最爲簡單的介質，它表明了一個具備恆定σa和σs的空間，它使用了Henyey–Greenstein相位函數來表現這個介質中發生的散射現象，在這個相位函數中具備恆定的g值。類名爲HomogeneousMedium。

由於整個媒介σt是恆定的,因此根據Beer定律能夠用來計算沿光線的透射率。這裏須要注意浮點數運算偏差的問題。使用浮點數的正無窮來初始化Ray::tMax，它確保了光線足夠能夠取得最大範圍內的相交信息。然而，對於Ray::tMax使用Infinity的狀況下，在應用Beer定律時會產生了一個小問題。原則上，咱們只須要計算射線參數t範圍，乘以射線方向的長度,而後乘以σt，可是正無窮乘以0會獲得NaN，具體發生在當光線經過一個吸取率爲零的介質無限遠時，上面的代碼將嘗試執行0 *∞的乘法運算，最終獲得NaN，而不是0的預期透射率。因此這裏須要作容錯處理。

Spectrum HomogeneousMedium::Tr(const Ray &ray, Sampler &sampler) const {
    ProfilePhase _(Prof::MediumTr);
    return Exp(-sigma_t * std::min(ray.tMax * ray.d.Length(), MaxFloat));
}

3D網格

GridDensityMedium類講介質的密度存儲在3D網格中， 相似於ImageTexture類用2D網格採樣來表示圖像的方式。GridDensityMedium經過對這些樣本進行插值以計算他們之間的位置採樣點的介質密度。

GridDensityMedium::Tr()會調用Density()，利用提供的樣本重構給定點的體積密度函數。在這個過程當中它將會用WorldToMedium()把世界座標轉換成局部座標。σa和σs將會按照對應的位置進行插值運算。

首先調用Density()以周圍樣本座標與距離插值計算點的座標。（具體能夠參照7.1.7節，說到底就是個雙線性插值）

Float GridDensityMedium::Density(const Point3f &p) const {
    // Compute voxel coordinates and offsets for _p_
    Point3f pSamples(p.x * nx - .5f, p.y * ny - .5f, p.z * nz - .5f);
    Point3i pi = (Point3i)Floor(pSamples);
    Vector3f d = pSamples - (Point3f)pi;

    // Trilinearly interpolate density values to compute local density
    Float d00 = Lerp(d.x, D(pi), D(pi + Vector3i(1, 0, 0)));
    Float d10 = Lerp(d.x, D(pi + Vector3i(0, 1, 0)), D(pi + Vector3i(1, 1, 0)));
    Float d01 = Lerp(d.x, D(pi + Vector3i(0, 0, 1)), D(pi + Vector3i(1, 0, 1)));
    Float d11 = Lerp(d.x, D(pi + Vector3i(0, 1, 1)), D(pi + Vector3i(1, 1, 1)));
    Float d0 = Lerp(d.y, d00, d10);
    Float d1 = Lerp(d.y, d01, d11);
    return Lerp(d.z, d0, d1);
}

爲了插值計算樣本週圍的點

D()返回給定整數座標的樣本密度，它的惟一任務是處理超出邊界的樣本位置，併爲給定的樣本計算適當的數組偏移量。與MIPMaps不一樣，在區域外密度講返回爲0。

BSSRDF

雙向散射表面反射分佈函數：
\(L_o(p_o,\omega_o)=\int_A\int_{H^2(n)} S(p_o,\omega_o,p_i,\omega_i)L_i(p_i,\omega_i) \left|cos\theta_i \right|d\omega_idA\)

次表面光線傳輸使用11.1和11.2中介紹的體積散射過程與15.1中介紹的體積光傳輸方程來描述。BSSRDF中的S是對邊界上給定的一對點和方向之間發生的散射現象進行建模的一種歸納表示。

類名爲BSSRDF，位於core/bssrdf.h和bssrdf.cpp。
BSSDF必須實現傳遞當前相交表面以及散射介質折射率給基類的構造函數。所以，這裏假設在這個介質中的折射率是不變的，這也是BSSRDF模型中普遍使用的假設。

BSSRDF實現必須提供的關鍵方法是估算八維分佈函數S()。

爲了支持通常形狀圖元，咱們將引入一個更簡單、可分離的SeparableBSSRDF。

SeparableBSSRDF(const SurfaceInteraction &po, Float eta,
const Material *material, TransportMode mode)
: BSSRDF(po, eta), ns(po.shading.n), ss(Normalize(po.shading.dpdu)),
ts(Cross(ns, ss)), material(material), mode(mode) { }
const Normal3f ns;
const Vector3f ss, ts;
const Material *material;
const TransportMode mode;

SeparableBSSRDF接口將將BSSRDF轉換爲三個獨立的組成部分的可分離形式(一個空間和兩個方向):
\(S(p_o,\omega_o,p_i,\omega_i)\approx(1-F_r(cos\theta_o))S_p(p_o,p_i)S_\omega(\omega_i)\)

Spectrum S(const SurfaceInteraction &pi, const Vector3f &wi) {
    Float Ft = 1 - FrDielectric(Dot(po.wo, po.shading.n), 1, eta);
    return Ft * Sp(pi) * Sw(wi);
}

對於上面給出的可分離變量的表達式，由次表面散射引發的照度的積分可簡化爲：
\(L_o(p_o,\omega_o)=\int_A\int_{H^2(n)}S(p_o,\omega_o,p_i,\omega_i)L_i(p_i,\omega_i)\left|cos\theta_i\right|d\omega_idA(p_i)=(1-F_r(cos\theta_o))\int_AS_p(p_o,p_i)\int_{H^2(n)}S_\omega(\omega_i)L_i(p_i,\omega_i)\left|cos\theta_i\right|d\omega_idA(p_i)\)
咱們定義Sω(ωi)做爲擴展版本的菲涅耳透光率
\(S_\omega(\omega_i)=\frac{1-Fr(cos\theta_i)}{c\pi}\)

歸一化因子c是精心挑選的，因此Sw比上cos權重半球：
\(\int_{H^2}S_\omega(\omega)cos\theta d\omega=1\)
換句話說：
\(c=\int^{2\pi}_0\int^{\frac{\pi}{2}}_{0}\frac{1-Fr(\eta,cos\theta)}{\pi}sin\theta\cos\theta d\theta d\phi=1-2\int^{\frac{\pi}{2}}_{0}F_r(\eta,cos\theta)sin\theta cos\theta d\theta\)

這個積分稱爲菲涅爾反射函數的第一時刻。第i個菲涅耳時刻的通常定義是：

\(F_r,i(\eta)=\int^{\frac{\pi}{2}}_{0}F_r(\eta,cos\theta)sin\theta cos^i \theta d\theta\)

pbrt提供了FresnelMoment1()和FresnelMoment2()，用於擬合來計算相應的菲尼爾時刻。使用FresnelMoment1()能夠很方便地實現SeparableBSSRDF::Sw()

Spectrum Sw(const Vector3f &w) const {
    Float c = 1 - 2 * FresnelMoment1(1 / eta);
    return (1 - FrDielectric(CosTheta(w), 1, eta)) / (c * Pi);
}

將空間和方向參數解耦大大減小了S的維數,但沒有解決支持通常形狀實現較爲困難的問題。這裏咱們將引入第二個近似，它假設曲面不只是局部平面的，並且影響BSSRDF值的是點之間的距離，而不是它們的實際位置。這將Sp簡化爲一個只涉及兩點po和pi的距離的函數Sr:

\(S_p(p_o,p_i)\approx S_r(\left|\right|p_o-p_i\left|\right|)\)

Spectrum Sp(const SurfaceInteraction &pi) const {
    return Sr(Distance(po.p, pi.p));
}

SeparableBSSRDF在平面形狀的狀況下效果較好，隨着形狀逐漸偏離平面，偏差將會愈來愈大。

TABULATED BSSRDF

類名爲TabulatedBSSRDF，是SeparableBSSRDF的子類。TabulatedBSSRDF 實現了父類的接口，它提供表化的BSSRDF的表現方法。它能夠控制普遍的次表面參數，能夠用於模擬真實世界中的測算的BSSRDFs。

TabulatedBSSRDF使用與FourierBSDF 反射模型相同類型的自適應樣條插值方法。

次表面散射材質

SubsurfaceMaterial，定義於materials/subsurface.h與materials/subsurface.cpp。KdSubsurfaceMaterial,定義於materials/kdsubsurface.h與materials/kdsubsurface.cpp。二者的區別在於如何設定散射屬性。

const Float scale;
std::shared_ptr<Texture<Spectrum>> Kr, Kt, sigma_a, sigma_s;
std::shared_ptr<Texture<Float>> uRoughness, vRoughness;
std::shared_ptr<Texture<Float>> bumpMap;
const Float eta;
const bool remapRoughness;
BSSRDFTable table;

ComputeScatteringFunctions()使用貼圖計算指定座標的散射屬性值，使用成員變量scale來縮放吸取值與散射係數。

void SubsurfaceMaterial::ComputeScatteringFunctions(SurfaceInteraction *si, MemoryArena &arena, TransportMode mode,bool allowMultipleLobes) const {
    <Perform bump mapping with bumpMap, if present >
    <Initialize BSDF for SubsurfaceMaterial>
    Spectrum sig_a = scale * sigma_a->Evaluate(*si).Clamp();
    Spectrum sig_s = scale * sigma_s->Evaluate(*si).Clamp();
    si->bssrdf = ARENA_ALLOC(arena, TabulatedBSSRDF)(*si, this, mode, eta, sig_a, sig_s, table);
}