OpenGL 2.0 图形编程自学教程 00:导读——OpenGL 是什么

0 阅读18分钟

系列导航:本篇是第一篇 | 下一篇:01 环境搭建与第一个窗口

本篇目标:读完本篇,你将能够——

  • 说清楚计算机图形学研究什么,以及“渲染”在数学上要计算的是什么;
  • 理解光栅化和光线追踪这两种渲染思路的区别,并证明 Z-buffer 算法的正确性;
  • 证明“光栅化 + 透视校正插值”和“光线投射”算出的是同一张图;
  • 说出 OpenGL 的本质:它是一份规范、一个状态机,采用客户端-服务端模型异步执行;
  • 了解 OpenGL 从固定管线到可编程管线的演进,明白 2.0 版本的地位;
  • 知道什么是三维点云、它从哪里来、和三角形网格有什么不同,并亲手“扫描”出一份点云;
  • 知道用 Python 写 OpenGL 要注意什么,以及本系列 36 篇的学习路线。

前置知识:会写 Python,用过一点 NumPy。数学只要求高中水平(知道向量是什么更好)。线性代数在第 08、09 篇从头讲。

你在屏幕上看到的每一帧游戏画面、每一个三维建模软件的视图、每一张网页上的 3D 地图,背后都在回答同一个问题:给定一个用数字描述的三维场景,屏幕上每个像素应该是什么颜色?研究这个问题的学科叫计算机图形学(computer graphics),把场景变成图像的过程叫渲染(rendering)。

OpenGL 是最经典的一套渲染 API。它诞生于 1992 年,至今仍运行在几乎所有桌面电脑上,它的嵌入式版本 OpenGL ES 和网页版本 WebGL 运行在几十亿部手机和浏览器里。很多图形学课程和教材都用它做教学语言。

本系列用 Python 学 OpenGL 2.0。不过要先说清楚一点:OpenGL 只是一套函数接口,它告诉你“怎么调用”,不告诉你“为什么”。为什么要用 4×4 矩阵表示三维变换?为什么深度缓冲里存的值不是线性的?为什么光照公式长那个样子?这些答案都在图形学里。所以本系列以图形学原理为主线,每一篇都讲三件事:原理、数学推导、OpenGL 实现。

本篇是导读,还不写 OpenGL 代码,而是先回答几个根本问题:渲染到底在算什么?OpenGL 是什么样的东西?为什么学 2.0?然后给出整个系列的学习地图。

一、计算机图形学研究什么

1. 图形学的几个分支

按照经典教材《Fundamentals of Computer Graphics》(虎书)第 1 章的划分,图形学的主要领域有:

领域研究的问题例子
建模(modeling)怎样用数字描述形状和外观三角形网格、点云、曲面、材质参数
渲染(rendering)怎样把场景变成图像本系列的主题
动画(animation)怎样让物体随时间运动骨骼动画、物理模拟
交互(interaction)怎样让人操作虚拟场景相机控制、拾取物体
可视化、图像处理、虚拟现实……各种应用方向科学数据可视化、VR 头显

本系列关注实时渲染:每秒生成几十帧图像,让画面能随用户的操作立即变化。OpenGL 就是为它设计的。建模和交互会在用到时顺带讲(第 21 篇的模型加载、第 11 篇的相机控制、第 27 篇的拾取)。此外,第六部分专门讲三维点云的获取、存储和渲染,它是三维扫描、自动驾驶、机器人等领域最基本的数据形式。

2. 渲染在数学上要算什么

先把问题说精确。一个最简单的场景包括:

  • 几何:NN 个三角形 T1,…,TNT_1, \dots, T_N。每个三角形有 3 个顶点,每个顶点有位置,还可以带颜色、法线、纹理坐标等属性(attribute);
  • 相机:眼睛的位置、看的方向、视野的大小;
  • 图像:W×HW \times H 个像素组成的网格。

对每个像素,我们在它的中心取一个采样点 ss。从眼睛出发、穿过 ss 的射线叫视线(viewing ray)。沿着视线看出去,最先挡住视线的那个三角形,就是这个像素“看到”的东西。用公式写:

k∗(s)=arg⁡min⁡k : 视线(s) 与 Tk 相交  dk(s)k^*(s) = \mathop{\arg\min}_{k \,:\, \text{视线}(s) \text{ 与 } T_k \text{ 相交}} \; d_k(s)

其中 dk(s)d_k(s) 是视线与三角形 TkT_k 交点的深度(离眼睛的远近)。有了 k∗(s)k^*(s),像素颜色就是三角形 Tk∗T_{k^*} 在交点处的颜色;视线谁都没碰到时,就用背景色。

这就是渲染的核心问题之一,叫可见性问题(visibility problem)或隐藏面消除(hidden surface removal)。另一个核心问题是着色(shading):交点处的颜色怎么算,需要考虑光照、材质、纹理。第三、四部分专门讲着色。本篇先只讨论可见性。

3. 为什么是三角形

实时渲染几乎只画三角形。四边形、多边形、曲面,最后都会被切成三角形再交给 GPU。原因是三角形有几个其他形状没有的好性质:

  1. 三点必共面。三角形一定是平的,而 4 个点可能不在同一个平面上;
  2. 三角形一定是凸的,判断“点在不在里面”很简单;
  3. 任何多边形都能剖分成三角形,任何光滑曲面都能用足够多的三角形逼近;
  4. 三角形内部的每一点都有唯一的“重心坐标”,可以用它把顶点属性平滑地插值到内部。

第 4 条是后面所有插值(颜色、纹理坐标、法线)的基础,这里给出证明。

定理 1(重心坐标的存在唯一性):设平面上三点 a,b,c\mathbf{a}, \mathbf{b}, \mathbf{c} 不共线。则平面上任意一点 p\mathbf{p} 都可以唯一地写成 p=α a+β b+γ c,α+β+γ=1.\mathbf{p} = \alpha\,\mathbf{a} + \beta\,\mathbf{b} + \gamma\,\mathbf{c}, \qquad \alpha + \beta + \gamma = 1. 并且 p\mathbf{p} 在三角形 abc\mathbf{abc} 内部(含边界)当且仅当 α,β,γ≥0\alpha, \beta, \gamma \ge 0。

证明:把 α=1−β−γ\alpha = 1 - \beta - \gamma 代入,条件等价于

p−a=β (b−a)+γ (c−a).\mathbf{p} - \mathbf{a} = \beta\,(\mathbf{b} - \mathbf{a}) + \gamma\,(\mathbf{c} - \mathbf{a}).

三点不共线,说明向量 b−a\mathbf{b} - \mathbf{a} 和 c−a\mathbf{c} - \mathbf{a} 线性无关。平面是二维的,两个线性无关的向量构成它的一组基,所以 p−a\mathbf{p} - \mathbf{a} 在这组基下的坐标 (β,γ)(\beta, \gamma) 存在而且唯一,α=1−β−γ\alpha = 1 - \beta - \gamma 也随之唯一确定。

再证第二部分。三角形就是三个顶点的凸包,即所有形如 αa+βb+γc\alpha\mathbf{a} + \beta\mathbf{b} + \gamma\mathbf{c}(α,β,γ≥0\alpha, \beta, \gamma \ge 0,α+β+γ=1\alpha + \beta + \gamma = 1)的点组成的集合。由唯一性,p\mathbf{p} 只有一组重心坐标,所以 p\mathbf{p} 属于凸包,当且仅当这唯一的一组坐标全部非负。■\blacksquare

(α,β,γ)(\alpha, \beta, \gamma) 叫 p\mathbf{p} 的重心坐标(barycentric coordinates)。若三个顶点分别带有属性值 fa,fb,fcf_a, f_b, f_c(比如颜色),三角形内一点的属性就取 f(p)=αfa+βfb+γfcf(\mathbf{p}) = \alpha f_a + \beta f_b + \gamma f_c。它在顶点处等于顶点的值,在内部线性变化。第 07 篇会推导重心坐标的计算公式,并给出它的几何意义(面积比)。

除了三角形,OpenGL 还能画点(GL_POINTS)和线(GL_LINES)。点在本系列里有专门的用途:三维扫描设备得到的点云只有点、没有三角形,第八节会介绍它。

二、两种渲染思路:光栅化与光线追踪

1. 同一个问题,两种循环顺序

可见性问题要对所有“像素-三角形”组合 (s,k)(s, k) 做判断:视线 ss 碰没碰到 TkT_k?碰到的话有多深?然后对每个像素取最近的一个。这是一个二重循环,而二重循环有两种嵌套顺序:

fig1_paradigms.png

图 1:两种渲染思路的侧视图。眼睛在左边,竖着的一列小方块是屏幕上的 8 个像素,两条粗线是从侧面看过去的两个三角形。左:光栅化把三角形的顶点投影到屏幕上,看它覆盖了哪些像素;右:光线投射从每个像素发出一条射线,找最近的交点。两边最终给像素涂的颜色完全相同。

图像顺序(image order):外层循环像素,内层循环三角形。对每个像素发出一条视线,和所有三角形求交,留下最近的交点。这就是光线投射(ray casting),1968 年由 Appel 提出。如果在交点处继续发出反射、折射、阴影射线,就是 1980 年 Whitted 提出的光线追踪(ray tracing)。

对象顺序(object order):外层循环三角形,内层循环像素。把每个三角形投影到屏幕上,找出它覆盖的像素,逐个处理。这就是光栅化(rasterization)。“光栅”(raster)指的就是像素网格。

两种顺序计算的是同一个 arg⁡min⁡\arg\min,区别只在于“先固定谁”。光线投射在内层循环里就能直接取最小值;光栅化在处理某个三角形时,并不知道后面还有没有更近的三角形,所以需要为每个像素记住“到目前为止最近的深度”。这个记忆就是下面要讲的 Z-buffer。

OpenGL 是光栅化 API。本系列讲的整条渲染管线,本质上就是图 1 左边那个循环在硬件里的实现。

2. Z-buffer:光栅化怎么解决遮挡

