pbrt将反射分为以下四类: 漫反射, 光滑镜面反射, 完美镜面反射, 回溯反射. 反射分布方程可以被分为各向同性与各向异性的, 若将物体旋转时在各个角度上具有相同的反射结果即位各向同性, 反之则为各向异性.
BSDF表示
#
pbrt中BxDF接口代表特定种类表面散射的实现, BSDF则是围绕BxDF指针做一层封装并提供额外功能.
几何设置与约定
#
pbrt中法线计算发生在由表面的切线, 副切线与法线组成的坐标系中, 分别与x , y , z x,y,z x , y , z
BxDF接口
#
单独的BRDF与BTDF方法都定义在BxDF接口中.
class BxDF
: public TaggedPointer<DiffuseTransmissionBxDF, DiffuseBxDF, CoatedDiffuseBxDF,
CoatedConductorBxDF, DielectricBxDF, ThinDielectricBxDF,
HairBxDF, MeasuredBxDF, ConductorBxDF, NormalizedFresnelBxDF> {
public:
// BxDF Interface
PBRT_CPU_GPU inline BxDFFlags Flags() const;
using TaggedPointer::TaggedPointer;
std::string ToString() const;
PBRT_CPU_GPU inline SampledSpectrum f(Vector3f wo, Vector3f wi,
TransportMode mode) const;
PBRT_CPU_GPU inline pstd::optional<BSDFSample> Sample_f(
Vector3f wo, Float uc, Point2f u, TransportMode mode = TransportMode::Radiance,
BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const;
PBRT_CPU_GPU inline Float PDF(
Vector3f wo, Vector3f wi, TransportMode mode,
BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const;
PBRT_CPU_GPU
SampledSpectrum rho(Vector3f wo, pstd::span<const Float> uc,
pstd::span<const Point2f> u2) const;
SampledSpectrum rho(pstd::span<const Point2f> u1, pstd::span<const Float> uc2,
pstd::span<const Point2f> u2) const;
PBRT_CPU_GPU inline void Regularize();
};
Flags方法用于返回上文所述的材质类型以及区分反射与透射, 部分光线传播算法会通过该返回值决定行为.
enum BxDFFlags {
Unset = 0,
Reflection = 1 << 0,
Transmission = 1 << 1,
Diffuse = 1 << 2,
Glossy = 1 << 3,
Specular = 1 << 4,
// Composite _BxDFFlags_ definitions
DiffuseReflection = Diffuse | Reflection,
DiffuseTransmission = Diffuse | Transmission,
GlossyReflection = Glossy | Reflection,
GlossyTransmission = Glossy | Transmission,
SpecularReflection = Specular | Reflection,
SpecularTransmission = Specular | Transmission,
All = Diffuse | Glossy | Specular | Reflection | Transmission
};
BxDF的关键方法是f, 它根据方向返回BxDF的值, 方向需要位于前文所述的材质的本地空间中. BxDF接口认为不同波长的光是被解耦的, 某个波长上的能量不会被反射到不同的波长上, 因此f的返回值通过SampledSpectrum来表达. 荧光材质会在不同的波长间重新分布能量, 需要返回n × n n \times n n × n n n n BxDF的构造函数以及各种方法都没有SampledSpectrum中具体某个波长的信息, 这是不需要的. TransportMode用于表示出射方向是对着相机还是对着光源, 不对称散射会用到这一特性.
PBRT_CPU_GPU inline SampledSpectrum BxDF::f(Vector3f wo, Vector3f wi, TransportMode mode) const {
auto f = [&](auto ptr) -> SampledSpectrum { return ptr->f(wo, wi, mode); };
return Dispatch(f);
}
BxDF还需要提供用于重要性抽样的方法, Sample_f用于实现这一任务. 由于光线路径是逆向构造的, 因此这里出射方向作为参数, 采样得到的入射方向. 参数中的uc与u用于实现一维与二维的采样, 通常uc用于选择反射或透射, u用于选择方向.
PBRT_CPU_GPU inline pstd::optional<BSDFSample> BxDF::Sample_f(Vector3f wo, Float uc, Point2f u,
TransportMode mode,
BxDFReflTransFlags sampleFlags) const {
auto sample_f = [&](auto ptr) -> pstd::optional<BSDFSample> {
return ptr->Sample_f(wo, uc, u, mode, sampleFlags);
};
return Dispatch(sample_f);
}
用户可以通过sampleFlags参数限制采样结果为透射或反射, 像在不透明表面采样透射样本的不合法采样会导致采样失效.
enum class BxDFReflTransFlags {
Unset = 0,
Reflection = 1 << 0,
Transmission = 1 << 1,
All = Reflection | Transmission
};
若采样成功, BSDF的值, 入射方向, 分布密度以及对应的flags会被返回, 这里返回的wi位于本地空间中, 在BSDF中会转到渲染空间.
struct BSDFSample {
// BSDFSample Public Methods
BSDFSample() = default;
PBRT_CPU_GPU
BSDFSample(SampledSpectrum f, Vector3f wi, Float pdf, BxDFFlags flags, Float eta = 1,
bool pdfIsProportional = false)
: f(f),
wi(wi),
pdf(pdf),
flags(flags),
eta(eta),
pdfIsProportional(pdfIsProportional) {}
PBRT_CPU_GPU
bool IsReflection() const { return pbrt::IsReflective(flags); }
PBRT_CPU_GPU
bool IsTransmission() const { return pbrt::IsTransmissive(flags); }
PBRT_CPU_GPU
bool IsDiffuse() const { return pbrt::IsDiffuse(flags); }
PBRT_CPU_GPU
bool IsGlossy() const { return pbrt::IsGlossy(flags); }
PBRT_CPU_GPU
bool IsSpecular() const { return pbrt::IsSpecular(flags); }
std::string ToString() const;
SampledSpectrum f;
Vector3f wi;
Float pdf = 0;
BxDFFlags flags;
Float eta = 1;
bool pdfIsProportional = false;
};
半球反射
#
rho方法的两种重载分别用于计算第四章所介绍的半球-方向反射量与半球-半球反射量, 样本数量由调用者决定.
PBRT_CPU_GPU SampledSpectrum BxDF::rho(Vector3f wo, pstd::span<const Float> uc,
pstd::span<const Point2f> u2) const {
if (wo.z == 0)
return {};
SampledSpectrum r(0.);
DCHECK_EQ(uc.size(), u2.size());
for (size_t i = 0; i < uc.size(); ++i) {
// Compute estimate of $\rho_\roman{hd}$
pstd::optional<BSDFSample> bs = Sample_f(wo, uc[i], u2[i]);
if (bs && bs->pdf > 0)
r += bs->f * AbsCosTheta(bs->wi) / bs->pdf;
}
return r / uc.size();
}
SampledSpectrum BxDF::rho(pstd::span<const Point2f> u1, pstd::span<const Float> uc,
pstd::span<const Point2f> u2) const {
DCHECK_EQ(uc.size(), u1.size());
DCHECK_EQ(u1.size(), u2.size());
SampledSpectrum r(0.f);
for (size_t i = 0; i < uc.size(); ++i) {
// Compute estimate of $\rho_\roman{hh}$
Vector3f wo = SampleUniformHemisphere(u1[i]);
if (wo.z == 0)
continue;
Float pdfo = UniformHemispherePDF();
pstd::optional<BSDFSample> bs = Sample_f(wo, uc[i], u2[i]);
if (bs && bs->pdf > 0)
r += bs->f * AbsCosTheta(bs->wi) * AbsCosTheta(wo) / (pdfo * bs->pdf);
}
return r / (Pi * uc.size());
}
BSDF中的Delta分布
#
Delta分布主要用于完美镜面反射, 但是对于Sample_f与PDF方法, 返回Delta函数对应的无穷大密度在c++中是无法得到正确的渲染结果的. 对于Sample_f方法, 将Delta函数从PDF与BSDF分离, 可以发现Delta函数会被抵消, 由于PDF与Delta函数一致, 因此PDF设置为1即可. 对于PDF方法, 它的返回值为0, 因为恰好位于Delta函数对应的反射方向的概率过小.
f r ( p , ω o , ω i ) p ( ω i ) = δ ( ω ′ − ω i ) f r rem ( p , ω o , ω i ) δ ( ω ′ − ω i ) p rem ( ω i ) = f r rem ( p , ω o , ω i ) p rem ( ω i )
\begin{equation}
\frac{f_r(p, \omega_o, \omega_i)}{p(\omega_i)} = \frac{\delta(\omega' - \omega_i)f_r^{\text{rem}}(p, \omega_o, \omega_i)}{\delta(\omega' - \omega_i)p^{\text{rem}(\omega_i)}} = \frac{f_r^{\text{rem}}(p, \omega_o, \omega_i)}{p^{\text{rem}(\omega_i)}}
\end{equation}
p ( ω i ) f r ( p , ω o , ω i ) = δ ( ω ′ − ω i ) p rem ( ω i ) δ ( ω ′ − ω i ) f r rem ( p , ω o , ω i ) = p rem ( ω i ) f r rem ( p , ω o , ω i ) BSDF类
#
BSDF主要用于处理从材质本地空间到渲染空间的转换. BSDF的构造函数会传入几何法线, 着色法线以及切线, shadingFrame用于存储着色所用的坐标系以及处理与渲染空间之间的转换.
BSDF() = default;
PBRT_CPU_GPU
BSDF(Normal3f ns, Vector3f dpdus, BxDF bxdf)
: bxdf(bxdf), shadingFrame(Frame::FromXZ(Normalize(dpdus), Vector3f(ns))) {}
BSDF的f方法用于封装BxDF中的对应方法, 其中的模板版本用于在用户已知内部BxDF类型的情况下直接调用它的f方法, 这用于在GPU上避免pbrt中的动态方法分发.
PBRT_CPU_GPU
SampledSpectrum f(Vector3f woRender, Vector3f wiRender,
TransportMode mode = TransportMode::Radiance) const {
Vector3f wi = RenderToLocal(wiRender), wo = RenderToLocal(woRender);
if (wo.z == 0)
return {};
return bxdf.f(wo, wi, mode);
}
template <typename BxDF>
PBRT_CPU_GPU SampledSpectrum f(Vector3f woRender, Vector3f wiRender,
TransportMode mode = TransportMode::Radiance) const {
Vector3f wi = RenderToLocal(wiRender), wo = RenderToLocal(woRender);
if (wo.z == 0)
return {};
const BxDF *specificBxDF = bxdf.CastOrNullptr<BxDF>();
return specificBxDF->f(wo, wi, mode);
}
BSDF的Sample_f会通过std::optional处理采样失效的情况.
PBRT_CPU_GPU
pstd::optional<BSDFSample> Sample_f(
Vector3f woRender, Float u, Point2f u2,
TransportMode mode = TransportMode::Radiance,
BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const {
Vector3f wo = RenderToLocal(woRender);
if (wo.z == 0 || !(bxdf.Flags() & sampleFlags))
return {};
// Sample _bxdf_ and return _BSDFSample_
pstd::optional<BSDFSample> bs = bxdf.Sample_f(wo, u, u2, mode, sampleFlags);
if (bs)
DCHECK_GE(bs->pdf, 0);
if (!bs || !bs->f || bs->pdf == 0 || bs->wi.z == 0)
return {};
PBRT_DBG("For wo = (%f, %f, %f), ns %f %f %f sampled f = %f %f %f %f, pdf = %f, "
"ratio[0] = %f wi = (%f, %f, %f)\n",
wo.x, wo.y, wo.z, shadingFrame.z.x, shadingFrame.z.y, shadingFrame.z.z,
bs->f[0], bs->f[1], bs->f[2], bs->f[3], bs->pdf,
(bs->pdf > 0) ? (bs->f[0] / bs->pdf) : 0, bs->wi.x, bs->wi.y, bs->wi.z);
bs->wi = LocalToRender(bs->wi);
return bs;
}
漫反射
#
Lambertian反射模型是最简单的BRDF之一, 它表达了将入射光线均匀的反射到各个方向的完美漫反射. 大部分角度下Lambertian模型都可以很好的表达漫反射, 但是在与法线接近垂直的掠射角下镜面反射会导致明显的偏差.
Lambertian反射的表达式如下, 其中的π \pi π π \pi π
∫ ω f r ( p , ω o , ω i ) cos θ d ω = ∫ 0 2 π ∫ 0 π 2 R π cos θ sin θ d θ d ϕ = R
\begin{equation}
\int_\omega f_r(p, \omega_o, \omega_i) \cos\theta \mathrm{d}\omega = \int_0^{2\pi}\int_0^{\frac{\pi}{2}} \frac{R}{\pi} \cos\theta \sin\theta \mathrm{d}\theta \mathrm{d}\phi = R
\end{equation}
∫ ω f r ( p , ω o , ω i ) cos θ d ω = ∫ 0 2 π ∫ 0 2 π π R cos θ sin θ d θ d ϕ = R pbrt在采样Lambertian时会将cos \cos cos
// Sample cosine-weighted hemisphere to compute _wi_ and _pdf_
Vector3f wi = SampleCosineHemisphere(u);
if (wo.z < 0)
wi.z *= -1;
Float pdf = CosineHemispherePDF(AbsCosTheta(wi));
圆盘采样
#
对极坐标r , θ r,\theta r , θ
为了保证样本点在面积上均匀分布, 需要满足p ( x , y ) = 1 π p(x,y)=\frac{1}{\pi} p ( x , y ) = π 1 p ( r , θ ) = r π p(r,\theta)=\frac{r}{\pi} p ( r , θ ) = π r
p ( r ) = ∫ 0 2 π p ( r , θ ) d θ = 2 r p ( θ ∣ r ) = p ( r , θ ) p ( r ) = 1 2 π
\begin{equation}
\begin{aligned}
p(r)&=\int_0^{2\pi}p(r,\theta)\mathrm{d}\theta=2r\\
p(\theta|r)&=\frac{p(r,\theta)}{p(r)}=\frac{1}{2\pi}
\end{aligned}
\end{equation}
p ( r ) p ( θ ∣ r ) = ∫ 0 2 π p ( r , θ ) d θ = 2 r = p ( r ) p ( r , θ ) = 2 π 1 根据重要性抽样得到的结果如下, 此时根据Jacobian可得d r d θ = π x d x d y drd\theta=\frac{\pi}{\sqrt{x}}dxdy d r d θ = x π d x d y
r = ϵ 1 θ = 2 π ϵ 2
\begin{equation}
\begin{aligned}
r&=\sqrt{\epsilon_1}\\
\theta&=2\pi\epsilon_2
\end{aligned}
\end{equation}
r θ = ϵ 1 = 2 π ϵ 2 Shirley-Chiu方法可以解决这一问题, 根据Cartesian轴及y = x , y = − x y=x,y=-x y = x , y = − x [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] max ( x , y \max(x,y max ( x , y a a a b b b a b 2 \frac{ab}{2} 2 ab θ r 2 2 \frac{\theta r^2}{2} 2 θ r 2 π a b 8 \frac{\pi ab}{8} 8 πab
r = { x ∣ x ∣ > ∣ y ∣ y ∣ x ∣ ≤ ∣ y ∣ θ = { π 4 y x ∣ x ∣ > ∣ y ∣ π 2 − x y ∣ x ∣ ≤ ∣ y ∣
\begin{equation}
\begin{aligned}
r&=
\begin{cases}
x & |x|>|y|\\
y & |x| \le |y|
\end{cases}\\
\theta&=
\begin{cases}
\frac{\pi}{4}\frac{y}{x} & |x| > |y|\\
\frac{\pi}{2} - \frac{x}{y} & |x| \le |y|
\end{cases}
\end{aligned}
\end{equation}
r θ = { x y ∣ x ∣ > ∣ y ∣ ∣ x ∣ ≤ ∣ y ∣ = { 4 π x y 2 π − y x ∣ x ∣ > ∣ y ∣ ∣ x ∣ ≤ ∣ y ∣ 此时Jacobian行列式满足如下关系, 因此p ( r , θ ) = ∣ J ∣ − 1 p ( u , v ) = ∣ J ∣ − 1 1 4 = r π p(r,\theta)=|J|^{-1}p(u,v)=|J|^{-1}\frac{1}{4}=\frac{r}{\pi} p ( r , θ ) = ∣ J ∣ − 1 p ( u , v ) = ∣ J ∣ − 1 4 1 = π r
∣ J ∣ = { π 4 x = π 4 r ∣ x ∣ > ∣ y ∣ π 4 y = π 4 r ∣ x ∣ ≤ ∣ y ∣
\begin{equation}
\begin{aligned}
|J|&=
\begin{cases}
\frac{\pi}{4x}=\frac{\pi}{4r} & |x|>|y|\\
\frac{\pi}{4y}=\frac{\pi}{4r} & |x|\le|y|
\end{cases}
\end{aligned}
\end{equation}
∣ J ∣ = { 4 x π = 4 r π 4 y π = 4 r π ∣ x ∣ > ∣ y ∣ ∣ x ∣ ≤ ∣ y ∣ 余弦加权半球采样
#
渲染方程中有余弦项, 将它与BSDF一起做重要性抽样可以有效的考虑余弦的影响, 例如减小与法线夹角较大的光线被采样到的概率. 令p ( ω ) = c cos θ p(\omega)=c\cos\theta p ( ω ) = c cos θ
∫ Θ p ( ω ) d ω = 1 ∫ 0 2 π ∫ 0 π 2 c cos θ sin θ d θ d ϕ = 1 c = 1 π p ( θ , ϕ ) = 1 π cos θ sin θ
\begin{equation}
\begin{aligned}
&\int_\Theta p(\omega)\mathrm{d}\omega=1\\
&\int_0^{2\pi}\int_0^{\frac{\pi}{2}}c\cos\theta\sin\theta \mathrm{d}\theta \mathrm{d}\phi=1\\\
&c=\frac{1}{\pi}\\
&p(\theta,\phi)=\frac{1}{\pi}\cos\theta\sin\theta
\end{aligned}
\end{equation}
∫ Θ p ( ω ) d ω = 1 ∫ 0 2 π ∫ 0 2 π c cos θ sin θ d θ d ϕ = 1 c = π 1 p ( θ , ϕ ) = π 1 cos θ sin θ Malley方法通过将均匀圆盘采样的结果投影到半球上实现重要性抽样. 由于此时r = sin θ r=\sin\theta r = sin θ r , ϕ r,\phi r , ϕ θ , ϕ \theta,\phi θ , ϕ p ( θ , ϕ ) = p ( r , ϕ ) cos θ = sin θ cos θ π p(\theta,\phi)=p(r,\phi)\cos\theta=\frac{\sin\theta\cos\theta}{\pi} p ( θ , ϕ ) = p ( r , ϕ ) cos θ = π s i n θ c o s θ
∣ J ∣ = ∣ ∂ r ∂ θ ∂ r ∂ ϕ ∂ ϕ ∂ θ ∂ ϕ ∂ ϕ ∣ = ∣ cos θ 0 0 1 ∣ = cos θ
\begin{equation}
\begin{aligned}
|J|=\begin{vmatrix}
\frac{\partial r}{\partial \theta} & \frac{\partial r}{\partial \phi}\\
\frac{\partial \phi}{\partial \theta} & \frac{\partial \phi}{\partial \phi}
\end{vmatrix}
=\begin{vmatrix}
\cos\theta & 0\\
0 & 1
\end{vmatrix}=\cos\theta
\end{aligned}
\end{equation}
∣ J ∣ = ∂ θ ∂ r ∂ θ ∂ ϕ ∂ ϕ ∂ r ∂ ϕ ∂ ϕ = cos θ 0 0 1 = cos θ 镜面反射与透射
#
镜面反射只会将光线反射到某个方向集合中, 本节主要关注完美镜面反射, 微表面理论对反射结果的影响会在后面讨论.
物理原理
#
pbrt主要关注几何光学, 在波长远小于物体尺寸时这是可行的, 但对光的本质的探究也是必要的. 根据电磁学的定义光是电磁场的振荡, pbrt只关注其中电场的部分即电子的运动, 当光与表面接触时会刺激表面上的电子使其剧烈振荡, 移动的电子会导致电场中振荡的叠加.
依据电磁学理论, 物质可以分为绝缘体, 导体与半导体. 绝缘体中电子不会脱离原子; 导体中电子可以自由移动, 但移动过程中会有能量衰减, 通常完全吸收会发生在距离表面小于0.1微米处; 半导体兼顾二者的特性, 例如硅在可见光内具有金属特性, 在红外光中则具有透明特性.
折射率
#
由光所激发的电子产生的振荡通常会比原本的振荡的速度要小, 速度减小的程度被称为折射率(index of refraction, IOR), 它通过真空光速与当前材质中的光速的比值来定义, 通常在1.0~2.5之间, 同时光的波长也会影响IOR. IOR的突变会导致反射率的增加, IOR的变化是我们能观察到物体的原因.
镜面反射法则
#
令入射光线方向为ω i \omega_i ω i
ω r = − ω i + 2 ( n ⋅ ω i ) n
\begin{equation}
\omega_r = -\omega_i + 2(\mathbf{n} \cdot \omega_i)\mathbf{n}
\end{equation}
ω r = − ω i + 2 ( n ⋅ ω i ) n Snell定律
#
Snell定律描述了入射光线与折射光线在方向上的关系, 这可以通过Fermat原理推导得到, 即光的传播路径是耗时最小的路径, 通过求极值即可证明. 从Snell定律可以看出折射方向取决于IOR的比值即η = η i η t \eta=\frac{\eta_i}{\eta_t} η = η t η i
η i sin θ i = η t sin θ t ϕ t = ϕ i + π
\begin{equation}
\begin{aligned}
\eta_i \sin\theta_i &= \eta_t \sin\theta_t\\
\phi_t &= \phi_i + \pi
\end{aligned}
\end{equation}
η i sin θ i ϕ t = η t sin θ t = ϕ i + π 经过不同介质后, 由于光的不同波长的部分具有不同的IOR, 这导致传播方向不同, 这被称为散射, 通常会表现为经过介质后产生彩虹状的光锥.
利用Snell定律可以按如下方式表示折射光线. 当光线从物体内部发出时, 由于法线朝外, ω i ⋅ n \omega_i \cdot \mathbf{n} ω i ⋅ n
ω t = ω i η + ( cos θ i η − cos θ t ) n
\begin{equation}
\omega_t = \omega_i\eta + (\cos\theta_i\eta - \cos\theta_t) \mathbf{n}
\end{equation}
ω t = ω i η + ( cos θ i η − cos θ t ) n 当光线穿过光学密度更小(折射后的介质的IOR较小)时, 当入射天顶角超过θ c = sin − 1 ( η − 1 ) \theta_c = \sin^{-1}(\eta^{-1}) θ c = sin − 1 ( η − 1 ) 90 ∘ 90^{\circ} 9 0 ∘
Fresnel方程
#
Fresnel方程描述了光线折射与反射的量, 它是Maxwell方程在光滑平面上的解.
将光线分解为相对于表面的垂直与水平偏振, 它们因反射而产生的振幅的变化是不同的. 若θ t > π 2 \theta_t>\frac{\pi}{2} θ t > 2 π 1 1 1
r ∣ ∣ = E r ∣ ∣ E i ∣ ∣ = η t cos θ i − η i cos θ t η t cos θ i + η i cos θ t = cos θ i − η cos θ t cos θ i + η cos θ t r ⊥ = E r ⊥ E i ⊥ = η i cos θ i − η t cos θ t η i cos θ i + η t cos θ t = η cos θ i − cos θ t η cos θ i + cos θ t
\begin{equation}
\begin{aligned}
r_{||} &= \frac{E_r^{||}}{E_i^{||}} &= \frac{\eta_t \cos\theta_i - \eta_i \cos\theta_t}{\eta_t \cos\theta_i + \eta_i \cos\theta_t} &= \frac{\cos\theta_i - \eta \cos\theta_t}{\cos\theta_i + \eta \cos\theta_t}\\
r_{\perp} &= \frac{E_r^{\perp}}{E_i^{\perp}} &= \frac{\eta_i \cos\theta_i - \eta_t \cos\theta_t}{\eta_i \cos\theta_i + \eta_t \cos\theta_t} &= \frac{\eta \cos\theta_i - \cos\theta_t}{\eta \cos\theta_i + \cos\theta_t}
\end{aligned}
\end{equation}
r ∣∣ r ⊥ = E i ∣∣ E r ∣∣ = E i ⊥ E r ⊥ = η t cos θ i + η i cos θ t η t cos θ i − η i cos θ t = η i cos θ i + η t cos θ t η i cos θ i − η t cos θ t = cos θ i + η cos θ t cos θ i − η cos θ t = η cos θ i + cos θ t η cos θ i − cos θ t 电磁学关注的是反射波的振幅与相位, 而pbrt的几何光学中更关注光线所携带的功率, 这可以通过振幅的平方来表示. 若入射光是无偏振的即含有等量的水平与垂直偏振, 此时反射功率为水平与垂直反射功率的平均值, 即位Fresnel反射率.
F r = 1 2 ( r ∣ ∣ 2 + r ⊥ 2 )
\begin{equation}
F_r = \frac{1}{2}(r_{||}^2 + r_{\perp}^2)
\end{equation}
F r = 2 1 ( r ∣∣ 2 + r ⊥ 2 ) 实时渲染中通常使用Schlick近似.
F r = F r 0 + ( 1 − F r 0 ) ( 1 − cos θ i ) 5 F r 0 = ( η i − η t η i + η t ) 2
\begin{equation}
\begin{aligned}
F_r &= F_{r_0}+(1-F_{r_0})(1-\cos\theta_i)^5\\
F_{r_0} &= \left(\frac{\eta_i - \eta_t}{\eta_i + \eta_t}\right)^2
\end{aligned}
\end{equation}
F r F r 0 = F r 0 + ( 1 − F r 0 ) ( 1 − cos θ i ) 5 = ( η i + η t η i − η t ) 2 导体Fresnel方程
#
导体的IOR需要用复数表示, 实部描述光线速度的减小, 虚部描述光线在材质内传播时的衰减.
虚部为0 0 0 z z z α \alpha α
E ( z ) = e i α η z R ( E ( z ) ) = cos ( α η z )
\begin{equation}
\begin{aligned}
E(z)&=e^{i\alpha\eta z}\\
\mathcal{R}(E(z))&=\cos(\alpha\eta z)
\end{aligned}
\end{equation}
E ( z ) R ( E ( z )) = e i α ηz = cos ( α ηz ) 当使用复数η c = η + i k \eta_c=\eta+ik η c = η + ik
R ( E ( z ) ) = e − α z k cos ( α η z )
\begin{equation}
\begin{aligned}
\mathcal{R}(E(z))=e^{-\alpha zk}\cos(\alpha\eta z)
\end{aligned}
\end{equation}
R ( E ( z )) = e − α z k cos ( α ηz ) 复数的2范数为实部与虚部和的平方, 此时可泛化上述Fresnel反射率的计算.
∣ ∣ a + b i ∣ ∣ 2 = a 2 + b 2
\begin{equation}
\begin{aligned}
||a+bi||_2=\sqrt{a^2+b^2}
\end{aligned}
\end{equation}
∣∣ a + bi ∣ ∣ 2 = a 2 + b 2 导体BRDF
#
本节描述的是从光从绝缘体介质中接触导体表面时的BRDF, 只考虑光滑的情况. 光滑导体只将光线返回到一个反向, 反射光线的辐亮度由Fresnel反射率决定, 因此BRDF如下.
f r ( p , ω o , ω i ) = F r ( ω r ) δ ( ω i − ω r ) ∣ cos θ r ∣
\begin{equation}
f_r(p, \omega_o, \omega_i) = F_r(\omega_r) \frac{\delta(\omega_i - \omega_r)}{|\cos\theta_r|}
\end{equation}
f r ( p , ω o , ω i ) = F r ( ω r ) ∣ cos θ r ∣ δ ( ω i − ω r ) 绝缘体BSDF
#
pbrt在选择采样反射或透射时, 根据Fresnel反射率以及用户设置的Flags来决定. 由于只考虑光滑表面, BRDF和BTDF与上一节的导体BRDF都是相同的.
薄绝缘体BSDF
#
pbrt通过ThinDielectricBxDF模拟光线穿过薄绝缘体平行的内外表面的现象, 在薄导体内部多次反射后的总反射率如下.
R ′ = R + T R T + T R R R T + ⋯ = R + T 2 R 1 − R 2
\begin{equation}
R' = R + TRT + TRRRT + \cdots = R + \frac{T^2R}{1-R^2}
\end{equation}
R ′ = R + TRT + TRRRT + ⋯ = R + 1 − R 2 T 2 R 非对称散射与折射
#
BRDF都是对称的, 即入射与反射方向可以交换, 但BTDF不是, 因为IOR的变化会导致折射后光线范围的变化, 为保持能量守恒光线的辐亮度不再保持一致.
根据能量守恒我们可以得到如下等式.
d 2 ϕ o = d 2 ϕ i
\begin{equation}
d^2\phi_o = d^2\phi_i
\end{equation}
d 2 ϕ o = d 2 ϕ i 依据L = d 2 ϕ d ω d A ⊥ L=\frac{d^2\phi}{\mathrm{d}\omega \mathrm{d}A^{\perp}} L = d ω d A ⊥ d 2 ϕ
L o cos θ o d A d ω o = L i cos θ i d A d ω i
\begin{equation}
L_o \cos\theta_o \mathrm{d}A \mathrm{d}\omega_o = L_i \cos\theta_i \mathrm{d}A \mathrm{d}\omega_i
\end{equation}
L o cos θ o d A d ω o = L i cos θ i d A d ω i 微分Snell定律可得如下关系, 其中ϕ \phi ϕ d ϕ i = d ϕ o \mathrm{d}\phi_i=\mathrm{d}\phi_o d ϕ i = d ϕ o sin θ \sin\theta sin θ η \eta η
η o cos θ o d θ o = η i cos θ i d θ i η o cos θ o sin θ o d ω o d ϕ o = η i cos θ i sin θ i d ω i d ϕ i η o 2 cos θ o d ω o = η i 2 cos θ i d ω i
\begin{equation}
\begin{aligned}
\eta_o \cos\theta_o \mathrm{d}\theta_o
&= \eta_i \cos\theta_i \mathrm{d}\theta_i\\
\eta_o \frac{\cos\theta_o}{\sin\theta_o} \frac{\mathrm{d}\omega_o}{\mathrm{d}\phi_o}
&= \eta_i \frac{\cos\theta_i}{\sin\theta_i} \frac{\mathrm{d}\omega_i}{\mathrm{d}\phi_i}\\
\eta_o^2 \cos\theta_o \mathrm{d}\omega_o
&= \eta_i^2 \cos\theta_i \mathrm{d}\omega_i
\end{aligned}
\end{equation}
η o cos θ o d θ o η o sin θ o cos θ o d ϕ o d ω o η o 2 cos θ o d ω o = η i cos θ i d θ i = η i sin θ i cos θ i d ϕ i d ω i = η i 2 cos θ i d ω i 将微分后的Snell定律代入可得如下关系.
L o = η o 2 η i 2 L i
\begin{equation}
L_o = \frac{\eta_o^2}{\eta_i^2} L_i
\end{equation}
L o = η i 2 η o 2 L i 由于路径追踪实际上为光线传播的逆过程, 在折射时需要考虑非对称折射的现象.
f t ( p , ω i , ω o ) = L o 0 L i 0 cos θ i d ω o = η o 2 η i 2 cos θ i d ω i f t ( p , ω o , ω i ) = L i 1 L o 1 cos θ o d ω i = η i 2 η o 2 cos θ o d ω o f t ( p , ω i , ω o ) = f t ( p , ω o , ω i ) η 2
\begin{equation}
\begin{aligned}
f_t(p, \omega_i, \omega_o) &= \frac{L_o^0}{L_i^0 \cos\theta_i \mathrm{d}\omega_o} = \frac{\eta_o^2}{\eta_i^2 \cos\theta_i \mathrm{d}\omega_i}\\
f_t(p, \omega_o, \omega_i) &= \frac{L_i^1}{L_o^1 \cos\theta_o \mathrm{d}\omega_i} = \frac{\eta_i^2}{\eta_o^2 \cos\theta_o \mathrm{d}\omega_o}\\
f_t(p, \omega_i, \omega_o) &= \frac{f_t(p, \omega_o, \omega_i)}{\eta^2}
\end{aligned}
\end{equation}
f t ( p , ω i , ω o ) f t ( p , ω o , ω i ) f t ( p , ω i , ω o ) = L i 0 cos θ i d ω o L o 0 = η i 2 cos θ i d ω i η o 2 = L o 1 cos θ o d ω i L i 1 = η o 2 cos θ o d ω o η i 2 = η 2 f t ( p , ω o , ω i ) 微表面理论
#
大部分物体表面在微观上是粗糙的, 这可以通过微表面理论模拟, 这避免了对建模精度的要求, 且更适用于Monte Carlo.
法线分布函数
#
微表面分布主要定义了法线的分布, 由于光源与观察者都与表面相距较远, 这种概率表示是可行的. 宏表面与微表面所占的面积关系如下, 即微表面沿宏观法线投影的面积与宏表面面积相等.
∫ d A μ ( ω m ( p ) ⋅ n ) d p = ∫ d A d p
\begin{equation}
\int_{\mathrm{d}A_\mu}(\omega_m(p) \cdot \mathbf{n}) dp = \int_{\mathrm{d}A} dp
\end{equation}
∫ d A μ ( ω m ( p ) ⋅ n ) d p = ∫ d A d p 将所有具有某个方向的法线的微表面在宏表面上的相对投影面积定义为法线分布函数, 此时可以得到法线分布函数的归一化性质.
∫ Ω D ( ω m ) ( ω m ⋅ n ) d ω m = ∫ Ω D ( ω m ) cos θ m d ω m = 1
\begin{equation}
\int_\Omega D(\omega_m)(\omega_m \cdot \mathbf{n}) \mathrm{d}\omega_m = \int_\Omega D(\omega_m) \cos\theta_m \mathrm{d}\omega_m = 1
\end{equation}
∫ Ω D ( ω m ) ( ω m ⋅ n ) d ω m = ∫ Ω D ( ω m ) cos θ m d ω m = 1 若转为斜率分布, 可以得到如下关系.
∫ − ∞ ∞ ∫ − ∞ ∞ P 22 ( x m ~ , y m ~ ) d x m ~ d y m ~ = 1
\begin{equation}
\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} P^{22}(x_{\tilde{m}},y_{\tilde{m}}) \mathrm{d}x_{\tilde{m}} \mathrm{d}y_{\tilde{m}}=1
\end{equation}
∫ − ∞ ∞ ∫ − ∞ ∞ P 22 ( x m ~ , y m ~ ) d x m ~ d y m ~ = 1 根据平面的定义x m x + y m y + z m z + c = 0 x_mx+y_my+z_mz+c=0 x m x + y m y + z m z + c = 0
( x m ~ , y m ~ ) = ( − x m z m , − y m z m ) = − tan θ m ( cos ϕ m , sin ϕ m ) ( x m , y m , z m ) = ( − x m ~ , − y m ~ , 1 ) x m ~ 2 + y m ~ 2 + 1 = ( − x m ~ , − y m ~ , 1 ) cos θ m
\begin{equation}
\begin{aligned}
(x_{\tilde{m}},y_{\tilde{m}})&=(-\frac{x_m}{z_m},-\frac{y_m}{z_m})=-\tan\theta_m(\cos\phi_m,\sin\phi_m)\\
(x_m, y_m, z_m)&=\frac{(-x_{\tilde{m}},-y_{\tilde{m}},1)}{\sqrt{x_{\tilde{m}}^2+y_{\tilde{m}}^2+1}}=(-x_{\tilde{m}},-y_{\tilde{m}},1)\cos\theta_m
\end{aligned}
\end{equation}
( x m ~ , y m ~ ) ( x m , y m , z m ) = ( − z m x m , − z m y m ) = − tan θ m ( cos ϕ m , sin ϕ m ) = x m ~ 2 + y m ~ 2 + 1 ( − x m ~ , − y m ~ , 1 ) = ( − x m ~ , − y m ~ , 1 ) cos θ m Jacobi行列式的计算结果如下.
d x m ~ d y m ~ = sin θ cos 3 θ d θ d ϕ
\begin{equation}
\mathrm{d}x_{\tilde{m}}\mathrm{d}y_{\tilde{m}}=\frac{\sin\theta}{\cos^3\theta}\mathrm{d}\theta \mathrm{d}\phi
\end{equation}
d x m ~ d y m ~ = cos 3 θ sin θ d θ d ϕ 此时可以得到法线分布函数与斜率分布函数的转换关系.
D ( ω m ) = P 22 ( x m ~ , y m ~ ) d x m ~ d y m ~ cos θ sin θ d θ d ϕ = P 22 ( x m ~ , y m ~ ) cos 4 θ
\begin{equation}
D(\omega_m) = \frac{P^{22}(x_{\tilde{m}},y_{\tilde{m}}) \mathrm{d}x_{\tilde{m}} \mathrm{d}y_{\tilde{m}}}{\cos\theta \sin\theta \mathrm{d}\theta \mathrm{d}\phi} = \frac{P^{22}(x_{\tilde{m}},y_{\tilde{m}})}{\cos^4\theta}
\end{equation}
D ( ω m ) = cos θ sin θ d θ d ϕ P 22 ( x m ~ , y m ~ ) d x m ~ d y m ~ = cos 4 θ P 22 ( x m ~ , y m ~ ) 若法线分布函数为形状无关的, 即修改粗糙度等价于对法线分布函数进行拉伸, 则其斜率分布需要满足如下形式.
P 22 ( x m ~ , y m ~ ) = 1 α x α y f ( x m ~ 2 α x 2 + y m ~ 2 α y 2 ) = 1 α x α y f ( tan θ m cos 2 ϕ m α x 2 + sin 2 ϕ m α y 2 )
\begin{equation}
\begin{aligned}
P^{22}(x_{\tilde{m}},y_{\tilde{m}})
&=\frac{1}{\alpha_x\alpha_y}f(\sqrt{\frac{x_{\tilde{m}^2}}{\alpha_x^2} + \frac{y_{\tilde{m}^2}}{\alpha_y^2}})\\
&=\frac{1}{\alpha_x\alpha_y}f(\tan\theta_m\sqrt{\frac{\cos^2\phi_m}{\alpha_x^2} + \frac{\sin^2\phi_m}{\alpha_y^2}})
\end{aligned}
\end{equation}
P 22 ( x m ~ , y m ~ ) = α x α y 1 f ( α x 2 x m ~ 2 + α y 2 y m ~ 2 ) = α x α y 1 f ( tan θ m α x 2 cos 2 ϕ m + α y 2 sin 2 ϕ m ) 对于渲染任务最常用的Trowbridge-Reitz法线分布函数定义如下.
P 22 ( x m ~ , y m ~ ) = 1 π α x α y ( 1 + x m ~ 2 α x 2 + y m ~ 2 α y 2 ) 2 D ( ω m ) = 1 π α x α y cos 4 θ m ( 1 + tan 2 θ m ( cos 2 ϕ m α x 2 + sin 2 ϕ m α y 2 ) ) 2
\begin{equation}
\begin{aligned}
P^{22}(x_{\tilde{m}},y_{\tilde{m}}) &= \frac{1}{\pi\alpha_x\alpha_y(1+\frac{x_{\tilde{m}^2}}{\alpha_x^2} + \frac{y_{\tilde{m}^2}}{\alpha_y^2})^2}\\
D(\omega_m) &= \frac{1}{\pi\alpha_x\alpha_y\cos^4\theta_m(1+\tan^2\theta_m(\frac{\cos^2\phi_m}{\alpha_x^2} + \frac{\sin^2\phi_m}{\alpha_y^2}))^2}
\end{aligned}
\end{equation}
P 22 ( x m ~ , y m ~ ) D ( ω m ) = π α x α y ( 1 + α x 2 x m ~ 2 + α y 2 y m ~ 2 ) 2 1 = π α x α y cos 4 θ m ( 1 + tan 2 θ m ( α x 2 c o s 2 ϕ m + α y 2 s i n 2 ϕ m ) ) 2 1 遮蔽函数
#
微表面可能会被遮挡, 需要通过遮蔽函数来表示某个方向上被遮蔽的比率.
∫ Ω D ( ω m ) G 1 ( ω , ω m ) ( ω ⋅ ω m ) d ω m = ω ⋅ n = cos θ
\begin{equation}
\int_\Omega D(\omega_m)G_1(\omega,\omega_m)(\omega\cdot\omega_m)\mathrm{d}\omega_m=\omega\cdot\mathbf{n}=\cos\theta
\end{equation}
∫ Ω D ( ω m ) G 1 ( ω , ω m ) ( ω ⋅ ω m ) d ω m = ω ⋅ n = cos θ Smith遮蔽函数认为法线与高度是独立的, 即微表面类似于一个一个浮动的小片段, 此时可以将遮蔽函数分解为局部遮蔽与全部遮蔽, 局部遮蔽只需要法线可见, 即ω ⋅ ω m > 0 \omega\cdot\omega_m>0 ω ⋅ ω m > 0 ( 0 , 0 , 1 ) (0,0,1) ( 0 , 0 , 1 ) x x x ( sin θ , 0 , cos θ ) (\sin\theta, 0, \cos\theta) ( sin θ , 0 , cos θ ) ∫ − ∞ ∞ x m ~ P 2 − ( x m ~ ) d x m ~ = 0 \int_{-\infty}^{\infty}x_{\tilde{m}}P^{2-}(x_{\tilde{m}})\mathrm{d}x_{\tilde{m}}=0 ∫ − ∞ ∞ x m ~ P 2 − ( x m ~ ) d x m ~ = 0 ∫ − ∞ ∞ P 2 − ( x m ~ ) d x m ~ = 1 \int_{-\infty}^{\infty}P^{2-}(x_{\tilde{m}})\mathrm{d}x_{\tilde{m}}=1 ∫ − ∞ ∞ P 2 − ( x m ~ ) d x m ~ = 1
G 1 ( ω ) = cos θ ∫ Ω D ( ω m ) max ( 0 , ω ⋅ ω m ) d ω m = cos θ ∫ Ω D ( ω m ) cos θ m max ( 0 , ω ⋅ ω m ) cos θ m d ω m = cos θ ∫ − ∞ ∞ ∫ − ∞ ∞ P 22 ( x m ~ , y m ~ ) max ( 0 , ( − x m ~ sin θ + cos θ ) cos θ m ) cos θ m d x m ~ d y m ~ = cos θ ∫ − ∞ ∞ ∫ − ∞ cot θ P 22 ( x m ~ , y m ~ ) ( − x m ~ sin θ + cos θ ) d x m ~ d y m ~ = cot θ ∫ − ∞ cot θ P 2 − ( x m ~ ) ( − x m ~ + cot θ ) d x m ~ = cot θ ∫ − ∞ ∞ P 2 − ( x m ~ ) x m ~ d x m ~ − cot θ ∫ − ∞ ∞ P 2 − ( x m ~ ) d x m ~ + cot θ + ∫ − ∞ cot θ P 2 − ( x m ~ ) ( − x m ~ + cot θ ) d x m ~ = 1 1 + 1 cot θ ∫ cot θ ∞ P 2 − ( x m ~ ) ( x m ~ − cot θ ) d x m ~ = 1 1 + Λ ( ω )
\begin{equation}
\begin{aligned}
G_1(\omega)
&=\frac{\cos\theta}{\int_\Omega D(\omega_m)\max(0,\omega\cdot\omega_m)\mathrm{d}\omega_m}\\
&=\frac{\cos\theta}{\int_\Omega D(\omega_m)\cos\theta_m\frac{\max(0,\omega\cdot\omega_m)}{\cos\theta_m}\mathrm{d}\omega_m}\\
&=\frac{\cos\theta}{\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}P^{22}(x_{\tilde{m}},y_{\tilde{m}}) \frac{\max(0, (-x_{\tilde{m}}\sin\theta + \cos\theta)\cos\theta_m)}{\cos\theta_m} \mathrm{d}x_{\tilde{m}} \mathrm{d}y_{\tilde{m}}}\\
&=\frac{\cos\theta}{\int_{-\infty}^{\infty}\int_{-\infty}^{\cot\theta}P^{22}(x_{\tilde{m}},y_{\tilde{m}}) (-x_{\tilde{m}}\sin\theta + \cos\theta) \mathrm{d}x_{\tilde{m}} \mathrm{d}y_{\tilde{m}}}\\
&=\frac{\cot\theta}{\int_{-\infty}^{\cot\theta}P^{2-}(x_{\tilde{m}}) (-x_{\tilde{m}} + \cot\theta) \mathrm{d}x_{\tilde{m}}}\\
&=\frac{\cot\theta}{\int_{-\infty}^{\infty}P^{2-}(x_{\tilde{m}}) x_{\tilde{m}} \mathrm{d}x_{\tilde{m}}-\cot\theta\int_{-\infty}^{\infty}P^{2-}(x_{\tilde{m}}) \mathrm{d}x_{\tilde{m}}+\cot\theta+\int_{-\infty}^{\cot\theta}P^{2-}(x_{\tilde{m}}) (-x_{\tilde{m}} + \cot\theta) \mathrm{d}x_{\tilde{m}}}\\
&=\frac{1}{1+\frac{1}{\cot\theta}\int_{\cot\theta}^{\infty}P^{2-}(x_{\tilde{m}}) (x_{\tilde{m}} - \cot\theta) \mathrm{d}x_{\tilde{m}}}\\
&=\frac{1}{1+\Lambda(\omega)}
\end{aligned}
\end{equation}
G 1 ( ω ) = ∫ Ω D ( ω m ) max ( 0 , ω ⋅ ω m ) d ω m cos θ = ∫ Ω D ( ω m ) cos θ m c o s θ m m a x ( 0 , ω ⋅ ω m ) d ω m cos θ = ∫ − ∞ ∞ ∫ − ∞ ∞ P 22 ( x m ~ , y m ~ ) c o s θ m m a x ( 0 , ( − x m ~ s i n θ + c o s θ ) c o s θ m ) d x m ~ d y m ~ cos θ = ∫ − ∞ ∞ ∫ − ∞ c o t θ P 22 ( x m ~ , y m ~ ) ( − x m ~ sin θ + cos θ ) d x m ~ d y m ~ cos θ = ∫ − ∞ c o t θ P 2 − ( x m ~ ) ( − x m ~ + cot θ ) d x m ~ cot θ = ∫ − ∞ ∞ P 2 − ( x m ~ ) x m ~ d x m ~ − cot θ ∫ − ∞ ∞ P 2 − ( x m ~ ) d x m ~ + cot θ + ∫ − ∞ c o t θ P 2 − ( x m ~ ) ( − x m ~ + cot θ ) d x m ~ cot θ = 1 + c o t θ 1 ∫ c o t θ ∞ P 2 − ( x m ~ ) ( x m ~ − cot θ ) d x m ~ 1 = 1 + Λ ( ω ) 1 各向同性的Trowbridge-Reitz的Λ ( ω ) \Lambda(\omega) Λ ( ω ) α = α x 2 cos 2 ϕ + α y 2 sin 2 ϕ \alpha=\sqrt{\alpha_x^2\cos^2\phi+\alpha_y^2\sin^2\phi} α = α x 2 cos 2 ϕ + α y 2 sin 2 ϕ
Λ ( ω ) = 1 cot θ ∫ − ∞ ∞ ∫ cot θ ∞ 1 π α 2 ( 1 + x m ~ 2 α 2 + y m ~ 2 α 2 ) 2 ( x m ~ − cot θ ) d x m ~ d y m ~ = 1 cot θ ∫ − ∞ ∞ ∫ cot θ α ∞ α π ( 1 + ( x m ~ α ) 2 + ( y m ~ α ) 2 ) 2 ( x m ~ α − cot θ α ) d x m ~ α d y m ~ α = 1 π x 0 ∫ − ∞ ∞ ∫ x 0 ∞ 1 ( 1 + x 2 + y 2 ) 2 ( x − x 0 ) d x d y = 1 π x 0 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 3 2 d x ∫ − ∞ ∞ 1 ( 1 + ( y 1 + x 2 ) 2 ) 2 d y 1 + x 2 = 1 π x 0 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 3 2 d x ∫ − ∞ ∞ 1 ( 1 + y 2 ) 2 d y = 1 π x 0 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 3 2 d x ∫ − ∞ ∞ 1 ( 1 + tan 2 θ ) 2 d tan θ = 1 π x 0 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 3 2 d x ∫ − π 2 π 2 cos 2 θ + 1 2 d θ = 1 2 x 0 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 3 2 d x = 1 2 x 0 ∫ x 0 ∞ ( x − x 0 ) ( − 1 x ) d ( 1 + x 2 ) − 1 2 = 1 2 x 0 ∫ ( 1 + x 0 2 ) − 1 2 0 x 0 1 x 2 − 1 − 1 d x = 1 2 ∫ 0 ( 1 + x 0 2 ) − 1 2 1 x 0 − x 1 − x 2 d x = 1 2 ( ( 1 + x 0 2 ) − 1 2 x 0 + ∫ 0 ( 1 + x 0 2 ) − 1 2 d 1 − x 2 ) = 1 2 ( ( 1 + x 0 ) − 1 2 x 0 + x 0 ( 1 + x 0 ) − 1 2 − 1 ) = 1 + 1 x 0 2 − 1 2 = 1 + α 2 tan 2 θ − 1 2
\begin{equation}
\begin{aligned}
\Lambda(\omega)
&=\frac{1}{\cot\theta}\int_{-\infty}^{\infty}\int_{\cot\theta}^{\infty}\frac{1}{\pi\alpha^2(1+\frac{x_{\tilde{m}^2}}{\alpha^2} + \frac{y_{\tilde{m}^2}}{\alpha^2})^2} (x_{\tilde{m}} - \cot\theta)\mathrm{d}x_{\tilde{m}}\mathrm{d}y_{\tilde{m}}\\
&=\frac{1}{\cot\theta}\int_{-\infty}^{\infty}\int_{\frac{\cot\theta}{\alpha}}^{\infty}\frac{\alpha}{\pi(1+(\frac{x_{\tilde{m}}}{\alpha})^2 + (\frac{y_{\tilde{m}}}{\alpha})^2)^2} (\frac{x_{\tilde{m}}}{\alpha} - \frac{\cot\theta}{\alpha})d\frac{x_{\tilde{m}}}{\alpha}d\frac{y_{\tilde{m}}}{\alpha}\\
&=\frac{1}{\pi x_0}\int_{-\infty}^{\infty}\int_{x_0}^{\infty}\frac{1}{(1+x^2+y^2)^2}(x-x_0)dxdy\\
&=\frac{1}{\pi x_0}\int_{x_0}^{\infty}(x-x_0)(1+x^2)^{-\frac{3}{2}}\mathrm{d}x\int_{-\infty}^{\infty}\frac{1}{(1+(\frac{y}{\sqrt{1+x^2}})^2)^2}d\frac{y}{\sqrt{1+x^2}}\\
&=\frac{1}{\pi x_0}\int_{x_0}^{\infty}(x-x_0)(1+x^2)^{-\frac{3}{2}}\mathrm{d}x\int_{-\infty}^{\infty}\frac{1}{(1+y^2)^2}\mathrm{d}y\\
&=\frac{1}{\pi x_0}\int_{x_0}^{\infty}(x-x_0)(1+x^2)^{-\frac{3}{2}}\mathrm{d}x\int_{-\infty}^{\infty}\frac{1}{(1+\tan^2\theta)^2}d\tan\theta\\
&=\frac{1}{\pi x_0}\int_{x_0}^{\infty}(x-x_0)(1+x^2)^{-\frac{3}{2}}\mathrm{d}x\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\frac{\cos2\theta+1}{2}\mathrm{d}\theta\\
&=\frac{1}{2x_0}\int_{x_0}^{\infty}(x-x_0)(1+x^2)^{-\frac{3}{2}}\mathrm{d}x\\
&=\frac{1}{2 x_0}\int_{x_0}^{\infty}(x-x_0)(-\frac{1}{x})d(1+x^2)^{-\frac{1}{2}}\\
&=\frac{1}{2 x_0}\int_{(1+x_0^2)^{-\frac{1}{2}}}^{0}\frac{x_0}{\sqrt{\frac{1}{x^2}-1}}-1\mathrm{d}x\\
&=\frac{1}{2}\int_{0}^{(1+x_0^2)^{-\frac{1}{2}}}\frac{1}{x_0}-\frac{x}{\sqrt{1-x^2}}\mathrm{d}x\\
&=\frac{1}{2}(\frac{(1+x_0^2)^{-\frac{1}{2}}}{x_0}+\int_{0}^{(1+x_0^2)^{-\frac{1}{2}}}d\sqrt{1-x^2})\\
&=\frac{1}{2}(\frac{(1+x_0)^{-\frac{1}{2}}}{x_0}+x_0(1+x_0)^{-\frac{1}{2}}-1)\\
&=\frac{\sqrt{1+\frac{1}{x_0^2}}-1}{2}\\
&=\frac{\sqrt{1+\alpha^2 \tan^2\theta}-1}{2}
\end{aligned}
\end{equation}
Λ ( ω ) = cot θ 1 ∫ − ∞ ∞ ∫ c o t θ ∞ π α 2 ( 1 + α 2 x m ~ 2 + α 2 y m ~ 2 ) 2 1 ( x m ~ − cot θ ) d x m ~ d y m ~ = cot θ 1 ∫ − ∞ ∞ ∫ α c o t θ ∞ π ( 1 + ( α x m ~ ) 2 + ( α y m ~ ) 2 ) 2 α ( α x m ~ − α cot θ ) d α x m ~ d α y m ~ = π x 0 1 ∫ − ∞ ∞ ∫ x 0 ∞ ( 1 + x 2 + y 2 ) 2 1 ( x − x 0 ) d x d y = π x 0 1 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 2 3 d x ∫ − ∞ ∞ ( 1 + ( 1 + x 2 y ) 2 ) 2 1 d 1 + x 2 y = π x 0 1 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 2 3 d x ∫ − ∞ ∞ ( 1 + y 2 ) 2 1 d y = π x 0 1 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 2 3 d x ∫ − ∞ ∞ ( 1 + tan 2 θ ) 2 1 d tan θ = π x 0 1 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 2 3 d x ∫ − 2 π 2 π 2 cos 2 θ + 1 d θ = 2 x 0 1 ∫ x 0 ∞ ( x − x 0 ) ( 1 + x 2 ) − 2 3 d x = 2 x 0 1 ∫ x 0 ∞ ( x − x 0 ) ( − x 1 ) d ( 1 + x 2 ) − 2 1 = 2 x 0 1 ∫ ( 1 + x 0 2 ) − 2 1 0 x 2 1 − 1 x 0 − 1 d x = 2 1 ∫ 0 ( 1 + x 0 2 ) − 2 1 x 0 1 − 1 − x 2 x d x = 2 1 ( x 0 ( 1 + x 0 2 ) − 2 1 + ∫ 0 ( 1 + x 0 2 ) − 2 1 d 1 − x 2 ) = 2 1 ( x 0 ( 1 + x 0 ) − 2 1 + x 0 ( 1 + x 0 ) − 2 1 − 1 ) = 2 1 + x 0 2 1 − 1 = 2 1 + α 2 tan 2 θ − 1 遮蔽阴影函数
#
除遮挡视线外, 光源也可能被微表面遮挡. 若采用G ( ω o , ω i ) = G 1 ( ω o ) G 1 ( ω i ) G(\omega_o, \omega_i)=G_1(\omega_o)G_1(\omega_i) G ( ω o , ω i ) = G 1 ( ω o ) G 1 ( ω i ) ω o = ω i \omega_o=\omega_i ω o = ω i
不失一般性, 令视线沿x轴, S ( x ) S(x) S ( x ) x = 0 x=0 x = 0 x x x g ( x ) d x g(x)\mathrm{d}x g ( x ) d x [ x , x + d x ] [x,x+\mathrm{d}x] [ x , x + d x ] S ( x + d x ) S(x + \mathrm{d}x) S ( x + d x )
S ( x + d x ) = S ( x ) ( 1 − g ( x ) d x )
\begin{equation}
S(x+\mathrm{d}x)=S(x)(1-g(x)\mathrm{d}x)
\end{equation}
S ( x + d x ) = S ( x ) ( 1 − g ( x ) d x ) 将该式视为一阶Taylor展开, 此时可以得到微分方程, S ( x ) S(x) S ( x )
d S ( x ) d x = − S ( x ) g ( x ) S ( x ) = S ( 0 ) e − ∫ 0 x g ( x ) d x
\begin{equation}
\begin{aligned}
\frac{\mathrm{d}S(x)}{\mathrm{d}x}&=-S(x)g(x)\\
S(x)&=S(0)e^{-\int_0^x g(x)\mathrm{d}x}
\end{aligned}
\end{equation}
d x d S ( x ) S ( x ) = − S ( x ) g ( x ) = S ( 0 ) e − ∫ 0 x g ( x ) d x 令视线在x x x μ v \mu_v μ v μ ( x ) \mu(x) μ ( x ) S ( 0 ) S(0) S ( 0 ) s ( μ v − μ ( 0 ) ) s(\mu_v-\mu(0)) s ( μ v − μ ( 0 ))
S ( 0 ) = s ( μ v − μ ( 0 ) ) = { 1 μ v ≥ μ ( 0 ) 0 μ v < μ ( 0 )
\begin{equation}
\begin{aligned}
S(0)=s(\mu_v-\mu(0))=
\begin{cases}
1 &\mu_v\ge\mu(0)\\
0 &\mu_v<\mu(0)
\end{cases}
\end{aligned}
\end{equation}
S ( 0 ) = s ( μ v − μ ( 0 )) = { 1 0 μ v ≥ μ ( 0 ) μ v < μ ( 0 ) 令表面高度为h ( x ) h(x) h ( x ) α \alpha α h ( x ) < h ( 0 ) + μ v x h(x)<h(0)+\mu_v x h ( x ) < h ( 0 ) + μ v x β \beta β h ( x + d x ) > h ( 0 ) + μ v ⋅ ( x + d x ) h(x+\mathrm{d}x)>h(0)+\mu_v\cdot(x+\mathrm{d}x) h ( x + d x ) > h ( 0 ) + μ v ⋅ ( x + d x ) P 3 ( h , μ ∣ x ) P_3(h,\mu|x) P 3 ( h , μ ∣ x ) x x x [ x , x + d x ] [x,x+\mathrm{d}x] [ x , x + d x ] μ ( x ) > μ v \mu(x)>\mu_v μ ( x ) > μ v g ( x ) g(x) g ( x ) P ( α , β ) P(\alpha,\beta) P ( α , β )
g ( x ) = 1 d x P ( β ∣ α ) = 1 d x P ( α , β ) P ( α ) = 1 d x ∫ μ v ∞ ∫ h ( 0 ) + μ v ( x + d x ) − μ d x h ( 0 ) + μ v x P 3 ( h , μ ∣ x ) d h d μ ∫ − ∞ ∞ ∫ − ∞ h ( 0 ) + μ v x P 3 ( h , μ ∣ x ) d h d μ = 1 d x ∫ μ v ∞ ( μ − μ v ) d x P 3 ( h ( 0 ) + μ v x , μ ∣ x ) d μ ∫ − ∞ ∞ ∫ − ∞ h ( 0 ) + μ v x P 3 ( h , μ ∣ x ) d h d μ = ∫ μ v ∞ ( μ − μ v ) P h ( h ( 0 ) + μ v x ) P 2 − ( μ ) d μ ∫ − ∞ h ( 0 ) + μ v x P h ( h ) d h = Λ ( μ v ) μ v P h ( h ( 0 ) + μ v x ) f ( h ( 0 ) + μ v x )
\begin{equation}
\begin{aligned}
g(x)
&=\frac{1}{\mathrm{d}x}P(\beta|\alpha)\\
&=\frac{1}{\mathrm{d}x}\frac{P(\alpha,\beta)}{P(\alpha)}\\
&=\frac{1}{\mathrm{d}x}\frac{\int_{\mu_v}^\infty\int_{h(0)+\mu_v(x+\mathrm{d}x)-\mu \mathrm{d}x}^{h(0)+\mu_vx}P_3(h,\mu|x)dh \mathrm{d}\mu}{\int_{-\infty}^{\infty}\int_{-\infty}^{h(0)+\mu_vx}P_3(h,\mu|x)dhd\mu}\\
&=\frac{1}{\mathrm{d}x}\frac{\int_{\mu_v}^\infty(\mu-\mu_v)dxP_3(h(0)+\mu_vx,\mu|x)\mathrm{d}\mu}{\int_{-\infty}^{\infty}\int_{-\infty}^{h(0)+\mu_vx}P_3(h,\mu|x)dhd\mu}\\
&=\frac{\int_{\mu_v}^\infty(\mu-\mu_v)P_h(h(0)+\mu_vx)P^{2-}(\mu)\mathrm{d}\mu}{\int_{-\infty}^{h(0)+\mu_vx}P_h(h)dh}\\
&=\Lambda(\mu_v)\frac{\mu_vP_h(h(0)+\mu_vx)}{f(h(0)+\mu_vx)}
\end{aligned}
\end{equation}
g ( x ) = d x 1 P ( β ∣ α ) = d x 1 P ( α ) P ( α , β ) = d x 1 ∫ − ∞ ∞ ∫ − ∞ h ( 0 ) + μ v x P 3 ( h , μ ∣ x ) d h d μ ∫ μ v ∞ ∫ h ( 0 ) + μ v ( x + d x ) − μ d x h ( 0 ) + μ v x P 3 ( h , μ ∣ x ) d h d μ = d x 1 ∫ − ∞ ∞ ∫ − ∞ h ( 0 ) + μ v x P 3 ( h , μ ∣ x ) d h d μ ∫ μ v ∞ ( μ − μ v ) d x P 3 ( h ( 0 ) + μ v x , μ ∣ x ) d μ = ∫ − ∞ h ( 0 ) + μ v x P h ( h ) d h ∫ μ v ∞ ( μ − μ v ) P h ( h ( 0 ) + μ v x ) P 2 − ( μ ) d μ = Λ ( μ v ) f ( h ( 0 ) + μ v x ) μ v P h ( h ( 0 ) + μ v x ) 在Smith遮蔽函数下s ( μ v − μ ( 0 ) ) s(\mu_v-\mu(0)) s ( μ v − μ ( 0 )) S ( h ) S(h) S ( h )
S ( h , ω ) = e − Λ ( ω ) ∫ 0 ∞ μ v P h ( h + μ v x ) f ( h + μ v x ) d x = e − Λ ( ω ) [ ln f ( h + μ v x ) ] 0 ∞ = e − Λ ( ω ) ( ln 1 − ln f ( h ) ) = e Λ ( ω ) ln f ( h ) = f ( h ) Λ ( ω )
\begin{equation}
\begin{aligned}
S(h, \omega)
&= e^{-\Lambda(\omega)\int_0^{\infty}\frac{\mu_vP_h(h+\mu_vx)}{f(h+\mu_vx)}\mathrm{d}x}\\
&= e^{-\Lambda(\omega)\left[\ln f(h+\mu_v x)\right]_0^{\infty}}\\
&= e^{-\Lambda(\omega)(\ln 1 - \ln f(h))}\\
&= e^{\Lambda(\omega)\ln f(h)}\\
&= f(h)^{\Lambda(\omega)}
\end{aligned}
\end{equation}
S ( h , ω ) = e − Λ ( ω ) ∫ 0 ∞ f ( h + μ v x ) μ v P h ( h + μ v x ) d x = e − Λ ( ω ) [ l n f ( h + μ v x ) ] 0 ∞ = e − Λ ( ω ) ( l n 1 − l n f ( h )) = e Λ ( ω ) l n f ( h ) = f ( h ) Λ ( ω ) 此时可以得到Smith遮蔽函数的全局项, 可以看到与上一节的证明结果是相同的.
G 1 ( ω ) = ∫ − ∞ ∞ f ( h ) Λ ( ω ) P h ( h ) d h = ∫ − ∞ ∞ f ( h ) Λ ( ω ) d f ( h ) = [ f ( h ) 1 + Λ ( ω ) 1 + Λ ( ω ) ] − ∞ ∞ = 1 1 + Λ ( ω )
\begin{equation}
\begin{aligned}
G_1(\omega)
&= \int_{-\infty}^{\infty} f(h)^{\Lambda(\omega)} P_h(h) dh\\
&= \int_{-\infty}^{\infty} f(h)^{\Lambda(\omega)} df(h)\\
&= \left[\frac{f(h)^{1+\Lambda(\omega)}}{1+\Lambda(\omega)}\right]_{-\infty}^{\infty}\\
&= \frac{1}{1+\Lambda(\omega)}
\end{aligned}
\end{equation}
G 1 ( ω ) = ∫ − ∞ ∞ f ( h ) Λ ( ω ) P h ( h ) d h = ∫ − ∞ ∞ f ( h ) Λ ( ω ) df ( h ) = [ 1 + Λ ( ω ) f ( h ) 1 + Λ ( ω ) ] − ∞ ∞ = 1 + Λ ( ω ) 1 此时可以推导出遮蔽阴影函数.
G ( ω o , ω i ) = ∫ − ∞ ∞ f ( h ) Λ ( ω o ) P h ( h ) d h + ∫ − ∞ ∞ f ( h ) Λ ( ω i ) P h ( h ) d h = ∫ − ∞ ∞ f ( h ) Λ ( ω o ) + Λ ( ω i ) P h ( h ) d h = 1 1 + Λ ( ω o ) + Λ ( ω i )
\begin{equation}
\begin{aligned}
G(\omega_o,\omega_i)
&= \int_{-\infty}^{\infty}f(h)^{\Lambda(\omega_o)}P_h(h)dh + \int_{-\infty}^{\infty}f(h)^{\Lambda(\omega_i)}P_h(h)dh\\
&= \int_{-\infty}^{\infty}f(h)^{\Lambda(\omega_o)+\Lambda(\omega_i)}P_h(h)dh\\
&= \frac{1}{1+\Lambda(\omega_o)+\Lambda(\omega_i)}
\end{aligned}
\end{equation}
G ( ω o , ω i ) = ∫ − ∞ ∞ f ( h ) Λ ( ω o ) P h ( h ) d h + ∫ − ∞ ∞ f ( h ) Λ ( ω i ) P h ( h ) d h = ∫ − ∞ ∞ f ( h ) Λ ( ω o ) + Λ ( ω i ) P h ( h ) d h = 1 + Λ ( ω o ) + Λ ( ω i ) 1 采样可见法线分布
#
根据法线分布函数与遮蔽阴影函数我们可以得到视线为ω \omega ω
D ω ( ω m ) = G 1 ( ω ) cos θ D ( ω m ) max ( 0 , ω ⋅ ω m )
\begin{equation}
D_\omega(\omega_m)=\frac{G_1(\omega)}{\cos\theta}D(\omega_m)\max(0,\omega\cdot\omega_m)
\end{equation}
D ω ( ω m ) = cos θ G 1 ( ω ) D ( ω m ) max ( 0 , ω ⋅ ω m ) 由于逆变换法没有解析形式, pbrt将Trowbridge-Reitz分布看作一个被缩放的半球, 投影到半球上为椭圆. 因此将法线变换到z z z
// Transform _w_ to hemispherical configuration
Vector3f wh = Normalize(Vector3f(alpha_x * w.x, alpha_y * w.y, w.z));
if (wh.z < 0)
wh = -wh;
构造以视线为z z z T2轴保证采样到缩放半圆的法线在该轴为负.
// Find orthonormal basis for visible normal sampling
Vector3f T1 = (wh.z < 0.99999f) ? Normalize(Cross(Vector3f(0, 0, 1), wh))
: Vector3f(1, 0, 0);
Vector3f T2 = Cross(wh, T1);
将半球投影到x y xy x y cos θ = n ⋅ ω \cos\theta=\mathbf{n}\cdot\omega cos θ = n ⋅ ω wh.z的半圆. 令h = 1 − x 2 h=\sqrt{1-x^2} h = 1 − x 2 x x x y y y ( − h , h ) (-h, h) ( − h , h ) ( − h cos θ , h ) (-h\cos\theta, h) ( − h cos θ , h ) y ′ = 1 + cos θ 2 y + h ( 1 − cos θ ) 2 y'=\frac{1+\cos\theta}{2}y+\frac{h(1-\cos\theta)}{2} y ′ = 2 1 + c o s θ y + 2 h ( 1 − c o s θ ) Lerp与之等价.
// Generate uniformly distributed points on the unit disk
Point2f p = SampleUniformDiskPolar(u);
// Warp hemispherical projection for visible normal sampling
Float h = std::sqrt(1 - Sqr(p.x));
p.y = Lerp((1 + wh.z) / 2, h, p.y);
最后将采样得到的法线变换回原本的空间, 由于法线的变换为逆变换的转置, 这里仍然乘上粗糙度.
// Reproject to hemisphere and transform normal to ellipsoid configuration
Float pz = std::sqrt(std::max<Float>(0, 1 - LengthSquared(Vector2f(p))));
Vector3f nh = p.x * T1 + p.y * T2 + pz * wh;
CHECK_RARE(1e-5f, nh.z == 0);
return Normalize(
Vector3f(alpha_x * nh.x, alpha_y * nh.y, std::max<Float>(1e-6f, nh.z)));
实时渲染IBL的重要性抽样中, 法线与视线方向是一致的, 因此可以对下式进行球面坐标上的重要性抽样.
∫ − π π ∫ 0 π 2 D ( θ , ϕ ) cos θ sin θ d θ d ϕ = 1
\begin{equation}
\int_{-\pi}^{\pi}\int_0^{\frac{\pi}{2}}D(\theta,\phi)\cos\theta\sin\theta \mathrm{d}\theta \mathrm{d}\phi = 1
\end{equation}
∫ − π π ∫ 0 2 π D ( θ , ϕ ) cos θ sin θ d θ d ϕ = 1 各向同性的Trowbridge-Reitz在θ , ϕ \theta,\phi θ , ϕ
p ( θ ) = ∫ 0 2 π D ( θ , ϕ ) cos θ sin θ d ϕ = ∫ 0 2 π 1 π α 2 cos 4 θ ( 1 + tan 2 θ α 2 ) 2 cos θ sin θ d ϕ = 2 α 2 cos θ sin θ cos 4 θ ( α 2 + tan 2 θ ) 2 = 2 α 2 cos θ sin θ ( α 2 cos 2 θ + sin 2 θ ) 2 = 2 α 2 cos θ sin θ ( ( α 2 − 1 ) cos 2 θ + 1 ) 2 p ( ϕ ) = p ( θ , ϕ ) p ( θ ) = 1 2 π
\begin{equation}
\begin{aligned}
p(\theta)
&=\int_{0}^{2\pi}D(\theta,\phi)\cos\theta\sin\theta \mathrm{d}\phi\\
&=\int_{0}^{2\pi}\frac{1}{\pi\alpha^2\cos^4\theta(1+\frac{\tan^2\theta}{\alpha^2})^2}\cos\theta\sin\theta \mathrm{d}\phi\\
&=\frac{2\alpha^2\cos\theta\sin\theta}{\cos^4\theta(\alpha^2+\tan^2\theta)^2}\\
&=\frac{2\alpha^2\cos\theta\sin\theta}{(\alpha^2\cos^2\theta+\sin^2\theta)^2}\\
&=\frac{2\alpha^2\cos\theta\sin\theta}{((\alpha^2-1)\cos^2\theta+1)^2}\\
p(\phi)&=\frac{p(\theta,\phi)}{p(\theta)}=\frac{1}{2\pi}
\end{aligned}
\end{equation}
p ( θ ) p ( ϕ ) = ∫ 0 2 π D ( θ , ϕ ) cos θ sin θ d ϕ = ∫ 0 2 π π α 2 cos 4 θ ( 1 + α 2 t a n 2 θ ) 2 1 cos θ sin θ d ϕ = cos 4 θ ( α 2 + tan 2 θ ) 2 2 α 2 cos θ sin θ = ( α 2 cos 2 θ + sin 2 θ ) 2 2 α 2 cos θ sin θ = (( α 2 − 1 ) cos 2 θ + 1 ) 2 2 α 2 cos θ sin θ = p ( θ ) p ( θ , ϕ ) = 2 π 1 由于逆变换法的需要, 计算二者的CDF.
P ( θ ) = ∫ 0 θ p ( θ 0 ) d θ 0 = ∫ 0 θ − 2 α 2 cos θ 0 ( ( α 2 − 1 ) cos 2 θ 0 + 1 ) 2 d cos θ 0 = ∫ 1 cos θ − 2 α 2 x ( ( α 2 − 1 ) x 2 + 1 ) 2 d x = ∫ 1 cos θ α 2 ( α 2 − 1 ) d 1 ( α 2 − 1 ) x 2 + 1 = α 2 ( α 2 − 1 ) ( 1 ( α 2 − 1 ) cos 2 θ + 1 − 1 α 2 ) = 1 − cos 2 θ ( ( α 2 − 1 ) cos 2 θ + 1 ) P ( ϕ ) = ∫ 0 ϕ p ( ϕ 0 ) d ϕ 0 = ϕ 2 π
\begin{equation}
\begin{aligned}
P(\theta)
&=\int_0^{\theta}p(\theta_0)\mathrm{d}\theta_0\\
&=\int_0^{\theta}\frac{-2\alpha^2\cos\theta_0}{((\alpha^2-1)\cos^2\theta_0+1)^2}d\cos\theta_0\\
&=\int_1^{\cos\theta}\frac{-2\alpha^2x}{((\alpha^2-1)x^2+1)^2}\mathrm{d}x\\
&=\int_1^{\cos\theta}\frac{\alpha^2}{(\alpha^2-1)}d\frac{1}{(\alpha^2-1)x^2+1}\\
&=\frac{\alpha^2}{(\alpha^2-1)}(\frac{1}{(\alpha^2-1)\cos^2\theta+1}-\frac{1}{\alpha^2})\\
&=\frac{1-\cos^2\theta}{((\alpha^2-1)\cos^2\theta+1)}\\
P(\phi)
&=\int_0^{\phi}p(\phi_0)\mathrm{d}\phi_0=\frac{\phi}{2\pi}
\end{aligned}
\end{equation}
P ( θ ) P ( ϕ ) = ∫ 0 θ p ( θ 0 ) d θ 0 = ∫ 0 θ (( α 2 − 1 ) cos 2 θ 0 + 1 ) 2 − 2 α 2 cos θ 0 d cos θ 0 = ∫ 1 c o s θ (( α 2 − 1 ) x 2 + 1 ) 2 − 2 α 2 x d x = ∫ 1 c o s θ ( α 2 − 1 ) α 2 d ( α 2 − 1 ) x 2 + 1 1 = ( α 2 − 1 ) α 2 ( ( α 2 − 1 ) cos 2 θ + 1 1 − α 2 1 ) = (( α 2 − 1 ) cos 2 θ + 1 ) 1 − cos 2 θ = ∫ 0 ϕ p ( ϕ 0 ) d ϕ 0 = 2 π ϕ 通过逆变换法得到最终变换.
θ = P − 1 ( u 0 ) = cos − 1 1 − u 0 ( α 2 − 1 ) u 0 + 1 ϕ = P − 1 ( u 1 ) = 2 π u 1
\begin{equation}
\begin{aligned}
\theta
&=P^{-1}(u_0)
=\cos^{-1}\sqrt{\frac{1-u_0}{(\alpha^2-1)u_0+1}}\\
\phi
&=P^{-1}(u_1)
=2\pi u_1
\end{aligned}
\end{equation}
θ ϕ = P − 1 ( u 0 ) = cos − 1 ( α 2 − 1 ) u 0 + 1 1 − u 0 = P − 1 ( u 1 ) = 2 π u 1 Torrance-Sparrow模型
#
综合以上结论我们可以推导出微表面上的Torrance-Sparrow BRDF并应用于材质中.
半向量变换
#
对于镜面反射表面法线ω m \omega_m ω m ω o \omega_o ω o ω o \omega_o ω o ω o \omega_o ω o θ o = 0 \theta_o=0 θ o = 0
d ω m d ω i = sin θ m d θ m d ϕ m sin θ i d θ i d ϕ i = sin θ m d θ m d ϕ m sin 2 θ m 2 d θ m d ϕ m = sin θ m 4 cos θ m sin θ m = 1 4 ( ω i ⋅ ω m ) = 1 4 ( ω o ⋅ ω m )
\begin{equation}
\begin{aligned}
\frac{\mathrm{d}\omega_m}{\mathrm{d}\omega_i}
&= \frac{\sin\theta_m \mathrm{d}\theta_m \mathrm{d}\phi_m}{\sin\theta_i \mathrm{d}\theta_i \mathrm{d}\phi_i}\\
&= \frac{\sin\theta_m \mathrm{d}\theta_m \mathrm{d}\phi_m}{\sin2\theta_m 2\mathrm{d}\theta_m \mathrm{d}\phi_m}\\
&= \frac{\sin\theta_m}{4\cos\theta_m\sin\theta_m}\\
&= \frac{1}{4(\omega_i\cdot\omega_m)}\\
&= \frac{1}{4(\omega_o\cdot\omega_m)}
\end{aligned}
\end{equation}
d ω i d ω m = sin θ i d θ i d ϕ i sin θ m d θ m d ϕ m = sin 2 θ m 2 d θ m d ϕ m sin θ m d θ m d ϕ m = 4 cos θ m sin θ m sin θ m = 4 ( ω i ⋅ ω m ) 1 = 4 ( ω o ⋅ ω m ) 1 Torrance-Sparrow PDF
#
ω i \omega_i ω i
p ( ω i ) = D ω o ( ω m ) d ω m d ω i = D ω o ( ω m ) 4 ( ω o ⋅ ω m )
\begin{equation}
p(\omega_i)=D_{\omega_o}(\omega_m)\frac{\mathrm{d}\omega_m}{\mathrm{d}\omega_i}=\frac{D_{\omega_o}(\omega_m)}{4(\omega_o\cdot\omega_m)}
\end{equation}
p ( ω i ) = D ω o ( ω m ) d ω i d ω m = 4 ( ω o ⋅ ω m ) D ω o ( ω m ) Torrance-Sparrow BRDF
#
在Monte Carlo下反射方程可以做如下近似.
L o ( p , ω o ) = ∫ Ω f r ( p , ω o , ω i ) L i ( p , ω i ) cos θ i d ω i ≈ f r ( p , ω o , ω i ) L i ( p , ω i ) cos θ i p ( ω i ) = F ( ω o ⋅ ω m ) G ( ω i , ω o ) L i ( p , ω i )
\begin{equation}
\begin{aligned}
L_o(p,\omega_o)
&=\int_\Omega f_r(p,\omega_o,\omega_i)L_i(p,\omega_i)\cos\theta_id\omega_i\\
&\approx \frac{f_r(p,\omega_o,\omega_i)L_i(p,\omega_i)\cos\theta_i}{p(\omega_i)}\\
&=F(\omega_o\cdot\omega_m)G(\omega_i,\omega_o)L_i(p,\omega_i)
\end{aligned}
\end{equation}
L o ( p , ω o ) = ∫ Ω f r ( p , ω o , ω i ) L i ( p , ω i ) cos θ i d ω i ≈ p ( ω i ) f r ( p , ω o , ω i ) L i ( p , ω i ) cos θ i = F ( ω o ⋅ ω m ) G ( ω i , ω o ) L i ( p , ω i ) 由此可得Torrance-Sparrow BRDF.
f r ( p , ω o , ω i ) = D ω o ( ω m ) F ( ω o , ω m ) G ( ω i ) 4 ( ω o ⋅ ω m ) cos θ i = D ( ω m ) F ( ω o , ω m ) G ( ω i ) G ( ω o ) 4 cos θ i cos θ o ≈ D ( ω m ) F ( ω o , ω m ) G ( ω i , ω o ) 4 cos θ i cos θ o
\begin{equation}
\begin{aligned}
f_r(p,\omega_o,\omega_i)
&=\frac{D_{\omega_o}(\omega_m)F(\omega_o,\omega_m)G(\omega_i)}{4(\omega_o\cdot\omega_m)\cos\theta_i}\\
&=\frac{D(\omega_m)F(\omega_o,\omega_m)G(\omega_i)G(\omega_o)}{4\cos\theta_i\cos\theta_o}\\
&\approx\frac{D(\omega_m)F(\omega_o,\omega_m)G(\omega_i, \omega_o)}{4\cos\theta_i\cos\theta_o}
\end{aligned}
\end{equation}
f r ( p , ω o , ω i ) = 4 ( ω o ⋅ ω m ) cos θ i D ω o ( ω m ) F ( ω o , ω m ) G ( ω i ) = 4 cos θ i cos θ o D ( ω m ) F ( ω o , ω m ) G ( ω i ) G ( ω o ) ≈ 4 cos θ i cos θ o D ( ω m ) F ( ω o , ω m ) G ( ω i , ω o ) Torrance-Sparrow采样
#
若Sample_wm得到的ω i \omega_i ω i
Vector3f wm = mfDistrib.Sample_wm(wo, u);
Vector3f wi = Reflect(wo, wm);
if (!SameHemisphere(wo, wi))
return {};
粗糙绝缘体BSDF
#
粗糙绝缘体PDF
#
根据Snell定律我们可以把投影到与半向量垂直的平面上的部分抵消, 由此可得以下关系. 此时得到的半向量总是指向光学密度较大的一侧, 需要后续修正. 对于反射由于η i = η o \eta_i=\eta_o η i = η o
ω m = η i ω i + η o ω o ∥ η i ω i + η o ω o ∥ = η ω i + ω o ∥ η ω i + ω o ∥
\begin{equation}
\omega_m=\frac{\eta_i\omega_i+\eta_o\omega_o}{\Vert\eta_i\omega_i+\eta_o\omega_o\Vert}=\frac{\eta\omega_i+\omega_o}{\Vert\eta\omega_i+\omega_o\Vert}
\end{equation}
ω m = ∥ η i ω i + η o ω o ∥ η i ω i + η o ω o = ∥ η ω i + ω o ∥ η ω i + ω o 此时Jacobi行列式结果如下, 这可以由几何关系推导得到, 具体见Water et al. 2007 .
d ω m d ω i = ∣ ω i ⋅ ω m ∣ η i 2 ( ( η i ω i + η o ω o ) ⋅ ω m ) 2 = ∣ ω i ⋅ ω m ∣ ( ( ω i ⋅ ω m ) + ( ω o ⋅ ω m ) η ) 2
\begin{equation}
\begin{aligned}
\frac{\mathrm{d}\omega_m}{\mathrm{d}\omega_i}
&=\frac{|\omega_i\cdot\omega_m|\eta_i^2}{((\eta_i\omega_i+\eta_o\omega_o)\cdot\omega_m)^2}\\
&=\frac{|\omega_i\cdot\omega_m|}{((\omega_i\cdot\omega_m)+\frac{(\omega_o\cdot\omega_m)}{\eta})^2}
\end{aligned}
\end{equation}
d ω i d ω m = (( η i ω i + η o ω o ) ⋅ ω m ) 2 ∣ ω i ⋅ ω m ∣ η i 2 = (( ω i ⋅ ω m ) + η ( ω o ⋅ ω m ) ) 2 ∣ ω i ⋅ ω m ∣
粗糙绝缘体BSDF
#
根据PDF可以将Torrance-Sparrow BRDF拓展到BTDF, 得到完整的BSDF.
f t ( p , ω o , ω i ) = D ω o ( ω m ) ( 1 − F ( ω o ⋅ ω m ) ) G ( ω i ) ∣ ω i ⋅ ω m ∣ ( ( ω i ⋅ ω m ) + ( ω o ⋅ ω m ) η ) 2 ∣ cos θ i ∣ = D ( ω m ) ( 1 − F ( ω o ⋅ ω m ) ) G ( ω i , ω o ) ( ( ω i ⋅ ω m ) + ( ω o ⋅ ω m ) η ) 2 ∣ ω i ⋅ ω m ∣ ∣ ω o ⋅ ω m ∣ ∣ cos θ i ∣ ∣ cos θ o ∣
\begin{equation}
\begin{aligned}
f_t(p,\omega_o,\omega_i)
&=\frac{D_{\omega_o}(\omega_m)(1-F(\omega_o\cdot\omega_m))G(\omega_i)|\omega_i\cdot\omega_m|}{((\omega_i\cdot\omega_m)+\frac{(\omega_o\cdot\omega_m)}{\eta})^2|\cos\theta_i|}\\
&=\frac{D(\omega_m)(1-F(\omega_o\cdot\omega_m))G(\omega_i,\omega_o)}{((\omega_i\cdot\omega_m)+\frac{(\omega_o\cdot\omega_m)}{\eta})^2}\frac{|\omega_i\cdot\omega_m||\omega_o\cdot\omega_m|}{|\cos\theta_i||\cos\theta_o|}
\end{aligned}
\end{equation}
f t ( p , ω o , ω i ) = (( ω i ⋅ ω m ) + η ( ω o ⋅ ω m ) ) 2 ∣ cos θ i ∣ D ω o ( ω m ) ( 1 − F ( ω o ⋅ ω m )) G ( ω i ) ∣ ω i ⋅ ω m ∣ = (( ω i ⋅ ω m ) + η ( ω o ⋅ ω m ) ) 2 D ( ω m ) ( 1 − F ( ω o ⋅ ω m )) G ( ω i , ω o ) ∣ cos θ i ∣∣ cos θ o ∣ ∣ ω i ⋅ ω m ∣∣ ω o ⋅ ω m ∣ 粗糙绝缘体采样
#
与之前章节描述的一样, 对法线进行重要性抽样并随机选取反射或折射.
塑料BSDF
#
塑料BSDF被建模为漫反射基础表面涂上较薄的绝缘体层, 与薄绝缘体模型类似, 塑料BSDF需要求解所有可能的内部反射情况之和.
光滑塑料
#
假设涂层为光滑表面, 对于只在内部漫反射一次的情况, 由于表面光滑, 漫反射的出射方向是确定的, 使用Dirac delta函数表示, 使用微分Snell定律计算Jacobian, 此时结果如下, 其中θ i t \theta_{it} θ i t θ o t \theta_{ot} θ o t
L o = ∫ Ω ( 1 − F ( ω ⋅ n ) ) L ω η 2 cos ( ω ⋅ n ) δ ( ω o t ) cos ( ω o t ⋅ n ) d ω = ( 1 − F ( ω o ⋅ n ) ) L o t η 2 = ( 1 − F ( ω o ⋅ n ) ) η 2 ∫ Ω i t R π L ω i t cos ω i t d ω i t = ( 1 − F ( ω o ⋅ n ) ) η 2 ∫ Ω i ( 1 − F ( ω i ⋅ n ) ) R π L ω i 1 η 2 cos ω i η 2 d ω i = ( 1 − F ( ω o ⋅ n ) ) ( 1 − F ( ω i ⋅ n ) ) η 2 ∫ Ω i R π L ω i cos ω i d ω i
\begin{equation}
\begin{aligned}
L_o
&=\int_{\Omega} (1 - F(\omega\cdot\mathbf{n})) L_\omega \eta^2 \cos(\omega\cdot\mathbf{n}) \frac{\delta(\omega_{ot})}{\cos(\omega_{ot}\cdot\mathbf{n})} \mathrm{d}\omega\\
&=(1 - F(\omega_{o}\cdot\mathbf{n})) L_{ot} \eta^2\\
&=(1 - F(\omega_{o}\cdot\mathbf{n})) \eta^2 \int_{\Omega_{it}} \frac{R}{\pi} L_{\omega_{it}} \cos\omega_{it} \mathrm{d}\omega_{it}\\
&=(1 - F(\omega_{o}\cdot\mathbf{n})) \eta^2 \int_{\Omega_{i}} (1 - F(\omega_{i}\cdot\mathbf{n})) \frac{R}{\pi} L_{\omega_{i}} \frac{1}{\eta^2} \cos\omega_{i} \eta^2 \mathrm{d}\omega_{i}\\
&=(1 - F(\omega_{o}\cdot\mathbf{n})) (1 - F(\omega_{i}\cdot\mathbf{n})) \eta^2 \int_{\Omega_{i}} \frac{R}{\pi} L_{\omega_{i}} \cos\omega_{i} \mathrm{d}\omega_{i}\\
\end{aligned}
\end{equation}
L o = ∫ Ω ( 1 − F ( ω ⋅ n )) L ω η 2 cos ( ω ⋅ n ) cos ( ω o t ⋅ n ) δ ( ω o t ) d ω = ( 1 − F ( ω o ⋅ n )) L o t η 2 = ( 1 − F ( ω o ⋅ n )) η 2 ∫ Ω i t π R L ω i t cos ω i t d ω i t = ( 1 − F ( ω o ⋅ n )) η 2 ∫ Ω i ( 1 − F ( ω i ⋅ n )) π R L ω i η 2 1 cos ω i η 2 d ω i = ( 1 − F ( ω o ⋅ n )) ( 1 − F ( ω i ⋅ n )) η 2 ∫ Ω i π R L ω i cos ω i d ω i 整理为BSDF可得如下结果.
f r ( p , ω o , ω i ) = ( 1 − F ( ω o ⋅ n ) ) ( 1 − F ( ω i ⋅ n ) ) η 2 R π = F t R π
\begin{equation}
\begin{aligned}
f_r(p, \omega_o, \omega_i)
&=(1 - F(\omega_{o}\cdot\mathbf{n})) (1 - F(\omega_{i}\cdot\mathbf{n})) \eta^2 \frac{R}{\pi}\\
&=F_t \frac{R}{\pi}\\
\end{aligned}
\end{equation}
f r ( p , ω o , ω i ) = ( 1 − F ( ω o ⋅ n )) ( 1 − F ( ω i ⋅ n )) η 2 π R = F t π R 对于漫反射+镜面反射+漫反射的情况, 令ω \omega ω ω ⋅ n \omega\cdot\mathbf{n} ω ⋅ n
f r ( p , ω o , ω i ) = F t R π ∫ ω R π F ( ω ⋅ n ) cos θ d ω = F t R π ∫ 0 2 π ∫ 0 π 2 R π F ( cos θ ) sin θ cos θ d θ d ϕ = F t 2 R 2 π ∫ 0 π 2 F ( cos θ ) sin θ cos θ d θ = F t 2 R 2 π ∫ 0 π 2 − F ( cos θ ) 2 d cos 2 θ = F t R 2 π ∫ 0 1 F ( x ) d x = F t F r R 2 π
\begin{equation}
\begin{aligned}
f_r(p, \omega_o, \omega_i)
&= F_t \frac{R}{\pi} \int_\omega \frac{R}{\pi} F(\omega \cdot \mathbf{n}) \cos\theta \mathrm{d}\omega\\
&= F_t \frac{R}{\pi} \int_0^{2\pi} \int_0^{\frac{\pi}{2}} \frac{R}{\pi} F(\cos\theta) \sin\theta \cos\theta \mathrm{d}\theta \mathrm{d}\phi\\
&= F_t \frac{2R^2}{\pi} \int_0^{\frac{\pi}{2}} F(\cos\theta) \sin\theta \cos\theta \mathrm{d}\theta\\
&= F_t \frac{2R^2}{\pi} \int_0^{\frac{\pi}{2}} -\frac{F(\cos\theta)}{2} d\cos^2\theta\\
&= F_t \frac{R^2}{\pi} \int_0^1 F(\sqrt{x}) \mathrm{d}x\\
&= F_t \frac{F_r R^2}{\pi}
\end{aligned}
\end{equation}
f r ( p , ω o , ω i ) = F t π R ∫ ω π R F ( ω ⋅ n ) cos θ d ω = F t π R ∫ 0 2 π ∫ 0 2 π π R F ( cos θ ) sin θ cos θ d θ d ϕ = F t π 2 R 2 ∫ 0 2 π F ( cos θ ) sin θ cos θ d θ = F t π 2 R 2 ∫ 0 2 π − 2 F ( cos θ ) d cos 2 θ = F t π R 2 ∫ 0 1 F ( x ) d x = F t π F r R 2 递推任意次反射, 不考虑表面反射, 结果如下. 重要性抽样可以使用余弦加权半球采样, 因为1 − F ( ω ⋅ n ) 1-F(\omega\cdot\mathbf{n}) 1 − F ( ω ⋅ n ) cos θ \cos\theta cos θ
f r ( p , ω o , ω i ) = ∑ n = 0 ∞ F t R n + 1 F r n π = F t π R 1 − R F r
\begin{equation}
\begin{aligned}
f_r(p, \omega_o, \omega_i)
&= \sum_{n=0}^\infty F_t \frac{R^{n+1} F_r^n}{\pi}\\
&= \frac{F_t}{\pi} \frac{R}{1-RF_r}
\end{aligned}
\end{equation}
f r ( p , ω o , ω i ) = n = 0 ∑ ∞ F t π R n + 1 F r n = π F t 1 − R F r R 粗糙塑料
#
对与带有微表面的涂层, 我们仍然首先证明单次漫反射, 由于BSDF可逆和积分顺序可交换, 最终结果如下.
L o = ∫ Ω o t f o ( p , ω o , ω o t ) L ω o t η 2 cos ( ω o t ⋅ n ) d ω o t = ∫ Ω o t f o ( p , ω o , ω o t ) cos ( ω o t ⋅ n ) η 2 d ω o t ∫ Ω i t R π L ω i t cos ( ω i t ⋅ n ) d ω i t = η 2 R π T ( ω o ) ∫ Ω i t ∫ Ω i f i ( p , ω i t , ω i ) L ω i cos ( ω i ⋅ n ) cos ( ω i t ⋅ n ) d ω i d ω i t = η 2 R π T ( ω o ) ∫ Ω i t f i ( p , ω i t , ω i ) cos ( ω i t ⋅ n ) d ω i t ∫ Ω i L ω i cos ( ω i ⋅ n ) d ω i = η 2 R π T ( ω o ) ∫ Ω i T ( ω i ) L ω i cos ( ω i ⋅ n ) d ω i
\begin{equation}
\begin{aligned}
L_o
&=\int_{\Omega_{ot}} f_o(p, \omega_o, \omega_{ot}) L_{\omega_{ot}} \eta^2 \cos(\omega_{ot}\cdot\mathbf{n}) \mathrm{d}\omega_{ot}\\
&=\int_{\Omega_{ot}} f_o(p, \omega_o, \omega_{ot}) \cos(\omega_{ot}\cdot\mathbf{n}) \eta^2 \mathrm{d}\omega_{ot} \int_{\Omega_{it}} \frac{R}{\pi} L_{\omega_{it}} \cos(\omega_{it}\cdot\mathbf{n}) \mathrm{d}\omega_{it}\\
&=\eta^2 \frac{R}{\pi} T(\omega_o) \int_{\Omega_{it}} \int_{\Omega_{i}} f_i(p, \omega_{it}, \omega{i}) L_{\omega_{i}} \cos(\omega_{i}\cdot\mathbf{n}) \cos(\omega_{it}\cdot\mathbf{n}) \mathrm{d}\omega_{i} \mathrm{d}\omega_{it}\\
&=\eta^2 \frac{R}{\pi} T(\omega_o) \int_{\Omega_{it}} f_i(p, \omega_{it}, \omega{i}) \cos(\omega_{it}\cdot\mathbf{n}) \mathrm{d}\omega_{it} \int_{\Omega_{i}} L_{\omega_{i}} \cos(\omega_{i}\cdot\mathbf{n}) \mathrm{d}\omega_{i}\\
&=\eta^2 \frac{R}{\pi} T(\omega_o) \int_{\Omega_{i}} T(\omega_i) L_{\omega_{i}} \cos(\omega_{i}\cdot\mathbf{n}) \mathrm{d}\omega_{i}\\
\end{aligned}
\end{equation}
L o = ∫ Ω o t f o ( p , ω o , ω o t ) L ω o t η 2 cos ( ω o t ⋅ n ) d ω o t = ∫ Ω o t f o ( p , ω o , ω o t ) cos ( ω o t ⋅ n ) η 2 d ω o t ∫ Ω i t π R L ω i t cos ( ω i t ⋅ n ) d ω i t = η 2 π R T ( ω o ) ∫ Ω i t ∫ Ω i f i ( p , ω i t , ω i ) L ω i cos ( ω i ⋅ n ) cos ( ω i t ⋅ n ) d ω i d ω i t = η 2 π R T ( ω o ) ∫ Ω i t f i ( p , ω i t , ω i ) cos ( ω i t ⋅ n ) d ω i t ∫ Ω i L ω i cos ( ω i ⋅ n ) d ω i = η 2 π R T ( ω o ) ∫ Ω i T ( ω i ) L ω i cos ( ω i ⋅ n ) d ω i 整理为BSDF可得如下结果, 后续过程不再证明, 与光滑情况下结果类似. 透射率与反射率可做预计算, 比如mitsuba3 的实现. 若材质处理计算量大, 可用Fresnel替换, 比如tungsten 的实现.
f r ( p , ω o , ω i ) = η 2 R π T ( ω i ) T ( ω o ) = F t R π
\begin{equation}
\begin{aligned}
f_r(p, \omega_o, \omega_i)
&= \eta^2 \frac{R}{\pi} T(\omega_i) T(\omega_o)
&= F_t \frac{R}{\pi}
\end{aligned}
\end{equation}
f r ( p , ω o , ω i ) = η 2 π R T ( ω i ) T ( ω o ) = F t π R BSDF测量值
#
现实世界中的很多材质无法通过微表面模型表达, pbrt采用测量值来实现数据驱动的反射模型, 通过采样估计的粗糙度对应的法线分布来提高测量效率.
测量过程
#
参数化测量
#
由于测量出的BRDF在量级上可能差距较大, 不适合线性插值, pbrt将采样值变换为Monte Carlo权重的形式再插值.
M ( k ) = f r ( ω o , ω i ( k ) ) cos θ i ( k ) p ( ω i ( k ) )
\begin{equation}
M^{(k)} = \frac{f_r(\omega_o, \omega_i^{(k)})\cos\theta_i^{(k)}}{p(\omega_i^{(k)})}
\end{equation}
M ( k ) = p ( ω i ( k ) ) f r ( ω o , ω i ( k ) ) cos θ i ( k ) BRDF求解
#
令M M M R ( ω o , u ) = ω i R(\omega_o, u)=\omega_i R ( ω o , u ) = ω i
f r ( ω o , ω i ) = M ( R − 1 ( ω o , ω i ) ) p ( ω i ) cos θ i
\begin{equation}
f_r(\omega_o,\omega_i)=\frac{M(R^{-1}(\omega_o,\omega_i))p(\omega_i)}{\cos\theta_i}
\end{equation}
f r ( ω o , ω i ) = cos θ i M ( R − 1 ( ω o , ω i )) p ( ω i ) 泛化微表面模型
#
对于测量值BRDF, 可以使用最基本的Smith遮蔽函数的定义来得到PDF, 因为可以通过数值积分计算这里不采用Smith遮蔽函数的解析形式.
p ( ω i ) = D ω o ( ω m ) 4 ( ω o ⋅ ω m ) = D ( ω m ) max ( 0 , ω o ⋅ ω m ) 4 ( ω o ⋅ ω m ) ∫ Ω D ( ω ) max ( 0 , ω o , ω ) d ω = D ( ω m ) 4 σ ( ω o )
\begin{equation}
\begin{aligned}
p(\omega_i)
&= \frac{D_{\omega_o}(\omega_m)}{4(\omega_o\cdot\omega_m)}\\
&= \frac{D(\omega_m)\max(0,\omega_o\cdot\omega_m)}{4(\omega_o\cdot\omega_m)\int_\Omega D(\omega)\max(0,\omega_o,\omega)\mathrm{d}\omega}\\
&= \frac{D(\omega_m)}{4\sigma(\omega_o)}
\end{aligned}
\end{equation}
p ( ω i ) = 4 ( ω o ⋅ ω m ) D ω o ( ω m ) = 4 ( ω o ⋅ ω m ) ∫ Ω D ( ω ) max ( 0 , ω o , ω ) d ω D ( ω m ) max ( 0 , ω o ⋅ ω m ) = 4 σ ( ω o ) D ( ω m ) 此时可以求解BRDF.
f r ( ω o , ω i ) = M ( R − 1 ( ω o , ω i ) ) D ( ω m ) 4 σ ( ω o ) cos θ i
\begin{equation}
f_r(\omega_o,\omega_i)=\frac{M(R^{-1}(\omega_o,\omega_i))D(\omega_m)}{4\sigma(\omega_o)\cos\theta_i}
\end{equation}
f r ( ω o , ω i ) = 4 σ ( ω o ) cos θ i M ( R − 1 ( ω o , ω i )) D ( ω m ) 初始微表面模型
#
法线分布函数只有两个维度, 不会有维度组合数量的问题, 因此可以通过估计D ( ω ) D(\omega) D ( ω ) f r ( p , ω j , ω j ) f_r(p, \omega_j, \omega_j) f r ( p , ω j , ω j ) D ( ω j ) D(\omega_j) D ( ω j )
f r ( p , ω , ω ) = D ( ω ) F ( ω ⋅ ω ) G ( ω , ω ) 4 cos 2 θ ∝ D ( ω ) G 1 ( ω ) cos 2 θ D ( ω ) ∝ f r ( p , ω , ω ) cos θ ∫ Ω D ( ω m ) max ( 0 , ω ⋅ ω m ) d ω m
\begin{equation}
\begin{aligned}
f_r(p,\omega,\omega)&=\frac{D(\omega)F(\omega\cdot\omega)G(\omega,\omega)}{4\cos^2\theta}\propto\frac{D(\omega)G_1(\omega)}{\cos^2\theta}\\
D(\omega)&\propto f_r(p,\omega,\omega)\cos\theta\int_\Omega D(\omega_m)\max(0,\omega\cdot\omega_m)\mathrm{d}\omega_m
\end{aligned}
\end{equation}
f r ( p , ω , ω ) D ( ω ) = 4 cos 2 θ D ( ω ) F ( ω ⋅ ω ) G ( ω , ω ) ∝ cos 2 θ D ( ω ) G 1 ( ω ) ∝ f r ( p , ω , ω ) cos θ ∫ Ω D ( ω m ) max ( 0 , ω ⋅ ω m ) d ω m 基础数据结构
#
通过PiecewiseLinear2D实现插值与采样.
struct MeasuredBxDFData {
// MeasuredBxDFData Public Members
pstd::vector<float> wavelengths;
PiecewiseLinear2D<3> spectra;
PiecewiseLinear2D<0> ndf;
PiecewiseLinear2D<2> vndf;
PiecewiseLinear2D<0> sigma;
bool isotropic;
PiecewiseLinear2D<2> luminance;
MeasuredBxDFData(Allocator alloc)
: ndf(alloc),
sigma(alloc),
vndf(alloc),
luminance(alloc),
spectra(alloc),
wavelengths(alloc) {}
static MeasuredBxDFData *Create(const std::string &filename, Allocator alloc);
std::string ToString() const {
return StringPrintf("[ MeasuredBxDFData filename: %s ]", filename);
}
std::string filename;
};
求解
#
MeasuredBxDF中存储MeasuredBxDFData指针并实现相应的功能.
毛发散射
#
各个波长在毛发上的反射率是类似的, 因此毛发的颜色主要由吸收的折射光的量决定. 这可以通过DielectricBxDF与体渲染结合实现, 但这种方式计算量较大且无法实现重要性抽样. 特化的毛发BSDF可以解决该问题.
几何
#
毛发通过Curve实现, pbrt假设毛发为朝向入射光的带状曲面, 毛发BSDF基于圆柱形曲面, 在毛发不相互交织且不比像素宽很多的情况下, 这种转换是可行的.
对于光线与圆柱曲面的交点u , v u,v u , v u u u ∂ p ∂ u \frac{\partial p}{\partial u} ∂ u ∂ p x x x ∂ p ∂ v \frac{\partial p}{\partial v} ∂ v ∂ p y y y x x x θ \theta θ y z yz yz y y y ϕ \phi ϕ y z yz yz
θ = π 2 − cos − 1 ( ω ⋅ x ) ϕ = π − sign ( ω ⋅ z ) ( π − cos − 1 ( ω ⋅ y ) )
\begin{equation}
\begin{aligned}
\theta &= \frac{\pi}{2}-\cos^{-1}(\omega\cdot x)\\
\phi &= \pi-\text{sign}(\omega\cdot z)(\pi - \cos^{-1}(\omega\cdot y))
\end{aligned}
\end{equation}
θ ϕ = 2 π − cos − 1 ( ω ⋅ x ) = π − sign ( ω ⋅ z ) ( π − cos − 1 ( ω ⋅ y )) 根据v v v − 1 + 2 v -1+2v − 1 + 2 v γ \gamma γ
γ = sin − 1 ( − 1 + 2 v )
\begin{equation}
\gamma = \sin^{-1}(-1+2v)
\end{equation}
γ = sin − 1 ( − 1 + 2 v ) 毛发散射
#
毛发包含发鳞, 发皮质, 发髓三个部分, 发髓由一些列表面的凸起构成, 发皮质与发髓包含吸收光的色素, 以发髓的颜色为主. pbrt假设毛发表面是与圆柱夹角为α \alpha α
光线会在毛发内部多次反射与折射, 令p p p p = 0 p=0 p = 0 R R R p = 1 p=1 p = 1 T T TT TT p = 3 p=3 p = 3 T R T TRT TRT M p M_p M p N p N_p N p A p A_p A p p ≤ 3 p\le 3 p ≤ 3
f ( ω o , ω i ) = ∑ p = 0 ∞ f p ( ω o , ω i ) = ∑ p = 0 ∞ M p ( θ o , θ i ) A p ( ω o ) N p ( ϕ ) cos θ i
\begin{equation}
\begin{aligned}
f(\omega_o,\omega_i)
&=\sum_{p=0}^{\infty} f_p(\omega_o,\omega_i)\\
&=\sum_{p=0}^{\infty} \frac{M_p(\theta_o,\theta_i)A_p(\omega_o)N_p(\phi)}{\cos\theta_i}\\
\end{aligned}
\end{equation}
f ( ω o , ω i ) = p = 0 ∑ ∞ f p ( ω o , ω i ) = p = 0 ∑ ∞ cos θ i M p ( θ o , θ i ) A p ( ω o ) N p ( ϕ ) HairBxDF构造函数如下, 其中h h h η \eta η σ a \sigma_a σ a β m \beta_m β m β n \beta_n β n α \alpha α
PBRT_CPU_GPU
HairBxDF(Float h, Float eta, const SampledSpectrum &sigma_a, Float beta_m,
Float beta_n, Float alpha);
纵向散射
#
纵向散射方程如下, 其中I 0 I_0 I 0 v v v
M p ( θ o , θ i ) = e − sin θ i sin θ o v 2 v sinh ( 1 v ) I 0 ( cos θ o cos θ i v )
\begin{equation}
M_p(\theta_o,\theta_i)=\frac{e^{-\frac{\sin\theta_i\sin\theta_o}{v}}}{2v\text{sinh}(\frac{1}{v})}I_0(\frac{\cos\theta_o\cos\theta_i}{v})
\end{equation}
M p ( θ o , θ i ) = 2 v sinh ( v 1 ) e − v s i n θ i s i n θ o I 0 ( v cos θ o cos θ i ) pbrt根据β m \beta_m β m p p p
v 0 = ( 0.726 β m + 0.812 β m 2 + 3.7 β m 2 0 ) 2 v 1 = 0.25 v 0 v 2 = v 3 = ⋯ = v max = 4 v 0
\begin{equation}
\begin{aligned}
v_0 &= (0.726\beta_m+0.812\beta_m^2+3.7\beta_m^20)^2\\
v_1 &= 0.25v_0\\
v_2 &= v_3 = \cdots = v_{\max} = 4v_0
\end{aligned}
\end{equation}
v 0 v 1 v 2 = ( 0.726 β m + 0.812 β m 2 + 3.7 β m 2 0 ) 2 = 0.25 v 0 = v 3 = ⋯ = v m a x = 4 v 0 纤维吸收
#
吸收率受光线传播距离影响, 这需要计算折射光线的纵向角与方位角. 纵向角可以通过Snell定律获取, 方位角可以通过将ω t \omega_t ω t y z yz yz
η ′ = η 2 − sin 2 θ o cos θ o
\begin{equation}
\eta'=\frac{\sqrt{\eta^2-\sin^2\theta_o}}{\cos\theta_o}
\end{equation}
η ′ = cos θ o η 2 − sin 2 θ o 根据直线与圆相交的性质可以得到在y z yz yz
l = 2 cos γ t cos θ t
\begin{equation}
l=\frac{2\cos\gamma_t}{\cos\theta_t}
\end{equation}
l = cos θ t 2 cos γ t 利用后面体渲染章节介绍的Beer定律可以得到经过吸收后的透射率.
T r = e − σ a l
\begin{equation}
T_r=e^{-\sigma_a l}
\end{equation}
T r = e − σ a l 此时可以得到A p A_p A p
A p = { A p − 1 T f = ( 1 − f ) 2 T p f p − 1 p > 0 f p = 0
\begin{equation}
\begin{aligned}
A_p&=
\begin{cases}
A_{p-1}Tf=(1-f)^2 T^p f^{p-1} & p>0\\
f & p=0
\end{cases}
\end{aligned}
\end{equation}
A p = { A p − 1 T f = ( 1 − f ) 2 T p f p − 1 f p > 0 p = 0 ∑ p = p max ∞ ( 1 − f ) 2 T p f p − 1 = ( 1 − f ) 2 T p max f p max − 1 1 − T f
\begin{equation}
\sum_{p=p_{\max}}^{\infty}(1-f)^2 T^p f^{p-1}=\frac{(1-f)^2 T^{p_{\max}} f^{p_{\max}-1}}{1-Tf}
\end{equation}
p = p m a x ∑ ∞ ( 1 − f ) 2 T p f p − 1 = 1 − T f ( 1 − f ) 2 T p m a x f p m a x − 1 方位散射
#
方位散射只需要考虑ϕ \phi ϕ y z yz yz
ϕ ( p , h ) = 2 p γ t − 2 γ o + p π
\begin{equation}
\phi(p,h)=2p\gamma_t-2\gamma_o+p\pi
\end{equation}
ϕ ( p , h ) = 2 p γ t − 2 γ o + p π 对于粗糙度, pbrt选用截尾Logistic分布, 与Gaussian分布相比Logistic具有相似的形状且具有积分解析式.
l t ( x , s , [ a , b ] ) = l ( x , s ) ∫ a b l ( x ′ , s ) d x ′
\begin{equation}
l_t(x,s,[a,b])=\frac{l(x,s)}{\int_a^b l(x',s)\mathrm{d}x'}
\end{equation}
l t ( x , s , [ a , b ]) = ∫ a b l ( x ′ , s ) d x ′ l ( x , s ) pbrt根据β n \beta_n β n s s s
s = π 8 ( 0.265 β n + 1.194 β n 2 + 5.372 β n 22 )
\begin{equation}
s=\sqrt{\frac{\pi}{8}}(0.265\beta_n+1.194\beta_n^2+5.372\beta_n^{22})
\end{equation}
s = 8 π ( 0.265 β n + 1.194 β n 2 + 5.372 β n 22 ) 散射模型求解
#
头发表面的凸起角度α \alpha α θ \theta θ R R R 2 α 2\alpha 2 α T T TT TT − α -\alpha − α T R T TRT TRT − 4 α -4\alpha − 4 α
可逆性
#
毛发BSDF是不可逆的, 这是由凸起方向以及A p A_p A p N p N_p N p
采样
#
首先以A p A_p A p p p p M p M_p M p N p N_p N p
毛发吸收系数
#
毛发的颜色由发皮质中的色素对光的吸收能力决定, 但毛发吸收系数与毛发颜色并不直接相关. 人类毛发中黑色素的浓度决定由金色到黑色的变化, 黄黑色素决定由橙色到红色的变化, pbrt中可以通过指定二者的浓度获取期望的发色.
PBRT_CPU_GPU RGBUnboundedSpectrum HairBxDF::SigmaAFromConcentration(Float ce, Float cp) {
RGB eumelaninSigma_a(0.419f, 0.697f, 1.37f);
RGB pheomelaninSigma_a(0.187f, 0.4f, 1.05f);
RGB sigma_a = ce * eumelaninSigma_a + cp * pheomelaninSigma_a;
#ifdef PBRT_IS_GPU_CODE
return RGBUnboundedSpectrum(*RGBColorSpace_sRGB, sigma_a);
#else
return RGBUnboundedSpectrum(*RGBColorSpace::sRGB, sigma_a);
#endif
}
