GAMES101(现代计算机图形学入门)笔记

GAMES101

课程链接:https://www.bilibili.com/video/BV1X7411F744

作业仓库:https://github.com/Tanyuu/GAMES101-Homework1-8

Lecture 02 Review of Linear Algebra(线性代数复习)

叉乘矩阵:

向量叉乘不同于点乘(ab=aTb=[xayaza][xayaza]\vec{a}\cdot\vec{b}=\vec{a}^T\vec{b}=\begin{bmatrix}x_a&y_a&z_a\end{bmatrix}\begin{bmatrix}x_a\\y_a\\z_a\end{bmatrix}),不是两个向量的矩阵运算。为了解决这个问题,我们引入了叉乘矩阵(Dual Matrix),对于三维向量的叉乘:

a×b=[a]×b=[0zayaza0xayaxa0][xayaza]\vec a\times\vec b=[a]_\times b= \begin{bmatrix} 0&-z_a&y_a\\ z_a&0&-x_a\\ -y_a&x_a&0 \end{bmatrix} \begin{bmatrix}x_a\\y_a\\z_a\end{bmatrix}

Lecture 03 Transformation(变换)

变换包含若干细分,如Viewing(视角变换)和Modeling(模型变换);

二维模型变换

仿射变换(Affine Transformations)与齐次坐标(Homogeneous Coordinates)

仿射变换形式如下:

[xy]=[abcd][xy]+[ef]\begin{bmatrix} x'\\y' \end{bmatrix} = \begin{bmatrix} a&b\\c&d \end{bmatrix} \begin{bmatrix} x\\y \end{bmatrix} + \begin{bmatrix} e\\f \end{bmatrix}

比如平移变换(Translation Transformation)(T(tx,ty)T{(t_x,t_y)})就是一种仿射变换:

[xy]=[1001][xy]+[txty]\begin{bmatrix} x'\\y' \end{bmatrix} = \begin{bmatrix} 1&0\\0&1 \end{bmatrix} \begin{bmatrix} x\\y \end{bmatrix} + \begin{bmatrix} t_x\\t_y \end{bmatrix}

对于这种形式,仿射变换并不能像其他变换一样用单纯的用矩阵乘法进行表示。为了方便表示,我们引入了齐次坐标(我更愿意翻译为“同坐标”),形式如下:

[xyω]\begin{bmatrix} x\\y\\\omega \end{bmatrix}

ω\omega 为0,则说明其表示的是向量,否则则为点。此外,若齐次坐标表示点时, [xyω]\begin{bmatrix} x\\y\\\omega \end{bmatrix} 等价于 ,[xωyω1]\begin{bmatrix}\frac x \omega\\\frac y\omega\\1\end{bmatrix}即代表点 (xω,yω)(\frac x\omega,\frac y\omega)

齐次坐标还有一些神奇的性质,它使得点与向量间的运算封闭:q

向量+向量=向量

点+向量=点

向量-向量=向量

点-向量=点+(-向量)=点

向量-点=点

点+点=中点

齐次坐标表示下,平移变换(T(tx,ty)T{(t_x,t_y)})即可用矩阵乘法表示为:

[xyω]=[10tx01ty001][xyω]=[x+ωtxy+ωtyω]\begin{bmatrix} x'\\y'\\\omega' \end{bmatrix} = \begin{bmatrix} 1&0&t_x\\ 0&1&t_y\\ 0&0&1 \end{bmatrix} \begin{bmatrix} x\\y\\\omega \end{bmatrix} = \begin{bmatrix} x+\omega t_x\\y+\omega t_y\\\omega \end{bmatrix}

解决了前述问题。

缩放变换(Scale Transformation)与反射变换(Reflection Transformation)

image.png

横轴缩放 sxs_x 倍,纵轴缩放 sys_y 倍的缩放变换 S(sx,sy)S{(s_x,s_y)} 表示为:

[xyω]=[sx000sy0001][xyω]=[xsxysyω]\begin{bmatrix} x'\\y'\\\omega' \end{bmatrix} = \begin{bmatrix} s_x&0&0\\ 0&s_y&0\\ 0&0&1 \end{bmatrix} \begin{bmatrix} x\\y\\\omega \end{bmatrix}=\begin{bmatrix} xs_x\\ys_y\\\omega \end{bmatrix}

s<0s_*<0 则认为是反射变换;

切变变换(Shear Transformation)

image.png

符合 x=x+axy,y=y+ayxx'=x+a_xy,y'=y+a_yx 的切变变换表示为:

[xyω]=[1ax0ay10001][xyω]=[x+axyy+ayxω]\begin{bmatrix} x'\\y'\\\omega' \end{bmatrix} = \begin{bmatrix} 1&a_x&0\\ a_y&1&0\\ 0&0&1 \end{bmatrix} \begin{bmatrix} x\\y\\\omega \end{bmatrix}=\begin{bmatrix} x+a_xy\\y+a_yx\\\omega \end{bmatrix}

旋转变换(Rotate Transformation)

旋转变换默认绕 (0,0)(0,0) 进行,以 x+x^+ 为起始,逆时针为正;

image.png

旋转 θ\theta 的旋转变换 RθR_\theta 表示为:

[xyω]=[cosθsinθ0sinθcosθ0001][xyω]\begin{bmatrix} x'\\y'\\\omega' \end{bmatrix} = \begin{bmatrix} \cos\theta&-\sin\theta&0\\ \sin\theta&\cos\theta&0\\ 0&0&1 \end{bmatrix} \begin{bmatrix} x\\y\\\omega \end{bmatrix}

结合下述的逆变换,我们可以发现:

R(θ)=[cosθsinθsinθcosθ]=R(θ)TR(θ)=R(θ)1R{(-\theta)}= \begin{bmatrix} \cos\theta&\sin\theta\\ -\sin\theta&\cos\theta \end{bmatrix} =R{(\theta)}^T\\ R{(-\theta)}=R{(\theta)}^{-1}

逆变换(Inverse)

AA'AAMM 变换变换得到的(A=MAA'=MA),即 AA 可由 AA' 进行 M1M^{-1} 变换得到(A=M1AA=M^{-1}A'

变换顺序(Transform Ordering)

若先进行 R(45°)R{(45\degree)} ,再进行 T(1,0)T{(1,0)} ,则过程表示为:

AR(45°)R(45°)AT(1,0)T(1,0)R(45°)AA\xrightarrow{R{(45\degree)}}R{(45\degree)}A\xrightarrow{T{(1,0)}}T{(1,0)}R{(45\degree)}A

最终的式子 T(1,0)R(45°)AT{(1,0)}R{(45\degree)}A 实际上运算顺序是从右至左,但是由于矩阵乘法具有结合律,可以先计算 T(1,0)R(45°)T{(1,0)}R{(45\degree)} ,再右乘上 AA

一般地,如果向量 AA 依此经过变换 M1,M2,,MnM_1,M_2,\dots,M_n ,则可以表示为:

Mn(M2(M1(A)))=(MnM2M1)AM_n(\dots M_2(M_1(A)))=(M_n\dots M_2M_1)A

考虑平移矩阵(T(tx,ty)T_{(t_x,t_y)}):

[xyω]=[10tx01ty001][xyω]\begin{bmatrix} x'\\y'\\\omega' \end{bmatrix} = \begin{bmatrix} 1&0&t_x\\ 0&1&t_y\\ 0&0&1 \end{bmatrix} \begin{bmatrix} x\\y\\\omega \end{bmatrix}

其左上的 (2,2)(2,2) 矩阵是单位阵,可以换成其他几种变换矩阵的左上 (2,2)(2,2) 矩阵,实际上的运算顺序为先做线性变换,再做平移变换;

上述的旋转变换是绕原点进行的,如果要绕 cc 点进行,可以先通过平移将待旋转图形平移 c-c ,即将旋转中心平移至原点,再进行旋转,再平移回原位,以 cc 为旋转中心旋转 α\alpha 表示如下:

AT(c)T(c)AR(α)R(α)T(c)AT(c)T(c)R(α)T(c)AA\xrightarrow{T{(-c)}}T{(-c)}A\xrightarrow{R{(\alpha)}}R{(\alpha)}T{(c)}A\xrightarrow{T{(-c)}}T{(c)}R{(\alpha)}T{(-c)}A

Lecture 04 Transformation Cont.(变换 续)

三维模型变换

平移变换

T(tx,ty,tz)=[100tx010ty001tz0001]T{(t_x,t_y,t_z)}= \begin{bmatrix} 1&0&0&t_x\\ 0&1&0&t_y\\ 0&0&1&t_z\\ 0&0&0&1 \end{bmatrix}

缩放变换与反射变换

S(sx,sy,sz)=[sx0000sy0000sz00001]S{(s_x,s_y,s_z)}= \begin{bmatrix} s_x&0&0&0\\ 0&s_y&0&0\\ 0&0&s_z&0\\ 0&0&0&1 \end{bmatrix}

旋转变换

绕坐标轴旋转

Rx(α)=[10000cosαsinα00sinαcosα00001]Ry(α)=[cosα0sinα00100sinα0cosα00001]Rz(α)=[cosαsinα00sinαcosα0000100001]R_x(\alpha)= \begin{bmatrix} 1&0&0&0\\ 0&\cos\alpha&-\sin\alpha&0\\ 0&\sin\alpha&\cos\alpha&0\\ 0&0&0&1 \end{bmatrix}\\ R_y(\alpha)= \begin{bmatrix} \cos\alpha&0&\sin\alpha&0\\ 0&1&0&0\\ -\sin\alpha&0&\cos\alpha&0\\ 0&0&0&1 \end{bmatrix}\\ R_z(\alpha)= \begin{bmatrix} \cos\alpha&-\sin\alpha&0&0\\ \sin\alpha&\cos\alpha&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{bmatrix}\\

由于循环对称性,x+=y+×z+,z+=x+×y+x^+=y^+\times z^+,z^+=x^+\times y^+ ,而 y+=z+×x+y^+=z^+\times x^+ ,所以 Ry(α)R_y(\alpha) 形式不同。直观来说,单位正方体绕坐标轴旋转一个锐角时,绕 x+,z+x^+,z^+ 的情况都是从较小轴转向较大轴,而 y+y^+ 的情况由于循环对称而相反;

欧拉角(Euler angles)旋转变换

将任意旋转分解为 x,y,zx,y,z 轴上的旋转 α,β,γ\alpha,\beta,\gamma 并分别单独变换;

Rxyz(α,β,γ)=Rx(α)Ry(β)Rz(γ)R_{xyz}(\alpha,\beta,\gamma)=R_x(\alpha)R_y(\beta)R_z(\gamma)

罗德里格斯旋转公式(Rodrigues’ Rotation Formula)

绕旋转轴 n\vec n (默认直线过原点,且 n\vec n 为三维单位列向量),\bigodot 方向逆时针旋转 α\alpha 角度的旋转矩阵:

R(n,α)=[0cos(α)I000001]+(1cos(α))[n0][n0]T+sin(α)[0nzny0nz0nx0nynx000000]R(\vec n,\alpha)= \begin{bmatrix} &&&0\\ &\cos(\alpha)I&&0\\ &&&0\\ 0&0&0&1 \end{bmatrix} + (1-\cos(\alpha)) \begin{bmatrix} \vec n\\ 0 \end{bmatrix} \begin{bmatrix} \vec n\\ 0 \end{bmatrix} ^T+\sin(\alpha) \begin{bmatrix} 0&-n_z&n_y&0\\ n_z&0&-n_x&0\\ -n_y&n_x&0&0\\ 0&0&0&0 \end{bmatrix}

其中 II 为单位阵,第三项的左上 3×33\times 3sin(α)[n]×\sin(\alpha)[\vec n]_\times

