0 学习视频
B站GAMES101-现代计算机图形学入门-闫令琪
※ 接上文【CG】计算机图形学(Computer Graphics)基础(其壹)
7 光线追踪
7.1 为什么需要光线追踪?
-
光栅化无法妥善处理全局效果
-
(软)阴影
-
尤其是当光线反射不止一次时
-
-
光栅化很快(近似效果),但质量相对低
-
光线追踪很精准,但很慢
-
光栅化:实时,光线追踪:离线(电影)
-
一帧需要~10K处理器核心(CPU)小时渲染
-
7.2 基础光线追踪算法
7.2.1 光线Light Rays
三个关于光线的想法
- 光线延直线传播(虽然这是错的,波动性)
- 如果它们交叉,它们不会彼此“碰撞”(虽然这仍是错的,波动性)
- 光线从光源发出,不断反射到眼睛(但是物理学中是路径反转的——光路的可逆性)
“如果你凝视着深渊时,深渊也会凝视着你。”——Friedrich Wilhelm Nietzsche(尼采)
7.3 光线投射Ray Casting
Appel 1968
-
通过每个像素投射一个光线来生成图像
-
通过向光源发送光线检查阴影
7.3.1 光线投射——生成眼睛光线Eye Rays
针孔Pinhole相机模型
7.3.2 递归(Whitted-Style)光线追踪
“一个升级的用于着色显示照明模型” T.Whitted,CACM 1980
渲染时间:
- VAX 11/780(1979)74m
- PC(2006) 6s
- GPU(2012) 1/30s
左边球折射,右边球反射
7.4 射线表面交点
7.4.1 射线方程Ray Equation
射线由其起点和方向矢量定义
例子:
射线方程: r ( t ) = o + t d \textbf{r}(t)=\textbf{o}+t\textbf{d} r(t)=o+td 0 ≤ t < ∞ 0≤t<∞ 0≤t<∞
r \textbf{r} r:ray
t t t:time
o \textbf{o} o:origin
d \textbf{d} d:(normalized)direction
7.4.2 光线与球体的交点
射线: r ( t ) = o + t d \textbf{r}(t)=\textbf{o}+t\textbf{d} r(t)=o+td 0 ≤ t < ∞ 0≤t<∞ 0≤t<∞
球: p : ( p − c ) 2 − R 2 = 0 \textbf{p}:(\textbf{p}-\textbf{c})^2-R^2=0 p:(p−c)2−R2=0
什么是一个交点?
交点 P P P必须同时满足射线和球的方程
解方程:
( o + t d − c ) 2 − R 2 = 0 (\textbf{o}+t\textbf{d}-\textbf{c})^2-R^2=0 (o+td−c)2−R2=0
a t 2 + b t + c = 0 at^2+bt+c=0 at2+bt+c=0
当 a = d ⋅ d a=\textbf{d}\cdot \textbf{d} a=d⋅d, b = 2 ( o − c ) ⋅ d b=2(\textbf{o}-\textbf{c})\cdot\textbf{d} b=2(o−c)⋅d, c = ( o − c ) ⋅ ( o − c ) − R 2 c=(\textbf{o}-\textbf{c})\cdot(\textbf{o}-\textbf{c})-R^2 c=(o−c)⋅(o−c)−R2
t = − b ± b 2 − 4 a c 2 a t=\frac{-b±\sqrt{b^2-4ac}}{2a} t=2a−b±b2−4ac
7.4.3 光线与隐式表面求交
射线: r ( t ) = o + t d \textbf{r}(t)=\textbf{o}+t\textbf{d} r(t)=o+td 0 ≤ t < ∞ 0≤t<∞ 0≤t<∞
一般隐式表面: p : f ( p ) = 0 \textbf{p}:f(\textbf{p})=0 p:f(p)=0
替代射线方程: f ( o + t d ) = 0 f(\textbf{o}+t\textbf{d})=0 f(o+td)=0
解得正实数根
7.4.4 光线与三角形面(显式表面)求交
为什么?
- 渲染:可视化,阴影,光照…
- 几何:内/外部测试
- 任意一点与光源相连(射线)
- 奇数个交点:该点在内部
- 偶数个交点:该点在外部
- 任意一点与光源相连(射线)
怎么计算?
解决以下问题:
- 简单的想法:对每个三角形判断是否和光线有交点
- 简单,但是慢(如何加速?)
- 注意:具有0或1个交点(忽略多个交点)
7.4.5 光线和三角形的交点
三角形肯定处于一个平面内
- 射线与平面的交点
- 判断交点是否在三角形内部
7.4.6 平面方程Plane Equation
平面由法向量和平面上的一个点定义
例子:
平面方程(如果 p p p满足该方程,那么 p p p在该平面上)
p : ( p − p ′ ) ⋅ N = 0 \textbf{p}:(\textbf{p}-\textbf{p}^\prime)\cdot\textbf{N}=0 p:(p−p′)⋅N=0 a x + b y + c z + d = 0 ax+by+cz+d=0 ax+by+cz+d=0
7.4.7 光线与平面的交点
光线方程:
r ( t ) = o + t d \textbf{r}(t)=\textbf{o}+t\textbf{d} r(t)=o+td 0 ≤ t < ∞ 0≤t<∞ 0≤t<∞
平面方程:
p : ( p − p ′ ) ⋅ N = 0 \textbf{p}:(\textbf{p}-\textbf{p}^\prime)\cdot\textbf{N}=0 p:(p−p′)⋅N=0
求得交点
设 p = r ( t ) \textbf{p}=\textbf{r}(t) p=r(t)并解得 t t t
( p − p ′ ) ⋅ N = ( o + t d − p ′ ) ⋅ N = 0 (\textbf{p}-\textbf{p}^\prime)\cdot\textbf{N}=(\textbf{o}+t\textbf{d}-\textbf{p}^\prime)\cdot\textbf{N}=0 (p−p′)⋅N=(o+td−p′)⋅N=0
t = ( p ′ − o ) ⋅ N d ⋅ N ( 0 ≤ t < ∞ ) t=\frac{(\textbf{p}^\prime-\textbf{o})\cdot\textbf{N}}{\textbf{d}\cdot\textbf{N}}(0≤t<∞) t=d⋅N(p′−o)⋅N(0≤t<∞)
7.4.8 MÖller Trumbore算法
一种更快速的直接给出重心坐标的方法
讨论其中的部分推导!
7.5 加速光线与表面求交
7.5.1 光线追踪——性能挑战
7.5.1.1 简单的光线与场景的交点
- 用每个物体彻底测试与射线的交点
- 找到最近的交点(最小的t值)
7.5.1.2 问题
- 天真算法 = # p i x e l s × # o j e c t s ( × # b o u n c e s ) =\#pixels \times \#ojects(\times \#bounces) =#pixels×#ojects(×#bounces)
- 非常慢!
7.6 包围体积(盒)Bounding Volumes
7.6.1 包围盒Bounding Volumes
避免交点的快速方法:使用一个简单的盒子包围一个复杂的物体
- 物体被盒子完全地包含
- 如果它没有命中盒子,它就不会命中物体
- 所以先测试BVol,然后测试是否命中物体
7.6.2 用盒子Box的光线交点(3D)
理解:盒子Box是三个不同的对面的形成的交集
具体地:
通常使用一个轴对齐包围盒Axis-Aligned Bounding Box(AABB)
即任意包围盒的边都是沿着x,y或z轴的
7.6.3 光线与轴对齐盒子求交
2D的例子,3D中是一样的!计算光线与面的交点并取得 t m i n / t m a x t_{min}/t_{max} tmin/tmax间隔的交点
对x平面和y平面的线段求交集得到最终结果
-
回顾:一个盒子(3D)=三对无限大的平面
-
关键想法
- 只有当光线进入所有对面时光线才进入了盒子
- 只要光线离开任何一个对面,光线就离开的盒子
-
对于每一对对面,都计算计算tmin和tmax(负数也没关系)
-
对于3D盒子, t e n t e r = m a x { t m i n } , t e x i t = m i n { t m a x } t_{enter}=max\{t_{min}\},t_{exit}=min\{t_{max}\} tenter=max{tmin},texit=min{tmax}
-
如果 t e n t e r < t e x i t t_{enter}<t_{exit} tenter<texit,那么说光线在盒子中停留了一会儿(所以它们必然相交!)
-
然而,光线不是一条直线(是一条射线)
- 应该检查时间 t t t 是否是负数以满足物理上的正确性!
-
如果 t e x i t < 0 t_{exit}<0 texit<0呢?
- 盒子一定在光线的”背后“——不可能有交点!
-
如果 t e x i t ≥ 0 t_{exit}≥0 texit≥0 且 t e n t e r < 0 t_{enter}<0 tenter<0呢?
- 光线的起点在盒子中——一定有交点!
-
总结,光线与AABB有交点当且仅当(iff)
- t e n t e r < t e x i t & & t e x i t ≥ 0 t_{enter}<t_{exit}\&\&t_{exit}≥0 tenter<texit&&texit≥0
7.6.4 为什么使用轴对齐Axis-Aligned?
计算更简洁
7.7 使用AABBs加速光线追踪
7.7.1 均匀的格子Uniform grids
7.7.1.1 预处理-构建加速网格
-
找到包围盒
-
划分格子
-
判定哪些格子中可能有物体,存储起来
7.7.1.2 光线与场景求交(光线追踪)
逐步顺序地遍历光线经过的格子
对于每个格子
判断之前存储起来的格子(有物体)中的物体是否与光线相交
7.7.1.3 格子的分辨率?
一个格子
-
没有加速
很多格子
-
多次遍历,效率低下
启发式(技巧):
-
# c e l l s = C ∗ # o b j s \#cells=C*\#objs #cells=C∗#objs
-
C ≈ 27 i n 3 D C≈27\ in\ 3D C≈27 in 3D
7.7.1.4 均匀的格子——什么时候效果好?
大量物体在尺寸和空间上均匀分布时均匀的格子效果好
7.7.1.5 均匀的格子——什么时候效果不好?
“体育场中的茶壶”问题(空的地方多,不均匀)
7.7.2 空间划分Spatial partitions
7.7.2.1 空间划分的例子
7.7.2.2 KD树预处理
7.7.2.3 KD树的数据结构
内部结点存储
- 划分轴:x-,y-,或z-轴
- 划分位置:沿轴划分开的平面坐标
- 孩子:指向孩子结点
- 没有物体存储在内部结点中
叶子结点
- 物体列表
7.7.2.4 遍历一颗KD树
假设左边的蓝色区域是叶子结点:对所有物体求交
假设右上边的绿色区域是叶子结点:对所有物体求交
KD树的问题:
- 一个物体可能存在多个不同的格子里
- KD树的建立并不简单,需要考虑三角形与盒子求交
7.7.3 物体划分Object Partitions与包围盒层级Bounding Volume Hierarchy(BVH)
7.7.3.1 包围盒层级Bounding Volume Hierarchy(BVH)
- 找到一个包围盒
- 递归地把任何一个包围盒拆成两个部分
- 对拆成的两个部分都重新计算包围盒
- 必要时停下
- 把物体存储于各个叶子结点中
BVH解决了KD树的问题:
-
一个物体只可能出现在一个格子里
-
省去了包围盒与三角形求交的问题
BVH存在的问题:
BVH对空间的划分不是严格的划分开
怎么划分,很有讲究
7.7.3.2 构建BVHs
如何划分一个结点?
- 选择一个维度去划分
- 技巧1:总是选择沿着节点中最长的轴划分
- 技巧2:在中间物体的位置划分结点(涉及排序问题)
终止标准?
- 技巧:当结点中包含几个元素时停止(例如5个)
7.7.3.2.1 任意一个无序序列,如何在O(n)时间快速找到其中位数(或第i大的树)?
快速选择算法,受快速排序算法启发
7.7.3.3 BVHs的数据结构
内部结点存储
- 包围盒
- 孩子:指向孩子结点
叶子结点存储
- 包围盒
- 物体列表
结点代表场景中基元的子集
- 子树中的所有物体
7.7.3.4 BVH遍历算法
7.7.3.5 空间划分vs物体划分
空间划分(如KD树)
- 把空间划分为两个不重叠的区域
- 一个物体可能处于多个区域中
物体划分(如BVH)
- 把物体划分到不相交的子集中
- 每个子集的包围盒在空间中可能相交
7.8 辐射度量学Basic radiometry
7.8.1 度量学Radiometry
- 是一个测量系统,也是一个光照单位
- 准确测量光的空间特性
- 辐射通量Radiant flux,强度intensity,辐照度irradiance,辐亮度radiance
- 准确定义计算光照的物理表达方法
学习的方式
WHY,WHAT,then HOW
7.8.2 辐射能量和通量Radiant Energy and Flux(Power)
定义:辐射能量(Radiant energy)是电磁能量辐射。以焦耳为能量测量,并表示为符号 Q Q Q.
Q [ J = J o u l e ] Q\ [J=Joule] Q [J=Joule]
定义:辐射通量(Radiant Flux(Power))是单位时间发出、反射、传输或接收的能量。
Φ = d Q d t [ W = W a t t ] [ l m = l u m e n ] ∗ \Phi=\frac{dQ}{dt}\ [W=Watt]\ [lm=lumen]^* Φ=dtdQ [W=Watt] [lm=lumen]∗
7.8.2.1 Flux-单位时间内流过传感器的#photons
7.8.2.2 感兴趣的重要光线度量
辐射强度Radiant Intensity
辐照度Irradiance
辐亮度Radiance
7.8.3 辐射强度Radiant Intensity
定义:辐射(发光)强度Radiant(luminous) Intensity是单位能量从一个点光源发出的立体角solid angle。
坎德拉(candela)是七个SI基本单位之一
7.8.3.1 角度和立体角
角:圆上一个角所对弧长与半径之比
-
θ = l r \theta=\frac{l}{r} θ=rl
-
圆有 2 π 2\pi 2π的弧度radians
立体角:球上一个角所对面积与球径平方之比
- Ω = A r 2 \Omega=\frac{A}{r^2} Ω=r2A
- 球有
4
π
4\pi
4π的立体角steradians
7.8.3.2 微分立体角Differential Solid Angles
d
A
=
(
r
d
θ
)
(
r
sin
θ
d
ϕ
)
=
r
2
sin
θ
d
θ
d
ϕ
dA=(r\ d\theta)(r\ \sin\theta\ d\phi)=r^2\ \sin\theta\ d\theta\ d\phi
dA=(r dθ)(r sinθ dϕ)=r2 sinθ dθ dϕ
d ω = d A r 2 = s i n θ d θ d ϕ d\omega=\frac{dA}{r^2}=sin\theta\ d\theta\ d\phi dω=r2dA=sinθ dθ dϕ
Sphere: S 2 \textbf{Sphere:}\ \ S^2 Sphere: S2
Ω = ∫ S 2 d ω = ∫ 0 2 π ∫ 0 π sin θ d θ d ϕ \Omega=\int_{S^{2}}d\omega=\int_0^{2\pi}\int_0^{\pi}\sin\theta\ d\theta\ d\phi Ω=∫S2dω=∫02π∫0πsinθ dθ dϕ
7.8.3.3 ω \omega ω 作为一个方向向量
将使用 ω \omega ω 表示一个方向向量(单位长度)
7.8.3.4 各向同性点源Isotropic Point Source
Φ = ∫ S 2 I d ω = 4 π I \Phi=\int_{S^2}I\ d\omega=4\pi I Φ=∫S2I dω=4πI
I = Φ 4 π I=\frac{\Phi}{4\pi} I=4πΦ
7.8.3.5 现代的LED灯
输出:815 lumens(11W的LED相当于60W的白炽灯)
Radiant intensity?
假设各向同性(每个方向均匀):
I n t e n s i t y = 815 l u m e n s / 4 p i s r = 65 c a n d e l a s Intensity=815\ lumens / 4pi\ sr=65\ candelas Intensity=815 lumens/4pi sr=65 candelas
7.8.4 辐照度Irradiance
定义:入射到表面点上的每(垂直/投影)单位面积的能量。
E ( x ) ≡ d Φ ( x ) d A E(\textbf{x})\equiv\frac{d\Phi(\textbf{x})}{dA} E(x)≡dAdΦ(x)
[ W m 2 ] = [ l m m 2 = l u x ] [\frac{W}{m^2}]=[\frac{lm}{m^2}=lux] [m2W]=[m2lm=lux]
7.8.4.1 为什么地球有四季之分?
地球自转轴:~23.5°偏轴
7.8.5 辐亮度Radiance
辐射亮度是描述环境中光线分布的基本场量。
- 辐亮度是与光线有关的量。
- 渲染就是计算辐亮度。
定义:辐亮度(亮度)是每单位立体角和每单位投影面积上,由表面反射、发射或接收的能量。
回顾:
- 辐照度Irradiance:每单位投影面积的能量
- 辐射强度Intensity:每立体角的能量
结论:
- 辐亮度Radiance:每立体角的辐照度Irradiance
- 辐亮度Radiance:每单位投影面积的辐射强度Intensity
7.8.5.1 入射辐亮度Incident Radiance
入射辐亮度(Incident Radiance)是到达表面的每单位立体角的辐照度(Irradiance)。例如,光沿着给定的光线到达表面(表面和入射方向)。
7.8.5.2 出射辐亮度Exiting Radiance
出射辐亮度(Incident Radiance)是离开表面的单位投影面积的辐射强度(Intensity)。
例如,对于一个面光(area light),它是沿着给定光线(表面和出射方向上的点)发出的光。
7.8.6 辐照度Irradiance vs. 辐亮度Radiance
辐照度Irradiance:面积 d A dA dA 接收到的总能量。
辐亮度Radiance:面积 d A dA dA 从“方向” d ω d\omega dω 接收到的能量。
7.8.7 双向反射分布函数Bidirectional Reflectance Distribution Function(BRDF)
7.8.7.1 一个点上的反射
来自方向
ω
i
\omega_i
ωi 的光线(Radiance),转换为
d
A
dA
dA 接收的能量
E
E
E。而能量
E
E
E 将转换为对任何其他方向
ω
\omega
ω 的光线(Radiance)。
差分入射辐照度 Differential irradiance incoming:
d E ( ω i ) = L ( ω i ) cos θ i d ω i dE(\omega_i)=L(\omega_i)\cos\theta_id\omega_i dE(ωi)=L(ωi)cosθidωi
差分出射辐亮度 Differential radiance exiting(由 d E ( ω i ) dE(\omega_i) dE(ωi)产生):
d L r ( ω r ) dL_r(\omega_r) dLr(ωr)
7.8.7.2 BRDF
双向反射率分布函数(BRDF)表示有多少光从每个入射方向反射到每个出射方向 ω r \omega_r ωr。
f
r
(
ω
i
→
ω
r
)
=
d
L
r
(
ω
r
)
d
E
i
(
ω
r
)
=
d
L
r
(
ω
r
)
L
(
ω
i
)
cos
θ
i
d
ω
i
[
1
s
r
]
f_r(\omega_i→\omega_r)=\frac{dL_r(\omega_r)}{dE_i(\omega_r)}=\frac{dL_r(\omega_r)}{L(\omega_i)\cos\theta_id\omega_i}[\frac{1}{sr}]
fr(ωi→ωr)=dEi(ωr)dLr(ωr)=L(ωi)cosθidωidLr(ωr)[sr1]
7.8.7.3 反射方程The Reflection Equation
L r ( p , ω r ) = ∫ H 2 f r ( p , ω i → ω r ) L i ( p , ω i ) cos θ i d ω i L_r(p,\omega_r)=\int_{H^2}f_r(p,\omega_i→\omega_r)L_i(p,\omega_i)\cos\theta_id\omega_i Lr(p,ωr)=∫H2fr(p,ωi→ωr)Li(p,ωi)cosθidωi
7.8.7.4 困难:递归的表达Recursive Equation
反射方程中,反射光取决于入射光:
但是,入射光也取决于反射光(在场景中的另一点)(光线不断弹射)
7.8.7.5 渲染方程The Renaering Equation
重写反射方程:
L r ( p , ω r ) = ∫ H 2 f r ( p , ω i → ω r ) L i ( p , ω i ) cos θ i d ω i L_r(p,\omega_r)=\int_{H^2}f_r(p,\omega_i→\omega_r)L_i(p,\omega_i)\cos\theta_id\omega_i Lr(p,ωr)=∫H2fr(p,ωi→ωr)Li(p,ωi)cosθidωi
通过添加一个Emission term(自发光项)使其通用!(考虑到自发光的物体或是接受到其它光线的入射光)
L o ( p , ω o ) = L e ( p , ω o ) + ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i L_o(p,\omega_o)=L_e(p,\omega_o)+\int_{\Omega+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)d\omega_i Lo(p,ωo)=Le(p,ωo)+∫Ω+Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)dωi
注,假设所有方向都是朝外的!
7.8.7.6 理解渲染方程
- 反射方程
- 一个点光源
-
L r ( x , ω r ) = L e ( x , ω r ) + L i ( x , ω i ) f ( x , ω i , ω r ) ( ω i , n ) L_r(x,\omega_r)=L_e(x,\omega_r)+L_i(x,\omega_i)f(x,\omega_i,\omega_r)(\omega_i,n) Lr(x,ωr)=Le(x,ωr)+Li(x,ωi)f(x,ωi,ωr)(ωi,n)
- L r ( x , ω r ) L_r(x,\omega_r) Lr(x,ωr):反射光(输出图像)
- L e ( x , ω r ) L_e(x,\omega_r) Le(x,ωr):自发光的物体或是接受到其它光线的入射光
- L i ( x , ω i ) L_i(x,\omega_i) Li(x,ωi):入射光(来自光源)
- f ( x , ω i , ω r ) f(x,\omega_i,\omega_r) f(x,ωi,ωr):BRDF
- ( ω i , n ) (\omega_i,n) (ωi,n):入射角余弦值
-
多个点光源
L
r
(
x
,
ω
r
)
=
L
e
(
x
,
ω
r
)
+
Σ
L
i
(
x
,
ω
i
)
f
(
x
,
ω
i
,
ω
r
)
(
ω
i
,
n
)
L_r(x,\omega_r)=L_e(x,\omega_r)+\Sigma L_i(x,\omega_i)f(x,\omega_i,\omega_r)(\omega_i,n)
Lr(x,ωr)=Le(x,ωr)+ΣLi(x,ωi)f(x,ωi,ωr)(ωi,n)
- 多个点光源和一个面光源
L
r
(
x
,
ω
r
)
=
L
e
(
x
,
ω
r
)
+
∫
Ω
L
i
(
x
,
ω
i
)
f
(
x
,
ω
i
,
ω
r
)
cos
θ
i
d
ω
i
L_r(x,\omega_r)=L_e(x,\omega_r)+\int_{\Omega}L_i(x,\omega_i)f(x,\omega_i,\omega_r)\cos\theta_id\omega_i
Lr(x,ωr)=Le(x,ωr)+∫ΩLi(x,ωi)f(x,ωi,ωr)cosθidωi
- 存在其它互反射的表面
-
L r ( X , ω r ) = L e ( X , ω r ) + ∫ Ω L r ( X ′ , − ω i ) f ( X , ω i , ω r ) cos θ i d ω i L_r(X,\omega_r)=L_e(X,\omega_r)+\int_{\Omega}L_r(X^\prime,-\omega_i)f(X,\omega_i,\omega_r)\cos\theta_id\omega_i Lr(X,ωr)=Le(X,ωr)+∫ΩLr(X′,−ωi)f(X,ωi,ωr)cosθidωi
- L r ( X , ω r ) L_r(X,\omega_r) Lr(X,ωr):未知
- L e ( X , ω r ) L_e(X,\omega_r) Le(X,ωr):已知
- ∫ Ω L r ( X ′ , − ω i ) \int_{\Omega}L_r(X^\prime,-\omega_i) ∫ΩLr(X′,−ωi):未知
- f ( X , ω i , ω r ) f(X,\omega_i,\omega_r) f(X,ωi,ωr):已知
- cos θ i d ω i \cos\theta_id\omega_i cosθidωi:已知
简化上式:
-
I
(
u
)
=
e
(
u
)
+
∫
I
(
v
)
K
(
u
,
v
)
d
v
I(u)=e(u)+\int I(v)K(u,v)dv
I(u)=e(u)+∫I(v)K(u,v)dv
-
K
(
u
,
v
)
d
v
K(u,v)dv
K(u,v)dv:等式核心,光传输操作符
再简化:
L = E + K L L=E+KL L=E+KL
可以离散化为简单的矩阵方程(或联立线性方程组)(L、E为向量,K为光传输矩阵)
-
K
(
u
,
v
)
d
v
K(u,v)dv
K(u,v)dv:等式核心,光传输操作符
7.8.7.7 光线追踪和扩展
- 广义类数值蒙特卡罗方法
- 场景中所有光路近似集
L = E + K L L=E+KL L=E+KL
L
=
E
+
K
L
I
L
−
K
L
=
E
(
I
−
K
)
L
=
E
L
=
(
I
−
K
)
−
1
E
二项式定理
L
=
(
I
+
K
+
K
2
+
K
3
+
…
)
E
L
=
E
+
K
E
+
K
2
E
+
K
3
E
+
…
\begin{aligned} L & =E+KL \\ IL-KL & = E \\ (I-K)L&=E\\ L&=(I-K)^{-1}E\\ &二项式定理\\ L&=(I+K+K^2+K^3+…)E\\ L&=E+KE+K^2E+K^3E+… \end{aligned}
LIL−KL(I−K)LLLL=E+KL=E=E=(I−K)−1E二项式定理=(I+K+K2+K3+…)E=E+KE+K2E+K3E+…
其中
I
I
I为单位矩阵 。
将上式各项分开理解(光的多次弹射):
另一种理解方式(光栅化的着色):
7.9 概率论
7.9.1 随机变量
随机变量 X X X 表示随机实验结果的可能数值。
概率分布函数 X ∼ p ( x ) X\sim p(x) X∼p(x) 表示随机过程选择值 x x x 的相对概率。
例如,等概率分布函数,在取值范围内的所有值都具有相同的概率。
一个六面的骰子🎲, X X X的取值可能是1,2,3,4,5,6。
p ( 1 ) = p ( 2 ) = p ( 3 ) = p ( 4 ) = p ( 5 ) = p ( 6 ) = 1 6 p(1)=p(2)=p(3)=p(4)=p(5)=p(6)=\frac{1}{6} p(1)=p(2)=p(3)=p(4)=p(5)=p(6)=61
7.9.2 概率
x i x_i xi表示 n n n个离散值
p i p_i pi表示 x i x_i xi 出现的概率
概率分布的特点:
p i ≥ 0 p_i≥0 pi≥0
Σ i = 1 n p i = 1 \Sigma_{i=1}^np_i=1 Σi=1npi=1
对于等概率的六面骰子🎲,其每个面出现的概率是 p = 1 6 p=\frac{1}{6} p=61
7.9.3 随机变量的期望值
随机变量的期望值即从随机分布中反复抽取样本所得到的平均值。
x i x_i xi表示 n n n个离散值
p i p_i pi表示 x i x_i xi 出现的概率
随机变量 X X X的期望值:
E [ X ] = Σ i = 1 n x i p i E[X]=\Sigma_{i=1}^nx_ip_i E[X]=Σi=1nxipi
对于等概率的六面骰子🎲,其期望值是:
E [ X ] = Σ i = 1 n i 6 = ( 1 + 2 + 3 + 4 + 5 + 6 ) / 6 = 3.5 E[X]=\Sigma_{i=1}^n\frac{i}{6}=(1+2+3+4+5+6)/6=3.5 E[X]=Σi=1n6i=(1+2+3+4+5+6)/6=3.5
7.9.4 连续情况下的概率密度函数
X
∼
p
(
x
)
X\sim p(x)
X∼p(x)
对于一个随机变量
X
X
X,可以取一组连续值中的任意一个,其中一个特定值的相对概率由连续概率分布函数
p
(
x
)
p(x)
p(x) 给出。
p ( x ) p(x) p(x)的条件: p ( x ) > 0 p(x)>0 p(x)>0 且 ∫ p ( x ) d x = 1 \int p(x)dx=1 ∫p(x)dx=1
X X X 的期望值: E [ X ] = ∫ x p ( x ) d x E[X]=\int xp(x)dx E[X]=∫xp(x)dx
7.9.5 随机变量函数
随机变量 X X X 的函数 Y Y Y 也是随机变量
X ∼ p ( x ) X\sim p(x) X∼p(x)
Y = f ( X ) Y=f(X) Y=f(X)
Y Y Y的期望值:
E [ Y ] = E [ f ( X ) ] = ∫ f ( x ) p ( x ) d x E[Y]=E[f(X)]=\int f(x)p(x)dx E[Y]=E[f(X)]=∫f(x)p(x)dx
7.10 蒙特卡洛积分Monte Carlo Integration
7.10.1 蒙特卡洛积分Monte Carlo Integration
W
h
y
?
Why?
Why?我们想解一个积分,但是用分析的方法解它可能太难了
W h a t & H o w ? What\&How? What&How?通过对函数值的随机样本取平均值来估计函数的积分
接下来定义给定函数 f ( x ) f(x) f(x) 的定积分的蒙特卡洛估计量:
定积分: ∫ a b f ( x ) d x \int_a^bf(x)dx ∫abf(x)dx
随机变量: X i ∼ p ( x ) X_i\sim p(x) Xi∼p(x)
蒙特卡洛估计量: F N = 1 N Σ i = 1 N f ( X i ) p ( X i ) F_N=\frac{1}{N}\Sigma_{i=1}^N\frac{f(X_i)}{p(X_i)} FN=N1Σi=1Np(Xi)f(Xi)
一些注意事项:
- 样本越多,变化就越小。
- 在x上采样,在x上积分。
7.10.2 例子:等值蒙特卡洛估计量
等值随机变量: X i ∼ p ( x ) = C ( c o n s t a n t ) X_i\sim p(x)=C(constant) Xi∼p(x)=C(constant)
∫ a b p ( x ) d x = 1 \int_a^bp(x)dx=1 ∫abp(x)dx=1
⇒ ∫ a b C d x = 1 ⇒\int_a^bCdx=1 ⇒∫abCdx=1
⇒ C = 1 b − a ⇒C=\frac{1}{b-a} ⇒C=b−a1
接下来定义给定函数 f ( x ) f(x) f(x) 的定积分的蒙特卡洛估计量:
定积分: ∫ a b f ( x ) d x \int_a^bf(x)dx ∫abf(x)dx
等值随机变量: X i ∼ p ( x ) = 1 b − a X_i\sim p(x)=\frac{1}{b-a} Xi∼p(x)=b−a1
基础蒙特卡洛估计量: F N = b − a N Σ i = 1 N f ( X i ) F_N=\frac{b-a}{N}\Sigma_{i=1}^Nf(X_i) FN=Nb−aΣi=1Nf(Xi)
7.11 路径追踪Path Tracing
用来改进Whitted-Style 光线追踪
7.11.1 Whitted-Style 光线追踪的问题1
对于有光泽的材料,光线应该反射到哪里?
7.11.2 Whitted-Style 光线追踪的问题2
漫反射材质之间没有反射?
7.11.3 Whitted-Style 光线追踪是错的
渲染方程才是正确的
L o ( p , ω o ) = L e ( p , ω o ) + ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i L_o(p,\omega_o)=L_e(p,\omega_o)+\int_{\Omega+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)d\omega_i Lo(p,ωo)=Le(p,ωo)+∫Ω+Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)dωi
但这牵涉到
- 求解半球上的积分
- 递归的式子
如何用数值方法解一个积分?
7.11.4 一个简单的蒙特卡洛解法
假设我们想在下面的场景中渲染一个像素(点),只需要求解该点的直接光照:
这里忽略物体的发光,即使用反射方程:
L o ( p , ω o ) = ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i L_o(p,\omega_o)=\int_{\Omega+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)d\omega_i Lo(p,ωo)=∫Ω+Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)dωi
可以用蒙特卡洛积分来解决!
我们要计算 p p p点反射到相机的辐亮度(Radiance)
L o ( p , ω o ) = ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i L_o(p,\omega_o)=\int_{\Omega+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)d\omega_i Lo(p,ωo)=∫Ω+Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)dωi
蒙特卡洛估计量: ∫ a b f ( x ) d x ≈ 1 N Σ k = 1 N f ( X k ) p ( X k ) X k ∼ p ( x ) \int_a^bf(x)dx≈\frac{1}{N}\Sigma_{k=1}^N\frac{f(X_k)}{p(X_k)}\ \ \ \ X_k\sim p(x) ∫abf(x)dx≈N1Σk=1Np(Xk)f(Xk) Xk∼p(x)
f ( x ) : L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) f(x):L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i) f(x):Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)
概率密度函数 ( P D F ) : p ( ω i ) = 1 / 2 π 概率密度函数(PDF):p(\omega_i)=1/2\pi 概率密度函数(PDF):p(ωi)=1/2π (假设均匀采样半球)
因此,可以得到:
$$
L o ( p , ω o ) = ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i ≈ 1 N Σ i = 1 N L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) p ( ω i ) \begin{aligned} L_o(p,\omega_o)&=\int_{\Omega+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)d\omega_i\\ &≈\frac{1}{N}\Sigma_{i=1}^N\frac{L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)}{p(\omega_i)} \end{aligned} Lo(p,ωo)=∫Ω+Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)dωi≈N1Σi=1Np(ωi)Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)
这意味着,可以得到一个正确的直接光照着色算法!
7.11.5 全局光照Global lllumination
如果光线击中物体怎么办?
Q点 也反射光线到 P点!是多少?即 Q点 的直接光照!
得到支持全局光照的路径追踪算法:
7.11.6 路径追踪存在的问题1
随着 # b o u n c e \#bounce #bounce(反弹)次数的增加, # r a y \#ray #ray(光线)数量会爆炸:
# r a y s = N # b o u n c e \#rays=N^{\#bounce} #rays=N#bounce
关键点:当且仅当 N = ? ? ? N=??? N=???时, # r a y \#ray #ray(光线)数量不会爆炸
因此,从现在开始,我们总是假设在每个着色点只有1条光线被跟踪:
以上算法,
N
=
1
N=1
N=1,即路径追踪。
而,当 N ! = 1 N!=1 N!=1时,即分布式光线追踪。
7.11.7 光线生成
N = 1 N=1 N=1会造成噪声很大的问题!
但这没问题,只是通过每个像素追踪更多的路径(path),并平均他们的辐照度(Radiance)!
光线生成非常类似于光线跟踪中的光线投射:
7.11.8 路径追踪存在的问题2
这是一个递归的算法,没有终止条件,则会一直进行下去:
真实世界中,光线确实就是不停得弹射的!
限制光线弹射次数,则会造成能量损失:
Cutting # b o u n c e s \#bounces #bounces == cutting energy!
7.11.9 解决方法:俄罗斯轮盘赌ussian Roulette(RR)
俄罗斯轮盘赌的输赢全靠概率:
- 如果概率为 0 < P < 1 0<P<1 0<P<1,你就没事
- 否则概率为 1 − P 1-P 1−P
左轮手枪里能装6枚子弹,如果在其中装入2枚子弹,那么你的存活概率就为 P = 4 / 6 P=4/6 P=4/6。
以前,我们总是在一个着色点发射光线,并得到着色结果 L o \textbf{L}_\textbf{o} Lo
假设我们手动设置一个概率 P ( 0 < P < 1 ) P(0<P<1) P(0<P<1)
- 在概率为 P P P 的情况下,发射一条光线,返回着色结果除以 P P P: L o / P \textbf{L}_\textbf{o}/P Lo/P
- 在概率为 1 − P 1-P 1−P 的情况下,不发射光线,你会得到 0 0 0
这样做,你仍然能够期望得到 L o \textbf{L}_\textbf{o} Lo:
E = P ∗ ( L o / P ) + ( 1 − P ) ∗ 0 = L o E=P*(\textbf{L}_\textbf{o}/P)+(1-P)*0=\textbf{L}_\textbf{o} E=P∗(Lo/P)+(1−P)∗0=Lo
最终得到的路径追踪算法:
7.11.10 正确的路径追踪
现在,我们终于得到了一个正确的路径追踪!
但实际上,它的效率并不高。
7.11.11 光线采样
理解路径追踪效率不高的原因
总会有1束光线打在灯光上。因此,如果我们在着色点均匀地对半球进行采样就会浪费大量的光线。
7.11.12 光线采样(纯数学)
蒙特卡洛方法允许任意采样方式,因此我们可以对光线进行采样(因此没有光线被"浪费")。
假设在光线上均匀采样:
p d f = 1 / A pdf=1/A pdf=1/A(因为 ∫ p f d A = 1 \int pfdA=1 ∫pfdA=1)
但是渲染方程是在立体角上的积分:
L o = ∫ L i f r cos d ω Lo=\int L_ifr \cos d\omega Lo=∫Lifrcosdω
回顾蒙特卡罗积分:
在x上采样,并在x上积分
因此,需要改写渲染方程成对 d A dA dA 的积分,而不是对 d ω d\omega dω的积分
需要 d ω d\omega dω 和 d A dA dA 之间的关系
回顾:
立体角定义:单位球面上的投影面积
d ω = d A cos θ ′ ∥ x ′ − x ∥ 2 d\omega=\frac{dA \cos\theta'}{\|x'-x\|^2} dω=∥x′−x∥2dAcosθ′
接着,对渲染方程重写:
L o ( x , ω o ) = ∫ Ω + L i ( x , ω i ) f r ( x , ω i , ω o ) cos θ d ω i = ∫ A L i ( x , ω i ) f r ( x , ω i , ω o ) cos θ cos θ ′ ∥ x ′ − x ∥ 2 d A \begin{aligned} L_o(x,\omega_o)& =\int_{\Omega^+}L_i(x,\omega_i)f_r(x,\omega_i,\omega_o)\cos\theta \mathrm{d}\omega_i \\ &=\int_{A}L_{i}(x,\omega_{i})f_{r}(x,\omega_{i},\omega_{o})\frac{\cos\theta\cos\theta^{\prime}}{\|x^{\prime}-x\|^{2}} \mathrm{d}A \end{aligned} Lo(x,ωo)=∫Ω+Li(x,ωi)fr(x,ωi,ωo)cosθdωi=∫ALi(x,ωi)fr(x,ωi,ωo)∥x′−x∥2cosθcosθ′dA
之前,我们假设光是通过均匀半球取样“意外”发射的。
现在我们考虑来自两个部分的radiance:
- 光源(直接贡献,不需要有RR)
- 其他非光源物体的反射(间接反射,RR)
最终算法如下:
最后一个问题:我们如何知道光上的采样点是否被遮挡?
路径追踪(Path tracing)是100%正确的,它可以做到照片级的真实感!
7.11.13 光线追踪:早期概念 vs. 现代概念
- 从前
光线跟踪(Ray tracing) == Whitted-style ray tracing - 现代
- 所有光线传播的大集合,包括
- (单向/双向)路径跟踪(Unidirectional & bidirectiona)path tracing
- 光子映射(Photon mapping)
- Metropolis光线传输( light transport)
- VCM/UPBP
- 所有光线传播的大集合,包括
7.11.14 未提及的问题
- 半球均匀采样
- 怎么做?一般情况下,如何对任何函数进行采样?(采样理论)
- 蒙特卡罗积分能够用在任意的 pdfs 上
- 怎么样是最好的选择?(重要性采样理论)
- 随机数是否有好坏之分?
- 是的(低差异化序列)
- 我可以采样半球和光线
- 我能把它们结合起来吗?YES!(multiple imp. sampling)
- 一个像素的radiance是通过它的所有路径上的radiance的平均值
- 为什么(像素重建滤波器 pixel reconstruction filter)
- 像素的radiance是像素的颜色吗?
- 不是(需要经过伽马校正、曲线、颜色空间)
8 材质与外观
8.1 材质
8.1.1 自然界不同材质的外观
不同的材质在不同的光照下就会呈现出不同的外观。
8.1.2 什么是计算机图形学中的材质?
8.1.3 材质==BRDF
8.1.3.1 漫反射(Diffuse/Lambertian)材质(BRDF)
漫反射:光在每个出射方向上被均匀地反射
假设入射光也是均匀的(能量守恒):
L o ( ω o ) = ∫ H 2 f r L i ( ω i ) cos θ i d ω i = f r L i ∫ H 2 ( ω i ) cos θ i d ω i = π f r L i f r = ρ π \begin{aligned} &L_o(\omega_o)&& =\int_{H^2}f_rL_i(\omega_i)\cos\theta_i \mathrm{d}\omega_i \\ &&&=f_rL_i\int_{H^2}\sout{(\omega_i)}\cos\theta_i \mathrm{d}\omega_i \\ &&&=\pi f_rL_i \\ &f_r=\frac\rho\pi \end{aligned} Lo(ωo)fr=πρ=∫H2frLi(ωi)cosθidωi=frLi∫H2(ωi)cosθidωi=πfrLi
这就是完全不吸收能量的BRDF。
其中 ρ \rho ρ为反射率(albedo)color
8.1.3.2 抛光的金属(Glossy)材质(BRDF)
8.1.3.3 理想反射/折射(ldeal reflective / refractive)材质(BSDF*)
8.1.3.4 完美镜面反射(Perfect Specular Reflection)
ω o + ω i = 2 cos θ n ⃗ = 2 ( ω i ⋅ n ⃗ ) n ⃗ ω o = − ω i + 2 ( ω i ⋅ n ⃗ ) n ⃗ \begin{aligned}&\omega_{o}+\omega_{i}=2\cos\theta \vec{\mathrm{n}}=2(\omega_{i}\cdot\vec{\mathrm{n}})\vec{\mathrm{n}}\\&\omega_{o}=-\omega_{i}+2(\omega_{i}\cdot\vec{\mathrm{n}})\vec{\mathrm{n}}\end{aligned} ωo+ωi=2cosθn=2(ωi⋅n)nωo=−ωi+2(ωi⋅n)n
8.1.3.5 完美镜面反射BRDF
8.1.3.6 完美折射Specular Refraction
除了从表面反射之外,光还可以透过表面。
光进入新的介质时会发生折射。
8.1.3.7 斯奈尔定律Snell’s Law
透射角取决于
- 入射光线的折射率(IOR)
- 出射光线的折射率(IOR)
η
i
sin
θ
i
=
η
t
sin
θ
t
\eta_i\sin\theta_i=\eta_t\sin\theta_t
ηisinθi=ηtsinθt
折射率
η
∗
\eta^*
η∗:
8.1.3.8 折射定律Law of Refraction
η
i
sin
θ
i
=
η
t
sin
θ
t
cos
θ
t
=
1
−
sin
2
θ
t
=
1
−
(
η
i
η
t
)
2
sin
2
θ
i
=
1
−
(
η
i
η
t
)
2
(
1
−
cos
2
θ
i
)
\begin{aligned} \eta_i\sin\theta_i& =\eta_t\sin\theta_t \\ \cos\theta_t& =\sqrt{1-\sin^2\theta_t} \\ &=\sqrt{1-\left(\frac{\eta_i}{\eta_t}\right)^2\sin^2\theta_i} \\ &=\sqrt{1-\left(\frac{\eta_i}{\eta_t}\right)^2(1-\cos^2\theta_i)} \end{aligned}
ηisinθicosθt=ηtsinθt=1−sin2θt=1−(ηtηi)2sin2θi=1−(ηtηi)2(1−cos2θi)
存在 1 − ( η i η t ) 2 ( 1 − cos 2 θ i ) < 0 1-\left(\frac{\eta_i}{\eta_t}\right)^2(1-\cos^2\theta_i)<0 1−(ηtηi)2(1−cos2θi)<0的情况
全反射Total internal reflection:
当光从光学密度较高的介质向光学密度较低的介质移动时:
η i η t > 1 \frac{\eta_i}{\eta_t}>1 ηtηi>1
从足够大的角度入射到边界上的光线不会离开介质(没有折射的情况)
8.1.3.9 斯奈尔之窗Snell’s Window /Circle
8.1.3.10 菲涅尔项Fresnel Reflection /Term
反射率取决于入射角(和光的偏振)
本例:反射率随掠射角(grazing angle)的增加而增加
(掠射角(grazing angle):光线几乎是平着打到表面上)
菲涅耳项(电介质, η \eta η = 1.5)
菲涅尔项(导体)
菲涅尔项——公式
精确计算,需要考虑光的极化
R s = ∣ n 1 cos θ 1 − n 2 cos θ t n 1 cos θ 1 + 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 ∣ 2 . R_{\mathrm{s}}=\left|\frac{n_1\cos\theta_1-n_2\cos\theta_\mathrm{t}}{n_1\cos\theta_1+n_2\cos\theta_\mathrm{t}}\right|^2=\left|\frac{n_1\cos\theta_\mathrm{i}-n_2\sqrt{1-\left(\frac{n_1}{n_2}\sin\theta_\mathrm{i}\right)^2}}{n_1\cos\theta_\mathrm{i}+n_2\sqrt{1-\left(\frac{n_1}{n_2}\sin\theta_\mathrm{i}\right)^2}}\right|^2,\\R_{\mathrm{p}}=\left|\frac{n_1\cos\theta_\mathrm{t}-n_2\cos\theta_\mathrm{i}}{n_1\cos\theta_\mathrm{t}+n_2\cos\theta_\mathrm{i}}\right|^2=\left|\frac{n_1\sqrt{1-\left(\frac{n_1}{n_2}\sin\theta_\mathrm{i}\right)^2}-n_2\cos\theta_\mathrm{i}}{n_1\sqrt{1-\left(\frac{n_1}{n_2}\sin\theta_\mathrm{i}\right)^2}+n_2\cos\theta_\mathrm{i}}\right|^2. Rs= n1cosθ1+n2cosθtn1cosθ1−n2cosθt 2= n1cosθi+n21−(n2n1sinθi)2n1cosθi−n21−(n2n1sinθi)2 2,Rp= n1cosθt+n2cosθin1cosθt−n2cosθi 2= n11−(n2n1sinθi)2+n2cosθin11−(n2n1sinθi)2−n2cosθi 2.
R e f f = 1 2 ( R s + R p ) R_{\mathrm{eff}}=\frac12\left(R_{\mathrm{s}}+R_{\mathrm{p}}\right) 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 \begin{aligned}R(\theta)&=R_0+(1-R_0)(1-\cos\theta)^5\\R_0&=\left(\frac{n_1-n_2}{n_1+n_2}\right)^2\end{aligned} R(θ)R0=R0+(1−R0)(1−cosθ)5=(n1+n2n1−n2)2
8.1.4 微表面材质Microfacet Material
8.1.4.1 微表面理论Microfacet Theory
-
观察一个粗糙的表面
- 宏观尺度:表面平坦 & 材质粗糙(材质&外观)
- 微观尺度:表面凹凸 & 材质镜像(几何)
-
表面的单个元素就像镜子一样
- 被称为微表面
- 每个微表面都有自己的朝向(法线)
8.1.4.2 微表面(Microfacet)材质(BRDF)
- 关键:微表面的法线分布
- 表面光滑,所有微平面的朝向都可看为是朝上的,可转换为抛光的金属(glossy)材质
- 表面粗糙,所有微平面的朝向四散而开,可转换为漫反射(diffuse)材质
- 表面光滑,所有微平面的朝向都可看为是朝上的,可转换为抛光的金属(glossy)材质
因此,通过微表面模型,可以把表面的粗糙程度用微表面的法线分布来表示。
- 什么样的微表面将
w
i
w_i
wi 反射到
w
o
w_o
wo?(提示:所有微表面都是镜子)
8.1.4.3 微表面材质例子
PBR中一定会使用到微表面模型
8.1.5 各向同性/各向异性(Isotropic/Anisotropic)材质(BRDFs)
- 关键:微表面的方向性
-
各项同性:微表面不存在一定的方向/方向性很弱
-
各向异性:微表面存在明确的方向性
-
8.1.5.1 各向异性BRDFs
反射取决于方位角 ϕ \phi ϕ
f r ( θ i , ϕ i ; θ r , ϕ r ) ≠ f r ( θ i , θ r , ϕ r − ϕ i ) f_r(\theta_i,\phi_i;\theta_r,\phi_r)\neq f_r(\theta_i,\theta_r,\phi_r-\phi_i) fr(θi,ϕi;θr,ϕr)=fr(θi,θr,ϕr−ϕi)
微表面定向的结果,例如,刷亮的金属(brushed metal)
8.1.5.2 各向异性材质:刷亮的金属
8.1.5.3 各向异性材质:尼龙
8.1.5.4 各向异性材质:天鹅绒
8.1.6 材质(BRDF)的性质
- 非负
f r ( ω i → ω r ) ≥ 0 f_r(\omega_i\to\omega_r)\geq0 fr(ωi→ωr)≥0 - 线性
L r ( p , ω r ) = ∫ H 2 f r ( p , ω i → ω r ) L i ( p , ω i ) cos θ i d ω i L_r(\mathrm{p},\omega_r)=\int_{H^2}f_r(\mathrm{p},\omega_i\to\omega_r) L_i(\mathrm{p},\omega_i) \cos\theta_i \mathrm{d}\omega_i Lr(p,ωr)=∫H2fr(p,ωi→ωr)Li(p,ωi)cosθidωi
- 可逆
f r ( ω r → ω i ) = f r ( ω i → ω r ) f_r(\omega_r\to\omega_i)=f_r(\omega_i\to\omega_r) fr(ωr→ωi)=fr(ωi→ωr)
- 能量守恒
∀ ω r ∫ H 2 f r ( ω i → ω r ) cos θ i d ω i ≤ 1 \forall\omega_r\int_{H^2}f_r(\omega_i\to\omega_r) \cos\theta_i \mathrm{d}\omega_i\leq1 ∀ωr∫H2fr(ωi→ωr)cosθidωi≤1 - 各向同性 vs. 各向异性
- 如果是各项同性,则 f r ( θ i , ϕ i ; θ r , ϕ r ) = f r ( θ i , θ r , ϕ r − ϕ i ) f_r(\theta_i,\phi_i;\theta_r,\phi_r)=f_r(\theta_i,\theta_r,\phi_r-\phi_i) fr(θi,ϕi;θr,ϕr)=fr(θi,θr,ϕr−ϕi)
- 接着根据可逆性,
f
r
(
θ
i
,
θ
r
,
ϕ
r
−
ϕ
i
)
=
f
r
(
θ
r
,
θ
i
,
ϕ
i
−
ϕ
r
)
=
f
r
(
θ
i
,
θ
r
,
∣
ϕ
r
−
ϕ
i
∣
)
f_r(\theta_i,\theta_r,\phi_r-\phi_i)=f_r(\theta_r,\theta_i,\phi_i-\phi_r)=f_r(\theta_i,\theta_r,|\phi_r-\phi_i|)
fr(θi,θr,ϕr−ϕi)=fr(θr,θi,ϕi−ϕr)=fr(θi,θr,∣ϕr−ϕi∣)
8.2 测量BRDFs
8.2.1 动机
- 避免需要开发/派生模型
- 自动包含所有呈现的散射效果
- 可以精确渲染真实世界的材质
- 对产品设计、特效很有用
8.2.2 基于图像的BRDF测量
8.2.3 使用gonioreflectometer测量BRDFs
通用的方法:
提高效率:
- 各向同性表面将维度从四维降低到三维
- 可逆性使测量次数减少了一半
- 巧妙的光学系统……
8.2.4 BRDFs测量的挑战
- 在掠射角的精确测量,对于菲涅尔项很重要
- 用足够密集的采样来测量,以捕获高频镜面反射
- 逆反射
- 空间变化的反射率……
8.2.5 测量后BRDFs的表示
存储要求:
- 紧凑的表现形式
- 测量数据的精确表示
- 任意出射入射方向的有效评估
- 适用于重要性抽样的良好分布
8.2.6 表格表示法
将有规则间隔的样本存储在 ( θ i , θ o , ∣ ϕ i − ϕ o ∣ ) (\theta_i,\theta_o,|\phi_i-\phi_o|) (θi,θo,∣ϕi−ϕo∣)
重新参数化角度以更好地匹配镜面反射
通常需要将测量值重新采样到表格中
非常高的存储要求
8.3 高级的光线传播
-
无偏的光线传播方式 Unbiased light transport methods
- 双向路径追踪 Bidirectional path tracing (BDPT)
- Metropolis light transport (MLT)
-
有偏的光线传播方式 Biased light transport methods
- Photon mapping
- Vertex connection and merging (VCM)
-
实时辐射度算法(Instant radiosity) (VPL/ many light methods)
8.3.1 有偏的 vs. 无偏的蒙特卡洛估计
-
无偏的蒙特卡罗估计没有任何系统误差
- 无偏估计量的期望值将永远是正确的值,无论使用多少样本
-
否则,都是有偏的
- 在一种特殊情况下,由于使用了无限个样本,期望值收敛到正确值——一致的
8.3.2 双向路径追踪 Bidirectional path tracing (BDPT)
- 回顾:一条连接相机和光线的路径
- BDPT
- 跟踪来自摄像机和光源的子路径
- 连接两个子路径的端点
- 如果光传输在光一侧很复杂,则适合
- 难以实现且相当缓慢
8.3.3 Metropolis光线传输 Metropolis light transport (MLT)
- 马尔可夫链的蒙特卡罗方法(MCMC)应用
- 通过PDF根据当前样本生成下一个样本
- 非常擅长在局部难以探索的光路
- 关键点
- 对现有路径进行局部扰动以获得新路径
- 对现有路径进行局部扰动以获得新路径
- 优势
- 适用于复杂的光路传播
- 也是无偏的
- 劣势
- 收敛速度难以估计
- 不能保证每个像素的收敛速度相等
- 所以,通常会产生“脏”的结果
- 因此,通常不用于渲染动画
8.3.4 光子映射 Photon Mapping
-
有偏的方法 & 两阶段方法
-
非常擅长处理镜面-漫反射-镜面(SDS)路径和生成caustics
其中一种实现方法: -
第一阶段——光子追踪
- 从光源发射光子,使其弹来弹去,接着在漫射表面记录光子
- 从光源发射光子,使其弹来弹去,接着在漫射表面记录光子
-
第二阶段——光子收集(最终收集)
- 从相机发射子路径,使它们来回跳动,直到它们碰到漫反射表面
-
计算——局部密度估计(local density estimation)
- 想法:有更多光子的区域应该更亮
- 对于每个着色点,找其最近的
N
N
N个光子,接着找到它们覆盖的面积
-
为什么光子映射是一个有偏的方法?
-
局部密度估计(local density estimation)
d N / d A ! = Δ N / Δ A dN/dA\ ! =\ \Delta N/\Delta A dN/dA != ΔN/ΔA -
但在有限的意义上
- 发射更多的光子→
- 相同数量的 N N N 个光子覆盖一个较小的 Δ A \Delta A ΔA→
- Δ A \Delta A ΔA更接近 d A dA dA
-
因此,有偏但一致!
-
在渲染中更容易理解的有偏
- 有偏 == 模糊
- 一致 == 样本足够多就能达到不模糊的效果
-
为什么不做一个“固定范围”的密度估计?
- 因为随着光子发射的越多,同一个范围内的光子会越来越多, Δ A \Delta A ΔA永远不会更接近 d A dA dA
8.3.5 Vertex connection and merging (VCM)
- BDPT与光子映射的结合
- 关键思想
- 如果BDPT中的子路径不能连接,但可以合并,那么就不要浪费它们
- 利用光子映射处理附近"光子"的合并
8.3.6 实时辐射度算法 Instant Radiosity(IR)
- 有时也被称为多光源方法(many-light approaches)
- 关键思想
- 已被照亮的表面可被当作是光源,再来照亮其它物体
- 方法
- 从光源发射许多光的子路径,并假设每个子路径的端点(停在的地方)是虚拟点光源(Virtual Point Light,VPL)
- 像往常一样使用这些 VPL 渲染场景
- 优势:速度快,在漫反射场景中通常能得到很好的效果
- 劣势
- 当 VPL 接近着色点时会出现奇怪的亮点
- 不能处理有光泽的材质(glossy materials)
8.4 高级的外观建模
- 非表面模型
- 散射介质
- 毛发/毛皮/纤维(BCSDF)
- 粒状材质
- 表面模型
- 半透明材质(BSSRDF)
- 布料
- 有细节的材质(non-statistical BRDF)
- 程序化生成的模型
8.4.1 非表面模型 Non-Surface Models
8.4.1.1 散射介质 Participating Media
雾 Fog
云 Cloud
- 在光通过散射介质的任何一点上,它都可能(部分)被吸收和散射。
- 使用相位函数(Phase Function)来描述散射介质中的任意点x的光散射角分布。
渲染: - 随机选择一个弹射的方向
- 随机选择一段直线距离
- 在每个着色点上,连接到光源
应用:
8.4.1.2 毛发表现
8.4.1.2.1 Kajiya-Kay模型
8.4.1.2.2 Marschner模型
- 将头发看作是类玻璃状圆柱体
- 3种类型的光相互作用:R, TT, TRT(R:反射,T:透射)
8.4.1.3 动物毛皮表现
8.4.1.3.1 视作人的头发
不能代表扩散和饱和的外观
8.4.1.3.2 人类头发 vs. 动物毛发
8.4.1.3.3 髓质的重要性
随着髓质大小的增加:
8.4.1.3.3 双层圆柱模型 Double Cylinder Model
Marschner模型:
Double Cylinder Model:
8.4.1.4 粒状材质 Granular Material
- 什么是粒状材质?
- 我们能避免所有颗粒的显式建模吗?
- 是的,有程序定义
- 是的,有程序定义
8.4.2 表面模型 Surface Models
8.4.2.1 半透明材质 Translucent Material
玉石
水母
8.4.2.1.1 次表面散射 Subsurface Scattering
由于光在不同的点出射而引起的许多表面的视觉特性
- 违反了BRDF的一个基本假设
- BSSRDF:BRDF的延伸,由于在另一点上的入射差分辐亮度而在一点上产生的附加辐亮度:
S ( x i , ω i , x o , ω o ) S(x_i,\omega_i,x_o,\omega_o) S(xi,ωi,xo,ωo) - 渲染方程的延伸:对曲面上的所有点和所有方向进行积分(!)
L ( x o , ω o ) = ∫ A ∫ H 2 S ( x i , ω i , x o , ω o ) L i ( x i , ω i ) cos θ i d ω i d A L(x_o,\omega_o)=\int_A\int_{H^2}S(x_i,\omega_i,x_o,\omega_o) L_i(x_i,\omega_i)\cos\theta_i \mathrm{d}\omega_i \mathrm{d}A L(xo,ωo)=∫A∫H2S(xi,ωi,xo,ωo)Li(xi,ωi)cosθidωidA
8.4.2.1.2 两极近似 Dipole Approximation
- 通过引入两个点光源来近似光的漫反射
8.4.2.1.3 BRDF vs. BSSRDF
8.4.2.1.4 BSSRDF的应用
8.4.2.2 布料 Cloth
-
一系列缠绕的纤维!
-
两级的缠绕
-
机织或针织
8.4.2.2.1 布料:渲染成表面
- 给定编织模式,计算整体行为
- 使用BRDF渲染
8.4.2.2.2 布料:渲染成表面——局限性
8.4.2.2.3 布料:渲染成散射介质
- 单根纤维的特性及其分布 → 散射参数
- 呈现为散射介质
8.4.2.2.4 布料:渲染成实际纤维
- 显式渲染每一根纤维!
8.4.2.3 有细节的材质 Detaiied Appearance
- 看起来并不真实,为什么?
- 过于完美
- 过于完美
- 真实世界更加复杂多变
- 细节感带来真实感
8.4.2.3.1 综述:微表面材质 Microfacet BRDF
S
u
r
f
a
c
e
=
S
p
e
c
u
l
a
r
m
i
c
r
o
f
a
c
e
t
s
+
s
t
a
t
i
s
t
i
c
a
l
n
o
r
m
a
l
s
Surface=Specular\ microfacets + statistical\ normals
Surface=Specular microfacets+statistical normals
8.4.2.3.2 统计NDF vs. 实际NDF
8.4.2.3.3 细节定义
8.4.2.3.3 不同的细节
8.4.2.3.4 渲染?太难了!
8.4.2.3.5 路径采样困难问题
8.4.2.3.6 解决方案:像素上的BRDF
8.4.2.3.7 p-NDFs具有明显的特征
8.4.2.3.8 p-NDFs的形状
8.4.2.3.9 应用
8.4.2.3.10 近期趋势:波动光学
8.4.2.3.11 波动光学下的细节材质
8.4.2.4 程序化生成的模型 Procedural Appearance
- 我们可以定义没有纹理的细节吗?
- 是的!动态计算噪声函数
- 3D噪音→内部结构(如果被切割或损坏)
- 阈值化(噪声→二进制噪声)
- 是的!动态计算噪声函数
- 复杂的噪声函数可以非常强大
9 相机、镜头和光场
成像=合成+捕获
9.1 相机
9.1.1 相机里发生了什么?
9.1.1.1 针孔和透镜在传感器上形成图像
9.1.1.2 快门曝光传感器,实现精确的耐用度
9.1.1.3 传感器在曝光期间累积辐照度
9.1.1.4 为什么不是没有透镜的传感器?
每个传感器点将集成来自物体上所有点的光,因此所有像素值将是相似的,即传感器记录辐照度。
但有计算成像的研究正在进行…
9.2 针孔相机成像原理
9.2.1 针孔相机
9.2.2 最大针孔照片
9.3 视场 Field of View (FOV)
9.3.1 焦距(Focal Length)对视场(FOV)的影响
对于一个固定的传感器尺寸,减少焦距而增大视场。
F O V = 2 arctan ( h 2 f ) \mathrm{FOV}=2\arctan\left(\frac{h}{2f}\right) FOV=2arctan(2fh)
9.3.2 焦距(Focal Length) vs. 视场(FOV)
- 由于历史原因,通常用 35mm格式的胶片(36x24mm)上使用的镜头的焦距来表示视场角。
- 35mm格式的焦距示例:
- 17mm是广角镜头104°
- 50mm是"正常"镜头47°
- 200mm是长焦镜头12°
- 注意!当我们说当前的手机大约有28毫米“等效”焦距时,这是使用了上述惯例
9.3.3 传感器尺寸对视场的影响
9.3.4 传感器尺寸
9.3.5 在较小的传感器上保持视场
为保持视场,应按传感器的宽度/高度比例减小镜头焦距
9.4 曝光 Exposure
9.4.1 曝光是什么
- H = T × E H =T\times E H=T×E
- 曝光 = 时间 × 辐照度 曝光=时间\times 辐照度 曝光=时间×辐照度
- 曝光时间(T)
- 由快门控制
- 辐照度(E)
- 落在传感器单位面积上的光能量
- 由镜头光圈和焦距控制
9.4.2 摄影中的曝光控制
- 光圈大小
- 通过打开/关闭光圈来更改光圈(如果相机有光圈控制)
- 快门速度
- 更改传感器像素集成光的持续时间
- ISO增益(感光度)
- 更改传感器值和数字图像信号之间的放大(模拟 and/or 数字)
9.4.3 曝光:光圈,快门,增益(ISO)
9.4.3.1 ISO(增益)
曝光的第三个变量
胶片:以纹理灵敏度换取
数码:以噪声灵敏度换取
- 模数转换前的乘法信号
- 线性效果(ISO 200需要的光是ISO 100的一半)
9.4.3.2 佳能T2i中的ISO增益 vs 噪声
9.4.3.3 F-数(F-Stop):曝光水平
写作 F N FN FN或 F / N F/N F/N。 N N N是 f f f-数。
非正式的理解:光圈圆形孔径的直径分之一
即 N N N越大,直径越小。
9.4.3.4 物理快门(1/25秒曝光)
9.4.3.5 快门速度的副作用
运动模糊:握手、主体移动
双倍快门时间加倍运动模糊
注意:动态模糊不总是不好的!
提示:考虑抗锯齿
滚动快门:在不同时间拍摄的照片的不同部分,会造成扭曲(螺旋桨)
9.4.3.6 恒定曝光:F-stop vs. 快门速度
例如:这些光圈和快门速度的组合提供相等的曝光。
如果曝光太亮/太暗,可能需要调整光圈 and/or 快门向上/向下。
- 摄影师必须在景深和动态模糊之间进行权衡,以适应移动的拍摄对象
9.4.4 高速与低速摄影
9.4.4.1 高速摄影
正常曝光 = 极快快门速度×(大光圈 and/or 高ISO)
9.4.4.2 长曝光摄影
又称“延时摄影”,拉丝。
9.5 薄透镜近似
9.5.1 真实镜头设计高度复杂
9.5.2 真正的透镜元件并不理想-像差
真正的平凸透镜(球面形状)。透镜不会把光线汇聚到任何地方。
9.5.3 理想薄透镜-焦点
- 所有进入透镜的平行光线都要经过它的焦点。
- 通过焦点的所有射线在经过透镜后将平行。
- 焦距可以任意改变(在现实中,是的!)。
9.5.4 薄透镜方程
1
f
=
1
z
i
+
1
z
o
\frac1f=\frac1{z_i}+\frac1{z_o}
f1=zi1+zo1
- z o z_o zo:物距
- z i z_i zi:像距
- f f f:焦距
从任一方向穿过透镜中心的光线,其方向都不会发生改变。
9.5.5 高斯射线图
9.5.6 高斯射线追踪结构
共轭深度 z o z_o zo、 z i z_i zi之间的关系是什么?
9.6 景深模糊 Defocus Blur
9.6.1 计算混沌环(CoC)尺寸
C A = d ′ z i = ∣ z s − z i ∣ z i \frac CA=\frac{d^{\prime}}{z_i}=\frac{|z_s-z_i|}{z_i} AC=zid′=zi∣zs−zi∣
CoC与透镜光圈大小成正比
9.6.2 CoC vs. 透镜光圈大小
左侧:大光圈,模糊效果
右侧:小光圈,清晰效果
9.6.3 再看F-数(又称F-Stop)
形式定义:镜头f-数的定义为焦距除以光圈的直径。
真实镜头上常见的f-数:1.4,2,2.8,4.0,5.6,8,11,16,22,32。
2的f-stop有时写成f/2,反映了绝对孔径(A)可以通过将焦距(f)除以相对孔径(N)来计算。
9.6.4 F-Stop计算示例
9.6.5 CoC的大小与F-Stop成反比
C
=
A
∣
z
s
−
z
i
∣
z
i
=
f
N
∣
z
s
−
z
i
∣
z
i
C=A\frac{|z_s-z_i|}{z_i}=\frac fN\frac{|z_s-z_i|}{z_i}
C=Azi∣zs−zi∣=Nfzi∣zs−zi∣
9.7 薄透镜光线追踪
9.7.1 用镜头聚焦的静物片示例
9.7.2 景深模糊的光线跟踪(薄透镜)
(一个可能的)步骤:
- 选择传感器尺寸、镜头焦距和光圈尺寸
- 选择感兴趣物体的深度 z o z_o zo
- 根据薄透镜方程(聚焦)计算相应的传感器深度 z i z_i zi
渲染:
- 对于传感器上的每个像素 x ′ x^\prime x′(实际上是胶片)
- 取透镜平面上的随机点 x ′ ′ x^{\prime\prime} x′′
- 你知道穿过透镜的光线会照射到 x ′ ′ ′ x^{\prime\prime\prime} x′′′(使用薄透镜公式)
- 估计从 x ′ ′ x^{\prime\prime} x′′到 x ′ ′ ′ x^{\prime\prime\prime} x′′′的辐亮度(radiance)
9.8 景深 Depth of Field
将CoC设置为图像平面上允许的最大模糊点,该点在最终观看条件下会显得很尖锐
9.8.1 Depth of Field的CoC
例如,在场景中,相应的CoC被认为足够小的深度范围
9.8.2 Depth of Field
d
N
−
d
S
d
N
=
C
A
d
S
−
d
F
d
F
=
C
A
N
=
f
A
1
D
F
+
1
d
F
=
1
f
1
D
S
+
1
d
S
=
1
f
1
D
N
+
1
d
N
=
1
f
\begin{aligned}\frac{d_N-d_S}{d_N}&=\frac CA\\\frac{d_S-d_F}{d_F}&=\frac CA\\N&=\frac fA\\\frac1{D_{F}}+\frac1{d_{F}} &=\frac{1}{f} \\\frac1{D_{S}}+\frac1{d_{S}} &=\frac{1}{f}\\\frac1{D_N}+\frac1{d_N} &=\frac{1}{f} \end{aligned}
dNdN−dSdFdS−dFNDF1+dF1DS1+dS1DN1+dN1=AC=AC=Af=f1=f1=f1
D
O
F
=
D
F
−
D
N
\mathrm{DOF}=D_F-D_N
DOF=DF−DN
D
F
=
D
S
f
2
f
2
−
N
C
(
D
S
−
f
)
D
N
=
D
S
f
2
f
2
+
N
C
(
D
S
−
f
)
D_F=\frac{D_Sf^2}{f^2-NC(D_S-f)}\quad D_N=\frac{D_Sf^2}{f^2+NC(D_S-f)}
DF=f2−NC(DS−f)DSf2DN=f2+NC(DS−f)DSf2
9.9 光场/Lumigraph
9.9.1 我们看到了什么
9.9.2 全光函数 The Plenoptic Function
Q:我们所能看到的所有事物的集合是什么?
A:全光函数(Adelson & Bergen)
让我们从一个静止的人开始,试着把他能看到的一切都叙述出来……
9.9.3 灰度快照
P ( θ , ϕ ) P(\theta,\phi) P(θ,ϕ)
- 是光的强度
- 从单一角度看
- 在同一时间
- 在可见光谱波长上的平均值
也作 P ( x , y ) P(x,y) P(x,y),但球坐标更好
9.9.4 彩色快照
P
(
θ
,
ϕ
,
λ
)
{P}(\theta,\phi,\lambda)
P(θ,ϕ,λ)
- 是光的强度
- 从单一角度看
- 在同一时间
- 作为波长的函数
9.9.5 一部电影
P
(
θ
,
ϕ
,
λ
,
t
)
P(\theta,\phi,\lambda,t)
P(θ,ϕ,λ,t)
- 是光的强度
- 从单一角度看
- 随着时间的推移
- 作为波长的函数
9.9.6 全息电影
P
(
θ
,
ϕ
,
λ
,
t
,
V
x
,
V
y
,
V
z
)
P(\theta,\phi,\lambda,t,V_x,V_y,V_z)
P(θ,ϕ,λ,t,Vx,Vy,Vz)
- 是光的强度
- 从任何角度看
- 随着时间的推移
- 作为波长的函数
9.9.7 全光函数
P
(
θ
,
ϕ
,
λ
,
t
,
V
x
,
V
y
,
V
z
)
P(\theta,\phi,\lambda,t,V_x,V_y,V_z)
P(θ,ϕ,λ,t,Vx,Vy,Vz)
- 可以重建每一个可能的视图,在每一刻,从每一个位置,在每一个波长
- 包含了每张照片,每部电影,任何人看过的一切!它完全捕捉到了我们的视觉现实!对一个函数来说还不错……
9.9.8 采样全光函数(俯视图)
只需查找即可—— Quicktime VR
9.9.9 光线
不要担心时间和颜色:
9.9.10 光线重定义
无限线
假设光是恒定的(真空)
9.9.11 只需要全光表面
9.9.12 综合新颖的视角
光场是一个四维的函数,能向我们提供任意观测方向所能看到的结果。
9.9.13 光场
外凸空间
9.9.14 光场——组织
9.9.15 斯坦福摄像机矩阵
9.9.16 全景图像(“苍蝇眼”镜头)
使用透镜的空间多路复用光场捕获:
- 在空间分辨率和角度分辨率之间进行固定的权衡
9.10 光场照相机
9.10.1 Lytro光场照相机
Lytro:由伦(Ren Ng)教授(加州大学伯克利分校)创立的微透镜设计
最重要的功能
- 计算重聚焦(实际改变焦距和光圈大小等,在拍完照片后)
9.10.2 光场照相机
理解
- 每个像素(辐照度irradiance)现在存储为一个像素块(辐亮度radiance)
- 拍摄照片的特写镜头
如何从光场照片中获得一张“普通”照片? - 一个简单的情况——始终选择每个块底部的像素
- 然后是中央和顶部的
- 本质上就是“移动摄像头的位置’
计算/数字重定向
- 同样的想法:在视觉上改变焦距,相应地选择重新聚焦的光线方向
总而言之,所有这些功能都是可用的,因为 - 光场包含了一切
光场摄像机有问题吗?
- 空间分辨率不足(空间和方向信息均使用相同的胶片)
- 高成本(微透镜的复杂设计)
- 计算机图形学讲究权衡取舍
10 颜色感知
10.1 色彩的物理基础
10.1.1 光的基本组成
- 牛顿证明用棱镜可以把阳光细分成彩虹
- 产生的光不能用第二个棱镜进一步细分
10.1.2 光的可见光谱
电磁辐射
- 不同频率(波长)的振荡
10.1.3 谱功率密度 Spectral Power Distribution (SPD)
测量光的显著性
- 每种波长的光量
- 单位:
- 辐射计量单位/纳米(如瓦特/纳米)
- 也可以是无单位的
- 当绝对单位并不重要时,经常使用缩放到最大波长的“相对单位”来进行不同波长的比较
10.1.4 日光的谱功率密度各不相同
10.1.5 光源的谱功率密度
描述能量按波长分布
10.1.6 谱功率密度的线性性质
10.1.7 什么是颜色?
- 颜色是人类感知的一种现象,它不是光的普遍属性
- 不同波长的光不是“颜色
10.2 颜色的生物学基础
10.2.1 人眼的解剖学
10.2.2 视网膜感光细胞:视杆细胞和视锥细胞
视杆细胞是在非常暗的光线(“暗视”条件下)的主要受体,例如昏暗的月光。
- 大约有1.2亿视杆细胞存在眼中
- 只感知到灰色的阴影,没有颜色
视锥细胞是典型光照水平下的主要受体(“明视”)
- 大约有600万至700万个视锥细胞
- 三种锥体,每种锥体的光谱灵敏度不同
- 提供色彩感觉
10.2.3 人类视锥细胞的光谱响应
三种类型的视锥细胞:S、M和L(对应于短、中、长波长的峰值响应)
10.2.4 三种视锥细胞类型的比例差异很大
12例正常色觉者黄斑中心凹边缘视锥细胞的分布。注意不同视锥细胞类型的百分比变化很大。(伪彩色图像)
10.3 色彩的三刺激理论
10.3.1 人类视锥细胞的光谱响应
现在我们有三个探测器(S、M、L视锥细胞),每一个都有不同的光谱响应曲线。
S
=
∫
r
S
(
λ
)
s
(
λ
)
d
λ
M
=
∫
r
M
(
λ
)
s
(
λ
)
d
λ
L
=
∫
r
L
(
λ
)
s
(
λ
)
d
λ
\begin{gathered} \text{S}=\int r_{S}(\lambda)s(\lambda) d\lambda \\ \text{M} =\int r_M(\lambda)s(\lambda) d\lambda \\ \text{L} =\int r_{L}(\lambda)s(\lambda) d\lambda \end{gathered}
S=∫rS(λ)s(λ)dλM=∫rM(λ)s(λ)dλL=∫rL(λ)s(λ)dλ
10.3.2 人类的视觉系统
- 人的眼睛不能测量,大脑不能接收每种波长的光的信息
- 相反,眼睛"看到"只有三个响应值(S,M,L),这是唯一提供给大脑的信息
10.4 同色异谱
10.4.1 不同光谱能量分布的同色光 Metamers
Metamers是投射到相同(S,M,L)(3-dim)响应的两种不同的光谱(∞-dim)
- 对人类来说它们的颜色是一样的
Metamers的存在对色彩再现至关重要
- 不需要重现真实世界场景的全貌
- 例如:Metamers可以在只有三种颜色像素的显示器上再现现实世界场景的感知颜色
10.4.2 同色异谱是一种效应
10.4.3 同色异谱
色彩搭配背后的理论
10.5 色彩再现/匹配
10.5.1 加色
给定一组主光,每个主光具有自己的光谱分布(例如R、G、B显示像素):
s
R
(
λ
)
,
s
G
(
λ
)
,
s
B
(
λ
)
s_R(\lambda), s_G(\lambda), s_B(\lambda)
sR(λ),sG(λ),sB(λ)
调整这些灯光的亮度,并将它们添加到一起:
R
s
R
(
λ
)
+
G
s
G
(
λ
)
+
B
s
B
(
λ
)
R s_R(\lambda)+G s_G(\lambda)+B s_B(\lambda)
RsR(λ)+GsG(λ)+BsB(λ)
颜色现在由标量值描述:
R
,
G
,
B
R,G,B
R,G,B
10.5.2 加色配色实验
10.5.3 示例实验
10.5.4 实验2
10.5.5 CIE RGB颜色匹配实验
与之前的加色匹配相同的设置,但三原色是单色光(单波长)
测试光也是单色光
10.5.6 CIE RGB颜色匹配函数
图标示了每种CIE RGB主光的组合量,以匹配x轴上给出的波长的单色光。
10.5.7 使用匹配函数再现色彩
对于任何光谱,感知到的颜色通过以下公式匹配,用于缩放CIE RGB三原色
R
C
I
E
R
G
B
=
∫
λ
s
(
λ
)
r
ˉ
(
λ
)
d
λ
G
C
I
E
R
G
B
=
∫
λ
s
(
λ
)
g
ˉ
(
λ
)
d
λ
B
C
I
E
R
G
B
=
∫
λ
s
(
λ
)
b
ˉ
(
λ
)
d
λ
\begin{gathered} R_{\mathrm{CIE\ RGB}} =\int_{\lambda}s(\lambda)\bar{r}(\lambda) d\lambda \\ G_{\mathrm{CIE\ RGB}} =\int_{\lambda}s(\lambda)\bar{g}(\lambda) d\lambda \\ B_{\mathrm{CIE\ RGB}} =\int_{\lambda}s(\lambda)\bar{b}(\lambda) d\lambda \end{gathered}
RCIE RGB=∫λs(λ)rˉ(λ)dλGCIE RGB=∫λs(λ)gˉ(λ)dλBCIE RGB=∫λs(λ)bˉ(λ)dλ
10.5.8 标准色彩空间
标准化RGB(sRGB)
- 使特定的显示器RGB标准化
- 其他彩色设备通过校准模拟显示器
- 至今广泛采用的
- 色域(?)是有限的
10.5.9 一个通用的色彩空间:CIE XYZ
标准色原色X、Y、Z的虚设集
- 不存在具有这些匹配函数的原色
- Y是亮度(不分颜色的亮度)
这样设计
- 匹配函数是严格正的
- 覆盖所有可观察到的颜色
10.5.10 分离亮度、色度
亮度:Y
色度:x,y,z,定义为
x = X X + Y + Z y = Y X + Y + Z z = Z X + Y + Z \begin{gathered} x=\frac{X}{X+Y+Z} \\ y=\frac{Y}{X+Y+Z} \\ z=\frac{Z}{X+Y+Z} \end{gathered} x=X+Y+ZXy=X+Y+ZYz=X+Y+ZZ
- 由于x+y+z=1,我们只需要记录三个中的两个
- 通常选择x和y,得出特定亮度Y的(x,y)坐标
10.5.11 CIE彩度图
弯曲边界
- 称为谱轨迹
- 对应于单色光(每一点代表单一波长的纯色)
其中的任何颜色都不那么纯净
- 即混合的
10.5.12 色域
色域是由一组基色产生的一组色度。
不同的色彩空间代表不同的颜色范围。
所以它们有不同的色域,即,它们覆盖了色度图上的不同区域。
10.6 感知组织的颜色空间
10.6.1 HSV色彩空间(Hue-Saturation-Value)
轴线对应色彩的艺术特征
广泛应用于“取色器”中
10.6.2 色彩的感知维度
- Hue(色调)
- 颜色的“种类”,不考虑属性
- 比色相关:主导波长
- 艺术家的相关性:所选颜料颜色
- Saturation(饱和度)
- “色彩斑斓’
- 比色相关:纯度
- 艺术家的相关性:彩色管中的颜料分数
- Lightness (or value)(亮度)
- 总光量
- 色度相关:亮度
- 艺术家的相关性:色调较浅,色调较深
10.6.3 CIE LAB空间(又名L*a*b*)
一个常用的色彩空间,力求在感知上的一致性
- L*是亮度
- a*和b*是色对
- a*为红绿色
- b*为蓝黄色
10.6.4 互补色理论
CIE LAB的颜色空间维度有很好的神经学基础
- 大脑似乎在早期使用三个轴来编码颜色:
- 白色-黑色,红色-绿色,黄色-蓝色
- 白色-黑色轴是亮度,其他轴决定色调和饱和度
10.6.5 一切都是相对的
10.6.6 CMYK:一种减色的色彩空间
减色模型
- 混合得越多,颜色就越深
青色、洋红色、黄色和黑色广泛应用于印刷
问题:
如果混合C、M和Y得到K,为什么需要K?
11 动画与模拟
11.1 动画
“将事物带入生活”
- 交流的工具
- 审美问题往往主导技术问题
建模的一种扩展
- 将场景模型表示为时间的函数
输出:连续观看时提供动作感觉的图像序列
- 电影:24 fps(帧每秒)
- 视频(一般):30 fps
- 虚拟现实:90 fps
11.2 动画中的历史节点
11.2.1 最早的动画
11.2.2 动画的历史
11.2.3 最早的电影
最初用作科学工具而非娱乐
加速动画发展的关键技术
11.2.4 第一部手绘完整长度(>40分钟)动画
11.2.5 第一部数字计算机生成的动画
11.2.6 早期的计算机动画
11.2.7 数字恐龙!
11.2.8 第一部CG长片
11.2.9 计算机动画——十年前
11.2.9 计算机动画——五年前
11.3 关键帧动画
11.3.1 关键帧动画
动画师(如首席动画师)创建关键帧
助手(人或计算机)在帧之间创建帧增强
11.3.2 关键帧插值
将每个帧看作参数值的向量
11.3.3 各参数的关键帧插值
线性插值通常不够好
用于平滑/可控插值的召回曲线
11.4 物理模拟
11.4.1 牛顿定律
11.4.2 基于物理的动画
使用数值模拟生成物体的运动
11.4.3 例子:布料模拟
11.4.4 例子:流体
11.5 质点弹簧系统:动态系统建模示例
11.5.1 例子:质点弹簧绳
11.5.2 例子:毛发
11.5.3 例子:质点弹簧网格
11.5.4 一个简单的弹簧
理想的弹簧
f
a
→
b
=
k
S
(
b
−
a
)
f
b
→
a
=
−
f
a
→
b
\begin{aligned}&f_{a\to b}=k_S(\boldsymbol{b}-\boldsymbol{a})\\&f_{b\to a}=-\boldsymbol{f}_{a\to b}\end{aligned}
fa→b=kS(b−a)fb→a=−fa→b
力将各点拉在一起
强度与位移成正比(胡克定律)
k s k_s ks是一个弹性系数:硬度
问题:这个弹簧希望长度为零
11.5.5 非零长度弹簧
非零长度(rest length)的弹簧
11.5.6 导数的点表示法
如果 x x x 是感兴趣点位置的矢量,我们将使用点符号表示速度和加速度:
x x ˙ = v x ¨ = a \begin{aligned}&\boldsymbol{x}\\&\boldsymbol{\dot{x}}=\boldsymbol{v}\\&\boldsymbol{\ddot{x}}=\boldsymbol{a}\end{aligned} xx˙=vx¨=a
11.5.7 引入能量损失
简单运动阻尼
f
=
−
k
d
b
˙
f=-k_d\dot{\boldsymbol{b}}
f=−kdb˙
- 在运动中表现得像粘性拖曳物
- 减慢速度方向上的运动速度
- k d k_d kd是一个阻尼系数
问题:减慢所有运动
- 想要一个生锈的弹簧的振动减慢,但它是否也应该更慢地落到地上呢?
11.5.8 弹簧的内部阻尼
仅阻尼内部弹簧驱动的运动
- 粘滞阻力只影响弹簧长度的变化
- 不会减缓弹簧系统的群组运动(例如群组的全局平移或旋转)
- 注:这只是一种特定类型的阻尼
11.5.9 弹簧结构
行为是由结构联系决定的
这种结构不能抵抗切变
这种结构无法抵抗平面外弯曲
这种结构能抵抗切变,但具有各向异性的偏差
这种结构同样不能抵抗面外弯曲
这种结构能抗切变,方向性偏差小
这种结构同样不能抵抗面外弯曲
这种结构能抵抗切变,方向偏差较小
这种结构能抵御平面外弯曲
红弹簧应该弱得多
11.5.10 其它模拟方式:FEM(有限元法)代替弹簧
11.6 粒子系统 Particle Systems
11.6.1 粒子系统 Particle Systems
将动力系统模型化为大量粒子的集合体
每个粒子的运动由一组物理力(或非物理力)来定义
图形和游戏中的流行技术
- 易于理解、实现
- 可扩展性:粒子越少,速度越快,粒子越多,复杂性越高
挑战
- 可能需要许多颗粒(如流体)
- 可能需要加速度结构(例如为相互作用寻找最近的粒子)
11.6.2 粒子系统动画
对于动画中的每个帧
- 必要时创造新的粒子
- 计算每个粒子的受力
- 更新每个粒子的位置和速度
- 必要时去除死粒子
- 渲染粒子
11.6.3 粒子系统作用力
- 吸引力和排斥力
- 重力、电磁力……
- 弹簧、推进力……
- 阻尼力
- 摩擦力、空气阻力、粘滞力……
- 碰撞
- 墙壁、容器、固定的物体……
- 动态对象、角色身体部位……
11.6.4 引力
牛顿万有引力定律
- 粒子间的引力
F g = G m 1 m 2 d 2 G = 6.67428 × 1 0 − 11 N m 2 k g − 2 \begin{aligned}&F_{g}=G\frac{m_{1}m_{2}}{d^{2}}\\&G=6.67428\times10^{-11} \mathrm{Nm}^{2}\mathrm{kg}^{-2}\end{aligned} Fg=Gd2m1m2G=6.67428×10−11Nm2kg−2
11.6.5 例子:银河模拟
11.6.6 例子:基于粒子的流体
11.6.7 作为ODE的模拟群集
把每只鸟模型化为粒子
受制于非常简单的力:
- 对邻居中心的吸引力
- 对个别邻居的排斥
- 朝邻近区域平均轨迹方向的排列
用数值模拟大粒子系统的演化突发复杂行为(也见于鱼类、蜜蜂、…)
11.6.8 例子:分子动力学
11.6.9 例子:人群+“岩石”动力学
11.7 正向运动学 Forward Kinematics
11.7.1 正向运动学
- 关节骨架
- 拓扑结构(什么连接到什么)
- 节点的几何关系
- 树结构(在没有循环的情况下)
- 节类型 Joint types
-
销 Pin (一维旋转)
-
球 Ball (二维旋转)
-
棱柱接头 Prismatic joint (移动)
-
例子:二维简易双节臂
动画机提供角度,计算机确定最终效果器的位置p
11.7.2 例子:步行循环
11.7.3 运动学的利弊
- 优势
- 直接控制方便
- 实现起来很简单
- 劣势
- 动画可能与物理不一致
- 对艺术家来说很耗时
11.8 逆运动学 Inverse Kinematics
11.8.1 逆运动学
直接逆运动学:对于两段臂,可分析求解参数
为什么这个问题很难?
- 在空间中存在多种解决方案
- 可能不存在所谓的解决方案
一般N链IK问题的数值解
- 选择初始配置
- 定义误差度量(例如,目标与当前位置之间距离的平方)
- 计算误差梯度作为配置的函数
- 应用梯度下降法(或牛顿法或其他优化程序)
11.8.2 基于风格的IK
11.9 绑定 Rigging
11.9.1 绑定
绑定是角色上一套更高级别的控制,允许更快速直观地修改姿势、变形、表情等。
重点
- 就像木偶上的线
- 捕获所有有意义的人物变化
- 因人物而异
创造成本高
- 手工劳动
- 需要艺术和技术培训
11.9.2 绑定的例子
11.9.3 融合形状 Blend Shapes
直接在曲面之间插值,而不是骨架
例如,建模一组面部表情:
最简单的方案:采取顶点位置的线性组合
样条用于随着时间的推移控制权重选择
11.10 动作捕捉 Motion Capture
11.10.1 动作捕捉
数据驱动的动画序列创建方法
- 记录真实世界的表现(例如,执行某项活动的人)
- 从收集到的数据中提取姿势作为时间的函数
11.10.2 动作捕捉的利弊
- 优势
- 可以快速捕获大量真实数据
- 真实性很高
- 弱点
- 复杂而昂贵的设置
- 捕捉到的动画可能不符合艺术需要,需要修改
11.10.3 动作捕捉的装备
11.10.4 光学运动捕捉
- 物体标记
- 通过多个摄像头的三角测量定位
- 8+摄像机,240 Hz,遮挡时很困难
11.10.5 捕捉到的数据
11.10.6 面部动画的挑战
恐怖谷效应
- 机器人和图形学
- 随着人造人物外表接近人类的现实主义,我们的情绪反应变得消极,直到它在表达上达到足够令人信服的现实主义水平。
11.10.7 面捕动作捕捉
11.10.8 生产流水线
11.11 简单粒子模拟
首次研究单个粒子的运动
- 接着,概括为多种粒子
首先,假设粒子的运动是由速度矢量场决定的,速度矢量场是一个位置和时间的函数:
v
(
x
,
t
)
v(x,t)
v(x,t)
11.11.1 常微分方程 Ordinary Differential Equation(ODE)
计算质点随时间的位置需要解一阶常微分方程:
d x d t = x ˙ = v ( x , t ) \frac{dx}{dt}=\boldsymbol{\dot{x}}=v(x,t) dtdx=x˙=v(x,t)
“一阶”(First-order),是指所取的第一个导数。
“常”(Ordinary)意味着没有“部分”导数,即 x x x 只是 t t t 的一个函数。
11.11.2 求解粒子的位置
对于给定的粒子初始位置
x
0
x_0
x0,我们可以用正向数值积分法求解ODE。
11.11.3 欧拉方法
欧拉方法(又称前向欧拉、显式欧拉)
- 简单迭代法
- 常用的
- 非常不准确
- 通常情况不稳定
x t + Δ t = x t + Δ t x ˙ t x ˙ t + Δ t = x ˙ t + Δ t x ¨ t \begin{aligned}\boldsymbol{x}^{t+\Delta t}&=\boldsymbol{x}^{t}+\Delta t\boldsymbol{\dot{x}}^{t}\\\boldsymbol{\dot{x}}^{t+\Delta t}&=\boldsymbol{\dot{x}}^{t}+\Delta t\boldsymbol{\ddot{x}}^{t}\end{aligned} xt+Δtx˙t+Δt=xt+Δtx˙t=x˙t+Δtx¨t
11.11.4 欧拉方法-误差
在数值积分中,误差会累积
欧拉积分尤为糟糕
可通过减小 Δ t \Delta t Δt减小误差
11.11.5 欧拉方法的不稳定性
欧拉方法(显式/前向)
x t + Δ t = x t + Δ t v ( x , t ) x^{t+\Delta t}=x^t+\Delta t \boldsymbol{v}(\boldsymbol{x},t) xt+Δt=xt+Δtv(x,t)
两个关键问题:
不准确度随着时间的推移而增加,随着时间
Δ
t
\Delta t
Δt增加
不稳定性是一个常见的严重问题,它会导致模拟偏差(diverge)
11.11.6 误差和不稳定性
通过有限差分数值积分来解决会导致两个问题
- 误差
- 每次步骤的误差都会累积,随着模拟的进行,精度会降低
- 在图形应用中,精度可能不是关键
- 不稳定性
- 误差会加剧,即使底层系统不存在,模拟也会发生偏差
- 稳定性不足是模拟中的一个根本问题,不容忽视
11.12 对抗不稳定性
11.12.1 对抗不稳定性的几种方法
- 中点法/修正欧拉
- 平均起点和终点的速度
- 自适应步长
- 将一个步长和两个半步长进行递归比较,直到误差可以接受为止
- 隐式方法
- 在下一个步长中使用速度(硬性)
- 基于位置/Verlet合成
- 约束时间步后粒子的位置和速度
11.12.2 中点法 Midpoint Method
- 中点法
- 计算欧拉步长(a)
- 计算欧拉步中点的导数(b)
- 使用中点导数更新位置(c)
x m i d = x ( t ) + Δ t / 2 ⋅ v ( x ( t ) , t ) x ( t + Δ t ) = x ( t ) + Δ t ⋅ v ( x mid , t ) \begin{aligned} x_{\mathrm{mid}}& =x(t)+\Delta t/2\cdot v(x(t),t) \\ x(t+\Delta t)& =x(t)+\Delta t\cdot v(x_\text{mid},t) \end{aligned} xmidx(t+Δt)=x(t)+Δt/2⋅v(x(t),t)=x(t)+Δt⋅v(xmid,t)
11.12.3 修正欧拉 Modified Euler
- 修正欧拉
- 平均起点和终点的速度
- 更好的效果
x t + Δ t = x t + Δ t 2 ( x ˙ t + x ˙ t + Δ t ) x ˙ t + Δ t = x ˙ t + Δ t x ¨ t x t + Δ t = x t + Δ t x ˙ t + ( Δ t ) 2 2 x ¨ t \begin{gathered} \boldsymbol{x}^{t+\Delta t}=\boldsymbol{x}^t+\frac{\Delta t}2 (\dot{\boldsymbol{x}}^t+\dot{\boldsymbol{x}}^{t+\Delta t}) \\ \dot{\boldsymbol{x}}^{t+\Delta t}=\dot{\boldsymbol{x}}^t+\Delta t \ddot{\boldsymbol{x}}^t \\ \boldsymbol{x}^{t+\Delta t}=\boldsymbol{x}^t+\Delta t \dot{\boldsymbol{x}}^t+\frac{(\Delta t)^2}2 \ddot{\boldsymbol{x}}^t \end{gathered} xt+Δt=xt+2Δt(x˙t+x˙t+Δt)x˙t+Δt=x˙t+Δtx¨txt+Δt=xt+Δtx˙t+2(Δt)2x¨t
11.12.4 自适应步长 Adaptive Step Size
- 自适应步长
- 基于误差估计选择步长的技术
- 非常实用的技术
- 但可能需要很小的步骤!
- 重复,直到误差低于阈值:
- 计算 X T X_T XT 欧拉步长,size T T T
- 计算 X T / 2 X_{T/2} XT/2 两个欧拉步长,size T / 2 T/2 T/2
- 计算误差 ∣ ∣ X T − X T / 2 ∣ ∣ ||X_T-X_{T/2}|| ∣∣XT−XT/2∣∣
- 如果(误差>阈值),减小步长大小,然后重试
11.12.5 隐式欧拉方法 Implicit Euler Method
- 隐式方法
- 非正式地称为后向方法
- 使用下一步的速度和加速度,用于当前步
x t + Δ t = x t + Δ t x ˙ t + Δ t x ˙ t + Δ t = x ˙ t + Δ t x ¨ t + Δ t \begin{aligned}\boldsymbol{x}^{t+\Delta t}&=\boldsymbol{x}^{t}+\Delta t\boldsymbol{\dot{x}}^{t+\Delta t}\\\boldsymbol{\dot{x}}^{t+\Delta t}&=\boldsymbol{\dot{x}}^{t}+\Delta t\boldsymbol{\ddot{x}}^{t+\Delta t}\end{aligned} xt+Δtx˙t+Δt=xt+Δtx˙t+Δt=x˙t+Δtx¨t+Δt - 求解 x t + Δ t \boldsymbol{x}^{t+\Delta t} xt+Δt 和 x ˙ t + Δ t \boldsymbol{\dot{x}}^{t+\Delta t} x˙t+Δt 的非线性问题
- 使用找根算法,例如牛顿法
- 提供了更好的稳定性
如何确定/量化“稳定性”?
- 我们使用局部截断误差(每步)/总累积误差(综合)
- 绝对值无关紧要,但顺序是重要的
- 隐式欧拉是1阶的,这意味着
- 局部截断误差: O ( h 2 ) O(h^2) O(h2)
- 全局截断误差: O ( h ) O(h) O(h) ( h h h 是步数,即 Δ t \Delta t Δt)
- 对
O
(
h
)
O(h)
O(h) 的理解
如果我们把 h h h 减半,我们可以预期误差也减半
11.12.6 龙格-库塔法 Runge-kutta Families
一组解决ODE(常微分方程)的高级方法
- 特别擅长处理非线性问题
- 它的四阶版本是最广泛使用的,即 RK4
初始条件:
d
y
d
t
=
f
(
t
,
y
)
,
y
(
t
0
)
=
y
0
.
\frac{dy}{dt}=f(t,y),\quad y(t_0)=y_0.
dtdy=f(t,y),y(t0)=y0.
RK4 方法:
y
n
+
1
=
y
n
+
1
6
h
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
,
t
n
+
1
=
t
n
+
h
\begin{aligned}&y_{n+1}=y_n+\frac16h\left(k_1+2k_2+2k_3+k_4\right),\\&t_{n+1}=t_n+h\end{aligned}
yn+1=yn+61h(k1+2k2+2k3+k4),tn+1=tn+h
- k 1 = f ( t n , y n ) k_1= f(t_n,y_n) k1=f(tn,yn)
- k 2 = f ( t n + h 2 , y n + h k 1 2 ) k_2 = f\left(t_n + \frac{h}{2},y_n +h\frac{k_1}{2}\right) k2=f(tn+2h,yn+h2k1)
- k 3 = f ( t n + h 2 , y n + h k 2 2 ) k_3= f\left(t_n+\frac{h}{2},y_n+h\frac{k_2}{2}\right) k3=f(tn+2h,yn+h2k2)
- k 4 = f ( t n + h , y n + h k 3 ) k_4= f\left(t_n+h,y_n+hk_3\right) k4=f(tn+h,yn+hk3)
11.12.7 基于位置/Verlet集成
想法:
- 在修正欧拉前向步后,约束粒子的位置以防止发散和不稳定行为
- 使用受限位置来计算速度
- 这两种想法都会消耗能量,稳定
优缺点:
- 快速且简单
- 不以物理为基础,消耗能量(误差)
11.12.8 刚体模拟
简单的案例
- 类似于模拟粒子
- 考虑一下更多属性吧
d d t ( X θ X ˙ ω ) = ( X ˙ ω F / M Γ / I ) \left.\frac d{dt}\left(\begin{array}{c}\mathrm{X}\\\theta\\\mathrm{\dot{X}}\\\omega\end{array}\right.\right)=\left(\begin{array}{c}\mathrm{\dot{X}}\\\omega\\\mathrm{F}/M\\\Gamma/I\end{array}\right) dtd XθX˙ω = X˙ωF/MΓ/I
- X \mathrm{X} X:位置
- θ \theta θ:旋转角
- ω \omega ω:角速度
- F \mathrm{F} F:力
- Γ \Gamma Γ:转矩
- I I I:惯性动量
11.13 流体模拟
11.13.1 一种简单的基于位置的方法
关键思想
- 假设水是由小刚体球组成的
- 假设水不能被压缩(即恒定密度)
- 所以,只要密度在某个地方发生变化,就应该通过改变粒子的位置来“修正”
- 你需要知道与每个粒子位置相关的任何地方的密度梯度
- 更新?就是梯度下降!
11.13.2 欧拉方法 vs. 拉格朗日方法
模拟大量物质集合的两种不同观点
-
(“质点法””)拉格朗日方法
- 摄影师全程跟踪同一只鸟。
- 摄影师全程跟踪同一只鸟。
-
(“网格法”)欧拉方法
- 摄影师是静止不动的,只能拍摄所有经过一帧的鸟(比如说时间)。
- 摄影师是静止不动的,只能拍摄所有经过一帧的鸟(比如说时间)。
11.13.3 物点法 Material Point Method(MPM)
混合型,结合了欧拉和拉格朗日的观点
- 拉格朗日:考虑带有材料特性的粒子
- 欧拉:用网格做数值积分
- 交互作用:粒子向网格传递属性,网格进行更新,最后间接返回粒子