Github TinyRenderer渲染器课程实践记录: 三角形光栅化
Abstract
上一节：Bresemham直线绘制
Reference :
经过第一部分学习Bresenham直线绘制后，可以用三条线来画一个三角形：
void triangle(Vec2i t0, Vec2i t1, Vec2i t2, TGAImage& image, TGAColor color) {
line(t0, t1, image, color);
line(t1, t2, image, color);
line(t2, t0, image, color);
}
而本部分关注如何在三角形内部填充颜色。
扫描线算法 Line Sweeping
扫描线算法基本思想为，根据顶点的 \(y\) 值从下到上地，逐步从三角形左边界到右边界绘制直线，直到到达三角形最高点绘制完毕：
算法步骤采用三线性插值法，绘制要分别考虑被顶点 \(v1\) 所在水平线分割的上下部分，因为插值与三个向量代表的微分方向有关：
-
\(\vec{v_{0}v_{2}}\)
-
\(\vec{v_{0}v_{1}}\)
-
\(\vec{v_{1}v_{2}}\)
每绘制一条直线时，需从它左端点一个个画点画到其右端点，因此每条直线的左、右端点的采样是同时进行的。
在上图例中，左边界 \(\vec{v_{0}v_{2}}\) 上端点的采样范围是 \([v_{0}, v_{2}]\) ，插值比例系数来自于 \(v_0\) 和 \(v_2\) 之间竖直方向位移 \(h1\) 的插值结果。
同理，右边界分为 \(\vec{v_{0}v_{1}}\) 与 \(\vec{v_{1}v_{2}}\) 两部分。第一部分和第二部分端点的采样范围分别是 \([v_{0}, v_{1}]\) 与 \([v_{1}, v_{2}]\)，插值比例系数分别来自于竖直方向位移 \(h2\) 与 \(h3\) 以及单条扫描线左端到右端 \(x\) 的插值结果：
void triangle(Vec2i t0, Vec2i t1, Vec2i t2, TGAImage& image, TGAColor color) {
if(t0.y > t1.y) std::swap(t0, t1);
if(t0.y > t2.y) std::swap(t0, t2);
if(t1.y > t2.y) std::swap(t1, t2);
int h1 = t2.y - t0.y;
int h2 = t1.y - t0.y;
int h3 = t2.y - t1.y;
for(int i = 0; i <= h1; i++) {
bool halfPart = i > h2 || t1.y == t0.y;
// segmentHeight的值一定不为0，不会发生除0错误
int segmentHeight = halfPart ? h3 : h2;
float alpha = (float)i / h1;
float beta = (float)(i - (halfPart ? h2 : 0)) / segmentHeight;
Vec2i startPoint = t0 + (t2-t0)*alpha;
Vec2i endPoint = halfPart ? t1 + (t2-t1)*beta : t0 + (t1-t0)*beta;
if(startPoint.x > endPoint.x) std::swap(startPoint, endPoint);
for(int j = startPoint.x; j <= endPoint.x; j++) {
image.set(j, t0.y+i, color);
}
}
}
-
绘制循环中，变量 \(i\) 为 \(\vec{v_{0}v_{2}}\) 方向的插值步长，影响左端点采样。
-
halfPart 表征当前是否已经到达第二部分，也即 \(\vec{v_{1}v_{2}}\) 的插值方向；\(v_0\) 与 \(v_1\) 处同一水平线的情形同样被包含于 halfPart 的语境中。
-
halfPart 的值决定使用位移 \(h2\) 还是 \(h3\) 来进行右端点插值采样，并将值赋予变量 segmentHeight。
-
三线性插值的第一次插值：变量 \(alpha\) 代表 \(h1\) 位移上的比例系数；\(beta\) 代表 \(h2\) / \(h3\) 位移上的比例系数（\(h2\) 还是 \(h3\) 由 halfPart 决定）。
-
三线性插值的第二次插值（采样）：startPoint和endPoint分别为当前直线左右端点。此次插值（采样）受第一次插值算出的比例系数影响。
-
三线性插值的第三次插值（采样）：用一个循环计算startPoint和endPoint之间的点。
绘制效果：
可以看到各种情况都能正确绘制，左下角蓝色三角形就是 \(v_0\) 与 \(v_1\) 处同一水平线的情形。
重心坐标法 Barycentric Coordinates
假设我们能将三角形用一个包围盒BoundingBox框起来，然后依次遍历包围盒中的所有点，若点落在三角形内就进行着色，否则不进行任何操作。伪代码看起来像这样：
填充三角形(三角形) {
包围盒 = 找到三角形的包围盒(三角形的顶点);
for (包围盒里每一个点) {
if (这个点在三角形内) {
着色(此点);
}
}
}
这个算法重点是判定点是否在三角形内，对此采用重心坐标法。
给定一个2D三角形和一个点P，目标是计算相对于此三角形点P的重心坐标 \(\alpha\) 、\(\beta\)、\(\gamma\)，且 \(\alpha + \beta + \gamma = 1\)。实际上，只要计算出 \(\beta\) 与 \(\gamma\)，\(\alpha\) 用 \(1 - \beta - \gamma\) 表示即可，然后我们可以确定点P位置：
\(P = \alpha A + \beta B + \gamma C\) \((1)\)
可想象为，仅有质心的三个重物 \((1-\beta-\gamma), \beta, \gamma\) 分别置于三角形三点，则此时点P为此三角形的重心。（仔细想想，这是不是另一种版本的线性插值呢？）
我们稍微将思路调整一下，建立一个新坐标系。原点为A，基为 \(\vec{AC}\) 、\(\vec{AB}\)，则P可表示为：
\(P = A + \beta \vec{AB} + \gamma \vec{AC}\) \((2)\)
将A移到左边，则有：
\(\vec{AP} = \beta \vec{AB} + \gamma \vec{AC}\) \((3)\)
这样，我们便将情形转化到了向量表示上：
变换 \((3)\) 式：
\(\beta \vec{AB} + \gamma \vec{AC} + \vec{PA} = 0\) \((4)\)
将向量拆成 \((x,y)\) 笛卡尔坐标形式，便能得到线性方程：
\( \begin{equation} \begin{cases} \beta \cdot \vec{AB}_{x} + \gamma \cdot \vec{AC}_{x} + \vec{PA}_{x} = 0\\ \beta \cdot \vec{AB}_{y} + \gamma \cdot \vec{AC}_{y} + \vec{PA}_{y} = 0 \end{cases} \end{equation} \)
将其转为线性代数方程：
\( \begin{equation} \begin{cases} \left[\begin{array}{ccc} \beta & \gamma & 1 \end{array}\right]\left[\begin{array}{c} \vec{AB}_{x} \\ \vec{AC}_{x} \\ \vec{PA}_{x} \end{array}\right] = 0\\\\ \left[\begin{array}{ccc} \beta & \gamma & 1 \end{array}\right]\left[\begin{array}{c} \vec{AB}_{y} \\ \vec{AC}_{y} \\ \vec{PA}_{y} \end{array}\right] = 0 \end{cases} \end{equation} \)
现在设向量
-
\(\vec{n} = \left[\begin{array}{ccc} \beta & \gamma & 1 \end{array}\right]\)
-
\(\vec{i} = \left[\begin{array}{ccc} \vec{AB}_{x} & \vec{AC}_{x} & \vec{PA}_{x} \end{array}\right]\)
-
\(\vec{j} = \left[\begin{array}{ccc} \vec{AB}_{y} & \vec{AC}_{y} & \vec{PA}_{y} \end{array}\right]\)
神奇的事情出现了。没错，这个方程组表示，向量 \(\vec{n}\) 同时正交（垂直）于向量 \(\vec{i}\) 和 \(\vec{j}\) 。
感觉熟悉吗？换句话说，\(\vec{n} = \left[\begin{array}{ccc} \beta & \gamma & 1 \end{array}\right]\) 是 \(\vec{i} \times \vec{j}\) （叉积/外积）得到的法向量。
目前为止，你或许会感觉这些计算非常反直觉：\(\beta\) 与 \(\gamma\) 怎么就被包含在一个向量里了？向量 \(\vec{i}, \vec{j}\) 又是什么？我非常赞同，但谁叫数学如此巧妙呢？
有了上述结论，我们便可用 \(\vec{i} \times \vec{j}\) 来求我们需要的 \(\beta\)、\(\gamma\)，而 \(\vec{i}\)、\(\vec{j}\) 是可通过已知变量算出来的！然而\(\vec{i}\) 与 \(\vec{j}\)
叉积算出来的向量的形式不一定是 \(\left[\begin{array}{ccc} \beta & \gamma & 1 \end{array}\right]\)，而是：
\(k \cdot \left[\begin{array}{ccc} \beta & \gamma & 1 \end{array}\right] = \left[\begin{array}{ccc} k\beta & k\gamma & k \end{array}\right]\) ， \((k \in \mathbb{R}^*)\)
这是因为叉积算出来的本质上是 \(\vec{n}\) 的共线向量，恢复到 \(\vec{n}\) 只需将该向量的每个分量同时除以 \(k\) 即可，而 \(k\) 其实就是该向量的 \(z\) 值：
提前消除一个疑虑：叉积算出的向量 \(z\) 值为负是叉积顺序导致的结果，回忆一下叉积的右手螺旋定则，这会影响结果法向量的方向。
有了以上结论，我们可以着手计算了。得到的结果除以其 \(z\) 值之后的 \(\vec{n}\)，会有三种情况，分别是：
1. 三角形退化情况
此时 \(\vec{n}\) 的 \(z\) 值为 \(0\) 。回想一下，从叉积几何意义的角度来看，此时参与运算的 \(\vec{i}, \vec{j}\) 平行；
从另一个角度看，求重心坐标的“三角形”一定三点共线甚至两两重合，此时“三角形”至少退化degenerate为一条线段了，甚至变为一个点！不信的话我用叉积的坐标运算定义来证明：
\(\vec{n}_z = \vec{i}_x*\vec{j}_y-\vec{j}_x*\vec{i}_y = \vec{AB}_{x}*\vec{AC}_{y}-\vec{AC}_{x}*\vec{AB}_{y} = 0\)
\(\Rightarrow \vec{AB}_{x}*\vec{AC}_{y} = \vec{AC}_{x}*\vec{AB}_{y}\)
\(\Rightarrow \frac{\vec{AB}_{x}}{\vec{AB}_{y}} = \frac{\vec{AC}_{x}}{\vec{AC}_{y}}\)
也就是说，此时向量 \(\vec{AB}\) 与 \(\vec{AC}\) 共线（回忆我们在上面推导的向量共线等式），“三角形”看起来是这样的：
2. 点 \(P\) 在三角形外的情况 --- \(\beta\)、\(\gamma\) 至少有一个为负
光是看几何表示，我们就能确定，点 \(P\) 不在三角形内部时， \(\beta\)、\(\gamma\) 至少有一个为负。但是在计算过程中我们会发现，计算完 \(\vec{i} \times \vec{j}\) 得到的向量在还没有除以 \(z\) 之前，其 \(\beta\) 与 \(\gamma\) 有可能是正的。
但如果你脑袋转得够快，会明白这是由于叉积的顺序不正确，导致得出的向量与我们要的 \(\vec{n}\) 的方向相反。这时，只需将向量除以其 \(z\) 值，结果便会恢复正常。
3. 点 \(P\) 在三角形内的情况 --- 这是我们想要的！
此时便不用推导任何公式了，因为情形3正好与情形2相反：\(\beta\)、\(\gamma\) 均为正数。
呼，没想到一个看似简单的判定算法竟然花费我们这么多时间。结束了在公式上的挣扎，终于可以将结论付诸于代码了：
// 若返回的重心坐标有含负号的成员，则P被判定为无效点
Vec3f barycentric(Vec2i *pts, Vec2i P) {
// 叉积
Vec3f result = cross(Vec3f(pts[2][0]-pts[0][0], pts[1][0]-pts[0][0], pts[0][0]-P[0]), Vec3f(pts[2][1]-pts[0][1], pts[1][1]-pts[0][1], pts[0][1]-P[1]));
// 三角形退化情况
// result除以自身z之前，其z值可能为负
// z值为int型，它的绝对值小于1，说明它等于0，此时返回带负号的重心坐标
if (std::abs(u[2])<1) return Vec3f(-1,1,1);
// result除以z恢复为我们要的n
return Vec3f(1.f-(u.x+u.y)/u.z, u.y/u.z, u.x/u.z);
}
重心坐标算法写出来后，便可以完成整个包围盒光栅化算法了：
void triangle(Vec2i *pts, TGAImage &image, TGAColor color) {
Vec2i bboxmin(image.get_width()-1, image.get_height()-1);
Vec2i bboxmax(0, 0);
Vec2i clamp(image.get_width()-1, image.get_height()-1);
for (int i=0; i<3; i++) {
for (int j=0; j<2; j++) {
// 外层的std::max确保将坐标为负的点剔除
bboxmin[j] = std::max(0, std::min(bboxmin[j], pts[i][j]));
// 外层的std::min确保将超出屏幕外的点剔除
bboxmax[j] = std::min(clamp[j], std::max(bboxmax[j], pts[i][j]));
}
}
Vec2i P;
for (P.x=bboxmin.x; P.x<=bboxmax.x; P.x++) {
for (P.y=bboxmin.y; P.y<=bboxmax.y; P.y++) {
Vec3f bc_screen = barycentric(pts, P);
// 回忆重心坐标判定点的方法
if (bc_screen.x<0 || bc_screen.y<0 || bc_screen.z<0) continue;
image.set(P.x, P.y, color);
}
}
}
原作者写的代码非常优雅巧妙。我相信构造包围盒的代码中，
bboxmin[j] = std::max(0, std::min(bboxmin[j], pts[i][j]));
bboxmax[j] = std::min(clamp[j], std::max(bboxmax[j], pts[i][j]));
的外层是用来处理用户不当调用的边界判定。
想象一下，三角形顶点在由世界坐标变为屏幕坐标后才会被传到triangle函数进行光栅化，此时顶点坐标一定为非负且被限制在屏幕大小内，这时即使不加这两个外层判定，一样不会发生程序异常；但若用户未进行视区变换便将顶点的世界坐标直接传入triangle，那么这两层判定便能很好把握住边界处理。
扩展：简单的平面着色渲染
下图是对一个模型的二维三角形面片使用我们的代码进行随机上色的例子（其wireframe表示由右图展示，它由Bresenham直线算法绘制）：
但是这样的着色太过粗暴了，我们来考虑一个更接近现实的例子。
假设场景里有一个光源，多边形与光线正交时其接收光照强度最高；而若光线与多边形平行时，其接收的光照强度几乎为0。
设光线矢量为 \(\vec{l}\) ，平面法线为 \(\vec{n}\) 。换言之，光照强度与 \(\vec{l}\) 和 \(\vec{n}\) 之间的夹角 \(\theta\) 成负相关，这和 \(cos\theta\) 在区间 \([0, \frac{\pi}{2}]\) 单调递减的性质相契合。我们试试计算光线与平面法向量的内积：
\(\vec{l} \cdot \vec{n} = |\vec{l}| |\vec{n}|cos\theta\)
由于它们会被正规化为单位向量，也即 \(|\vec{l}| |\vec{n}| = 1\)，那么上式变为：\(\vec{l} \cdot \vec{n} = cos\theta\)，这正是我们想要的。
你或许会思考这种情形：光线从平面的背面照过来，或者计算出的法向量为负。此时光线与平面法向量夹角 \(\theta \gt \frac{\pi}{2}\) ，\(cos\theta \lt 0\)。那么只需在代码判断时，若内积为负则不进行着色。这在图形学中叫做 Back-face culling背面剔除
Vec3f light_dir(0,0,-1);
for (int i=0; i<model->nfaces(); i++) {
std::vector<int> face = model->face(i);
// 屏幕坐标，由世界坐标经视区变换得来
Vec2i screen_coords[3];
// 世界坐标
Vec3f world_coords[3];
for (int j=0; j<3; j++) {
Vec3f v = model->vert(face[j]);
// 视区变换
screen_coords[j] = Vec2i((v.x+1.)*width/2., (v.y+1.)*height/2.);
world_coords[j] = v;
}
// 法向量需用三角形顶点的世界坐标来计算
Vec3f n = (world_coords[2]-world_coords[0])^(world_coords[1]-world_coords[0]);
n.normalize();
float intensity = n*light_dir;
// Back-face culling
if (intensity>0){
triangle(screen_coords[0], screen_coords[1], screen_coords[2], image, TGAColor(intensity*255, intensity*255, intensity*255, 255));
}
}
效果如图：
你会发现口腔内的点竟然凸出来了，这是这种算法的一个缺点。下一节会用深度测试来修复这种问题。