如果旋转轴不经过原点,可以再通过两次平移完成旋转过程;

视角变换

视图变换(相机变换)(View / Camera Transformation)

将三维的点集由世界坐标变换为相机坐标;

相机有如下几种属性:位置(Position)e\vec e ,视线方向(Look-at / Gaze direction)g^\hat g ,向上方向(Up direction)t^\hat t (假设与视线方向垂直)

我们将原坐标系变换为相机坐标系,并符合如下规则:相机在原点、向上方向为 y+y^+ ,视线方向为 zz^-x+ is g×t,y+ is t,z+ is gx^+\text{ is }g\times t,y^+\text{ is }t,z^+\text{ is }-g);

假设变换矩阵为 MviewM_{view} ,考虑其的形式(Mview=RviewTviewM_{view}=R_{view}T_{view}):首先需要将相机点移到原点(Tview=T(e)T_{view}=T_{(-\vec e)}),再进行坐标轴变换;

考虑 Rview1R_{view}^{-1} 即从相机坐标系变换回原坐标系:参考

Rview1=[xg^×t^xt^xg^0yg^×t^yt^yg^0zg^×t^zt^zg^00001]R_{view}^{-1}= \begin{bmatrix} x_{\hat g\times\hat t}&x_{\hat t}&x_{-\hat g}&0\\ y_{\hat g\times\hat t}&y_{\hat t}&y_{-\hat g}&0\\ z_{\hat g\times\hat t}&z_{\hat t}&z_{-\hat g}&0\\ 0&0&0&1 \end{bmatrix}

得:

Rview=[xg^×t^yg^×t^zg^×t^0xt^yt^zt^0xg^yg^zg^00001]R_{view}= \begin{bmatrix} x_{\hat g\times\hat t}&y_{\hat g\times\hat t}&z_{\hat g\times\hat t}&0\\ x_{\hat t}&y_{\hat t}&z_{\hat t}&0\\ x_{-\hat g}&y_{-\hat g}&z_{-\hat g}&0\\ 0&0&0&1 \end{bmatrix}

投影变换(Projection Transformation)

正交投影(Orthographic Projection)

简略来说,进行相机变换之后,将 zz 坐标写为0, x,yx,y 坐标的 (,+)2(-\infty,+\infty)^2 映射到 [1,1]2[-1,1]^2 中;

但是正式来说,首先定义相机变换后空间中的一个立方体(视景体),具有在 xx 轴上的左右 [l,r][l,r] ,在 yy 轴上的上下 [b,t][b,t] ,在 zz 轴上的远近 [f,n][f,n] ,然后通过平移将此立方体几何中心移至原点,再通过缩放将其变换为 [1,1]3[-1,1]^3 的标准正方体(Canonical Cube);

Mortho=[2rl00002tb00002nf00001][100r+l2010t+b2001n+f20001]M_{ortho}= \begin{bmatrix} \frac 2 {r-l}&0&0&0\\ 0&\frac 2 {t-b}&0&0\\ 0&0&\frac 2 {n-f}&0\\ 0&0&0&1 \end{bmatrix} \begin{bmatrix} 1&0&0&-\frac {r+l}2\\ 0&1&0&-\frac {t+b}2\\ 0&0&1&-\frac {n+f}2\\ 0&0&0&1 \end{bmatrix}

透视投影(Perspective Projection)

透视投影是最常见的投影方式,符合近大远小的规律,平行线不再平行;

透视投影的景视体是一个以原点为锥点的截锥体(Frustum),我们首先将截锥体变换成立方体,再进行正交投影,就完成了透视投影的过程;变换过程要保证几条规律:将远平面挤压至与近平面大小相同,挤压过程中 ff 值不变,远平面中心点空间位置不变;

有:

[xyz1]=[nxny?z]\begin{bmatrix} x'\\ y'\\ z'\\ 1 \end{bmatrix} = \begin{bmatrix} n x\\ n y\\ ? \\ z \end{bmatrix}

考虑 ?? 的取值:

对于远近平面上的两点,需要保证:

[xyn1][nxnyn2n],[xyf1][nxnyf2f]\begin{bmatrix} x\\ y\\ n\\ 1 \end{bmatrix} \to \begin{bmatrix} n x\\ n y\\ n^2 \\ n \end{bmatrix} ,\\ \begin{bmatrix} x\\ y\\ f\\ 1 \end{bmatrix} \to \begin{bmatrix} nx\\ ny\\ f^2 \\ f \end{bmatrix}

得到了方程组:

{Ax+By+Cn+D=n2Ax+By+Cf+D=f2\begin{cases} Ax+By+Cn+D=n^2\\ Ax+By+Cf+D=f^2 \end{cases}

其中 x,y,n,fx,y,n,f 为已知量,A,B,C,DA,B,C,D 为未知量,可以解得 A=0,B=0,C=n+f,D=nfA=0,B=0,C=n+f,D=-nf

则:

Mpersportho=[n0000n0000n+fnf0010]M_{persp\to ortho}= \begin{bmatrix} n&0&0&0\\ 0&n&0&0\\ 0&0&n+f&-nf\\ 0&0&1&0 \end{bmatrix}

Mpersp=MorthoMpersporthoM_{persp}=M_{ortho}M_{persp\to ortho} (上式要求 n>f ,即 n, f <0 时的情况);

实际在这种模型中,zz 轴的映射关系为:

image.png

Lecture 05 Rasterization 1 (Triangles)(光栅化1-三角形)

视角变换之后,所有的物体均已经变换至 [1,1]3{[-1,1]}^3 的空间中,接下来应该考虑如何将其绘制出来;

宽高比 aspect_ratio=width/height ;

视角(field-of-view),fovY为垂直视场角,fovX为水平视场角,此课程中fov默认为垂直视场角;

在视角变换之后,相机位于原点, y+y^{+} 轴为上, zz^{-} 轴为朝向,即有 l=r,b=tl=-r,b=-t

考虑 fovY 、aspect_ratio 与 l, r, b, t 的关系,如图所示:

KX7XK9_UX9ZJIO_G_NDAH.jpg

tanfovY2=tnaspect=rt\tan\frac{fovY}2=\frac t{|n|}\\ aspect=\frac r t

我们将屏幕视为二维数组,光栅化(rasterize)即为在屏幕上绘制的过程;

此课程中规定,屏幕空间的像素坐标视作分别代表X、Y轴位置的整数数对,如下图蓝色位置为 (2,1) ;

image.png

对于 width×height\text {width}\times \text{height} 的屏幕,其像素坐标范围为 (0, 0) 到 (width-1, height-1) ;

虽然像素坐标为整数坐标,但是坐标为 (x, y) 的像素中心在 (x+0.5, y+0.5) ;

屏幕覆盖范围为 (0, 0) 到 (width, height) ;

我们需要完成从 [1,1]3{[-1,1]}^{3}[0,width]×[0,height][0,width]\times[0,height] 的映射;

我们先不考虑 z 坐标,将 xy 坐标拉伸至屏幕范围,有视口变换矩阵:

Mviewport=[width200width20height20height200100001]M_{viewport}= \begin{bmatrix} \frac {width}2&0&0&\frac {width}2\\ 0&\frac {height}2&0&\frac {height}2\\ 0&0&1&0\\ 0&0&0&1 \end{bmatrix}

三角形存在很多独特的性质:

  1. 保证共面(guaranteed to be planar);
  2. 保证凸形,便于判断内外(well-defined interior);
  3. 利用重心坐标系便于插值以实现渐变(barycentric interpolation);

接下来我们首先考虑三角形的光栅化;

采样(Sample)

定义布尔函数 inside(tri, x, y) ,其中xy不要求整数,含义为:

inside(tri,x,y)={1(x,y) in triangle tri0otherwise\text{inside}(tri,x,y)= \begin{cases} 1&(x,y) \text{ in triangle }tri\\ 0&\text{otherwise} \end{cases}

1
2
3
for(int x = 0; x < xmax; ++x)
for(int y = 0; y < ymax; ++y)
image[x][y] = inside(tri, x + 0.5, y + 0.5);

可以用同侧法判断点是否在三角形内;

关于点在边界上的情况,按统一标准处理即可;

考虑枚举所有像素的低效,我们可以只扫描三角形轴向包围盒(axis-aligned bounding box,AABB)中的像素([Xmin,Xmax]×[Ymin,Ymax][X_{min},X_{max}]\times[Y_{min},Y_{max}]);

对于细长且旋转的三角形,AABB 策略会被降低效率,可以使用增量遍历策略:

image.png

Lecture 06 Rasterization 2 (Antialiasing and Z-Buffering)(光栅化2-反走样和深度缓存)

反走样

走样:CG中不符合预期的采样(Sampling Artifacts (Errors/Mistakes/Inaccuracies) in Computer Graphics),包括锯齿、摩尔纹、车轮效应等;

对于余弦函数,我们称满足 cos2πfx\cos2\pi fx 形式的函数为频率为 ff 的函数;在采样过程中,高频部分会产生失真;

同样采样方法采样两种不同频率的信号出现相同的信息称之为走样;

滤波即为去掉某些频率的内容;

时域上的乘积等于频域上的卷积,所以采样就是重复原始信号的频谱;

image.png

如果采样率不足(fs较小)或高频丰富(f0较大)时,便会发生频谱混叠(aliasing),即为走样;

image.png

进一步说,除了增加采样率外,还可以通过先低通滤波再采样的方法;

image.png

在实际操作中,我们通过对inside函数进行卷积来实现低通滤波,理想情况下,inside函数可以返回像素面积中被三角形覆盖的比例,实际实现上有以下几种途径:

  1. 超采样(Super-Sample Anti-Aliasing,SSAA)

    我们将每个像素划分为更多小的像素,分别对小像素进行采样后进行平均;

  2. 多重采样(Multi-Sample AA,MSAA)

    不同于SSAA,SSAA对于所有子采样点着色,而MSAA只对当前模型所更新的子采样点进行着色;

    SSAA与MSAA都需要进行解析(Resolve),信号必须重新采样到指定的分辨率,这样我们才可以显示它;

    SSAA和MSAA会大量消耗算力,实际实现上可能通过不规则分布和像素复用等方法降低消耗;

  3. 快速近似(Fast Approximate AA,FXAA)

    在光栅化后带锯齿的图像上寻找边界,将边界替换为没有锯齿的边界;

  4. 时域抗锯齿(Temporal AA,TAA)

    类似于在时域上进行多重采样;

  5. 深度学习超分辨率(Deep Learning Super Sampling,DLSS)

深度绘制

画家算法(Painters’ Algorithm)

从远到近绘制,近的物体覆盖远的物体;

同距离物体的绘制顺序并不随意,会出现遮挡问题:

image.png

也可能出现无法解决的深度组合:

image.png

所以实际解决时无法使用此算法;

深度缓存

为了简单起见,此部分我们认为 z 总是正的,z 越小物体离相机越近

在渲染图像时同时生成两个结果:对于想要的图像存储在 frame buffer ,对于每个像素的深度信息存储在 depth buffer;

深度缓存对于每个像素,只维护此帧中最小的 z 值;

1
2
3
4
5
6
7
8
9
lnitialize depth buffer to ∞;
During rasterization:
foreach (triangle T)
foreach (sample (x,y,z) in T)
if (z < zbuffer[x,y]) // closest sample so far
framebuffer[x,y] = rgb; // update color
zbuffer[×,y] = z; // update depth
else
; // do nothing, this sample is occluded

由画家算法的 O(nlogn)O(n\log n) 复杂度降低为了 O(n)O(n) 复杂度,并且在不存在z-fighting假设下,此算法对三角形顺序不敏感;

在 MSAA 策略下,z-buffer 粒度会上升至每一个采样点;

z-buffer 不能处理透明物体的深度问题;

Lecture 07 Shading 1 (Illumination, Shading and Graphics Pipeline)(着色1-照明、布林-冯反射模型和阴影)

着色:对不同物体应用不同材质(material)的过程;

布林-冯反射模型(Blinn-Phong Reflectance Model)

将从一个点接受的光分为三部分:镜面反射(specular highlights)、漫反射(diffuse reflection)、环境光(ambient lighting);环境光实际上是由环境和周边物体的反射混合形成的,在此模型中被简化为一个常量;

image.png

有如下输入量:相机方向单位向量 v^\hat v ,表面法向单位向量 n^\hat n ,对于每个光源的光源方向单位向量 l^\hat l ,表面参数 kk (与颜色、材质、反射方式等有关的反射率);

模型是一种局部光照模型,指只考虑着色点、相机和光源之间的关系,不考虑着色点是否被其他物体遮挡;

对于漫反射项 LdL_d 、镜面反射项 LsL_s 、环境光项 LaL_a ,布林-冯模型有

L=Ld+Ls+LaL=L_d+L_s+L_a

漫反射项

漫反射指对于一条打向着色点的光线,其被均匀地反射到各个方向;

image.png

朗伯余弦定律(Lambert’s cosine law):单位面积接受的辐射能量与法向和光向的夹角有关,夹角有 cosθ=l^n^\cos\theta =\hat l\cdot \hat n

image.png

光衰减(平方反比定律):与点光源距离为1的单位面积接受的能量为 I ,则距离为 r 的单位面积接受的能量为 I/r2I/r^2

综上有漫反射光公式

Ld=kd(I/r2)max(0,n^l^)L_d=k_d(I/r^2)\max(0,\hat n\cdot \hat l)

漫反射光强 = 漫反射系数 * 原始光强/(光源距离^2 )* cos法光角;

高光项

image.png

v^\hat v 和 R 足够接近时得到的高光;v^\hat v 和 R 足够接近也就意味着 h^\hat hv^\hat vl^\hat l 的半程向量)和 n^\hat n 足够接近;

image.png

h^=bisector(v,l)=v^+l^v^+l^Ls=ks(I/r2)max(0,cosα)p=ks(I/r2)max(0,n^h^)p\hat h=\text{bisector}(v,l)\\ =\frac{\hat v+\hat l}{||\hat v+\hat l||}\\ L_s=k_s(I/r^2)\max(0,\cos\alpha)^p\\ =k_s(I/r^2)\max(0,\hat n\cdot \hat h)^p

镜面反射光强 = 镜面反射系数 * 原始光强/(光源距离^2 )* cos法半角;

通常认为高光为白色,镜面反射系数取白色光反射系数;

朗伯余弦定律在镜面反射项中被简化;

如果通过 R 与 v^\hat v 的夹角计算,会导致在计算 R 上耗费时间,所以实践中并不采用;

上式中的 pp 为镜面指数,用于加快余弦项的衰减速度,通常在100~200之间;镜面指数越大,高光面积越小;

image.png

环境光项

La=kaIaL_a=k_aI_a

此模型中环境光与光向、法向、视向无关,视作一常数;

阴影映射(Shadow Mapping)

如果一个点不在阴影中,那么说明这个点对于相机和光源都是可见的;

阴影映射用于处理单个点光源的二值硬阴影(软阴影越靠近物体边界越浅,来自于光源非点产生的半影);

在点光源处放置一个摄像机,对场景光栅化,存储每个可见点(与光源)的深度信息(Z-buffer);

从相机出发再次光栅化,对于可见点将其投影回光源判断其深度与存储值是否一致,一致即说明该点无阴影,不一致即说明该点在阴影下;

阴影映射本身受制于点光源处光栅化的分辨率,生成的阴影会包含锯齿;

Lecture 08 Shading 2 (Shading, Pipeline and Texture Mapping)(着色2-着色频率和图形管线)

着色频率(Shading Frequencies)

着色频率是指空间上进行着色的密度,有以下几种主要方法:

  1. 平面着色(flat shading)

    对于每个三角形进行着色计算,将结果应用于整个三角形上;

  2. 顶点着色(gouraud shading)

    对于每个顶点进行着色计算,三角形内的着色由插值得出;

  3. 像素着色(Phong shading)与布林-冯模型的冯是同一个人,与着色模型无关

    对于每个像素进行着色计算;

在面数较少时,像素着色频率越高效果越逼真,但在高面数情况下效果相近;

获取顶点/像素的法线

  1. 如果已知模型的原始几何属性(e.g. 球、平面等),可以直接通过几何方法计算法线;

  2. 一般情况,可以认为顶点法线是其相邻面法线的平均值
    image.png

    Nv=NiNiN_v=\frac {\sum N_i}{||\sum N_i||}

  3. 为了获取更准确的法线,可以按三角形面积等进行加权平均;

  4. 像素的法线可以通过重心插值获取;

图形(实时渲染)管线(Graphics (Realtime Rendering) Pipeline)

实时渲染管线指从场景到图像的过程;

image.png

http://shadertoy.com/view/Id3Gz2

Lecture 09 Shading 3 (Texture Mapping Cont.)(着色3-纹理映射)

三维物体的表面是二维的,表面可以与一张图形成映射关系;我们可以定义一张图,图上的每个三角形都映射到模型上的三角形,图上三角形的纹理就代表了模型上对应三角形的纹理;这张图就是uv贴图,得名于其两个坐标轴一般用 u 和 v 表示,无论其是方形或是矩形,通常认为范围有 u,v[0,1]u,v\in[0,1]

如果一块纹理在映射时,与上下左右的自己衔接不会产生缝隙,则被称为四方连续(tilable texture);

我们假设已经掌握了这种映射关系,在已知三角形顶点坐标的情况下,获取三角形内部点的坐标需要通过重心坐标(Barycentric Coordinates)来完成;

重心坐标(Barycentric coordinates)

重心坐标可以用于三角形内部的插值,在顶点处处理出的值(纹理映射、颜色、法向等)可以平滑地过渡到三角形内的任何一点;

image.png

ABC\triangle ABC 所在平面内任何一点 (x,y)(x,y) ,都可以表示成上图的形式;

由于和等于1的限制,实际上坐标由两个参数即可描述;

如果点在三角形内,则有 α,β,γ0\alpha ,\beta, \gamma\geqslant0

对于平面内任意一点 P ,考虑 PAC,PCB,PBA\triangle PAC,\triangle PCB, \triangle PBA 的有向面积(一一对应) AB,AA,ACA_B, A_A, A_C ,有

αP=AAAA+AB+ACβP=ABAA+AB+ACγP=ACAA+AB+AC\alpha_P=\frac {A_A}{A_A+A_B+A_C}\\ \beta_P=\frac {A_B}{A_A+A_B+A_C}\\ \gamma_P=\frac {A_C}{A_A+A_B+A_C}

从坐标角度考虑,对于平面内任意一点 P(x,y)P(x,y) ,其重心坐标有

αP=(xxB)(yCyB)+(yyB)(xCxB)(xAxB)(yCyB)+(yAyB)(xCxB)βP=(xxC)(yAyC)+(yyC)(xAxC)(xBxC)(yAyC)+(yByC)(xAxC)γP=1αPβP\alpha_P=\frac{-(x-x_B)(y_C-y_B)+(y-y_B)(x_C-x_B)}{-(x_A-x_B)(y_C-y_B)+(y_A-y_B)(x_C-x_B)}\\ \beta_P=\frac{-(x-x_C)(y_A-y_C)+(y-y_C)(x_A-x_C)}{-(x_B-x_C)(y_A-y_C)+(y_B-y_C)(x_A-x_C)}\\ \gamma_P=1-\alpha_P-\beta_P

ABC\triangle ABC 的重心有坐标 (13,13,13)(\frac 13,\frac 13,\frac 13)

对于三角形内一点 P(α,β,γ)P(\alpha,\beta,\gamma) 和三角形顶点上的属性 VA,VB,VCV_A,V_B,V_C ,有插值 VP=αVA+βVB+γVCV_P=\alpha V_A+\beta V_B+\gamma V_C

在投影变换下,不能够保证重心坐标不变;需要在三维空间进行插值,或对插值结果需要进行矫正;参考资料

1z0=α1z1+β1z2+γ1z3b0z0=αb1z1+βb2z2+γb3z3\frac 1{z_0}=\alpha\frac 1{z_1}+\beta\frac 1{z_2}+\gamma\frac 1{z_3}\\ \frac {b_0}{z_0}=\alpha\frac {b_1}{z_1}+\beta\frac {b_2}{z_2}+\gamma\frac {b_3}{z_3}

其中 b0b_0 即做透视校正插值后的正确顶点属性,b1,b2,b3b_1,b_2,b_3 即该三角形三个顶点的顶点属性,z0z_0 即做透视校正插值后的正确深度值,[α,β,γ][\alpha,\beta,\gamma] 即插值位置的重心坐标,z1,z2,z3z_1,z_2,z_3 即该三角形三个顶点的深度值。

纹理应用

1
2
3
4
for each rasterized screen sample (x,y):			// a pixel's center
(u,v) = evaluate texture coordinate at (x,y); // by barycentric coordinates
texcolor = texture.sample(u,v);
set sample's color to texcolor; // about kd in Blinn-Phone

image.png

一个像素与一个纹理元素(texel)(纹理像素)的大小关系并不固定,需要进行“纹理缩放”;

纹理放大

如果纹理分辨率低于渲染所需的分辨率,就会产生不符合预期的效果;

对于计算出来的浮点 (u,v)(u,v) ,我们需要选择一种合适的策略将其映射到纹素上;

双线性插值(Bilinear)

image.png

红点为计算出的 (u,v)(u,v) ,将其距离左下最近纹素的水平距离标记为 s ,竖直距离标记为 t ;

定义一个运算:线性插值(linear interpolation),对于 x[0,1]x\in[0,1] 和两个值 v1,v0v_1,v_0 ,有

lerp(x,v0,v1)=v0+x(v1v0)\text{lerp}(x,v_0,v_1)=v_0+x(v_1-v_0)

对双线性插值,定义有

u0=lerp(s,u00,u10)u1=lerp(s,u01,u11)uuv=lerp(t,u0,u1)u_0=\text{lerp}(s,u_{00},u_{10})\\ u_1=\text{lerp}(s,u_{01},u_{11})\\ u_{uv}=\text{lerp}(t,u_{0},u_{1})

双立方插值(Bicubic)区别在于考虑邻近的16个纹素进行立方插值;

纹理缩小

如果一个纹素的大小远小于一个像素的大小,那么便会产生走样现象,可以通过反走样方法进行处理,但是比较耗费资源,这里使用另一种策略;

Mipmap

我们将像素范围内的纹素近似平均值作为返回值;

Mipmap是一种一种快速的、近似的、正方形的范围查询;

我们需要存储纹理图的低分辨率版本(每级边长折半,多花费1/3的空间);

image.png

接下来计算单个像素的覆盖面积:考虑一个采样点和其水平、垂直方向的两个相邻采样点,这三个点映射在纹理图上的坐标;

image.png

L=max((dudx)2+(dvdx)2,(dudy)2+(dvdy)2)L=\max(\sqrt{(\frac {du}{dx})^2+(\frac {dv}{dx})^2},\sqrt{(\frac {du}{dy})^2+(\frac {dv}{dy})^2})

L 即可被近似为像素所在纹理上正方形的边长;

image.png

Mipmap级数 D 有 D=log2LD=\log_2L ,意味着对应区域在 Level D 的 Mipmap 上为一个像素,在该层进行查询即可;

三线性插值(Trilinear interpolation)

对于非整数的 D ,我们可以对 Mipmap 进行层间插值,这个过程被称为三线性插值(水平,垂直,层间);

image.png

Mipmap 会产生 Overblur (过度模糊)现象,由于在方形确定和三线性插值过程中使用了大量的近似;

image.png

各向异性过滤(Anisotropiic Filtering,Ripmap)

各向异性指在各个方向上其表现不同,在这里特指两个轴向方向上的表现不同;

各向异性过滤可以解决上图中的轴向矩形区域,对于斜向区域仍然存在问题;

image.png

不同于 Mipmap 只存储了上图中对角线上的图片,各向异性过滤通过存储不同长宽比的原始图像来实现轴向矩形区域的查询(多花费3倍的空间);

各向异性过滤的层数选项:对于4x各向异性过滤指取上图左上 3×33\times3 的图片区域进行计算,标记的是横纵压缩倍数;高倍数的各向异性过滤对计算能力消耗很小,但对显存消耗较大;

更进一步可以使用 EWA 过滤(EWA filtering)来实现真正意义上的各向同性,在此不作赘述;

环境贴图(环境光照)

在现代 GPU 中,纹理可以抽象地被表述成可以快速进行范围查询的一块数据;

如果我们通过纹理来描述环境光情况,再用这个纹理对物体进行渲染,便可以实现环境光照的效果;通过此种方法进行模拟,需要假设环境光来自无限远;

假设有一个镜子的球,将球面上反射的环境光记录为环境光信息,这就是球环境贴图(Spherical Environment Map);

image.png

接下来对球面进行展开存储,在类似于世界地图的投影方式中,两极部分会被扭曲变形;考虑使用立方体贴图(Cube Map)进行存储:

image.png

在球外添加一个正方体包围盒,从球心发出的射线打到包围盒上的位置即为该方向上的球面映射到立方体贴图上的位置;

凹凸贴图(法线贴图)(Bump / Normal Mapping)

纹理不止可以描述颜色,如果通过纹理记录表面相对高度,便可以在不增加面数的情况下来定义更精细的(假)模型;

image.png

假设有法线贴图 p=h(u,v)p=h(u,v) ,在一点 (u,v)(u,v) 处,原始法线为 (0,0,1)(0,0,1) ,经过法线贴图修饰后法线有

dpdu=c1[h(u+1,v)h(u,v)]dpdv=c1[h(u,v+1)h(u,v)]n^=(dpdu,dpdv,1).normalized()\frac {dp}{du}=c_1[h(u+1,v)-h(u,v)]\\ \frac {dp}{dv}=c_1[h(u,v+1)-h(u,v)]\\ \hat n=(-\frac {dp}{du},-\frac {dp}{dv},1).\text{normalized()}

c1,c2为在u,v两个方向上的影响系数,代表纹理定义的法线对于展示物体法线的影响系数;

由于原始法线 (0,0,1)(0,0,1) 假设,实际上计算出的法线是局部坐标系(切线空间)下的,需要再通过TBN矩阵变换回原坐标系;

切线空间与TBN矩阵

参考资料

TBN矩阵指的是由切线(Tangent)、副切线(Bitangent)(有些地方也叫副法线)、法线(Normal),组成的3x3的矩阵(这里的切线、副切线、法线全部指的是世界空间的向量,不是切线空间下的,TBN三向量正交且归一);

参考作业代码,有

T=[NxNyNx2+Nz2Nx2+Nz2NyNzNx2+Nz2]B=N×TT= \begin{bmatrix} \frac{N_xN_y}{\sqrt{N_x^2+N_z^2}}\\ \sqrt{N_x^2+N_z^2}\\ \frac{N_yN_z}{\sqrt{N_x^2+N_z^2}} \end{bmatrix}\\ B=N\times T

变换过程有

nWorld=[TxBxNxTyByNyTzBzNz]nTBNn_{World}= \begin{bmatrix} T_x & B_x & N_x\\ T_y & B_y & N_y\\ T_z & B_z & N_z \end{bmatrix} n_{TBN}

位移贴图(置换贴图)(Displacement Mapping)

相比于法线贴图只是修饰法线的结果,位移贴图是对顶点位置的移动;

观察边缘和阴影:

image.png

其代价为要求模型本身面数足够;为了解决这个问题 DirectX 使用动态曲面细分降低消耗;

Others

三维程序化噪声(3D Procedural Noise)

预计算着色(Precomputed Shading)、环境光遮蔽纹理(Ambient occlusion texture map)

体渲染(Volume Rendering)

Lecture 10 Geometry 1 (Introduction)(几何1-隐式描述)

隐式(Implicit)描述是通过描述点之间关系来规定点的属性,比如三维空间中满足 x^2+y^2+z^2=1 的点的集合 ;更一般地定义为满足 f(x,y,z)=0 的点的集合

隐式描述的缺点在于不直观,优点在于实现上可以快速判断给定点是否在指定集合中或是否在区域内/外;

显式(Explicit)描述是指所有点都是直接定义或者通过参数映射得到的(参数映射如下图);

image.png

相比于隐式描述,显式描述的优点在于能够快速确定有哪些点在指定集合中,缺点在于难以判断点在区域内/外;

隐式表述

代数曲面(Algebraic Surfaces)

表面是有关 x ,y ,z 的多项式的零点集合;

image.png

CSG(Constructive Solid Geometry)

通过基本几何体的布尔(集合?)运算来定义更复杂的几何体;

image.png

距离函数(有向距离场)(Distance Functions,Signed Distance Field)

通过定义距离函数描述空间上任何一点到几何体表面的最小(有向)距离来构造几何体;

image.png

距离函数可以通过函数运算来获得两个几何体的中值(混合);

image.png

水平集(Level Set Methods)

类似于距离函数,水平集是将距离函数定义域离散化到网格中;

image.png

分型(Fractals)

分型是一种自相似结构;

Lecture 11 Geometry 2 (Curves and Surfaces)(几何2-曲线和曲面的显式描述)

显式描述

点云(Point Cloud)

不考虑表面的连续性,将表面用表面上的足够多点表示;

通常用于大量数据的情况( point/pixel >>1 )

这种方式可以容易地表示任何类型的模型,但在样本密度低的部分难以绘制;

在处理时,点云经常被转换为多边形网格;

多边形面(Polygen Mesh)

易于处理;一般通过三角形或四边形网格实现;

由于需要处理图形间的连接关系,需要更复杂的数据结构;

是图形学中最常见的表现形式;

.obj 文件格式(The Wavefront Object File)

只是一个指定顶点、法线、文本坐标及其连接性的文本文件;‘

以下图为例,其中存储了一个正方体

image.png

其中 line 3-10 定义了8个顶点,line 12-25 定义了14个纹理坐标,line 27-34 定义了6个不同的平面法向量,line36-47定义了12个三角形,其中 4/3/1 含义为此处点坐标为4号,纹理坐标为3号,法线为1号 即为定义了使用第几个前述值;

曲线(Curves)

贝塞尔曲线(Bézier Curves)

用一系列的控制点去定义曲线;

对于 P0,..Pn+1P_0,..P_{n+1} 所定义的 n 阶贝塞尔曲线,需要满足以下性质:

  1. P0P_0 开始、在 Pn+1P_{n+1} 结束;
  2. 在起点处切向为 nP0P1n\overrightarrow{P_0P_1} 、在终点处切向为 nPnPn+1n\overrightarrow{P_{n}P_{n+1}}
  3. 对于曲线的仿射变换等于对控制点仿射变换后绘制曲线;(投影变换不属于仿射变换)
  4. 曲线在控制点形成的凸包内;
De Casteljau算法(de Casteljau Algorithm)

考虑已经给定的一系列控制点,如何绘制对应的贝塞尔曲线;

先将问题简化为三个点,此时生成的曲线为二阶贝塞尔曲线(quadratic Bezier);

image.png

对于曲线上参数为 t (t[0,1]t\in[0,1])的点,定义如上图;其中

b01(t)=(1t)b0+tb1b11(t)=(1t)b1+tb2b02(t)=(1t)b01+tb11b_0^1(t)=(1-t)b_0+tb_1\\ b_1^1(t)=(1-t)b_1+tb_2\\ b_0^2(t)=(1-t)b_0^1+tb_1^1

足够精度在定义域内枚举 t 即可计算曲线;

考虑四个点定义的贝塞尔曲线;

image.png

更广泛地,对于 n+1 个控制点(b0,...,bnb_0,...,b_n),n 阶贝塞尔曲线有代数形式

bn(t)=b0n(t)=j=0nbj(ni)ti(1t)nib^n(t)=b_0^n(t)=\sum_{j=0}^nb_j\begin{pmatrix}n\\i\end{pmatrix}t^i(1-t)^{n-i}

通过定义与参数有关的多项式对控制点进行插值;

上述例子都是二维的,对三维点同样有效;

逐段贝塞尔曲线(Piecewise Bézier Curves)

在贝塞尔曲线阶数较高(控制点较多)的时候,难以精细控制曲线;

考虑对贝塞尔曲线每四个控制点进行逐段定义;

在曲线连接处,保证另外两个控制点关于共用控制点反向、共线且等距即可保证曲线光滑;

曲线的连续性(Continuity)

C0C^0 连续:曲线首尾相接(值连续)(an=b0a_n=b_0);

C1C^1 连续:相接处切线大小相等、方向相同(导数连续)(an=b0=12(an1+b1)a_n=b_0=\frac 12(a_{n-1}+b_1));

C2C^2 连续即为二阶导数连续;

B-样条曲线(Basis Splines)

相比于贝塞尔曲线,B-样条曲线具有局部性;

B样条与非均匀有理B-样条(NURBS)暂不赘述;

曲面(Surfaces)

双立方贝塞尔曲面(Bicubic Bézier Surface Patch)

绘制过程可以认为是在两个方向上先后进行贝塞尔曲线绘制;

image.png

Lecture 12 Geometry 3(几何3-网格处理)

网格处理(Mesh Operation)

网格细分(Mesh Subdivision)

对网格进行加面,使其更光滑;

Loop 细分(Loop Subdivision)

Loop 细分(Loop 为人名,与循环无关)的三角形网格细分过程;

  1. 将三角形数量增多;连接三角形三边中点使三角形一分为四;

    image.png

  2. 调整三角形位置;

    1. 对于新顶点(边中点):

      image.png

    2. 对于旧顶点:

      image.png

Catmull-Clark细分(Catmull-Clark Subdivision)

对于一般情况(不全为三角形)可以使用此种方法;

首先将面定义为四边形面(quad face)和非四边形面(non-quad face),将度(相连边数)不为4的点定义为奇异点(extraordinary vertex);

对于每个面取点、边取中点,将面上点和其边中点相连(面取点策略可以取其重心);

在取点连线后,新奇异点数=旧奇异点数+旧非四边形面数,所有旧非四边形面不再存在;也就是说,只有第一次Catmull-Clark细分会增加奇异点数;

接下来进行对点位置的更新;

  1. 对于面上点:

    image.png

  2. 对于边上点:

    image.png

  3. 对于旧非奇异点:

    image.png

  4. 对于奇异点的位置更新课程没说,网络资料没有具体说明,猜测可以由3.的公式推广;

网格简化(Mesh Simplification)

在损失可接受的范围内对网格进行简化,节省时空负担;

边坍缩(Edge Collapse)

简单来说即为以点代边,将与边相邻的两个三角形网格坍缩掉,与边相连的其他边汇聚到一个顶点上;

坍缩边后点的位置需要通过二次误差度量(Quadric Error Metrics)得到;

二次误差度量即为将点放在与相关平面距离平方和最小的位置上;

image.png

可以选择二次度量误差最小的边来进行坍缩;

坍缩后会影响其他边的二次度量误差,需要更新受影响边的二次度量误差;

网格正规化(Mesh Regularization)

消除掉过于长的三角形,使得三角形趋近于正三角形;

Lecture 13 Ray Tracing 1(光线追踪1-Whitted风格光线追踪)

光线追踪的不同于光栅化的一种成像方式,解决了光栅化对全局效果的问题(软阴影、Glossy 反射(金属反射)、间接光照);光栅化可以认为是一种快速的、近似的成像方式;

光线可以有如下假设:

  1. 在均匀介质中直线传播,不考虑波动性;
  2. 光线间不考虑碰撞;
  3. 光线从光源发出,最终进入相机;也可以认为从相机发出光线,反向传播最终进入光源(光路可逆);

Whitted风格光线追踪(Whitted-Style Ray Tracing)

从相机出发,针对每个像素发出一条感知光线(camera ray / primary ray),直至与其相交的第一个三角面;在这个交点上判断对光源的可见性(阴影)(shadow ray),并考虑在此处发生的反射(和折射);沿反射方向继续发出感知光线(secondary ray),重复此过程;

image.png

在这个过程中,需要考虑求交、能量损失等若干技术问题;

射线-表面交点(Ray-surface Intersection)

光线认为是从起点 O ,向方向 d 传播的射线,有参数表示 r(t)=o+td(0t<)r(t)=o+td,(0\leqslant t<∞)

c 为球心,半径 R 的球有定义 p:(pc)2R2=0p:(p-c)^2-R^2=0

与光线参数表示联立有

(o+tdc)2R2=0(o+td-c)^2-R^2=0

只有 t 为未知量,解得

a=ddb=2(oc)dc=(oc)(oc)R2t=b±b24ac2a,tR,t0a=d\cdot d\\ b=2(o-c)\cdot d\\ c=(o-c)\cdot(o-c)-R^2\\ t={-b \pm \sqrt{b^2-4ac}\over 2a},t\in\R,t\geqslant0

隐式表面(Implicit Surface)

隐式表面有定义 p:f(p)=0p:f(p)=0

光线参数表示代入有 f(o+td)=0f(o+td)=0

解出正实根即可;

平面

平面可以由法线 N 和平面上一点 p’ 定义,定义有 p:(pp)N=0p:(p-p')\cdot N=0

联立有

(o+tdp)N=0t=(po)NdN(o+td-p')\cdot N=0\\ t={(p'-o)\cdot N \over d\cdot N}

解出正实根即可;

轴向平面

轴向平面只需要用与其垂直的坐标轴的坐标值即可定义;

以与 x 轴垂直的平面为例,有

t=pxoxdxt={p'_x-o_x\over d_x}

比一般平面计算速度更快;

三角形网格

光线与三角形的交点还可以判断光源是否在物体内;

先判断光线和三角形所在平面求交,再判断点是否在三角形内;也可以通过重心坐标直接得到(Möller Trumbore Algorithm);

对于方程

O+tD=(1b1b2)P0+b1P1+b2P2O+tD=(1-b_1-b_2)P_0+b_1P_1+b_2P_2

定义

E1=P1P0E2=P2P0S=OP0S1=D×E2S2=S×E1E_1=P_1-P_0\\ E_2=P_2-P_0\\ S=O-P_0\\ S_1=D\times E_2\\ S_2=S\times E_1\\

有解

[tb1b2]=1S1E1[S2E2S1SS2D]\begin{bmatrix} t\\ b_1\\ b_2 \end{bmatrix} = {\frac 1{S_1\cdot E_1}} \begin{bmatrix} S_2\cdot E_2\\ S_1\cdot S\\ S_2\cdot D \end{bmatrix}

(上述大写字母均为三位列向量)

如果 t,b1,b2,1b1b20t,b_1,b_2,1-b_1-b_2\geqslant0 即说明点在三角形内;

加速求交过程

对于光线和一个模型(场景)求交,朴素想法即为枚举面片进行计算,取 t 最小非负值,此种方法消耗很大,考虑加速这一过程;

加速结构的建立是在光线追踪之前;

包围盒(Bounding Volumes)

用简单集合体将模型包围,如果一个光线与包围盒没有交,那么其与包围盒内物体也没有交;

对于轴向对齐包围盒(AABB),可以认为其是三对平面所夹的空间;

考虑光线和包围盒求交过程;当光线与三组对面中的任何一个面均有交点时,认为光线进入包围盒,当光线穿过任何一个对面中的两个平面时,认为光线离开包围盒;

对于每一对面,分别计算 tmin,tmaxt_{min},t_{max} (无论正负),光线关于包围盒的进入点有 tenter=max{tmin}t_{enter}=\max\{t_{min}\} ,离开点有 texit=min{tmax}t_{exit}=\min\{t_{max}\}

如果 texit<0t_{exit}<0 ,我们认为包围盒在光线后,与光线没有交点;

如果 texit0,tenter<0t_{exit}\geqslant0,t_{enter}<0 ,我们认为光源在包围盒中;

总地来说,如果 tenter<texittexit0t_{enter}<t_{exit},t_{exit}\geqslant0 ,我们认为光线与包围盒有交点;

均匀网格(Uniform grids)

对场景求出包围盒后,将包围盒分为若干均等的小格子,对于包含物体表面的格子进行标记;

沿射线轨迹便理网格,对于遍历到的每个格子,测试射线与格子内对象的相交性,如果相交即停止,意味着找到了第一个交点。

image.png

格子划分密度上有经验 #cells=C#objs,C27 in 3D\#cells = C * \#objs, C\approx27\text{ in 3D}

在几何体分布较为均匀的场景优化效果较好;否则容易出现“‘Teapot in a stadium’ problem”;

空间划分(Spatial Partitions)

以 KD-Tree 为例详细说明;

image.png

八叉树(Oct-Tree)

对空间进行轴向划分,将每个节点进行八分割存储在子结点中,划分至子结点中有足够少的物体;

KD-Tree

不同于八叉树,每次只对节点进行一次轴向划分,不对于父子节点进行相同轴向的划分(三维上沿 x ,y ,z 方向上循环),同样直至子结点中有足够少的物体;

只需要在叶子节点存储模型;

在判断时,先判断与根节点是否存在交点,如果有交点则判断与两子节点是否存在交点;直至叶子节点,如果仍有交点则判断其与节点内部模型是否有交点;

KD-Tree 的主要问题是难以判断格点与哪些模型相交;此外,如果一个物体与多个叶子节点相交,其会被存储多份;

BSP-Tree

不同于KD-Tree进行轴向划分,对节点进行更自由的二分;

物体划分(Object Partitions)

BVH(Bounding Volume Hierarchy)树

BVH树对物体进行划分,如图所示,直至划分至子结点中有足够少的三角形;

image.png

BVH树解决了 KD-Tree 的两大问题;

其存在的问题是节点包围盒可能相交,影响有限;实际上会保证重叠部分尽量少;

关于如何划分节点:

  1. 类似 KD-Tree ,可以循环使用不同的轴向;也可以选择包围盒最长的一条轴进行划分;

  2. 对于如何划分两半,可以将三角形的重心按分割轴向排序,取中位三角形进行分隔,使得两部分三角形数量均匀,降低资源消耗;

    寻找中位数不需要排序,可以转化为top-k大数问题,通过快速选择算法在线性复杂度内解决;

在场景变化时,需要重新计算BVH树;

与 KD-Tree 类似,在中间节点只存储包围盒和子节点指针,只在叶子节点存储分割后的模型;

光线求交过程也与 KD-Tree 类似,有伪代码

1
2
3
4
5
6
7
8
9
10
11
12
13
Intersect (Ray ray, BVH node) {
if (ray misses node.bbox) return;

if (node is a leaf node)
{
test intersection with all objs;
return closest intersection;
}

hit1 = Intersect(ray, node.child1);
hit2 = Intersect(ray, node.child2);
return the closer of hit1, hit2;
}

Lecture 14 Ray Tracing 2(光线追踪2-辐射度量学(Radiometry))

辐射度量学是一种精确定义光与物体表面作用的工具,是路径追踪的基础;

辐射度量学同样着眼于光的几何属性,忽略波动性;

辐射度量学通过辐射通量(Radiant flux)、强度(Intensity)、照度(Intensity)、亮度(Radiance)等属性描述光,下述均使用英文;

Radiant Energy

Radiant Energy 是电磁辐射的能量,以焦耳(J)表示;

Q[J=Joule]Q[\text{J=Joule}]

Radiant Flux (Power)

Radiant Flux 是电磁辐射的功率,单位时间的能量,也可以被度量为单位时间内辐射出的光子数量,以瓦特(W)或流明(lm)表示;

Φ=dΦdt[W=Watt][lm=lumen]\Phi =\frac {\mathrm{d}\Phi}{\mathrm{d}t}[\text{W=Watt}][\text{lm=lumen}]

现代11W LED灯约815lm,等价于60W白炽灯;

立体角(Solid Angle)

类比于角度 θ=lr\theta =\frac lr ,圆上有 2π2\pi 角(radians),

image.png

立体角是角度的三维延伸,Ω=Ar2\Omega=\frac A{r^2} ,球上有 4π4\pi 立体角(steradians);

image.png

Radiant Intensity

(Radiant) Intensity 是对光源方向上功率的度量,每单位立体角的功率,以坎德拉(cd)表示;

I(ω)=dΦdω[Wsr][lmsr=cd=candela]I(\omega) =\frac {\mathrm{d}\Phi}{\mathrm{d}\omega}[\frac {\text{W}}{\text{sr}}][\frac {\text{lm}}{\text{sr}}\text{=cd=candela}]

单位立体角可以由天顶角 θ\theta (与 z 轴夹角)和方向角 ϕ\phi (绕 z 轴旋转的角度)定义;

image.png

dA=(rdθ)(rsinθdϕ)=r2sinθ dθdϕdω=dAr2=sinθ dθdϕI=Φ4π\mathrm{d}A=(r\mathrm{d}\theta)(r\sin\theta\mathrm{d}\phi)\\ =r^2\sin\theta\text{ }\mathrm{d}\theta\mathrm{d}\phi\\ \mathrm{d}\omega=\frac{\mathrm{d}A}{r^2}=\sin\theta\text{ }\mathrm{d}\theta\mathrm{d}\phi\\ I=\frac {\Phi}{4\pi}

ω\omega 可以作为方向向量;

Irradiance

Irradiance 是对物体表面单位面积接收到光的功率的度量;

E(x)dΦ(x)dA[Wm2][lmm2=lux]E(x)\equiv\frac{\mathrm{d}\Phi(x)}{\mathrm{d}A}[\frac W{m^2}][\frac{lm}{m^2}=\text{lux}]

要求单位平面与光线垂直;

Radiance

Radiance 是对光线传播中的度量,是每单位立体角单位面积上的功率;

image.png

L(p,ω)d2Φ(p,ω)dωdAcosθ[Wsr m2][cdm2=lmsr m2=nit]L(p,\omega)\equiv\frac{\mathrm{d^2}\Phi(p,\omega)}{\mathrm{d}\omega\mathrm{d}A\cos\theta} [\frac {\text{W}}{\text{sr }\text{m}^2}][\frac {\text{cd}}{\text{m}^2}=\frac {\text{lm}}{\text{sr }\text{m}^2}=\text{nit}]

dAcosθ\mathrm{d}A\cos\theta 为单位面积在传播方向上的有效面积;

联系前面的两个量,有

L(p,ω)=dE(p)dωcosθL(p,ω)=dI(p,ω)dAcosθL(p,\omega)=\frac{\mathrm{d}E(p)}{\mathrm{d}\omega\cos\theta}\\ L(p,\omega)=\frac{\mathrm{d}I(p,\omega)}{\mathrm{d}A\cos\theta}

图形学上常用物理量为 Irradiance 和 Radiance ,Irradiance可以理解为单位面积吸收的总能量,Radiance 为单位面积从 dω\mathrm{d}\omega 方向吸收的总能量,有

dE(p,ω)=Li(p,ω)cosθdωE(p)=H2Li(p,ω)cosθdω\mathrm{d}E(p,\omega)=L_i(p,\omega)\cos\theta\mathrm{d}\omega\\ E(p)=\int_{H^2}L_i(p,\omega)\cos\theta\mathrm{d}\omega

image.png

Lecture 15 Ray Tracing 3(光线追踪3-辐射度量学 续)

双向反射分布函数(Bidirectional Reflectance Distribution Function ,BRDF)

反射可以认为是从方向 ωi\omega_i 来的辐射转换为单位面积 dA\mathrm{d}A 上的功率 EE ,这个功率会被反射到任意其他方向 ωr\omega_r

image.png

BRDF描述了各个其他方向 ωr\omega_r 分配到了多少能量;

fr(ωiωr)=dLr(ωr)dEi(ωi)=dLr(ωr)Li(ωi)cosθidωi[1sr]f_r(\omega_i\to\omega_r)={\mathrm{d}L_r(\omega_r)\over\mathrm{d}E_i(\omega_i)}={\mathrm{d}L_r(\omega_r)\over L_i(\omega_i)\cos\theta_i\mathrm{d}\omega_i}[\frac 1{\text{sr}}]

那么某方向接收到的反射光即为

Lr(p,ωr)=H2fr(p,ωiωr)Li(p,ωi)cosθidωiL_r(p,\omega_r)=\int_{H^2}f_r(p,\omega_i\to\omega_r)L_i(p,\omega_i)\cos\theta_i\mathrm{d}\omega_i

image.png

如果物体自发光,那么就有了渲染方程

Lo(p,ωo)=Le(p,ωo)+Ω+fr(p,ωi,ωo)Li(p,ωi)(nωi)dωiL_o(p,\omega_o)=L_e(p,\omega_o)+\int_{\Omega^+}f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)(n\cdot\omega_i)\mathrm{d}\omega_i

LeL_e 部分为自发光,积分部分为反射光;ωi\omega_i 为指向光来向的向量;p 可以认为是反射点表面的性质;

记 E 为自发光,K为反射操作符,L 可以记为 L=E+KE+K2E+...L=E+KE+K^2E+... ,即为光源直射+一次反射+二次反射+……;

Lecture 16 Ray Tracing 4(光线追踪4-蒙特卡洛积分和路径追踪)

蒙特卡洛积分(Monte Carlo Integration)

计算难以解析求解的函数定积分的一种数值方法;

若求定积分 abf(x)dx\int_a^bf(x)\mathrm{d}x ,定义随机变量 Xip(x)X_i\sim p(x) ,有结果 FN=1Ni=1Nf(Xi)p(Xi)F_N=\frac 1N\sum_{i=1}^N\frac{f(X_i)}{p(X_i)} ;其中概率密度函数(probability density function ,PDF)要求 abp(x)dx=1\int_a^bp(x)\mathrm{d}x=1

如果随机变量均匀分布,则结果有 FN=baNi=1Nf(Xi)F_N=\frac {b-a}N\sum_{i=1}^Nf(X_i) ,此时即为黎曼积分;

N 越大,结果越准确;

在 p(x) 与 f(x) 形状一致时,误差最小;

路径追踪(Path Tracing)

Whitted风格光线追踪总是进行镜面反射 / 折射,停止在漫反射表面,存在难以处理金属反射(glossy reflection)质感、无漫反射等问题;

对于渲染方程,我们需要递归地通过蒙特卡洛解半球上的积分;

我们的概率密度函数可以取 p(ωi)=12πp(\omega_i)=\frac 1 {2\pi} ,有

Lo(p,ωo)Le(p,ωo)+2πNi=1Nfr(p,ωi,ωo)Li(p,ωi)(nωi)L_o(p,\omega_o)\approx L_e(p,\omega_o)+\frac {2\pi}N\sum_{i=1}^N{f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)(n\cdot\omega_i)}