最早的想法是画家算法(painter's algorithm):像油画家一样,先画远的,再画近的,近的自然盖住远的。它有两个根本缺陷:

  • 要先把三角形按深度排序,但一个三角形各处的深度不同,“谁更远”不一定有定义;
  • 三角形互相穿插(一个三角形的一半在另一个前面,另一半在后面),或者几个三角形循环遮挡(A 挡 B、B 挡 C、C 又挡 A)时,任何排序都是错的。

1974 年,Catmull 在博士论文中提出了 Z-buffer(深度缓冲)算法。它不比较三角形,而是在每个像素上比较深度:

  1. 准备两块和屏幕一样大的缓冲区:颜色缓冲 color[p] 初始化为背景色,深度缓冲 zbuf[p] 初始化为 +∞+\infty;
  2. 按任意顺序处理三角形。三角形覆盖像素 pp 时,产生一个片元(fragment),即“候选像素”,带有深度 zz 和颜色 cc;
  3. 深度测试:若 z<z < zbuf[p],就令 zbuf[p] =z= z,color[p] =c= c;否则丢弃这个片元。

下面证明它总能得到正确结果。

定理 2(Z-buffer 的正确性):固定一个像素 pp,设落在它上面的片元按处理顺序依次为 (z1,c1),(z2,c2),…,(zn,cn)(z_1, c_1), (z_2, c_2), \dots, (z_n, c_n)。处理完前 kk 个片元后: zbuf[p]=min⁡1≤i≤kzi,color[p]=cj,j=min⁡{ i≤k:zi=min⁡1≤l≤kzl }.\texttt{zbuf}[p] = \min_{1 \le i \le k} z_i,\qquad \texttt{color}[p] = c_j,\quad j = \min\{\, i \le k : z_i = \min_{1 \le l \le k} z_l \,\}. 即深度缓冲里是目前为止的最小深度,颜色缓冲里是最先到达的那个最近片元的颜色。

证明:对 kk 用数学归纳法。

  • k=1k = 1:初值 zbuf[p]=+∞\texttt{zbuf}[p] = +\infty,而 z1<+∞z_1 < +\infty,测试通过,zbuf[p]=z1\texttt{zbuf}[p] = z_1,color[p]=c1\texttt{color}[p] = c_1,结论成立。
  • 由 kk 推到 k+1k + 1:记 mk=min⁡i≤kzim_k = \min_{i \le k} z_i,由归纳假设 zbuf[p]=mk\texttt{zbuf}[p] = m_k。分两种情况:
    • 若 zk+1<mkz_{k+1} < m_k:测试通过,缓冲更新为 zk+1z_{k+1} 和 ck+1c_{k+1}。此时 mk+1=zk+1m_{k+1} = z_{k+1},并且前 kk 个片元的深度都 ≥mk>zk+1\ge m_k > z_{k+1},所以第一个达到最小值的就是第 k+1k+1 个,结论成立。
    • 若 zk+1≥mkz_{k+1} \ge m_k:测试失败,缓冲不变。此时 mk+1=mkm_{k+1} = m_k,最先达到最小值的片元仍然在前 kk 个里,没有变化,结论成立。■\blacksquare

取 k=nk = n 就得到最终结果。由此可以推出两个重要的结论:

推论 1(与顺序无关):如果像素 pp 上所有片元的深度互不相同,最终结果就与三角形的处理顺序无关。因为 min⁡\min 与元素的排列顺序无关,最小值又只在一个片元上取到。

推论 2(深度相等时与顺序有关):若有几个片元深度恰好相等并且都是最小的,留下的是最先到达的那个。如果把测试条件改成 z≤z \le zbuf[p](OpenGL 里的 GL_LEQUAL),留下的就变成最后到达的那个。第 13 篇会看到,两个面几乎重合时,深度在浮点精度下时而相等、时而大小颠倒,画面就会出现闪烁的条纹,这叫 z-fighting。

Z-buffer 每个像素只做一次比较,与场景复杂度无关;三角形可以按任意顺序、彼此独立地处理,互相穿插也没关系。因为它是逐像素比较深度,而在每一个像素上,“谁更近”总是有确定答案的。代价是需要一块额外的内存,这在 1974 年很贵,今天微不足道。正是这些性质,让它成了所有 GPU 的标准做法。OpenGL 里用 glEnable(GL_DEPTH_TEST) 打开深度测试,用 glDepthFunc(GL_LESS) 选择比较方式。

3. 两种方法算出的是同一张图

本篇附录代码用纯 Python 把两种方法各实现了一遍:场景里有 3 个带顶点颜色的三角形,其中前面两个互相穿插,渲染成 160×120 的图像。

fig2_two_renderers.png

图 2:附录代码的渲染结果。从左到右:光栅化 + Z-buffer 的结果;光线投射的结果;每个像素收到的片元数(overdraw);光栅化结束后深度缓冲里的值。前两张图完全一致;中间那两块橙色和白色的三角形互相穿插,交线处的遮挡关系是 Z-buffer 逐像素算出来的,画家算法处理不了这种情况。

两张图的对比结果(附录代码的输出):

覆盖情况不一致的像素:0 个
深度最大差异:2.66e-15
颜色最大差异:3.33e-16

差异只有 10−1510^{-15} 量级,是浮点舍入误差。两种方法在完全不同的空间里计算:光栅化在二维屏幕上判断“像素中心在不在投影后的三角形里”,光线投射在三维空间里判断“射线和三角形有没有交点”。结果为什么能一模一样?还有一个细节:附录代码里光栅化算深度用的是 1/(λ0/w0+λ1/w1+λ2/w2)1 / (\lambda_0/w_0 + \lambda_1/w_1 + \lambda_2/w_2),而不是简单地 λ0w0+λ1w1+λ2w2\lambda_0 w_0 + \lambda_1 w_1 + \lambda_2 w_2,这是为什么?下面的定理同时回答这两个问题。

先约定记号(与附录代码一致):眼睛在原点,看向 −z-z 方向。三维点 q=(qx,qy,qz)\mathbf{q} = (q_x, q_y, q_z) 的深度记为 w=−qzw = -q_z。只考虑位于眼睛前方的点(w>0w > 0)。透视投影把 q\mathbf{q} 映射到屏幕上的点

π(q)=(fx qxw, fy qyw)+(ox,oy),\pi(\mathbf{q}) = \left( \frac{f_x\, q_x}{w},\ \frac{f_y\, q_y}{w} \right) + (o_x, o_y),

其中 fx,fy,ox,oyf_x, f_y, o_x, o_y 是由视场角和窗口大小决定的常数(第 12 篇推导)。它的含义很直观:离得越远(ww 越大),投影到屏幕上越靠近中心。

定理 3(透视投影下的重心坐标):设三角形的三个顶点 v0,v1,v2\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2 的深度分别为 w0,w1,w2>0w_0, w_1, w_2 > 0。三角形上一点 q=b0v0+b1v1+b2v2\mathbf{q} = b_0\mathbf{v}_0 + b_1\mathbf{v}_1 + b_2\mathbf{v}_2(b0+b1+b2=1b_0+b_1+b_2 = 1)投影到屏幕上后,它相对于投影三角形 π(v0)π(v1)π(v2)\pi(\mathbf{v}_0)\pi(\mathbf{v}_1)\pi(\mathbf{v}_2) 的重心坐标是 λi=bi wib0w0+b1w1+b2w2,i=0,1,2.\lambda_i = \frac{b_i\, w_i}{b_0 w_0 + b_1 w_1 + b_2 w_2}, \qquad i = 0, 1, 2. 反过来,已知屏幕上的重心坐标 λi\lambda_i,三维空间里的重心坐标和深度是 bi=λi/wi∑jλj/wj,w=1∑jλj/wj.b_i = \frac{\lambda_i / w_i}{\sum_j \lambda_j / w_j}, \qquad w = \frac{1}{\sum_j \lambda_j / w_j}. 特别地,λi≥0\lambda_i \ge 0 对所有 ii 成立,当且仅当 bi≥0b_i \ge 0 对所有 ii 成立。

证明:深度 w=−qzw = -q_z 是坐标的线性函数,所以 q\mathbf{q} 的深度为 w=b0w0+b1w1+b2w2w = b_0 w_0 + b_1 w_1 + b_2 w_2。看投影的 xx 坐标(yy 完全相同)。去掉常数 fxf_x 和 oxo_x 不影响结论,只需考虑 qx/wq_x / w:

qxw=∑ibi vi,xw=∑ibi wiw⋅vi,xwi=∑iλi⋅vi,xwi.\frac{q_x}{w} = \frac{\sum_i b_i\, v_{i,x}}{w} = \sum_i \frac{b_i\, w_i}{w} \cdot \frac{v_{i,x}}{w_i} = \sum_i \lambda_i \cdot \frac{v_{i,x}}{w_i}.

第二步是在每一项上乘以再除以 wiw_i。这里 λi=biwi/w\lambda_i = b_i w_i / w,而 vi,x/wiv_{i,x}/w_i 正是第 ii 个顶点投影后的 xx 坐标(同样去掉了常数)。又因为

∑iλi=∑ibiwiw=ww=1,\sum_i \lambda_i = \frac{\sum_i b_i w_i}{w} = \frac{w}{w} = 1,

而对系数和为 1 的组合,仿射函数 u↦fxu+oxu \mapsto f_x u + o_x 满足 fx∑iλiui+ox=∑iλi (fxui+ox)f_x \sum_i \lambda_i u_i + o_x = \sum_i \lambda_i\,(f_x u_i + o_x),把常数加回去不影响等式,所以 π(q)=∑iλi π(vi)\pi(\mathbf{q}) = \sum_i \lambda_i\, \pi(\mathbf{v}_i)。由定理 1 的唯一性,λi\lambda_i 就是 π(q)\pi(\mathbf{q}) 在投影三角形中的重心坐标。

再求反变换。由 λi=biwi/w\lambda_i = b_i w_i / w 得 bi=λiw/wib_i = \lambda_i w / w_i。两边对 ii 求和,利用 ∑ibi=1\sum_i b_i = 1:

1=w∑iλiwi⟹w=1∑jλj/wj,bi=λi/wi∑jλj/wj.1 = w \sum_i \frac{\lambda_i}{w_i} \quad\Longrightarrow\quad w = \frac{1}{\sum_j \lambda_j / w_j}, \qquad b_i = \frac{\lambda_i / w_i}{\sum_j \lambda_j / w_j}.

最后,由于所有 wi>0w_i > 0、w>0w > 0,λi\lambda_i 与 bib_i 的符号总是相同。■\blacksquare

这个定理有三个直接的推论,正好解释了附录代码的结果:

  1. 覆盖判断等价。像素中心 ss 在投影三角形内(所有 λi≥0\lambda_i \ge 0),当且仅当三维三角形上存在投影到 ss 的点(所有 bi≥0b_i \ge 0),也就是穿过 ss 的视线与三角形相交。所以两种方法认定的“覆盖”完全一致。
  2. 深度要插值 1/w1/w,不能直接插值 ww。1/w=∑jλj⋅(1/wj)1/w = \sum_j \lambda_j \cdot (1/w_j),即 1/w1/w 在屏幕空间里是线性的,而 ww 本身不是。附录代码先插值 1/w1/w 再取倒数,得到的深度与光线投射的参数 tt 精确相等。
  3. 属性要做透视校正。插值顶点颜色、纹理坐标时,必须用三维空间的重心坐标 bib_i,而不是屏幕上的 λi\lambda_i,否则纹理会扭曲。这叫透视校正插值(perspective-correct interpolation),第 07 篇会专门讲,第 15 篇会展示不做校正时纹理变形的样子。GPU 的光栅化器会自动完成这一步。

4. 代价:为什么实时渲染选择了光栅化

两种方法结果相同,代价却很不一样。附录代码统计了两者做的基本测试次数:

光栅化光线投射
基本测试包围盒内“像素中心在不在三角形里”:25551 次射线-三角形求交:57600 次
测试次数与什么成正比每个三角形的包围盒面积之和像素数 × 三角形数

光线投射要测 19200×3=5760019200 \times 3 = 57600 次,因为每条射线都要和每个三角形求交。一般地,像素数为 PP、三角形数为 NN 时,朴素光线投射要做 P⋅NP \cdot N 次求交。对一个 1920×1080、一百万个三角形的场景,这是 2×10122 \times 10^{12} 次,每帧都要做。光线追踪必须借助 BVH 之类的加速结构,把每条射线的代价降到 O(log⁡N)O(\log N) 量级。

光栅化的代价大约是“每个三角形的固定开销 × NN”加上“所有三角形覆盖的像素总数”。后者等于 PP 乘以平均每个像素被覆盖的次数(叫深度复杂度或 overdraw)。图 2 第三张图就是 overdraw:本例中被覆盖的 7489 个像素共收到 11456 个片元,平均 1.53 个。光栅化不需要任何加速结构,而且有三个特别适合硬件实现的性质:

  1. 只需要局部信息。处理一个三角形时,只需要它自己的 3 个顶点,不需要知道场景里的其他东西。所以三角形可以像流水线一样一个接一个送进去,处理完就扔掉;
  2. 内层计算规则而简单。判断像素在不在三角形里、插值属性,都是几次乘加运算,而且相邻像素的计算几乎一样;
  3. 天然并行。不同三角形、不同像素之间互不依赖(只在写 Z-buffer 时需要协调),可以同时处理成千上万个。

这就是为什么从 1980 年代的 SGI 图形工作站到今天的显卡,实时渲染硬件都选择了光栅化。光线追踪的画面更真实(反射、折射、软阴影都很自然),长期以来主要用于电影等离线渲染。2018 年起,消费级显卡加入了光线追踪专用硬件,今天的游戏通常是光栅化为主、光线追踪为辅的混合渲染。第 27 篇会实现射线-三角形求交,用于鼠标拾取。

说明:附录代码里光栅化用时 0.09 秒,光线投射用时 3 秒,差距远大于测试次数之比。这主要是因为 ray_triangle 函数里用了 NumPy 的小数组运算,每次调用都有固定开销。不要把这当作两种算法的公平测速,测试次数才是有意义的比较。

三、OpenGL 是什么

1. 一份规范,而不是一个库

很多初学者以为 OpenGL 是一个可以下载安装的“库”,其实不是。OpenGL 是一份规范(specification),它是一份文档,规定了有哪些函数、每个函数的参数是什么、调用后应该产生什么结果。比如 OpenGL 2.0 规范是一份几百页的 PDF,由 Mark Segal 和 Kurt Akeley 执笔,可以在 Khronos 官网免费下载。

真正干活的是实现(implementation),它通常就是显卡驱动的一部分。NVIDIA、AMD、Intel、Apple 各自按照规范实现 OpenGL,随驱动一起安装。所以:

  • 同一段代码在不同显卡上,跑的是不同厂商写的程序,性能和偶尔的 bug 都可能不同;
  • 你的机器支持哪个版本的 OpenGL,取决于显卡和驱动,不取决于你装了什么软件。第 01 篇会教你查询;
  • 规范最初由 SGI 公司在 1992 年发布,之后由 ARB(架构评审委员会,Architecture Review Board)维护,2006 年起移交给 Khronos Group 这个行业联盟。Khronos 同时负责 OpenGL ES、WebGL 和 Vulkan。

在 Windows 上还有一个历史遗留问题:系统自带的 opengl32.dll 只直接导出 OpenGL 1.1 的函数。更高版本的函数要在运行时通过 wglGetProcAddress 向驱动查询函数地址再调用。C/C++ 程序通常用 GLEW、glad 这类“加载库”来做这件事。PyOpenGL 会自动处理,所以在 Python 里直接调用 glCreateShader 就行,不用操心这一层。

2. OpenGL 不管什么

规范只管“渲染”这一件事,下面这些它都不管:

OpenGL 不负责本系列用什么
创建窗口、处理键盘鼠标glfw(第 01 篇)
创建 OpenGL 上下文glfw(底层调用 WGL / GLX / CGL 等系统接口)
读取图片文件Pillow(第 15 篇)
向量、矩阵运算NumPy,矩阵自己写(第 08~12 篇)
读取三维模型文件自己写 OBJ 解析器(第 21 篇)

这里的上下文(context)是 OpenGL 的核心概念:它是一个容器,装着 OpenGL 的全部状态和你创建的所有对象(缓冲、纹理、着色器……),并且和一个窗口绑定。调用任何 OpenGL 函数之前,必须先有一个“当前上下文”。在 Python 里最常见的错误之一,就是还没创建窗口就调用了 GL 函数,结果程序直接崩溃。glfw 用两行代码就能完成这件事:glfw.create_window(...) 和 glfw.make_context_current(window)。

3. 客户端-服务端模型与异步执行

OpenGL 规范第 2 章开头就写明:OpenGL 采用客户端-服务端模型(client-server model)。你的程序是客户端,负责发出命令;OpenGL 实现(驱动 + GPU)是服务端,负责执行命令。

fig3_client_server.png 图 3:上:一次 GL 调用从 Python 程序经过 PyOpenGL、驱动,进入命令缓冲,再由 GPU 执行,结果写入帧缓冲并送到显示器。下:CPU 和 GPU 的时间线。CPU 发出命令后立刻返回,GPU 稍后才执行;只有 glFinish() 这类函数会让 CPU 停下来等 GPU。

这个模型最重要的推论是:OpenGL 命令是异步执行的。你调用 glDrawArrays 时,驱动只是把这条命令放进命令缓冲区就返回了,GPU 可能要过一会儿才真正开始画。CPU 和 GPU 像流水线的两道工序一样并行工作:GPU 画这一帧时,CPU 已经在准备下一帧了。这样两边都不闲着。

异步执行带来几个实际的后果,后面各篇都会遇到:

  1. 不能用 Python 的计时函数直接测某个 GL 调用的耗时。time.perf_counter() 只能测到“命令入队”的时间。要测 GPU 的真实耗时,需要用 glFinish() 强制等待(但这样会打断流水线),或者用 GPU 计时查询(第 34 篇);
  2. 从 GPU 读回数据会让 CPU 停下来等。比如 glReadPixels 读帧缓冲、glGetError 查询错误,驱动都必须先把前面排队的命令执行完。每帧大量调用这类函数,会让程序变慢;
  3. glFlush() 与 glFinish() 的区别:glFlush() 只要求驱动“尽快把命令发给 GPU”,立即返回;glFinish() 要等所有命令都执行完才返回。双缓冲程序每帧调用 glfw.swap_buffers(),它会隐含地刷新命令,所以一般两者都不用手动调用。

4. OpenGL 是一个状态机

OpenGL 的“红宝书”(《OpenGL Programming Guide》)第 1 章有一句名言:OpenGL 是一个状态机(state machine)。上下文里保存着大量状态变量:是否开启深度测试、当前清屏颜色、当前视口、当前使用哪个着色器程序、当前绑定了哪个缓冲……你设置一个状态后,它会一直保持,直到你再次改变它。绘制命令则按调用那一刻的全部状态来执行。

fig4_state_machine.png

图 4:上下文里保存着各种状态。左边的调用修改状态;有些状态是开关或数值,有些是“绑定点”,指向右边的某个对象。glDrawArrays 按此刻的全部状态执行:用 3 号程序、从 1 号缓冲读顶点、开启深度测试……

OpenGL 里的数据放在对象(object)里,比如缓冲对象存顶点数据,纹理对象存图像。每个对象用一个整数“名字”标识。但是 OpenGL 2.0 的大多数函数并不接受对象名作为参数,而是作用在“当前绑定到某个绑定点上的对象”上。这种模式叫先绑定、再操作(bind-to-edit):

# 示意代码(第 05 篇会完整运行)
vbo_a, vbo_b = glGenBuffers(2)                # 创建两个缓冲对象,得到两个整数名字
glBindBuffer(GL_ARRAY_BUFFER, vbo_a)          # 状态改变:GL_ARRAY_BUFFER 绑定点 → vbo_a
glBufferData(GL_ARRAY_BUFFER, data_a, GL_STATIC_DRAW)   # 数据写进“当前绑定的” vbo_a
glBindBuffer(GL_ARRAY_BUFFER, vbo_b)          # 绑定点改为 vbo_b
glBufferData(GL_ARRAY_BUFFER, data_b, GL_STATIC_DRAW)   # 这次写进 vbo_b

注意 glBufferData 的参数里根本没有 vbo_a 或 vbo_b。它写进哪个缓冲,完全取决于之前绑定了谁。

状态机设计的好处是函数参数少、驱动容易优化,不足是容易出现“状态泄漏”:某段代码改了一个状态忘记恢复,结果影响了后面完全不相关的绘制。这类 bug 很难查,因为出错的地方和原因相隔很远。本系列的代码会养成两个习惯:用到什么状态,就在用之前明确设置它;绘制完把绑定恢复为 0。OpenGL 4.5 引入的“直接状态访问”(DSA)允许直接按对象名操作,就是为了缓解这个问题(第 35 篇)。

5. 扩展机制

规范更新很慢,而硬件厂商希望尽快让新功能能用起来。OpenGL 的办法是扩展(extension):厂商可以在驱动里提供规范之外的新函数和新常量,名字带有前缀:

前缀含义例子
GL_NV_、GL_AMD_ 等某一家厂商的扩展GL_NV_register_combiners
GL_EXT_多家厂商共同支持GL_EXT_framebuffer_object
GL_ARB_ARB 正式批准,通常会进入下一版核心规范GL_ARB_vertex_buffer_object

一个功能常常走过“厂商扩展 → EXT → ARB → 核心”这条路。例如顶点缓冲先是 GL_ARB_vertex_buffer_object,在 1.5 版进入核心;帧缓冲对象(FBO)在 2.0 时代还是 GL_EXT_framebuffer_object 扩展,到 3.0 才进入核心。第 22 篇讲 FBO 时会用到这段历史。

四、从固定管线到可编程管线:OpenGL 版本简史

1. 时间线

fig5_timeline.png

图 5:桌面 OpenGL、OpenGL ES、WebGL 和新一代图形 API 的时间线。2.0 版(2004 年)开启了可编程管线时代;2007 年的 OpenGL ES 2.0 和 2011 年的 WebGL 1.0 都采用了同一套着色器模型。

版本年份主要新增内容
1.01992由 SGI 的 IRIS GL 发展而来;立即模式、矩阵栈、固定光照
1.11997纹理对象、顶点数组、多边形偏移
1.21998三维纹理、GL_CLAMP_TO_EDGE 等
1.32001多重纹理、立方体贴图、多重采样、压缩纹理
1.42002自动生成 mipmap、深度纹理、阴影比较
1.52003顶点缓冲对象(VBO)、遮挡查询
2.02004GLSL 1.10 着色器、多渲染目标、非 2 的幂纹理、点精灵、双面模板
2.12006GLSL 1.20、像素缓冲对象、sRGB 纹理
3.02008引入“弃用”机制;FBO、VAO 进入核心
3.12009删除被弃用的固定管线功能(除非驱动提供兼容扩展)
3.22009区分核心模式(core profile)与兼容模式(compatibility profile);几何着色器
3.32010GLSL 3.30;现代 OpenGL 教程最常用的基线版本
4.02010细分曲面着色器
4.32012计算着色器、调试输出
4.62017最后一个版本;支持 SPIR-V 着色器

2017 年的 4.6 是桌面 OpenGL 的最后一个版本,之后 Khronos 的工作重心转向了 2016 年发布的 Vulkan。苹果在 2018 年宣布弃用 OpenGL,macOS 上最高只支持到 4.1。不过 OpenGL 并没有消失:Windows 和 Linux 的驱动仍然完整支持它,大量科研、工业、教学软件仍在使用它,WebGL 也仍是网页 3D 的主力。

2. 固定管线是什么

OpenGL 1.x 的渲染流程是固定功能管线(fixed-function pipeline):每个阶段做什么都由硬件写死,程序员只能通过开关和参数来控制。例如光照:

# OpenGL 1.x 的固定光照(第 04 篇会演示)
glEnable(GL_LIGHTING)
glEnable(GL_LIGHT0)
glLightfv(GL_LIGHT0, GL_POSITION, (1.0, 1.0, 1.0, 0.0))
glMaterialfv(GL_FRONT, GL_DIFFUSE, (0.8, 0.2, 0.2, 1.0))

驱动会用一个固定的公式(大致是第 18 篇要推导的 Blinn-Phong 模型),在每个顶点上算光照,再把颜色插值到三角形内部。你想换一个光照公式?做不到。想在每个像素上算光照?也做不到。想做卡通渲染、法线贴图、阴影?只能靠各种扩展和巧妙的技巧拼凑。

1999 年 NVIDIA 的 GeForce 256 把“变换与光照”(T&L)放进了硬件,固定管线达到了顶峰,但它的局限也越来越明显。2001 年的 GeForce 3 开始支持可编程的顶点和像素处理,最初只能用类似汇编的语言编写(如 2002 年的 ARB_vertex_program)。同一时期,3Dlabs 公司牵头起草了 “OpenGL 2.0” 提案,核心就是一门类 C 的高级着色语言。它在 2003 年以 ARB 扩展的形式出现,2004 年正式成为 OpenGL 2.0 核心的一部分,这就是 GLSL(OpenGL Shading Language)。

3. OpenGL 2.0 带来了什么

fig6_pipeline.png

图 6:上:OpenGL 1.x 的固定管线。下:OpenGL 2.0 的可编程管线,“变换与光照”被顶点着色器取代,“纹理环境、颜色求和、雾”被片元着色器取代,其余阶段不变。

2.0 版最重要的变化,是让管线中的两个阶段可编程:

  • 顶点着色器(vertex shader):每个顶点运行一次。输入顶点的位置和属性,输出它在裁剪空间中的位置(第 10 篇),以及要传给片元的数据。它取代了固定管线的“变换与光照”;
  • 片元着色器(fragment shader):每个片元运行一次。输入插值后的数据,输出这个片元的颜色。它取代了固定管线的“纹理环境、颜色求和、雾”。

其他阶段(图元装配、裁剪、光栅化、深度测试、混合)仍是固定功能,只能配置,不能编程。一对最简单的 GLSL 1.10 着色器长这样(第 06 篇会逐行讲解):

// 顶点着色器
#version 110
attribute vec3 a_position;    // 每个顶点的输入
attribute vec3 a_color;
uniform mat4 u_mvp;           // 所有顶点共用的输入:变换矩阵
varying vec3 v_color;         // 传给片元着色器,光栅化时自动插值
void main() {
    v_color = a_color;
    gl_Position = u_mvp * vec4(a_position, 1.0);
}
// 片元着色器
#version 110
varying vec3 v_color;         // 插值后的颜色(GPU 已经做了第二节说的透视校正)
void main() {
    gl_FragColor = vec4(v_color, 1.0);
}

除了 GLSL,2.0 还把几个扩展纳入了核心:多渲染目标(一次绘制同时输出到多张图,glDrawBuffers)、非 2 的幂尺寸的纹理、点精灵(第 26 篇)、正反面分别设置的模板测试(第 25 篇)、分别设置 RGB 和 alpha 的混合方程(第 14 篇)。

4. 之后的故事

可编程管线证明了自己的价值以后,固定管线就成了累赘:驱动要同时维护两套机制,API 越来越臃肿。3.0 版(2008 年)引入了弃用机制,3.1 版删除了固定管线,3.2 版把规范分成两种模式:

  • 核心模式(core profile):只保留现代功能,删掉 glBegin、矩阵栈、固定光照等;
  • 兼容模式(compatibility profile):保留所有旧功能,老代码照常运行。

移动设备上的 OpenGL ES 2.0(2007 年)走得更彻底:它完全没有固定管线,必须用着色器。它的着色器模型几乎就是 OpenGL 2.0 的翻版。2011 年的 WebGL 1.0 又把 ES 2.0 原样搬进了浏览器。

五、为什么今天还要学 OpenGL 2.0 风格

OpenGL 2.0 是 2004 年的规范,今天的新项目通常用 OpenGL 3.3+、Vulkan、Metal、Direct3D 12 或 WebGPU。那为什么还要从 2.0 学起?下面是诚实的利弊对比:

说明
✅ 概念最少顶点缓冲、着色器、纹理、帧缓冲这些核心概念一个不少,但没有 VAO、UBO 等额外的“仪式性”代码,初学者可以把注意力放在原理上
✅ 可编程管线的全部思想都在第 03 篇的管线、第 10 篇的坐标变换、第 18 篇的光照,在 2.0 和 Vulkan 里本质完全相同
✅ 直接迁移到 ES 2.0 和 WebGL 1.0着色器几乎不用改,API 调用基本一致
✅ 到处都能运行Windows、Linux 的驱动默认创建兼容模式上下文,完全支持 2.0;macOS 的默认上下文就是 2.1
✅ 与 Python 配合简单PyOpenGL 对 2.0 的支持最成熟,网上资料也最多
❌ 有些写法已经过时attribute、varying、gl_FragColor 在 GLSL 3.30 里改了名字;没有 VAO
❌ 缺少现代功能没有几何着色器、实例化绘制、统一缓冲对象等;FBO 在 2.0 里还是扩展
❌ 不代表工业界的现状真实项目中的性能优化方法(多线程录制命令、显式内存管理)要到 Vulkan 才能学到

为了扬长避短,本系列做如下约定:

  1. 以着色器为中心。除第 04 篇用固定管线做对比演示外,从第 05、06 篇开始,所有绘制都通过 VBO 和 GLSL 着色器完成;
  2. 不使用 GLSL 的内置矩阵和内置顶点属性。2.0 允许在着色器里直接读 gl_ModelViewMatrix、gl_Vertex,但本系列一律用自定义的 attribute 和 uniform,矩阵在 Python 里用 NumPy 自己算。这样做和 ES 2.0、WebGL、现代 OpenGL 完全一致,而且矩阵的每一步数学都看得见;
  3. 着色器统一以 #version 110 开头,变量命名用前缀区分:a_ 表示 attribute,u_ 表示 uniform,v_ 表示 varying;
  4. 每一篇末尾都有“与现代 OpenGL 的差异”一节,说明本篇的代码在 3.3 core 里怎么写。学完全系列,你会同时掌握两种写法。第 35 篇做全面对照。

六、GPU 为什么能这么快

先算一笔账。一块 1920×1080 的屏幕有 2,073,600 个像素,每秒刷新 60 次,每秒就要输出 2073600×60=124,416,0002073600 \times 60 = 124{,}416{,}000 个像素。考虑到 overdraw(平均每个像素被覆盖两三次),每秒要处理三四亿个片元。如果每个片元的着色计算有一百次运算,每秒就需要几百亿次运算,而且还要同时变换几百万个顶点。一个 CPU 核心每秒大约只能做几十亿次标量运算,远远不够。

GPU 的解决办法是大规模并行。回忆第二节的结论:光栅化中,不同顶点之间、不同片元之间的计算互不依赖。顶点着色器对每个顶点运行同一段程序,片元着色器对每个片元运行同一段程序,只是输入数据不同。这种计算叫数据并行(data parallel)。

于是 GPU 采用了和 CPU 截然不同的设计:

CPUGPU
核心数量几个到几十个“大核”成千上万个简单的运算单元
设计目标让一条指令流跑得尽可能快(复杂的分支预测、乱序执行、大缓存)让海量相同的计算同时完成(总吞吐量最大)
执行方式每个核独立执行自己的程序一组线程(如 32 个)同时执行同一条指令,各自处理不同的数据(SIMT)
隐藏内存延迟的方法大缓存同时挂着大量线程,一组在等内存时切换到另一组

这对写着色器有直接的影响:

  1. 着色器只描述“一个”顶点或“一个”片元怎么算,不能写“遍历所有像素”的循环。外层循环是 GPU 硬件替你做的,就是第二节图 1 左边的那两层循环;
  2. 着色器之间不能互相通信。一个片元着色器读不到相邻片元的计算结果(第 22 篇的后处理用“先渲染到纹理、再读纹理”来绕过这个限制);
  3. 分支要谨慎。同一组线程执行同一条指令,如果 if 的两个分支各有线程要走,两个分支都要执行一遍,只是各自屏蔽掉不该走的线程。

七、用 Python 写 OpenGL:要注意什么

C/C++ 是 OpenGL 的“母语”,用 Python 写 OpenGL 有一个绕不开的问题:跨语言调用的开销。

PyOpenGL 是通过 Python 标准库的 ctypes 调用驱动里的 C 函数的。每次调用,都要把 Python 对象转换成 C 的参数,调用结束再把结果转回来。附录代码实测了在本机上通过 ctypes 调用一次 C 标准库函数 abs() 的耗时:

一次 ctypes 调用约 359 ns
100000 个三角形逐顶点提交需要 600000 次调用,至少 215 ms/帧,帧率不超过 5 fps

(你的机器上数字会不同,但量级相近。)这只是下限:PyOpenGL 还要做参数检查,默认每次调用后还会调用 glGetError 检查错误,实际开销更大,第 01 篇会实测。

这个数字说明了 OpenGL 1.x 立即模式(immediate mode)为什么不能用于正经的程序:

# 立即模式:每个顶点调用两次 GL 函数(第 04 篇演示,之后不再使用)
glBegin(GL_TRIANGLES)
for v, c in zip(vertices, colors):     # 30 万个顶点就要 60 万次调用
    glColor3f(*c)
    glVertex3f(*v)
glEnd()

在 C 语言里,一次函数调用只要几纳秒,这样写还能凑合;在 Python 里每帧就要两百多毫秒。正确的做法是:用 NumPy 数组一次性把所有顶点数据交给 GPU,存进顶点缓冲(VBO),然后每帧只用一个 glDrawArrays 调用画出全部三角形。这样 GL 调用的次数和顶点个数无关。这就是第 05 篇的内容,也是本系列贯穿始终的原则:

Python 端只做“指挥”,不做“体力活”。 逐顶点、逐像素的计算交给 GPU(着色器)或 NumPy(向量化运算),Python 循环只遍历“物体”这一级。

还要注意,GPU 处理的几乎都是 32 位浮点数,而 NumPy 默认用 64 位的 float64。传给 OpenGL 的数组必须明确指定 dtype=np.float32,否则数据会被错误解读,或者被悄悄转换一次。第 05 篇会讲清楚。

八、另一种三维数据:点云

到目前为止,场景都是由三角形组成的。三维世界还有另一种极其常见的数字形式:点云(point cloud)。自动驾驶汽车车顶的激光雷达、手机和机器人上的深度相机、无人机航测、文物和建筑的三维扫描,得到的原始数据都是点云。本系列用第六部分(第 29~32 篇)专门讲它,这里先建立基本概念,并用附录代码亲手“扫描”出一份点云。

1. 什么是点云

点云就是三维空间中一组点的集合:

P={p1,p2,…,pN},pi∈R3.\mathcal{P} = \{\mathbf{p}_1, \mathbf{p}_2, \dots, \mathbf{p}_N\}, \qquad \mathbf{p}_i \in \mathbb{R}^3 .

每个点除了坐标,还常常带有属性:颜色、反射强度(激光雷达测得的回波强弱)、法线、时间戳、语义标签等。在程序里,点云就是一个 N×KN \times K 的 NumPy 数组,每行一个点,前 3 列是坐标,后面是属性。这正好就是 OpenGL 顶点缓冲需要的格式。

它和三角形网格的区别如下:

三角形网格点云
基本元素顶点 + 三角形只有点
连接关系有:每个三角形记录由哪 3 个顶点组成没有:点与点之间互不相干
表面连续的,三角形内部可以用重心坐标插值离散的样本,点和点之间是空的
法线用两条边的叉积直接算出要从邻近的点估计(第 31 篇)
典型来源建模软件、从点云重建传感器直接测量
OpenGL 画法glDrawArrays(GL_TRIANGLES, ...)glDrawArrays(GL_POINTS, ...)

2. 点云从哪里来

激光雷达(LiDAR)。它向一个已知方向发出激光脉冲,测量回波返回的时间 Δt\Delta t。光走了一个来回,所以距离为

r=c Δt2,c≈3×108 m/s.r = \frac{c\,\Delta t}{2}, \qquad c \approx 3 \times 10^8\ \text{m/s}.

例如 Δt=66.7\Delta t = 66.7 纳秒,r=3×108×66.7×10−9/2≈10.0r = 3 \times 10^8 \times 66.7 \times 10^{-9} / 2 \approx 10.0 米。设激光方向的仰角为 θ\theta、方位角为 φ\varphi,在传感器自身的坐标系中(xx 向前、yy 向左、zz 向上),测到的点为

p=r (cos⁡θcos⁡φ, cos⁡θsin⁡φ, sin⁡θ)T.\mathbf{p} = r\,\big(\cos\theta\cos\varphi,\ \cos\theta\sin\varphi,\ \sin\theta\big)^\mathsf{T}.

括号里的向量长度为 1,因为 cos⁡2θ(cos⁡2φ+sin⁡2φ)+sin⁡2θ=1\cos^2\theta(\cos^2\varphi + \sin^2\varphi) + \sin^2\theta = 1。激光雷达旋转着向成千上万个方向发射,每转一圈就得到一帧点云。激光雷达就是物理世界中的光线投射:第二节图 1 右边的算法,用真实的激光束做了一遍。

深度相机(RGB-D 相机)。它输出一张“深度图”:每个像素 (x,y)(x, y) 上存的不是颜色,而是该像素看到的点的深度 ww。由第二节的透视投影公式 x=fxqx/w+oxx = f_x q_x / w + o_x,已知 xx 和 ww 就能解出 qxq_x,yy 同理:

qx=(x−ox) wfx,qy=(y−oy) wfy,qz=−w.q_x = \frac{(x - o_x)\, w}{f_x}, \qquad q_y = \frac{(y - o_y)\, w}{f_y}, \qquad q_z = -w .

这叫反投影(back-projection),是透视投影的逆运算。投影本身会丢掉深度,所以不可逆;但深度相机把深度测了出来,就能把每个像素还原成一个三维点。一张 640×480 的深度图最多给出 30 万个点。

多视图重建。从不同角度拍几十上百张照片,用摄影测量(structure from motion、多视图立体)算法恢复出场景的点云。

从网格表面采样。机器学习中常把 CAD 模型表面随机采样成点云(如 ModelNet 数据集)。采样一个三角形内的点,用的正是定理 1 的重心坐标。

3. 动手:模拟一次激光雷达扫描

附录代码的 lidar_scan 函数在第二节的场景里放了一台虚拟的激光雷达:位置在相机右前方,朝向场景中心,竖直方向 48 线、每线水平 160 个采样,竖直视场 50°、水平视场 90°。每束激光都用光线投射的 ray_triangle 函数求最近的交点,交点就是一个测量点。

【点云:模拟激光雷达扫描】
    发出 48×160 = 7680 束激光,命中 2919 束,得到 2919 个点
    距离范围 3.23 ~ 9.57,已保存到 scan.ply
    最近的 20%的点:平均距离 4.91,平均点间距 0.0674,点间距/距离 0.0142
    最远的 20%的点:平均距离 7.49,平均点间距 0.0900,点间距/距离 0.0120
    在相机视角下画出点云,被点覆盖的像素 5385 个

点云被保存成 PLY 格式。PLY 是斯坦福大学在 1990 年代为三维扫描数据设计的格式,文件开头是一段文字说明“有多少个点、每个点有哪些属性”,后面每行一个点:

ply
format ascii 1.0
element vertex 2919
property float x
property float y
property float z
property float range
end_header
0.7911 -2.7764 -6.0118 6.6843
0.8563 -2.7491 -6.0254 6.6439
……

这个文件可以直接用 MeshLab、CloudCompare 等免费软件打开查看。第 29 篇会讲 PLY 的二进制格式,以及 PCD(点云库 PCL 的格式)、KITTI 数据集的 .bin 格式、测绘行业的 LAS 格式。

fig7_point_cloud.png

图 7:左:俯视图。橙色方块是激光雷达,灰线是中间一行的激光束,彩色小点是测得的点(颜色表示到激光雷达的距离)。背景三角形上有大片空白,那是被前面两个三角形挡住、激光照不到的区域。中:第二节中相机看到的三角形。右:同一台相机看到的点云,每个点画成 2×2 个像素。

从图 7 可以看出点云的三个特点:

第一,点云只记录传感器看得见的表面。 激光雷达和相机不在同一个位置,从相机看过去,背景三角形上有几块区域本来应该能看到,点云里却是空的。那是激光雷达视角下的“阴影”。一台传感器只能得到场景的一个侧面,要得到完整的模型,需要从多个位置扫描再把它们拼起来(点云配准,常用 ICP 算法)。

第二,点的密度不均匀,离得越远越稀疏。 相邻两束激光的夹角为 Δφ\Delta\varphi。它们打到距离 rr 处、与激光垂直的平面上,两个落点之间的距离是以 rr 为半径、圆心角为 Δφ\Delta\varphi 的弦长:

s⊥=2rsin⁡Δφ2≈r Δφ(Δφ 很小时,sin⁡x≈x).s_\perp = 2r\sin\frac{\Delta\varphi}{2} \approx r\,\Delta\varphi \qquad (\Delta\varphi \text{ 很小时,} \sin x \approx x).

如果被打到的表面不与激光垂直,而是倾斜了一个入射角 α\alpha(激光方向与表面法线的夹角),落点在表面上的间距会被拉长为

s≈r Δφcos⁡α.s \approx \frac{r\,\Delta\varphi}{\cos\alpha}.

这是因为在很小的范围内,两束激光可以看成平行线,它们在垂直平面上相距 s⊥s_\perp,与倾斜平面相交时,交线与垂直平面的夹角正是 α\alpha,由直角三角形得 s⊥=scos⁡αs_\perp = s\cos\alpha。水平、竖直两个方向都是如此,所以单位面积上的点数(点密度)为

ρ≈cos⁡αr2 Δφ Δθ  ∝  cos⁡αr2.\rho \approx \frac{\cos\alpha}{r^2\,\Delta\varphi\,\Delta\theta} \;\propto\; \frac{\cos\alpha}{r^2}.

距离翻倍,点密度就降为四分之一。附录代码验证了这一点:本例中水平角间隔 Δφ=90°/159≈0.00988\Delta\varphi = 90° / 159 \approx 0.00988 弧度,最近的 20% 的点和最远的 20% 的点,平均距离相差 1.53 倍,点间距相差 1.34 倍;而“点间距 / 距离”两组分别是 0.0142 和 0.0120,比距离本身稳定得多,都在 Δφ\Delta\varphi 的 1.2~1.4 倍左右。多出来的倍数来自 1/cos⁡α1/\cos\alpha,粗略地说,相当于入射角在 35°~46° 之间。

第三,点和点之间是空的。 点云没有“面”。图 7 右边被点覆盖的像素只有 5385 个,而同样视角下三角形覆盖了 7489 个像素;离得近的地方,点与点之间还能看到背景。这就是为什么点云渲染时常把每个点画成好几个像素大小,甚至画成小圆盘。

4. 在 OpenGL 里画点云

OpenGL 画点云的方式非常直接:把 N×3N \times 3 的 float32 数组放进顶点缓冲,然后一次 glDrawArrays(GL_POINTS, 0, N)。每个顶点被光栅化成一个以其投影位置为中心、边长为“点大小”个像素的正方形。点大小可以用 glPointSize 统一设置,也可以在顶点着色器里给每个点单独写 gl_PointSize(2.0 里要先 glEnable(GL_VERTEX_PROGRAM_POINT_SIZE))。深度测试和三角形完全一样:附录代码的 render_points 正是这个过程的 Python 版本,它的 Z-buffer 逻辑与 rasterize 一模一样。

真实的点云有几个问题需要专门处理,这就是第六部分的内容:

问题解决方法篇目
从各种文件格式读取,以及没有真实数据时自己生成XYZ、PLY、PCD、KITTI .bin 的读写;可配置的模拟激光雷达第 29 篇
点之间有缝隙、方形点不好看、怎样着色才看得清随距离变化的点大小(推导 ∝1/w\propto 1/w)、圆形点、按高度/距离/强度着色第 30 篇
没有法线就没法打光,看不出起伏用邻近点的协方差矩阵估计法线(PCA),面元渲染,屏幕空间的 Eye-Dome Lighting第 31 篇
点太多:一次测绘扫描可达上亿个点体素下采样、八叉树、按屏幕尺寸选择细节层次(LOD)第 32 篇

本系列的点云数据:和三角形场景一样,点云数据也由代码生成,不用下载就能运行全部代码。第 29 篇会把本篇的 lidar_scan 扩展成一个可配置的模拟激光雷达,它生成的点云是第六部分的标准例子。同时,第 29 篇也会教你读取真实数据,并列出几个公开数据集供你自行下载:斯坦福三维扫描库(著名的 Stanford Bunny 兔子模型)、KITTI 自动驾驶数据集的激光雷达数据,以及机器学习常用的 ModelNet40。

九、本系列的符号与坐标约定

为了前后一致,全系列统一使用下面的符号和约定。第 08~12 篇会逐条给出依据。

符号

符号含义
v, p\mathbf{v},\ \mathbf{p}(粗体小写)向量、点
M\mathbf{M}(粗体大写)矩阵
n^\hat{\mathbf{n}}单位向量(本系列用“帽子”表示长度为 1)
a⋅b\mathbf{a} \cdot \mathbf{b},a×b\mathbf{a} \times \mathbf{b},∥v∥\lVert \mathbf{v} \rVert点积、叉积、长度
MT\mathbf{M}^\mathsf{T},M−1\mathbf{M}^{-1}转置、逆矩阵
(x,y,z,w)T(x, y, z, w)^\mathsf{T}齐次坐标:点取 w=1w = 1,方向取 w=0w = 0(第 09 篇)
M, V, P\mathbf{M},\ \mathbf{V},\ \mathbf{P}模型矩阵、观察矩阵、投影矩阵(第 10 篇)
pobj, pworld, peye, pclip, pndc, pwin\mathbf{p}_{\text{obj}},\ \mathbf{p}_{\text{world}},\ \mathbf{p}_{\text{eye}},\ \mathbf{p}_{\text{clip}},\ \mathbf{p}_{\text{ndc}},\ \mathbf{p}_{\text{win}}同一个点在物体、世界、观察、裁剪、NDC、窗口坐标系中的坐标
(λ0,λ1,λ2)(\lambda_0, \lambda_1, \lambda_2)屏幕空间的重心坐标
(b0,b1,b2)(b_0, b_1, b_2)三维空间(透视校正后)的重心坐标
ww观察空间中的深度 w=−zeyew = -z_{\text{eye}},在眼睛前方时为正
P={pi}i=1N\mathcal{P} = \{\mathbf{p}_i\}_{i=1}^{N}由 NN 个点组成的点云(第六部分)
rr,(θ,φ)(\theta, \varphi)激光雷达测得的距离,激光的仰角和方位角

约定

  1. 列向量,矩阵左乘:变换写成 p′=Mp\mathbf{p}' = \mathbf{M}\mathbf{p}。组合变换 M=TRS\mathbf{M} = \mathbf{T}\mathbf{R}\mathbf{S} 作用在点上时,从右往左生效:先缩放、再旋转、最后平移。这与 OpenGL 文档和绝大多数欧美教材一致;
  2. 观察空间是右手坐标系:眼睛在原点,xx 轴向右,yy 轴向上,看向 −z-z 方向;
  3. NDC(规范化设备坐标)是立方体 [−1,1]3[-1, 1]^3,zz 轴指向屏幕里面。注意这是一个左手坐标系,翻转发生在投影矩阵里(第 12 篇);
  4. 窗口坐标的原点在左下角,yy 轴向上,像素 (i,j)(i, j) 的中心在 (i+0.5, j+0.5)(i + 0.5,\ j + 0.5)。这和图像文件的惯例(原点在左上角、yy 轴向下)相反,加载纹理和截图时要上下翻转(第 02、15 篇)。本系列配图中的渲染结果也按左下角为原点显示;
  5. 颜色的每个分量都是 [0,1][0, 1] 之间的浮点数,RGBA 顺序;
  6. 数学推导中角度一律用弧度。旧 API(如 glRotatef、gluPerspective)的参数用角度,用到时会特别说明;
  7. 矩阵在 NumPy 中按数学写法存储(M[i, j] 是第 ii 行第 jj 列),而 OpenGL 期望列主序的数据。上传时怎样转换,第 08 篇专门讨论。

十、本系列的学习路线

fig8_roadmap.png

图 8:本系列学习路线,共 36 篇,分为七个部分。加粗的是里程碑篇目。

第一部分:起步与渲染管线

  • 00 导读:就是本篇。
  • 01 环境搭建与第一个窗口:安装 PyOpenGL 和 glfw,查询驱动支持的版本,理解事件循环和双缓冲,实测 Python 调用 GL 的开销。
  • 02 光栅显示与帧缓冲:像素、颜色、帧缓冲的组成;双缓冲为什么能防止画面撕裂。
  • 03 渲染管线全景:把图 6 的每个阶段讲透,每个阶段的输入是什么、输出是什么。
  • 04 第一个三角形:固定管线:用立即模式快速画出图形,建立直观印象,再分析它为什么被淘汰。
  • 05 顶点缓冲 VBO 与索引绘制:把 NumPy 数组送进显存,理解 stride、offset 和内存布局。
  • 06 GLSL 1.10 入门:编译、链接、调试着色器,读懂 attribute、uniform、varying。
  • 07 光栅化原理:边函数、重心坐标、透视校正插值、填充规则、走样。用 NumPy 写一个完整的软件光栅化器,和 GPU 的结果逐像素对照。本篇第二节就是它的预告。

第二部分:变换与坐标系(全系列的数学核心,建议精读)

  • 08 向量与矩阵:点积、叉积的几何意义,“矩阵就是基变换”,行主序与列主序。
  • 09 几何变换与齐次坐标:为什么三维变换要用 4×4 矩阵;绕任意轴旋转的 Rodrigues 公式;法线为什么要用逆转置矩阵变换。
  • 10 坐标系之旅:MVP:一个顶点从模型文件里的坐标,到屏幕上的像素,要经过的每一个坐标系和每一个矩阵。
  • 11 摄像机:LookAt 矩阵的推导,FPS 相机和轨道相机,欧拉角的万向锁,四元数。
  • 12 投影变换:正交投影和透视投影矩阵的完整推导,深度为什么是非线性的,near 和 far 怎么选。
  • 13 深度缓冲与可见性:Z-buffer 的硬件实现、深度精度与 z-fighting、背面剔除。

第三部分:颜色与纹理

  • 14 颜色与混合:混合方程、预乘 alpha、透明物体为什么要排序、Gamma 校正。
  • 15 纹理映射原理:纹理坐标、过滤、Mipmap;从采样定理理解走样和反走样。
  • 16 多纹理与程序化纹理:同时使用多张纹理,在着色器里用数学生成图案和噪声。

第四部分:光照与着色

  • 17 光照的物理直觉:光是怎样被测量的,朗伯余弦定律。
  • 18 Phong 与 Blinn-Phong 模型:经典局部光照模型的推导,各种光源。读完这篇,你会画出第一个有立体感的物体。
  • 19 着色频率:Flat、Gouraud、Phong 着色的区别,逐顶点和逐片元计算。
  • 20 材质、多光源与法线贴图:用切线空间和 TBN 矩阵在平面上“画”出凹凸。
  • 21 模型加载:解析 OBJ 文件,加载真实的三维模型。

第五部分:进阶技术

  • 22 FBO 与后处理:渲染到纹理,用卷积实现模糊和边缘检测。
  • 23 阴影映射:从光源的视角渲染深度,解决阴影痤疮和锯齿。
  • 24 立方体贴图:天空盒,镜面反射,根据斯涅尔定律计算折射。
  • 25 模板缓冲:物体描边,平面镜反射。
  • 26 粒子系统与点精灵:火焰、烟雾、雪花。
  • 27 拾取与射线求交:用鼠标点选三维物体,推导 Möller–Trumbore 算法。
  • 28 曲线与曲面(选修):Bézier 曲线、B 样条、细分曲面。

第六部分:三维点云

  • 29 点云数据与文件格式:激光雷达、深度相机、多视图重建的测量原理;XYZ、PLY(文本与二进制)、PCD、KITTI .bin 等格式的读写;把本篇的 lidar_scan 扩展成可配置的模拟激光雷达,生成第六部分的标准数据。
  • 30 点云渲染:GL_POINTS 与 gl_PointSize,推导透视下点大小应如何随距离变化,用点精灵画圆形点,按高度、距离、强度着色。读完这篇,你会有一个自己的点云查看器。
  • 31 法线估计与点云光照:k 近邻与 KD 树;用邻近点的协方差矩阵估计法线,并证明法线就是最小特征值对应的特征向量;法线朝向一致化;面元(surfel)渲染与 Eye-Dome Lighting。
  • 32 大规模点云:内存与带宽估算,体素下采样,八叉树,按屏幕尺寸选择细节层次(LOD),视锥剔除。

第七部分:综合与展望

  • 33 综合项目:二选一。一是地形漫游,包括高度图地形、相机、光照和阴影;二是点云查看器,加载真实点云,支持多种着色模式和鼠标测距。
  • 34 调试与性能:系统的调试方法,常见性能陷阱与优化。
  • 35 从 2.0 走向现代 OpenGL:全面对照 3.3 core,展望 Vulkan、WebGPU、基于物理的渲染(PBR)和光线追踪。

前置知识

部分需要的数学需要的编程
第一部分高中代数和几何Python、NumPy 基础
第二部分向量、矩阵乘法(第 08 篇从头讲)、三角函数同上
第三、四部分点积的几何意义;一点微积分(积分的直观含义)GLSL(第 06 篇讲)
第五部分无新增同上
第六部分协方差矩阵、特征值与特征向量(第 31 篇从头讲)第 31、32 篇另需 SciPy(KD 树)
第七部分无新增同上

运行环境

本篇只需要 NumPy 和 matplotlib:

pip install numpy matplotlib

从第 01 篇开始需要 PyOpenGL、glfw 和 Pillow:

pip install PyOpenGL glfw Pillow

显卡方面,近十几年的独立显卡和集成显卡都支持 OpenGL 2.0。第 01 篇会教你确认自己的机器支持哪个版本,以及在虚拟机、远程桌面等特殊环境下怎么办。

学习建议

  1. 一定要跑代码,并且改代码。图形学的好处是结果看得见:改一个参数,画面立刻就变。
  2. 推导至少亲手推一遍,尤其是第 07、11、12、18 篇。
  3. 遇到黑屏不要慌。OpenGL 程序出错时,最常见的表现就是一片黑,什么提示都没有。第 06 篇讲着色器的错误日志,第 34 篇给出系统的排查清单。
  4. 第 07 篇和第 12 篇是两个里程碑:读完第 07 篇,你会彻底明白 GPU 是怎么把三角形变成像素的;读完第 12 篇,你会明白三维坐标是怎么变到屏幕上的。
  5. 只关心点云的读者也请先读完第一、二部分,再读第 22 篇(FBO,Eye-Dome Lighting 要用)和第 26 篇(点精灵),然后直接进入第六部分。

十一、与现代 OpenGL 的差异

本系列每一篇的末尾都有这样一节,说明本篇内容在 OpenGL 3.3 核心模式里有什么不同。本篇没有代码,这里先给出全貌,后面各篇再分别细讲:

方面本系列(2.0 风格)OpenGL 3.3 core
创建上下文glfw 默认设置即可要用 window_hint 指定 3.3 版本和核心模式;macOS 上还要开启 forward-compatible
着色器版本声明#version 110#version 330 core
顶点着色器的输入attribute vec3 a_position;layout(location = 0) in vec3 a_position;
顶点到片元的数据varying顶点着色器里 out,片元着色器里 in
片元着色器的输出内置变量 gl_FragColor自己声明 out vec4 frag_color;
纹理采样函数texture2D(...)texture(...)
顶点数组对象(VAO)没有必须绑定一个 VAO 才能绘制
矩阵自己用 NumPy 算,通过 uniform 传入(本系列约定)相同(内置矩阵栈已删除)
glBegin/glEnd、固定光照、GL_QUADS只在第 04 篇演示已删除
帧缓冲对象(FBO)扩展核心功能
在着色器里设置点大小先 glEnable(GL_VERTEX_PROGRAM_POINT_SIZE);圆形点要再开 GL_POINT_SPRITE 才能用 gl_PointCoord开关改名为 GL_PROGRAM_POINT_SIZE;点精灵始终开启,gl_PointCoord 直接可用

创建 3.3 核心模式上下文的 glfw 代码只比默认情况多几行:

glfw.window_hint(glfw.CONTEXT_VERSION_MAJOR, 3)
glfw.window_hint(glfw.CONTEXT_VERSION_MINOR, 3)
glfw.window_hint(glfw.OPENGL_PROFILE, glfw.OPENGL_CORE_PROFILE)
glfw.window_hint(glfw.OPENGL_FORWARD_COMPAT, True)   # macOS 必需
window = glfw.create_window(800, 600, "OpenGL 3.3 core", None, None)

可以看到,两者的差异主要在“写法”上,背后的原理(管线、坐标变换、光照、纹理)完全相同。这也是本系列先学原理、再对照写法的原因。

十二、参考资料

本系列的讲法参考了以下教材、规范与论文(只借鉴思路与结构,内容均为重新组织和推导):

图形学教材

  1. S. Marschner, P. Shirley. Fundamentals of Computer Graphics, 4th ed. CRC Press, 2016. ——“虎书”,本系列图形学理论的主线。本篇对应第 1 章(导论)、第 2 章(三角形与重心坐标)、第 4 章(光线追踪)、第 8 章(图形管线)。
  2. E. Angel, D. Shreiner. Interactive Computer Graphics: A Top-Down Approach with Shader-Based OpenGL, 6th ed. Addison-Wesley, 2012. ——用着色器版 OpenGL 讲图形学,与本系列的思路最接近。本篇对应第 1 章。
  3. T. Akenine-Möller, E. Haines, N. Hoffman et al. Real-Time Rendering, 4th ed. CRC Press, 2018. ——实时渲染的百科全书。本篇对应第 2 章(渲染管线)、第 3 章(GPU)。

OpenGL 专著与规范

  1. D. Shreiner, M. Woo, J. Neider, T. Davis. OpenGL Programming Guide, 5th ed. Addison-Wesley, 2005. ——“红宝书”,对应 OpenGL 2.0。本篇对应第 1 章(状态机与管线概述)。
  2. R. J. Rost. OpenGL Shading Language, 2nd ed. Addison-Wesley, 2006. ——“橙宝书”,GLSL 设计者之一的著作,讲解 GLSL 1.10/1.20 及其设计动机。本篇对应第 1、2 章。
  3. M. Segal, K. Akeley. The OpenGL Graphics System: A Specification, Version 2.0. 2004. ——官方规范,第 2 章开头讲客户端-服务端模型。下载:registry.khronos.org/OpenGL/

数学

  1. F. Dunn, I. Parberry. 3D Math Primer for Graphics and Game Development, 2nd ed. CRC Press, 2011.
  2. E. Lengyel. Mathematics for 3D Game Programming and Computer Graphics, 3rd ed. Cengage, 2011.

论文

  1. A. Appel. Some Techniques for Shading Machine Renderings of Solids. AFIPS Spring Joint Computer Conference, 1968.(光线投射)
  2. E. Catmull. A Subdivision Algorithm for Computer Display of Curved Surfaces. PhD thesis, University of Utah, 1974.(Z-buffer)
  3. I. E. Sutherland, R. F. Sproull, R. A. Schumacker. A Characterization of Ten Hidden-Surface Algorithms. ACM Computing Surveys, 1974.(隐藏面消除算法综述,含画家算法的分析)
  4. T. Whitted. An Improved Illumination Model for Shaded Display. Communications of the ACM, 1980.(递归光线追踪)
  5. T. Möller, B. Trumbore. Fast, Minimum Storage Ray-Triangle Intersection. Journal of Graphics Tools, 1997.(附录代码中的射线-三角形求交,第 27 篇推导)

点云(第六部分会详细引用)

  1. H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, W. Stuetzle. Surface Reconstruction from Unorganized Points. SIGGRAPH 1992.(用 PCA 估计点云法线的经典方法,第 31 篇)
  2. H. Pfister, M. Zwicker, J. van Baar, M. Gross. Surfels: Surface Elements as Rendering Primitives. SIGGRAPH 2000.(面元渲染)
  3. R. B. Rusu, S. Cousins. 3D is here: Point Cloud Library (PCL). ICRA 2011.(PCL 库与 PCD 格式)
  4. A. Geiger, P. Lenz, R. Urtasun. Are We Ready for Autonomous Driving? The KITTI Vision Benchmark Suite. CVPR 2012.(KITTI 数据集)
  5. M. Schütz. Potree: Rendering Large Point Clouds in Web Browsers. Master's thesis, TU Wien, 2016.(大规模点云的八叉树 LOD 渲染,第 32 篇)

在线资源

  1. PyOpenGL 文档:pyopengl.sourceforge.net/;glfw 官网:www.glfw.org/
  2. LearnOpenGL(J. de Vries):learnopengl.com/ ——基于 3.3 core 的优秀教程,可与本系列的“现代差异”一节对照阅读。
  3. The Stanford 3D Scanning Repository:graphics.stanford.edu/data/3Dscan… ——Stanford Bunny 等经典扫描数据(PLY 格式)。
  4. 中文资源:闫令琪的 GAMES101《现代计算机图形学入门》课程,可作为图形学原理部分的补充。

小结

  • 计算机图形学研究怎样用数字描述、生成和操作图像。渲染要回答的问题是:每个像素沿视线看过去,最先看到哪个三角形,它在那里是什么颜色。
  • 实时渲染的主要图元是三角形,因为三角形一定是平的、凸的,而且内部每一点都有唯一的重心坐标,可以用来插值顶点属性。
  • 光线投射外层循环像素,光栅化外层循环三角形,二者计算的是同一个 arg⁡min⁡\arg\min。光栅化只需局部信息、计算规则、天然并行,所以成为 GPU 的选择,OpenGL 是光栅化 API。
  • Z-buffer 在每个像素上保留目前为止最小的深度,用归纳法可以证明它总是正确的;深度互不相同时结果与绘制顺序无关,深度相等时留下最先到达的片元。
  • 透视投影下,1/w1/w 在屏幕空间是线性的。屏幕重心坐标 λi\lambda_i 与三维重心坐标 bib_i 的关系是 bi=(λi/wi)/∑j(λj/wj)b_i = (\lambda_i/w_i) / \sum_j (\lambda_j / w_j)。正因为如此,光栅化与光线投射才能逐像素地得到相同结果。
  • OpenGL 是一份规范,由显卡驱动实现;它只管渲染,窗口、输入、图片、矩阵都要靠别的库。它是一个状态机,大多数函数作用在“当前绑定的对象”上;它采用客户端-服务端模型,命令异步执行。
  • OpenGL 2.0(2004 年)引入了 GLSL,用顶点着色器和片元着色器取代了固定管线的两个阶段,并直接奠定了 OpenGL ES 2.0 和 WebGL 1.0 的基础。
  • 用 Python 写 OpenGL,GL 调用的次数不能和顶点数成正比。批量数据交给 NumPy 和 VBO,逐顶点、逐像素的计算交给着色器。
  • 点云是没有连接关系的三维点集,来自激光雷达(r=cΔt/2r = c\Delta t/2)、深度相机(反投影)、多视图重建等。它只记录传感器看得见的表面,点密度约与 cos⁡α/r2\cos\alpha / r^2 成正比,点与点之间是空的。OpenGL 用 GL_POINTS 一次画出整个点云,深度测试与三角形完全相同。

练习题

  1. 概念题:同学 A 和同学 B 在各自的电脑上运行同一段 OpenGL 程序,打印出的 GL_VERSION 字符串不同,画面也有细微差别。这正常吗?请从“规范与实现”的角度解释。如果要让程序在两台机器上都能运行,应该注意什么?

  2. 计算题:一台 2560×1440 分辨率、144 Hz 刷新率的显示器,假设平均 overdraw 为 2.5,每秒要处理多少个片元?如果用一个主频 3 GHz、每个时钟周期完成一次运算的 CPU 核心来做,平均每个片元能分到几个时钟周期?由此说明为什么需要 GPU。

  3. 计算题:某个像素依次收到 4 个片元,深度和颜色分别为 (5.0, 红)(5.0,\ \text{红})、(3.0, 绿)(3.0,\ \text{绿})、(4.0, 蓝)(4.0,\ \text{蓝})、(3.0, 白)(3.0,\ \text{白})。深度缓冲初值为 +∞+\infty。 (a) 用 GL_LESS(z<z < zbuf)时,写出每处理一个片元后深度缓冲和颜色缓冲的值,最终颜色是什么? (b) 改用 GL_LEQUAL(z≤z \le zbuf)呢? (c) 用定理 2 和推论 2 解释 (a)、(b) 结果的差别。

  4. 证明题:画家算法要求把三角形按“远近”排序。请构造(画示意图并说明)三个三角形 A、B、C,使得 A 挡住 B 的一部分、B 挡住 C 的一部分、C 又挡住 A 的一部分。证明对这三个三角形,不存在任何一种绘制顺序能让画家算法得到正确结果,而 Z-buffer 可以。

  5. 推导题:一条线段的两个端点深度分别为 w0=2w_0 = 2、w1=6w_1 = 6。在屏幕上取线段投影的中点(屏幕空间重心坐标 λ0=λ1=0.5\lambda_0 = \lambda_1 = 0.5)。 (a) 用定理 3 求这一点的真实深度 ww,以及它在三维线段上的重心坐标 (b0,b1)(b_0, b_1); (b) 如果错误地直接在屏幕上线性插值深度,会得到多少? (c) 解释为什么屏幕上的中点不是三维线段的中点。(提示:远处的部分在屏幕上被“压缩”了。)

  6. 代码阅读题:下面的代码执行完后,纹理单元 0 和纹理单元 1 上分别绑定的是哪个纹理?glTexParameteri 修改的是哪个纹理的参数?

    tex_a, tex_b = glGenTextures(2)
    glActiveTexture(GL_TEXTURE1)
    glBindTexture(GL_TEXTURE_2D, tex_a)
    glActiveTexture(GL_TEXTURE0)
    glBindTexture(GL_TEXTURE_2D, tex_b)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    
  7. 思考题:有人这样测量绘制耗时:

    t0 = time.perf_counter()
    glDrawArrays(GL_TRIANGLES, 0, 3_000_000)
    print(time.perf_counter() - t0)
    

    结果显示只用了几微秒,于是得出“画 100 万个三角形只要几微秒”的结论。这个结论错在哪里?如果在 glDrawArrays 后面加一句 glFinish() 再计时,测到的又是什么?这样做有什么副作用?

  8. 估算题:假设通过 PyOpenGL 调用一次 GL 函数平均需要 2 微秒。用立即模式画一个 5 万个三角形的模型,每个顶点调用 glColor3f、glNormal3f、glVertex3f 各一次(三角形之间不共享顶点)。仅函数调用就要多少时间?帧率最高是多少?如果改用 VBO,每帧只需要大约 10 次 GL 调用,又是多少?

  9. 编程题:修改附录代码。 (a) 把场景列表反过来,即 rasterize(SCENE[::-1]),观察深度测试的通过和失败次数怎么变化,最终图像有没有变化。尝试所有 6 种绘制顺序,哪种顺序深度测试通过的次数最少?为什么最少的次数不等于被覆盖的像素数 7489? (b) 把分辨率改为 320×240,两种方法的测试次数各变成原来的多少倍?三角形的个数没变,为什么代价还是变大了?

  10. 估算题:KITTI 数据集的激光雷达每帧约有 12 万个点,每个点存 x,y,zx, y, z 和反射强度 4 个 float32。 (a) 一帧数据有多大? (b) 激光雷达每秒扫描 10 帧,录制 1 分钟有多少数据? (c) 如果加载时不小心用了 NumPy 默认的 float64,一帧变成多大?上传给 OpenGL 时如果仍声明为 GL_FLOAT,会发生什么?

  11. 推导题:一台深度相机的参数为 fx=fy=500f_x = f_y = 500,(ox,oy)=(320,240)(o_x, o_y) = (320, 240)(像素坐标按本系列约定,原点在左下角)。像素 (420,140)(420, 140) 测得深度 w=2.0w = 2.0。 (a) 用第八节的反投影公式求这个点在相机坐标系中的三维坐标; (b) 再用第二节的投影公式 x=fxqx/w+oxx = f_x q_x / w + o_x 把它投影回去,验证得到原来的像素; (c) 为什么只有一张普通照片(没有深度)时,不能把像素还原成三维点?

  12. 编程题:修改附录代码的点云部分。 (a) 把激光雷达放到相机的位置,即 SENSOR = np.array([0.0, 0.0, 0.0]),SENSOR_TARGET 不变。重新扫描并在相机视角下渲染,背景三角形上的“阴影”空洞还在吗?为什么? (b) 把激光雷达放回原处,把水平采样数 N_AZ 从 160 改为 320。点数和平均点间距(nn_spacing 的平均值)怎么变化?“阴影”空洞能不能靠提高分辨率消除?

参考答案

  1. 完全正常。OpenGL 是一份规范,GL_VERSION 报告的是当前驱动实现支持的版本,以及厂商和驱动版本信息。不同厂商(甚至同一厂商的不同驱动版本)各自实现规范,规范允许的范围内(如浮点精度、某些边界情况的处理)结果可以有细微差别,偶尔也会遇到驱动 bug。要在两台机器上都能运行,应当:只使用两台机器都支持的版本里的功能(本系列只依赖 2.0);用到扩展时先查询是否支持;不依赖规范没有明确规定的行为。

  2. 每秒片元数 =2560×1440×144×2.5=3,686,400×144×2.5=530,841,600×2.5≈1.33×109= 2560 \times 1440 \times 144 \times 2.5 = 3{,}686{,}400 \times 144 \times 2.5 = 530{,}841{,}600 \times 2.5 \approx 1.33 \times 10^9。每个片元能分到的时钟周期 =3×109/(1.33×109)≈2.26= 3 \times 10^9 / (1.33 \times 10^9) \approx 2.26 个。而一个片元的着色计算(插值、纹理采样、光照)通常需要几十到几百次运算,还要加上顶点处理和光栅化本身的开销,单个 CPU 核心远远不够。片元之间互相独立,正适合 GPU 的大规模并行。

  3. (a) GL_LESS:

    • 片元 1:5.0<∞5.0 < \infty,通过,zbuf = 5.0,颜色 = 红;
    • 片元 2:3.0<5.03.0 < 5.0,通过,zbuf = 3.0,颜色 = 绿;
    • 片元 3:4.0<3.04.0 < 3.0 不成立,丢弃,仍为 3.0、绿;
    • 片元 4:3.0<3.03.0 < 3.0 不成立,丢弃,仍为 3.0、绿。 最终颜色为绿。 (b) GL_LEQUAL:前三步相同;片元 4:3.0≤3.03.0 \le 3.0 成立,通过,颜色 = 白。最终颜色为白。 (c) 最小深度 3.0 在片元 2 和片元 4 上同时取到。由定理 2,GL_LESS 留下最先达到最小值的片元(片元 2,绿);而 GL_LEQUAL 允许“相等也覆盖”,留下最后一个达到最小值的片元(片元 4,白)。这正是推论 2:深度相等时,结果取决于比较方式和到达顺序。两种比较方式的最终深度都是 3.0,这部分与顺序无关。
  4. 构造:三个长条形三角形围成一个“风车”。从上往下看,A 的右端压在 B 的左端上面,B 的下端压在 C 的上端上面,C 的左端又压在 A 的下端上面(三个三角形都略微倾斜,所以每个三角形的一端高、一端低)。证明:画家算法要求被遮挡的先画。A 挡 B 部分区域,要求 B 在 A 之前画;同理 C 在 B 之前、A 在 C 之前。于是 B 在 A 之前、A 在 C 之前、C 在 B 之前,推出 B 在 B 之前,矛盾。所以不存在满足要求的顺序。Z-buffer 在每个像素上独立比较深度,在 A、B 重叠的像素上,A 更近,就留下 A;在 B、C 重叠的像素上留下 B;在 C、A 重叠的像素上留下 C。每个像素都有确定的最近片元,由定理 2 和推论 1,按任意顺序绘制都能得到正确结果。(互相穿插的三角形是同样的道理,附录代码里中间两个三角形就是例子。)

  5. (a) ∑jλj/wj=0.5/2+0.5/6=0.25+0.0833=0.3333\sum_j \lambda_j / w_j = 0.5/2 + 0.5/6 = 0.25 + 0.0833 = 0.3333,所以 w=1/0.3333=3w = 1/0.3333 = 3。b0=0.25/0.3333=0.75b_0 = 0.25 / 0.3333 = 0.75,b1=0.0833/0.3333=0.25b_1 = 0.0833/0.3333 = 0.25。验证:b0w0+b1w1=0.75×2+0.25×6=3b_0 w_0 + b_1 w_1 = 0.75 \times 2 + 0.25 \times 6 = 3,一致。 (b) 直接线性插值:0.5×2+0.5×6=40.5 \times 2 + 0.5 \times 6 = 4,比真实值大了 1。 (c) 透视投影中,屏幕坐标与 1/w1/w 成正比:越远的东西在屏幕上越小。线段远端(w=6w = 6)那一半在屏幕上被压缩得更厉害,所以屏幕上的中点对应三维线段上更靠近近端的点(b0=0.75b_0 = 0.75,离近端只有四分之一的长度)。

  6. 纹理单元 1 上绑定的是 tex_a,纹理单元 0 上绑定的是 tex_b。glTexParameteri 作用在“当前活动纹理单元的 GL_TEXTURE_2D 绑定点上的纹理”。最后一次 glActiveTexture 选的是单元 0,单元 0 上绑的是 tex_b,所以修改的是 tex_b 的参数。这是 bind-to-edit 模式的典型例子,glTexParameteri 的参数里并没有出现纹理的名字。

  7. 错在混淆了“命令入队”和“命令执行”。OpenGL 是异步执行的,glDrawArrays 只是把命令放进命令缓冲区就返回了,测到的几微秒只是 Python 调用和入队的时间,GPU 可能还没开始画。加上 glFinish() 后,CPU 会等到 GPU 执行完所有排队的命令,测到的是“入队 + GPU 执行完(包括之前排队的其他命令)”的总时间,比较接近真实的 GPU 耗时,但仍然不精确。副作用是打断了 CPU 和 GPU 的流水线:CPU 等 GPU 时,GPU 已经没有新命令可做,两边轮流空闲,整体帧率会下降。所以 glFinish() 只适合临时调试,精确测量应该用 GPU 计时查询(第 34 篇)。

  8. 顶点数 =50000×3=150000= 50000 \times 3 = 150000,调用次数 =150000×3+2=450002= 150000 \times 3 + 2 = 450002(加上 glBegin 和 glEnd),约 450000×2 μs=0.9450000 \times 2\ \mu\text{s} = 0.9 秒。每帧至少 0.9 秒,帧率最高约 1.1 fps。改用 VBO 后,每帧约 10×2 μs=20 μs10 \times 2\ \mu\text{s} = 20\ \mu\text{s},函数调用开销可以忽略不计,帧率只受 GPU 的处理能力限制。

  9. (a) 反过来以后,包围盒测试次数(25551)和片元数(11456)都不变,深度测试通过 7951 次、失败 3505 次(原来是 10994 次和 462 次),最终图像与原来完全相同(最大差异为 0)。原顺序是先画最远的大三角形,后画的近处三角形不断覆盖它,所以通过次数多;反过来以后,近处的先画,远处的片元大多在深度测试中被丢弃。这验证了推论 1。6 种顺序的通过次数分别为:(0,1,2)(0,1,2) 10994、(0,2,1)(0,2,1) 11087、(1,0,2)(1,0,2) 8671、(1,2,0)(1,2,0) 7858、(2,0,1)(2,0,1) 9443、(2,1,0)(2,1,0) 7951。按 (1,2,0)(1, 2, 0) 的顺序(先画两个前面的三角形,最后画背景)通过次数最少,为 7858 次。理论下界是被覆盖的像素数 7489,因为每个被覆盖像素收到的第一个片元一定通过测试,之后的片元只有更近时才通过。背景三角形 0 在所有重叠处都比三角形 1、2 远,最后画它就不会多出通过次数;多出来的 7858−7489=3697858 - 7489 = 369 次来自两个互相穿插的三角形 1 和 2。它们重叠 831 个像素(正是收到 3 个片元的那 831 个像素),其中三角形 1 更近的有 462 个,三角形 2 更近的有 369 个。先画 1 再画 2,2 会在那 369 个像素上再通过一次;先画 2 再画 1,1 会在 462 个像素上再通过一次。互相穿插的三角形不存在处处“从前往后”的顺序,所以最少也要多出 min⁡(462,369)=369\min(462, 369) = 369 次,(1,2,0)(1,2,0) 和 (2,1,0)(2,1,0) 的通过次数恰好是 7489+3697489 + 369 和 7489+4627489 + 462。实际渲染中,“大致从前往后画”是常用的优化,因为 GPU 可以在运行片元着色器之前就做深度测试,提前丢弃被挡住的片元(early-z,第 13 篇)。 (b) 320×240 时:光线投射的求交次数为 320×240×3=230400320 \times 240 \times 3 = 230400,正好是原来的 4 倍;光栅化的包围盒测试为 100474 次,是原来的 3.93 倍,片元数为 46085 个,是原来的 4.02 倍。(不是精确的 4 倍,是因为三角形边缘的像素是否被覆盖随分辨率变化。)两种方法的代价都和像素数成正比:光线投射的每条射线对应一个像素;光栅化的内层循环遍历的是三角形覆盖的像素。分辨率翻倍,像素数变为 4 倍,代价也大约变为 4 倍。只有光栅化的“每个三角形一次”的投影、求面积等准备工作与分辨率无关。

  10. (a) 120000×4×4=1,920,000120000 \times 4 \times 4 = 1{,}920{,}000 字节,约 1.92 MB。 (b) 1 分钟 =600= 600 帧,600×1.92 MB≈1.15 GB600 \times 1.92\ \text{MB} \approx 1.15\ \text{GB}。这就是为什么大规模点云必须考虑下采样和细节层次(第 32 篇)。 (c) float64 每个数 8 字节,一帧变成 3.84 MB,内存翻倍。更严重的是,如果告诉 OpenGL 数据是 GL_FLOAT(4 字节),它会把每个 8 字节的双精度数拆成两个 4 字节的单精度数来解读,得到的坐标全是错的,画面上会是一团乱七八糟的点。所以传给 OpenGL 的数组一定要先 astype(np.float32)(第 05 篇)。

  11. (a) qx=(420−320)×2.0/500=0.4q_x = (420 - 320) \times 2.0 / 500 = 0.4,qy=(140−240)×2.0/500=−0.4q_y = (140 - 240) \times 2.0 / 500 = -0.4,qz=−w=−2.0q_z = -w = -2.0,即 (0.4, −0.4, −2.0)(0.4,\ -0.4,\ -2.0)。 (b) x=500×0.4/2.0+320=420x = 500 \times 0.4 / 2.0 + 320 = 420,y=500×(−0.4)/2.0+240=140y = 500 \times (-0.4) / 2.0 + 240 = 140,回到了原来的像素。 (c) 投影公式 x=fxqx/w+oxx = f_x q_x / w + o_x 中,把 q\mathbf{q} 换成 kqk\mathbf{q}(k>0k > 0),qxq_x 和 ww 同时乘以 kk,xx 不变。也就是说,同一条视线上的所有点都投影到同一个像素,投影丢掉了深度。只有照片时,只知道点在哪条视线上,不知道在视线上的哪个位置。深度相机和激光雷达正是把这个丢掉的量测了出来。

  12. (a) 运行结果:得到 3238 个点,在相机视角下覆盖 7318 个像素;三角形覆盖、但没有被点覆盖的像素从原来的 2186 个降到 313 个,阴影空洞消失了,剩下的只是点与点之间的缝隙和边缘的零星像素。原因:激光雷达和相机在同一位置,激光束就是相机的视线。相机看得见的每一个表面点,都在某条从同一位置出发的射线上最先被碰到,所以激光雷达也看得见(只要采样够密)。阴影的根源是两个视点不同,这与第二节“两种方法看到同样的结果”是同一个道理。 (b) 运行结果:点数从 2919 变为 5859,几乎正好翻倍;平均点间距从 0.0706 降到 0.0390,约为原来的 0.55 倍(水平方向的点间距减半,竖直方向不变,最近邻距离取两者中较小的那个,所以接近减半)。被点覆盖的像素从 5385 增加到 6147,三角形覆盖、却没有点的像素从 2186 降到 1486。减少的那 700 个是点与点之间的缝隙被填上了;剩下的 1486 个里绝大部分是阴影:(a) 中视点相同、点的密度相近时,这个数字只有 313。阴影区域根本没有激光能照到,采样再密也没有用。要填补阴影,只能从别的位置再扫描一次。

附录:完整代码

下面的代码用纯 Python 把两种渲染方法各实现了一遍,并实测了一次 ctypes 调用的开销。它只依赖 NumPy 和 matplotlib,不需要 OpenGL。保存为 example.py 后运行 python example.py 即可。光线投射部分是纯 Python 循环,要运行几秒钟。

"""OpenGL 2.0 图形编程自学教程 00:导读 —— 附录代码

本篇还不调用 OpenGL,只用 Python 和 numpy 把“渲染”这件事的骨架走一遍:
  1. 搭一个小场景:3 个带顶点颜色的三角形,其中两个互相穿插;
  2. 用光栅化渲染一遍:外层循环三角形,内层循环它覆盖的像素,用 Z-buffer 解决遮挡;
  3. 用光线投射再渲染一遍:外层循环像素,内层循环所有三角形;
  4. 比较两张图,统计两种算法各做了多少次基本测试;
  5. 模拟一台激光雷达扫描这个场景,得到一份点云,存成 PLY 文件,再把点云画到屏幕上;
  6. 实测 Python 通过 ctypes 调用一次 C 函数的开销,估算“每个顶点调用一次 OpenGL”有多慢。
运行:python example.py
"""
import ctypes
import ctypes.util
import math
import pathlib
import sys
import time

import numpy as np
import matplotlib.pyplot as plt

W, H = 160, 120                  # 图像宽、高(像素)
FOV_Y = math.radians(60.0)       # 竖直方向的视场角
CLEAR_COLOR = (0.11, 0.13, 0.17) # 背景色,相当于 glClearColor

# ---------------- 1. 场景 ----------------
# 坐标直接给在“观察空间”里:眼睛在原点,看向 -z 方向,y 轴向上(OpenGL 的约定)。
# 每个三角形:3 个顶点的位置 + 3 个顶点的颜色(RGB,0~1)。
SCENE = [
    # 背后一块斜放的大三角形:底边近、顶点远
    (np.array([[-4.0, -2.8, -6.0], [4.0, -2.8, -6.0], [0.0, 3.2, -9.0]]),
     np.array([[0.16, 0.47, 0.84], [0.11, 0.69, 0.48], [0.38, 0.31, 0.84]])),
    # 前面两块互相穿插的三角形:从上往下看像一个“X”
    (np.array([[-2.5, -1.5, -4.0], [1.5, -1.0, -6.5], [-0.5, 2.0, -5.0]]),
     np.array([[0.92, 0.41, 0.20], [0.98, 0.75, 0.25], [0.92, 0.41, 0.20]])),
    (np.array([[-1.0, -1.8, -6.5], [2.5, -0.5, -4.0], [0.5, 2.2, -5.5]]),
     np.array([[0.95, 0.95, 0.93], [0.60, 0.58, 0.55], [0.95, 0.95, 0.93]])),
]


def project(v):
    """观察空间中的点 -> 窗口坐标 (x, y) 和深度 w(到眼睛的 -z 距离)。

    第 10、12 篇会把这一步拆成“投影矩阵 + 透视除法 + 视口变换”三步。
    """
    w = -v[2]
    t = math.tan(FOV_Y / 2)
    ndc_x = v[0] / (w * t * W / H)
    ndc_y = v[1] / (w * t)
    return (ndc_x + 1) / 2 * W, (ndc_y + 1) / 2 * H, w


def pixel_ray(x, y):
    """穿过像素 (x, y) 中心的视线方向。z 分量取 -1,于是参数 t 恰好等于深度。"""
    t = math.tan(FOV_Y / 2)
    ndc_x = 2 * (x + 0.5) / W - 1
    ndc_y = 2 * (y + 0.5) / H - 1
    return np.array([ndc_x * t * W / H, ndc_y * t, -1.0])


def edge(ax, ay, bx, by, px, py):
    """边函数:点 p 在有向边 a->b 左侧时为正,右侧为负,在边上为 0(第 07 篇详解)。"""
    return (bx - ax) * (py - ay) - (by - ay) * (px - ax)


# ---------------- 2. 光栅化:对象顺序 ----------------
def rasterize(scene):
    color = np.empty((H, W, 3))
    color[:] = CLEAR_COLOR
    zbuf = np.full((H, W), np.inf)          # 深度缓冲,初值为“无穷远”
    overdraw = np.zeros((H, W), dtype=int)  # 每个像素收到了几个片元
    stats = dict(inside_tests=0, fragments=0, depth_pass=0, depth_fail=0)

    for pos, col in scene:                                  # 外层:三角形
        (x0, y0, w0), (x1, y1, w1), (x2, y2, w2) = [project(v) for v in pos]
        area = edge(x0, y0, x1, y1, x2, y2)
        if area == 0:
            continue
        # 包围盒,截到屏幕范围内
        xmin = max(int(math.floor(min(x0, x1, x2))), 0)
        xmax = min(int(math.ceil(max(x0, x1, x2))), W - 1)
        ymin = max(int(math.floor(min(y0, y1, y2))), 0)
        ymax = min(int(math.ceil(max(y0, y1, y2))), H - 1)
        for y in range(ymin, ymax + 1):                     # 内层:包围盒里的像素
            for x in range(xmin, xmax + 1):
                px, py = x + 0.5, y + 0.5                   # 在像素中心采样
                stats["inside_tests"] += 1
                # 屏幕空间的重心坐标
                l0 = edge(x1, y1, x2, y2, px, py) / area
                l1 = edge(x2, y2, x0, y0, px, py) / area
                l2 = edge(x0, y0, x1, y1, px, py) / area
                if l0 < 0 or l1 < 0 or l2 < 0:
                    continue                                # 像素中心不在三角形内
                stats["fragments"] += 1
                overdraw[y, x] += 1
                # 1/w 在屏幕上是线性的,所以先插值 1/w 再取倒数得到深度(第 07 篇证明)
                inv_w = l0 / w0 + l1 / w1 + l2 / w2
                depth = 1.0 / inv_w
                if depth < zbuf[y, x]:                      # 深度测试,相当于 GL_LESS
                    stats["depth_pass"] += 1
                    zbuf[y, x] = depth
                    # 透视校正的重心坐标,用它插值顶点颜色
                    b = np.array([l0 / w0, l1 / w1, l2 / w2]) / inv_w
                    color[y, x] = b @ col
                else:
                    stats["depth_fail"] += 1
    return color, zbuf, overdraw, stats


# ---------------- 3. 光线投射:图像顺序 ----------------
def ray_triangle(o, d, a, b, c):
    """从 o 出发、方向为 d 的射线与三角形 abc 求交(Möller–Trumbore,第 27 篇推导)。

    命中时返回 (t, u, v),交点为 o + t d = (1-u-v)a + u b + v c;没命中返回 None。
    """
    e1, e2 = b - a, c - a
    p = np.cross(d, e2)
    det = e1 @ p
    if abs(det) < 1e-12:
        return None
    s = o - a
    u = (s @ p) / det
    q = np.cross(s, e1)
    v = (d @ q) / det
    if u < 0 or v < 0 or u + v > 1:
        return None
    t = (e2 @ q) / det
    return (t, u, v) if t > 0 else None


def raytrace(scene):
    color = np.empty((H, W, 3))
    color[:] = CLEAR_COLOR
    zbuf = np.full((H, W), np.inf)
    stats = dict(ray_tests=0, hits=0)
    eye = np.zeros(3)

    for y in range(H):                                      # 外层:像素
        for x in range(W):
            d = pixel_ray(x, y)
            for pos, col in scene:                          # 内层:所有三角形
                stats["ray_tests"] += 1
                hit = ray_triangle(eye, d, *pos)
                if hit is None:
                    continue
                stats["hits"] += 1
                t, u, v = hit
                if t < zbuf[y, x]:                          # 只保留最近的交点
                    zbuf[y, x] = t
                    color[y, x] = np.array([1 - u - v, u, v]) @ col
    return color, zbuf, stats


# ---------------- 4. 点云:模拟一次激光雷达扫描 ----------------
SENSOR = np.array([4.0, 1.5, -2.0])         # 激光雷达的位置(观察空间坐标),在相机右前方
SENSOR_TARGET = np.array([0.0, 0.0, -6.0])  # 激光雷达朝向这一点
N_EL, N_AZ = 48, 160                        # 竖直 48 线,每线水平 160 个采样
EL_FOV = math.radians(50.0)                 # 竖直视场
AZ_FOV = math.radians(90.0)                 # 水平视场


def lidar_scan(scene):
    """从 SENSOR 发出 N_EL × N_AZ 束激光,每束记录最近的交点,得到点云。

    返回 (N, 4) 数组,每行一个点:x, y, z(观察空间坐标)和到传感器的距离 r。
    """
    f = SENSOR_TARGET - SENSOR
    f /= np.linalg.norm(f)                                  # 正前方
    r = np.cross(f, [0.0, 1.0, 0.0])
    r /= np.linalg.norm(r)                                  # 右方
    u = np.cross(r, f)                                      # 上方(第 11 篇的 LookAt 也这样构造)
    points = []
    for el in np.linspace(-EL_FOV / 2, EL_FOV / 2, N_EL):
        for az in np.linspace(-AZ_FOV / 2, AZ_FOV / 2, N_AZ):
            # 球坐标 -> 方向向量;f、r、u 两两正交且长度为 1,所以 |d| = 1
            d = math.cos(el) * (math.cos(az) * f + math.sin(az) * r) + math.sin(el) * u
            best = math.inf
            for pos, _ in scene:
                hit = ray_triangle(SENSOR, d, *pos)
                if hit is not None and hit[0] < best:
                    best = hit[0]
            if best < math.inf:                             # |d| = 1,所以 t 就是距离
                points.append([*(SENSOR + best * d), best])
    return np.array(points)


def save_ply(path, points):
    """把点云存成 ASCII 格式的 PLY 文件(第 29 篇详细介绍 PLY 等点云格式)。"""
    with open(path, "w", encoding="ascii") as f:
        f.write("ply\nformat ascii 1.0\n")
        f.write(f"element vertex {len(points)}\n")
        for name in ("x", "y", "z", "range"):
            f.write(f"property float {name}\n")
        f.write("end_header\n")
        for p in points:
            f.write(" ".join(f"{v:.4f}" for v in p) + "\n")


def nn_spacing(xyz, chunk=500):
    """每个点到最近的另一个点的距离(暴力计算,分块以免占用太多内存)。"""
    out = np.empty(len(xyz))
    for i in range(0, len(xyz), chunk):
        d = np.linalg.norm(xyz[i:i + chunk, None, :] - xyz[None, :, :], axis=2)
        d[np.arange(len(d)), np.arange(i, i + len(d))] = np.inf   # 排除自己
        out[i:i + chunk] = d.min(axis=1)
    return out


def render_points(points, size=2):
    """在相机(原点)视角下画点云:每个点占 size×size 个像素,用 Z-buffer 处理遮挡。

    相当于 OpenGL 里 glPointSize(size) 之后调用 glDrawArrays(GL_POINTS, 0, N)。
    颜色按到传感器的距离从近(紫)到远(黄)。
    """
    color = np.empty((H, W, 3))
    color[:] = CLEAR_COLOR
    zbuf = np.full((H, W), np.inf)
    rng = points[:, 3]
    cols = plt.get_cmap("viridis")((rng - rng.min()) / (rng.max() - rng.min()))[:, :3]
    for p, c in zip(points[:, :3], cols):
        x, y, w = project(p)
        # 以 (x, y) 为中心、边长 size 的正方形,覆盖中心落在其中的像素
        x0 = math.ceil(x - size / 2 - 0.5)
        y0 = math.ceil(y - size / 2 - 0.5)
        for py in range(max(y0, 0), min(y0 + size, H)):
            for px in range(max(x0, 0), min(x0 + size, W)):
                if w < zbuf[py, px]:
                    zbuf[py, px] = w
                    color[py, px] = c
    return color


# ---------------- 5. 一次 ctypes 调用要多久 ----------------
def ctypes_call_cost(n=200_000):
    """调用 n 次 C 标准库的 abs(),返回平均每次调用的秒数。

    PyOpenGL 也是通过 ctypes 调用驱动里的 C 函数,所以这是一次 GL 调用开销的下限。
    """
    if sys.platform == "win32":
        libc = ctypes.cdll.msvcrt
    else:
        libc = ctypes.CDLL(ctypes.util.find_library("c"))
    f = libc.abs
    f.argtypes = [ctypes.c_int]
    f.restype = ctypes.c_int
    best = math.inf
    for _ in range(3):                                      # 测 3 次取最快的一次
        t0 = time.perf_counter()
        for i in range(n):
            f(i)
        best = min(best, (time.perf_counter() - t0) / n)
    return best


if __name__ == "__main__":
    print(f"图像 {W}×{H} = {W * H} 个像素,场景 {len(SCENE)} 个三角形\n")

    t0 = time.perf_counter()
    img_r, z_r, overdraw, st_r = rasterize(SCENE)
    t_r = time.perf_counter() - t0
    covered = int((overdraw > 0).sum())
    print("【光栅化】外层循环三角形")
    print(f"    包围盒内的“像素在不在三角形里”测试:{st_r['inside_tests']} 次")
    print(f"    产生片元 {st_r['fragments']} 个,覆盖 {covered} 个像素,"
          f"平均每个被覆盖像素 {st_r['fragments'] / covered:.2f} 个片元")
    print(f"    深度测试:通过 {st_r['depth_pass']} 次,被挡住 {st_r['depth_fail']} 次")
    print(f"    收到 1/2/3 个片元的像素数:{[int((overdraw == k).sum()) for k in (1, 2, 3)]}")
    print(f"    用时 {t_r:.2f} s(纯 Python 循环,只看相对大小)\n")

    t0 = time.perf_counter()
    img_t, z_t, st_t = raytrace(SCENE)
    t_t = time.perf_counter() - t0
    print("【光线投射】外层循环像素")
    print(f"    射线-三角形求交测试:{st_t['ray_tests']} 次 = {W * H} 像素 × {len(SCENE)} 三角形")
    print(f"    命中 {st_t['hits']} 次")
    print(f"    用时 {t_t:.2f} s\n")

    hit_r, hit_t = np.isfinite(z_r), np.isfinite(z_t)
    both = hit_r & hit_t
    print("【两种方法的结果对比】")
    print(f"    覆盖情况不一致的像素:{int((hit_r != hit_t).sum())} 个")
    print(f"    深度最大差异:{np.abs(z_r[both] - z_t[both]).max():.2e}")
    print(f"    颜色最大差异:{np.abs(img_r - img_t).max():.2e}\n")

    cloud = lidar_scan(SCENE)
    ply_path = pathlib.Path(__file__).with_name("scan.ply")
    save_ply(ply_path, cloud)
    print("【点云:模拟激光雷达扫描】")
    print(f"    发出 {N_EL}×{N_AZ} = {N_EL * N_AZ} 束激光,命中 {len(cloud)} 束,得到 {len(cloud)} 个点")
    print(f"    距离范围 {cloud[:, 3].min():.2f} ~ {cloud[:, 3].max():.2f},已保存到 {ply_path.name}")
    spacing = nn_spacing(cloud[:, :3])
    order = np.argsort(cloud[:, 3])
    k = len(cloud) // 5
    for name, idx in (("最近的 20%", order[:k]), ("最远的 20%", order[-k:])):
        print(f"    {name}的点:平均距离 {cloud[idx, 3].mean():.2f},"
              f"平均点间距 {spacing[idx].mean():.4f},"
              f"点间距/距离 {(spacing[idx] / cloud[idx, 3]).mean():.4f}")
    img_p = render_points(cloud)
    print(f"    在相机视角下画出点云,被点覆盖的像素 {int((img_p != CLEAR_COLOR).any(axis=2).sum())} 个\n")

    per_call = ctypes_call_cost()
    n_tri = 100_000
    n_calls = n_tri * 3 * 2          # 每个顶点一次 glColor3f + 一次 glVertex3f
    print("【Python 调用 C 函数的开销】")
    print(f"    一次 ctypes 调用约 {per_call * 1e9:.0f} ns")
    print(f"    {n_tri} 个三角形逐顶点提交需要 {n_calls} 次调用,"
          f"至少 {n_calls * per_call * 1e3:.0f} ms/帧,"
          f"帧率不超过 {1 / (n_calls * per_call):.0f} fps")

    # ---------------- 6. 画图 ----------------
    plt.rcParams["font.sans-serif"] = ["Microsoft YaHei", "SimHei", "DejaVu Sans"]
    fig, axes = plt.subplots(2, 2, figsize=(9, 7))
    shown = [(img_r, "光栅化"), (img_t, "光线投射"),
             (overdraw, "每个像素收到的片元数"), (img_p, "激光雷达点云")]
    for ax, (im, title) in zip(axes.flat, shown):
        ax.imshow(im, origin="lower", interpolation="nearest")   # 原点在左下角,同 OpenGL
        ax.set_title(title)
        ax.set_xticks([])
        ax.set_yticks([])
    plt.tight_layout()
    plt.show()

运行输出(用时和 ctypes 开销取决于机器,其余数字是确定的):

图像 160×120 = 19200 个像素,场景 3 个三角形

【光栅化】外层循环三角形
    包围盒内的“像素在不在三角形里”测试:25551 次
    产生片元 11456 个,覆盖 7489 个像素,平均每个被覆盖像素 1.53 个片元
    深度测试:通过 10994 次,被挡住 462 次
    收到 1/2/3 个片元的像素数:[4353, 2305, 831]
    用时 0.09 s(纯 Python 循环,只看相对大小)

【光线投射】外层循环像素
    射线-三角形求交测试:57600 次 = 19200 像素 × 3 三角形
    命中 11456 次
    用时 2.93 s

【两种方法的结果对比】
    覆盖情况不一致的像素:0 个
    深度最大差异:2.66e-15
    颜色最大差异:3.33e-16

【点云:模拟激光雷达扫描】
    发出 48×160 = 7680 束激光,命中 2919 束,得到 2919 个点
    距离范围 3.23 ~ 9.57,已保存到 scan.ply
    最近的 20%的点:平均距离 4.91,平均点间距 0.0674,点间距/距离 0.0142
    最远的 20%的点:平均距离 7.49,平均点间距 0.0900,点间距/距离 0.0120
    在相机视角下画出点云,被点覆盖的像素 5385 个

【Python 调用 C 函数的开销】
    一次 ctypes 调用约 359 ns
    100000 个三角形逐顶点提交需要 600000 次调用,至少 215 ms/帧,帧率不超过 5 fps

运行后文件夹里会多出一个 scan.ply,可以用 MeshLab 或 CloudCompare 打开。

代码要点:

  • project 和 pixel_ray 是同一个针孔相机的两种用法:前者把三维点投到屏幕上,后者从屏幕上的点反推出三维的视线。第 12 篇会把它们写成矩阵。
  • rasterize 的两层循环正是图 1 左边的伪代码,raytrace 的两层循环是图 1 右边的伪代码。
  • rasterize 里的 inv_w 和 b 就是定理 3 的两个公式。光线投射命中 11456 次,恰好等于光栅化产生的片元数,这也是定理 3 推论 1(覆盖判断等价)的体现。
  • lidar_scan 复用了光线投射的 ray_triangle,只是射线的起点从相机换成了激光雷达,方向从“穿过像素”换成了“按仰角和方位角均匀排列”。激光雷达就是物理世界的光线投射。
  • render_points 是 GL_POINTS 的 Python 版本:每个点投影后覆盖 size×size 个像素,深度测试与 rasterize 相同。第 30 篇会用真正的 OpenGL 画几十万个点。
  • 第 07 篇会把 rasterize 改写成 NumPy 向量化版本,速度快几百倍,并补上裁剪、填充规则等细节。