1
2
3
4
5
6
7
8
9
10
shade (p,wo) //忽略自发光
Randomly choose N directions wi~PDF
Lo = 0.0
For each wi
Trace a ray r(p, wi)
If ray r hit the light
Lo += ( 2 * pi / N ) * L_i * f_r * cosine
Else If ray r hit an object at q
Lo += ( 2 * pi / N ) * shade(q, -wi) * f_r * cosine
Return Lo

上面的伪代码复杂度不可接受,若使 N=1 则有

1
2
3
4
5
6
7
8
shade (p,wo) //忽略自发光
Randomly choose 1 direction wi~PDF
Lo = 0.0
Trace a ray r(p, wi)
If ray r hit the light
Return ( 2 * pi ) * L_i * f_r * cosine
Else If ray r hit an object at q
Return ( 2 * pi ) * shade(q, -wi) * f_r * cosine

N!=1 时即为分布式光线追踪(Distributed Ray Tracing),现很少使用;

对于降低 N 造成的偏差,我们可以通过对于一个像素多次进行路径追踪并平均结果得到;

image.png

对于射线生成(Ray Generation),有伪代码

1
2
3
4
5
6
7
8
ray_generation (camPos,pixel)
Uniformly choose N sample positions within the pixel
pixel_radiance = 0.0
For each sample in the pixel
shoot a ray r(camPos, cam_to_sample)
If ray r hit the scene at p
pixel_radiance += 1 / N * shade(p, sample_to_cam)
Return pixel_radiance

上述 shade 代码依然存在递归终止条件的问题,我们选择轮盘赌策略(Russian Roulette,RR),在满足一定条件时停止递归;

如果有一个概率 P ,我们有 P 的概率发出光线,对于得到的 LoL_oLo/PL_o/P ;有 1-P 的概率不发出光线,得到的结果为 0 ;得到的结果期望依然为 LoL_o

1
2
3
4
5
6
7
8
9
10
11
12
shade (p,wo) //忽略自发光
Manually specify a probability P_RR
Randomly select ksi in a uniform dist. in [0, 1]
If (ksi > P_RR) return 0.0

Randomly choose 1 direction wi~PDF
Lo = 0.0
Trace a ray r(p, wi)
If ray r hit the light
Return ( 2 * pi ) * L_i * f_r * cosine / P_RR
Else If ray r hit an object at q
Return ( 2 * pi ) * shade(q, -wi) * f_r * cosine / P_RR

也可以通过控制 P_RR 来控制期望弹射次数(递归深度)为 1/(1-P_RR) ;

此时即为路径追踪;

在低 SPP(samples per pixel)时,生成速度更快但会产生噪点;下面考虑加速;

如果光源很小,那么采样线打到光源的概率很低,浪费了很多的计算能力;

我们考虑对光源采样,将积分从 dω\mathrm{d}\omega 变换为 dA\mathrm{d}A

image.png

有关系

dω=dAcosθxx2\mathrm{d}\omega =\frac{\mathrm{d}A\cos\theta'}{||x'-x||^2}

渲染方程换元有

Lo(p,ωo)=Le(p,ωo)+Afr(p,ωi,ωo)Li(p,ωi)cosθcosθxp2dAL_o(p,\omega_o)=L_e(p,\omega_o)+\int_{A}f_r(p,\omega_i,\omega_o)L_i(p,\omega_i)\frac{\cos\theta\cos\theta'}{||x'-p||^2}\mathrm{d}A

考虑着色点的贡献来自于所有的光源和非光源,来自于光源的部分不需要 RR 策略来加速,来自于非光源的部分采用前述方法;

有伪代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
shade(p, wo)
# Contribution from the light source.
Uniformly sample the light at x' ( pdf_light = 1 / A)
Shoot a ray from p to x'
If the ray is not blocked in the middle
L_dir = L_i * f_r * cos theta * cos theta' / |x' - p|^2 / pdf_light

# Contribution from other reflectors.
L_indir = 0.0
Test Russian Roulette with probability P_RR
Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
Trace a ray r(p, wi)
If ray r hit a non-emitting object at q
L_indir = shade(q, -wi) * f_r * cos theta / pdf_hemi / P_RR

return L_dir + L_indir

对于点光源的处理比较困难,课程内暂不赘述;

事实上计算出的 radiance 并不是颜色,需要进行伽马矫正(gamma correction);

Lecture 17 Materials and Appearances(材质与外观)

对材质的定义等价于对 BRDF 的定义;

漫反射

对于漫反射模型,假设入射光被均匀地反射到半球上的各个方向(frf_r 为常数),渲染方程反射项有

Lo(ωo)=H2frLi(ωi)cosθidωi=frLiH2cosθidωi=πfrLiL_o(\omega_o)=\int_{H^2}f_rL_i(\omega_i)\cos\theta_i\mathrm{d}\omega_i\\ =f_rL_i\int_{H^2}\cos\theta_i\mathrm{d}\omega_i\\ =\pi f_rL_i

在能量守恒的假设下,有 fr=1πf_r=\frac 1 \pi ,更广泛地,如果我们定义反照率(albedo) ρ[0,1]\rho\in[0,1] ,即有 fr=ρπf_r=\frac \rho\pi

金属反射

对于不同金属,金属反射性质也不同;

image.png

镜面反射

由出射角等于入射角,有

ωo+ωi=2cosθn=2(ωin)nwo=wi+2(ωin)n\omega_o+\omega_i=2\cos\theta\vec n=2(\omega_i\cdot\vec n)\vec n\\ w_o=-w_i+2(\omega_i\cdot\vec n )\vec n

image.png

对于方位角,有

ϕo=(ϕi+π)mod2π\phi_o=(\phi_i+\pi)\mod 2\pi

image.png

镜面反射的 BRDF 值与 delta 函数有关,课程内暂不赘述;

折射

由折射定律(Snell’s Law)得

ηisinθi=ηtsinθt\eta_i\sin\theta_i=\eta_t\sin\theta_t

其中 ηi,ηt\eta_i,\eta_t 分别为入射介质和出射介质的折射率(真空为1);

image.png

cosθt=1(ηiηt)2(1cos2θi)\cos\theta_t=\sqrt{1-(\frac {\eta_i} {\eta_t})^2(1-\cos^2\theta_i)}

如果 1(ηiηt)2(1cos2θi)<01-(\frac {\eta_i} {\eta_t})^2(1-\cos^2\theta_i)<0 则折射不会发生,此时发生全反射;

image.png

其中可视锥体被称为 Snell’s Window ;

方位角同镜面反射;

描述折射性质的函数被称为 BTDF ,与 BRDF 统称为 BSDF ;

菲涅尔项(Fresnel Term)

菲涅尔项描述了反射和折射的比例;

image.png

对于绝缘体表面,有

image.png

S 和 P 描述了两种偏振光的表现,一般情况性质如红线描述;

对于金属(导体)表面,则有

image.png

对于两种偏振光的反射率有(n1,n2n_1,n_2 为两种介质的折射率)

image.png

对于无偏振的光,有

Reff=12(Rs+Rp)R_{eff}=\frac 12 (R_s+R_p)

由于上述计算比较复杂,有简化公式 Schlick’s approximation :

基准反射率即为法线入射时的反射率,为

R0=(n1n2n1+n2)2R_0=(\frac{n_1-n_2}{n_1+n_2})^2

关于入射角的反射率为

R(θ)=R0+(1R0)(1cosθ)5R(\theta)=R_0+(1-R_0)(1-\cos\theta)^5

对导体和绝缘体均可以近似;

微表面模型(Microfacet Material)

微表面模型假设物体表面是粗糙的,由很多微表面组成,每一个微表面只进行镜面反射,具有自己的法线;远看是外观,近看是几何;

通过描述微表面法线的分布来描述物体粗糙程度,例如

image.png

微表面模型的 BRDF 有

f(i,o)=F(i,h)D(h)G(i,o,h)4(n,i)(n,o)f(i,o)={F(i,h)D(h)G(i,o,h)\over4(n,i)(n,o)}

image.png

其中 h 为 i 和 o 的中间向量,F 为菲涅尔项,D 为微表面法向的概率密度函数,G 为几何项,定义了了微表面互相遮挡造成的影响;

image.png

在入射方向或观察方向和表面几乎平行时(Grazing Angle 掠射角度),表面遮挡比较严重,此时 G 项的值会比较小;

GGX模型下,对于粗糙度为 α\alpha 的各向同性材料,法线分布函数和遮挡通过函数有

D(n,h,α)=α2π((nh)2(α21)+1)2k=(α+1)28G(i,o,n)=ni(ni)(1k)+kno(no)(1k)+kD(n,h,\alpha)={\alpha^2\over\pi((n\cdot h)^2(\alpha^2-1)+1)^2}\\ k=\frac{(\alpha+1)^2}8\\ G(i,o,n)={n\cdot i\over(n\cdot i)(1-k)+k}{n\cdot o\over(n\cdot o)(1-k)+k}

各向同性/各向异性材料(Isotropic / Anisotropic Materials)

类似于电梯内表面的纵向抛光痕迹,一些材料表面是各向异性的;

image.png

各向同性即认为微表面不存在方向性或方向性很弱;

image.png

从 BRDF 考虑,高度角固定的情况下,如果对于固定方位角变化量,函数值不同,即为各向异性;

BRDF

BRDF 的性质

  1. 非负:fr(ωiωr)0f_r(\omega_i\to\omega_r)\geqslant0
  2. 线性;
  3. 可逆性:fr(ωiωr)=fr(ωrωi)f_r(\omega_i\to\omega_r)=f_r(\omega_r\to\omega_i)
  4. 能量守恒: ωr,H2fr(ωiωr)cosθidωi1\forall \omega_r,\int_{H^2}f_r(\omega_i\to\omega_r)\cos\theta_i\mathrm{d}\omega_i\leqslant1
  5. 各向同性与各向异性;

测量 BRDF

gonioreflectometer

image.png

BRDF 的存储

MERL BRDF Database 对各向同性材料进行了三个维度的测量 (θi,θo,ϕiϕo)(\theta_i,\theta_o,|\phi_i-\phi_o|) ,形成了 909018090*90*180 的数据;

Lecture 18 Advanced Topics in Rendering(渲染前沿)

无偏 / 有偏光线追踪

基于无偏蒙特卡洛算法(无论采用多少样本,期望均等于真实值)的光线追踪即为无偏光线追踪;

对于有偏蒙特卡洛算法,存在一种特殊情况,即为在样本趋近于无穷多时为无偏的,被称为一致的(consistent);

双向路径追踪(Bidirectional Path Tracing ,BDPT )(无偏)

BDPT 利用光向可逆性,从光源和相机分别发出一条半路径,两端点重合时即可成为一条路径;

在相同spp时,双向路径追踪效果更好,能够更有效地找到能量集中的路径;在光线传播集中在光源附近时表现更优;

Metropolis 光线传播(Metropolis Light Transport ,MLT)(无偏)

使用马尔可夫链(Markov Chain Monte Carlo ,MCMC)进行采样,可以基于一个样本生成相近的其他样本;

MLT 对于复杂光路(例如 SDS(Specular-Diffuse-Specular))效果更优,可以基于一条有效光路更快地找到有效的光路;

但是难以估计收敛速度,不能保证每个像素收敛速度相同(会形成不可控的脏点),因此不适合渲染动画;

光子映射(Photon Mapping)(有偏,一致)

适合渲染 caustics 现象和复杂光路(例如SDS);

image.png

光子映射的一种实现方法:

首先从光源出发模拟光子的运动,直至到达漫反射面;再从相机发出半路径,直至打到漫反射面;对相机出发的半路径终点进行局部密度估计(local density estimation)来计算像素的亮度;

局部密度估计即为对于某个点找周围的最近 n 个光子,结合这 n 个光子所占面积计算局部光子密度,进而计算亮度;

对于 n 更小的取值,结果噪声更多;对于 n 更大的取值,结果更模糊;

有偏体现在

dN/dAΔN/ΔA\mathrm{d}N/\mathrm{d}A\not=\Delta N/\Delta A

VCM(Vertex Connection and Merging)

结合了双向路径追踪和光子映射的过程;

假设双向路径终点再同一个表面上但是并不重合,就让相机端半路径作为光子进行光子映射;

实时辐射度(Instant Radiosity ,IR)

从光源打出半路径,将已经被照亮的表面当作光源(虚拟点光源,Virtual Point Light ,VPL)处理;

渲染时仅通过直接光照即可;

在接缝处会出现“漏光”,无法处理金属反射;

外观建模(Appearance Modeling)

非表面模型(Non-surface Models)

散射介质(Participating)

光线在穿过散射介质时会发生散射;

散射可以由相位函数(Phase Function)来定义

image.png

在散射云中,每次从上一个着色点随机选择一个方向和一段距离确定下一个着色点,在每一个着色点连接到光源进行着色;

不止云、雾和霾,很多材质都可以视为散射介质,比如说手指和胶体;

头发材质(Hair Appearance)

Marschner Model 将头发视为玻璃圆柱,表层和含颜色的填充层;

光线与单根头发的作用分为三个部分,反射(R),穿透(TT),内壁反射(TRT)

image.png

T 为 transmission(传递),R 为 reflection ;

通过多次散射达到对于头发表面的渲染;

基于毛发在解剖角度上存在髓质,有双圆柱模型,在上述模型的基础上添加了负责反射的髓质芯;此时光线和毛发的作用在上述基础上还增加了“发生散射并穿透(TTs)”和“发生散射并内部反射(TRTs)”;

image.png

颗粒材质(Granular Material)

如盐、糖、沙子;

通过对单元内组分的定义来进行渲染;

image.png

过于耗时,没有比较好的解决方案;

表面模型(Surface Models)

透射材质(Translucent Material)

光从一点进入物体,在物体内部发生次表面散射(Subsurface Scattering)后离开物体;比如手指、玉石、胶体;

次表面散射可以认为是对 BRDF 的一种延伸,BSSRDF 除了接受入射方向和出射方向,也需要接受入射位置和出射位置;

对于 S(xi,ωi,xo,ωo)S(x_i,\omega_i,x_o,\omega_o) ,有模型

Lo(xo,ωo)=AH2S(xi,ωi,xo,ωo)Li(xiωi)cosθidωidAL_o(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

image.png

模型可以简化为在真实光源相对于着色点的对称点添加一个虚拟光源,被称为 Dipole Approximation,暂不赘述;

布料(Cloth)

纤维(fibers)缠成股(ply),股缠绕成线(yarn),线织成布;

可以将织物视作空间中分布的体积,将其转化成对散射介质的渲染,这种方法速度较快;

也可以将纤维按毛发逻辑渲染,但是速度较慢;

程序化生成(Procedural Appearance)

对于三维物体内部纹理,如果直接存储会大量占用存储资源;我们选择对其进行程序化生成;通过一些 3D Noise 函数对于给定坐标返回纹理;

Lecture 19 Cameras, Lenses and Light Fields(相机、镜头和光场)

图形学主要有两种成像方式,即光栅化成像和光线追踪成像;现实成像是通过捕捉的方式进行的;

视场(Field of View ,FOV)

以垂直 FOV 为例

image.png

FOV=2arctan(h2f)FOV=2\arctan(\frac h{2f})

由于历史原因,以 36×24mm36\times24\text{mm} 为基准画幅,其对角线长为 h ,结合镜头的焦距来定义视场;

image.png

等效焦距的概念即为画幅下与标准画幅视场角相同的焦距;

曝光(Exposure)

Exposure=time×irradiance\text{Exposure}=\text{time}\times\text{irradiance}

time 为曝光时间,由快门控制;irradiance 由透镜孔径和焦距控制;

快门、光圈、ISO略;

光圈 f 数定义为焦距与透光直径的比值;

薄透镜近似(Thin Lens Approximation)

目前镜头都是由透镜组组成,课程内分析的是理想化的薄透镜,忽略其厚度;

image.png

特别地,假设薄透镜可以自由地改变焦距;

透镜的象距(ziz_i)和物距(zoz_o)有关系

1f=1zi+1zo\frac 1f=\frac1{z_i}+\frac 1{z_o}

image.png

可由相似三角形证明;

弥散圆(Circle of Confusion ,CoC)

景深外物体的一点会在感光元件上形成弥散圆(光晕),假设光圈直径为 A ,f 数为 N ,有

CA=zsziziC=fNzszizi\frac CA=\frac{|z_s-z_i|}{z_i}\\ C=\frac fN\frac{|z_s-z_i|}{z_i}

image.png

可以观察到,光圈直径越小(f数越大)弥散圆越小;

薄透镜射线追踪

定义传感器尺寸、光圈大小、焦平面位置、透镜焦距(ff)、透镜与焦平面的距离(zoz_o),后两者可以推导出透镜与传感器的距离(ziz_i);

选择传感器的一个像素中心和透镜上的点进行连线,根据透镜公式,可以得出这束光离开棱镜时的方向和位置,按光线追踪进一步处理即可;

image.png

景深(Depth of Field ,DOF)

CoC 足够小(小于等于像素大小)时我们便认为该物体是清晰的,物体清晰成像的范围便是景深;

image.png

DOF=DFDNDF=Dsf2f2NC(DSf)DN=Dsf2f2+NC(DSf)DOF=D_F-D_N\\ D_F={D_sf^2\over f^2-NC(D_S-f)}\\ D_N={D_sf^2\over f^2+NC(D_S-f)}

光场(Light Field / Lumigraph)

定义全光函数(The Plenoptic Function) P=(θ,ϕ,λ,t,VX,VY,VZ)P=(\theta,\phi,\lambda,t,V_X,V_Y,V_Z) 为 t 时刻在位置 VX,VY,VZV_X,V_Y,V_Z(θ,ϕ)(\theta,\phi) 方向看到的波长为 λ\lambda 的颜色亮度值;将现实世界定义为这七个维度的函数;

光场是对全光函数的一种简化:对于一个物体,在物体 / 包围盒表面上定义在任意位置向任意方向(外半球)的发光情况(4D);

光线除了点向定义,也可以由有序点对定义,所以光场也可以由两个平行的平面进行定义,光线便可以由两平面上的两点 (u,v),(s,t)(u,v),(s,t) 来定义;

昆虫的复眼实际上进行的就是光场成像,将一个像素替换为一个透镜,使我们能够将打在这个像素上不同方向上的光记录在不同位置;

光场照相机(Light Field Camera)

利用光场成像,我们可以在后期进行重新聚焦或移动相机的位置;

image.png

在每个透镜下存储的光线分别选择对应的光线即可进行重新成像;

光场相机存在分辨率不足、高成本等问题;

Lecture 20 Color and Perception(颜色和感知)

谱功率密度(Spectral Power Distribution ,SPD)

对于混合光,可以使用谱功率密度来描述不同波长光的分布;

SPD具有线性的性质,记光 A ,B 的谱功率密度分别为 SPDa ,SPDb ,将 A ,B 光混合得到的谱功率密度即为 SPDa + SPDb ;

颜色感知

颜色与 SPD 不完全相同,是人的感知;视网膜上主要有视锥细胞和视杆细胞两种感光细胞,其中视锥细胞主要感受光的颜色,视杆细胞主要感受光的强度;

视锥细胞分为 S ,M ,L 三类,分别对三种波长的光响应最强;不同的人这三种视锥细胞的分布偏差很大;

image.png

对于给定 SPD 的色光 s(λ)s(\lambda) ,三种视锥细胞感知的结果 S ,M ,L ,有

S=rS(λ)s(λ)dλM=rM(λ)s(λ)dλL=rL(λ)s(λ)dλS=\int r_S(\lambda)s(\lambda)\mathrm{d}\lambda\\ M=\int r_M(\lambda)s(\lambda)\mathrm{d}\lambda\\ L=\int r_L(\lambda)s(\lambda)\mathrm{d}\lambda

同色异谱(Metamerism)

对于两个不同 SPD 的色光,在我们的感知下可能相同;

颜色匹配过程正是利用率同色异谱现象;

颜色匹配(Color Reproduction / Matching)

加色系统是给定一组有自己频谱分布的主光(如 sR(λ),sG(λ),sB(λ)s_R(\lambda),s_G(\lambda),s_B(\lambda) ),调整这些主光的亮度并相加(RsR(λ)+GsG(λ)+BsB(λ)Rs_R(\lambda)+Gs_G(\lambda)+Bs_B(\lambda)),这个色光便可由 R ,G ,B 三个值描述;

image.png

认为左右颜色相同时即完成了颜色匹配;

CIE RGB 空间

CIE规定了三种波长作为 CIE RGB 的主光

image.png

对于任意波长的单色光,有匹配函数

image.png

对于给定 SPD 的色光 s(λ)s(\lambda) ,CIE RGB 主光亮度有

RCIE RGB=λs(λ)r(λ)dλGCIE RGB=λs(λ)g(λ)dλBCIE RGB=λs(λ)b(λ)dλR_{\text{CIE RGB}} = \int_\lambda s(\lambda)\overline r(\lambda)\mathrm d\lambda\\ G_{\text{CIE RGB}} = \int_\lambda s(\lambda)\overline g(\lambda)\mathrm d\lambda\\ B_{\text{CIE RGB}} = \int_\lambda s(\lambda)\overline b(\lambda)\mathrm d\lambda\\

此外,还有 CIE XYZ 颜色匹配系统,首先定义色彩匹配函数

image.png

我们观察到其中 y 在不同波长均有分布且对称,在一定程度上可以代表颜色的 luminance ;

我们将 XYZ 归一化,有

x=XX+Y+Zy=YX+Y+Zz=ZX+Y+Zx=\frac X{X+Y+Z}\\ y=\frac Y{X+Y+Z}\\ z=\frac Z{X+Y+Z}

我们选择 x ,y ,对应位置显示特定的 Y 值下的颜色,便形成了颜色玛蹄图

image.png

(13,13)(\frac 13,\frac 13) 处为色域的中心,是白色;

HSV 空间(Hue-Saturation-Value,色相-饱和度-明度)

image.png

CIE LAB 空间

image.png

L 轴为亮度轴,a 轴为红轴, b轴为黄轴;

其中红轴上的红绿两端互为互补色,黄蓝两端互为互补色;

CMYK-减色系统(Cyan-Magenta-Yellow-Key,蓝绿色-品红色-黄色-黑色)

印刷中常用的一种颜色空间,引入黑色是为了控制印刷成本;

Lecture 21 Animation(动画与模拟)

质点弹簧系统(Mass Spring System)

最基础的单元即是由弹簧(劲度系数 ksk_s ,长度 l)连接的两个质点 a ,b ,有

fab=ksbaba(bal)fba=fabf_{a\to b}=k_s\frac{b-a}{||b-a||}(||b-a||-l)\\ f_{b\to a} = -f_{a\to b}

上述模型会出现无限振动的问题,动能和势能永远相互转化,需要添加摩擦力;

x˙\dot x 为 x 的一阶导,x¨\ddot x 为 x 的二阶导;阻力项有

f=kdbaba(b˙a˙)babaf=-k_d\frac{b-a}{||b-a||}(\dot b-\dot a)\frac{b-a}{||b-a||}

其中 baba(b˙a˙)\frac{b-a}{||b-a||}(\dot b-\dot a) 规定了以 a 为参考系下 b 的速度在 ab\overline{ab} 方向上的投影(标量),baba\frac{b-a}{||b-a||} 规定了 ab\overrightarrow{ab} 方向的单位向量;

通过此种形式的模型即可对布料进行一定程度的仿真

image.png

其中蓝线比黑线更强,用于抵消斜向形变,红线比黑线更弱,用于抵消弯折形变;

除了质点弹簧系统,也可以使用有限元方法(Finite Element Method ,FEM)进行模拟;

粒子系统(Partice Systems)

粒子系统被用于模拟灰尘、雾、流体、鸟群等;

粒子受力比较复杂,可能会受到引力、电磁力、虚拟弹簧弹力、推进力、粒子间摩擦力。粒子间粘滞力、粒子间碰撞、空气阻力、与容器或固定物的碰撞、与动态物体的碰撞等等;

正运动学(Forward Kinematics)

通过定义不同的关节来模拟人体的结构;

有三种不同的关节:销(Pin)(一自由度旋转)、球(Ball)(二自由度旋转)、移动关节(Prismatic joint)(平移);

对于给定关节组的参数(类型、长度、旋转角度等),可以很轻松地算出关节末端的位置;

对于两个销关节,有参数 l1,θ1,l2,θ2l_1,\theta_1,l_2,\theta_2 ,可以计算关节末端的位置有

image.png

逆运动学(Inverse Kinematics)

通过定义关节组末端的位置和关节组的类型、长度等参数,反向确定关节组的旋转角度等其他自由度;

对于长度为 l1,l2l_1,l_2 的两个销关节,末端位置 (px,pz)(p_x,p_z) ,有角度参数

image.png

会涉及到多解和无解的问题;

Rigging

Rigging 是逆运动学的一个应用,类似操控布偶一样控制关节;对人物模型添加若干控制点,通过对控制点的控制来控制人物模型;

image.png

动作捕捉(Motion Capture)

通过在真人身上定义并捕捉控制点,将真人的动作映射到虚拟人物上;

真实感强、获取快速;但昂贵且设备复杂,真人演员可能无法满足动画需要;

其中光学动补应用较为广泛;

动画电影产品管线

image.png

Lecture 22 Animation Cont.(动画与模拟 续)

假设物体的运动由速度场决定,速度场可以定义为关于位置和时间的函数 v(x,t)v(x,t)

对于物体有

dxt=x˙=v(x,t){\mathrm dx\over\mathrm t}=\dot x=v(x,t)

这种写法被称为一阶常微分方程;

显式欧拉方法(Explicit Euler)

Δt\Delta t 时刻后的位置与速度可以由此时刻的位置和速度推导

xt+Δt=xt+Δtx˙tx˙t+Δt=x˙t+Δtx¨tx^{t+\Delta t}=x^t+\Delta t\dot x^t\\ \dot x^{t+\Delta t}=\dot x^t+\Delta t\ddot x^t

这种方法并不准确,不准确程度会随着 Δt\Delta t 的增加而增加,并且会导致模拟发散(不稳定性);

image.png

解决不稳定性的问题可以通过以下方法;

中点法(Midpoint Method)

用某个 Δt\Delta t 计算出一个 a 点,用 a 点的参数从起始点进行显式欧拉方法;

xmid=x(t)+Δt/2v(x(t),t)x(t+Δt)=x(t)+Δtv(xmid,t)x_{mid} = x(t)+\Delta t/2\cdot v(x(t),t)\\ x(t+\Delta t) = x(t)+\Delta t\cdot v(x_{mid},t)

自适应步长(Adaptive Step Size)

比较显示欧拉方法计算出的 x(t+Δt)x(t+\Delta t)x(t+Δt/2)(t+Δt/2+Δt/2)x(t+\Delta t/2)(t+\Delta t/2+\Delta t/2) ,如果误差可以接受即取后者结果,否则使 Δt=Δt/2\Delta t=\Delta t/2 递归处理;

image.png

隐式欧拉方法(Implicit Euler Method)

不同于显式欧拉方法,隐式欧拉方法使用下个时刻的参数进行模拟;

xt+Δt=xt+Δtx˙t+Δtx˙t+Δt=x˙t+Δtx¨t+Δtx^{t+\Delta t}=x^t+\Delta t\dot x^{t+\Delta t}\\ \dot x^{t+\Delta t}=\dot x^t+\Delta t\ddot x^{t+\Delta t}

实际上需要解出下一个时刻的参数;

Runge-Kutta Families —— RK4

Runge-Kutta 方法中有一个四阶稳定性的算法名为 RK4 ;

dydt=f(t,y)yn+1=yn+16+Δt(k1+2k2+2k+3+k4)tn+1=tn+Δtk1=f(tn,yn)k2=f(tn+Δt2,yn+Δtk12)k3=f(tn+Δt2,yn+Δtk22)k4=f(tn+Δt,yn+Δtk3){\mathrm dy\over\mathrm dt}=f(t,y)\\ y_{n+1}=y_n+\frac 16 + \Delta t(k_1+2k_2+2k+3+k_4)\\ t_{n+1}=t_n+\Delta t\\ k_1=f(t_n,y_n)\\ k_2=f(t_n+\frac {\Delta t}2,y_n+\Delta t\frac {k_1}2)\\ k_3=f(t_n+\frac {\Delta t}2,y_n+\Delta t\frac {k_2}2)\\ k_4=f(t_n+\Delta t,y_n+\Delta tk_3)

Position-Based / Verlet Intergration

此类方法并不是基于物理的,速度更快,不保证能量守恒;

稳定性(Stability)

对误差可以定义为每步的截断误差(local truncation error)和累积误差(total accumulated error);

通过描述误差与 Δt\Delta t 的关系来评价算法的稳定性;

其中隐式欧拉方法是一阶的,记 Δt=h\Delta t=h ,有截断误差为 O(h2)O(h^2) ,累积误差为 O(h)O(h)

刚体模拟(Rigid Body Simulation)

刚体不会发生形变;

刚体模拟在质点模拟的基础上,增加了角速度和角加速度

ddt(XθX˙ω)=(X˙ωF/MΓ/I){\mathrm d\over\mathrm dt} \begin{pmatrix} X\\ \theta\\ \dot X\\ \omega \end{pmatrix} = \begin{pmatrix} \dot X\\ \omega\\ F/M\\ \Gamma/I \end{pmatrix}

即位置、角度、速度、角速度对时间的导数分别为速度、角速度、加速度、角加速度(力矩/转动惯量);

流体模拟(Fluid Simulation)

认为水体是由不可压缩的刚体小球组成的;对于密度不正确的区域,通过移动小球来修正;为了计算密度,需要得知某点的密度关于所有粒子位置变化的导数;通过梯度下降来修正密度;

欧拉方法与拉格朗日方法(Euierian vs.Lagrangian)

此处欧拉方法与前述解拉格朗日方程的欧拉方法不同;

拉格朗日方法是“质点法”,通过对所有质点的正确模拟来达到对全局的模拟;

欧拉方法是“网格法”,针对于模拟空间网格本身的属性随时间变化的结果;

物质点方法(Material Point Method ,MPM)结合了这两种方法;


GAMES101(现代计算机图形学入门)笔记
https://tanyuu.github.io/2024.01-07/GAMES101笔记/
作者
F Juny
发布于
2024年2月1日
许可协议