光栅化的算法实现

最新推荐文章于 2025-05-05 19:30:04 发布

翻译最新推荐文章于 2025-05-05 19:30:04 发布 · 4.8k 阅读

光栅化的算法实现

一直用OpenGL绘制东西的时候，就会想到我们在写可编程管线的时候，都是使用gl_position去保存一个物体经过model(模型矩阵)矩阵，view(视口)矩阵以及projection(投影)矩阵变换后的位置。然后利用该位置信息在片段着色器中为其上色并显示在屏幕上。这便是我们OpenGL实现光栅化的过程。那到底OpenGL底层是如何做到这一点呢？我们尝试着脱离OpenGL如绘制一个图形。

光栅化算法框架

在了解光栅化问题时，我们需要清楚以往时用什么方法去解决可见性问题(即展示在我们屏幕上)。最常见的两种就是光栅化和光线追踪。

光线追踪解决可见性问题：

如上图，我们从摄像机的位置，向我们的屏幕(近平面)的每个像素中心发射一条光线，光线和场景中哪个对象相交，并且距离屏幕最近，则将该相交对象的颜色将保存在屏幕上。我们可以用以下伪代码表示该过程。

for (each pixel in image) { 
    Ray R = computeRayPassingThroughPixel(x,y); //获取光线的方向
    float tclosest = INFINITY; //记录最近的距离
    Triangle triangleClosest = NULL; //最近的三角形
    //寻找最近的三角形
    for (each triangle in scene) { 
        float thit; 
        if (intersect(R, object, thit)) {//判断相交成功返回1 
             if (thit < closest) { 
                 triangleClosest = triangle; //如果最近则保存该三角形
             } 
        } 
    } 
    if (triangleClosest) { 
        imageAtPixel(x,y) = triangleColorAtHitPoint(triangle, tclosest); //记录颜色
    } 
}

我们发现这里的代码只判断和三角形求交点，但是我们的图形是其他形状的呢？而我们的OpenGL，DirectX都是只处理三角形的，这又是为什么呢？这是因为三角形是最简单的图形，它们构成的片面是共面的而且最容易求交计算。我们知道光线追踪最昂贵的就是光线和对象的求交计算。所以使用三角形可以缓解我们计算的压力，因此我们在渲染其他图形时，也会使用三角剖分去将一个复杂几何图形转换为多个三角形的集合。

因为昂贵的求交运算，我们现在基本上都是使用光栅化解决可见性问题。但是在渲染高质量图形时，都会结合光栅化和光线跟踪，利用光栅化和z缓存计算眼睛第一个可见表面，反射和折射利用光线跟踪计算。这里还要注意区别光线跟踪和光线传输算法，前者是用来处理可见性问题，后者是用来计算某一点的颜色。

光栅化解决可见性问题：

栅格化采用相反的方法。为了解决可见性，它实际上是将三角形“投影”到屏幕上，换句话说，我们使用透视投影从该三角形的3D表示转换为2D表示。可以通过将组成三角形的顶点投影到屏幕上来轻松完成此操作。该算法的下一步是使用某种技术来填充该2D三角形覆盖的图像的所有像素。这两个步骤如下图所示：

从技术角度来看，它们非常容易执行。投影步骤只需要使用投影矩阵进行变换，并将坐标空间从观察空间重新映射到裁剪空间。找出结果三角形覆盖图像中的哪些像素也是非常简单的，将在后面进行描述。

与光线跟踪方法相比，该算法是什么样的？首先，请注意，我们需要对场景中的所有三角形进行迭代，而不是首先在光栅化中在外部循环中迭代图像中的所有像素。然后，在内部循环中，我们遍历图像中的所有像素，并找出当前像素是否“包含”在当前三角形的“投影图像”内（如上图）。换句话说，两种算法的内部和外部循环被交换。我们看一下光栅化算法的伪代码：

for (each triangle in scene) { //为每个三角形迭代
    //步骤1:利用投影变换去变化坐标点
    Vec2f v0 = perspectiveProject(triangle[i].v0); 
    Vec2f v1 = perspectiveProject(triangle[i].v1); 
    Vec2f v2 = perspectiveProject(triangle[i].v2); 
    for (each pixel in image) { //迭代所有像素
        //步骤2:当像素位于三角形内时，绘制该像素
        if (pixelContainedIn2DTriangle(v0, v1, v2, x, y)) { 
            image(x,y) = triangle[i].color; 
        } 
    } 
}

优化的三角形边界盒：

虽然我们的光栅化算法相比光线追踪有很大的性能优势，但是即使三角形只包含少数几个像素时，我们仍然需要循环迭代所有的像素。这是一种巨大的浪费，因此我们这里可以做一个优化，我们可以在投影三角形的时候，并计算它的2D边界框，这样只需要遍历边界框内的像素。这在渲染拥有数百万个三角形的复杂对象时，会提升巨大的性能，刚方法如下图：

我们对投影后三角形的三个点计算边框并取整，然后匹配像素只需要循环边框内的像素就行了。我们计算边框的伪代码如下：

Vec2f bbmin = INFINITY, bbmax = -INFINITY; 
Vec2f vproj[3]; //三角形的3个顶点
for (int i = 0; i < 3; ++i) { 
    vproj[i] = projectAndConvertToNDC(triangle[i].v[i]);//投影三角形并转换到标准设备空间
    // 将标准设备空间的点转换为屏幕像素
    vproj[i].x *= imageWidth; 
    vproj[i].y *= imageHeight; 
    //计算左上角和右下角的值
    if (vproj[i].x < bbmin.x) bbmin.x = vproj[i].x); 
    if (vproj[i].y < bbmin.y) bbmin.y = vproj[i].y); 
    if (vproj[i].x > bbmax.x) bbmax.x = vproj[i].x); 
    if (vproj[i].y > bbmax.y) bbmax.y = vproj[i].y); 
}

由于这里我们计算的NDC空间为[0,1](OpenGL中为[-1,1])，所以我们要把他转换到屏幕像素的坐标x $\in$ [0,imageWidth-1]和y $\in$ [0,imageHeight-1]。我们获取完2D边界框后可以进行如下循环：

uint xmin = std::max(0, std:min(imageWidth - 1, std::floor(min.x))); 
uint ymin = std::max(0, std:min(imageHeight - 1, std::floor(min.y))); 
uint xmax = std::max(0, std:min(imageWidth - 1, std::floor(max.x))); 
uint ymax = std::max(0, std:min(imageHeight - 1, std::floor(max.y))); 
for (y = ymin; y <= ymin; ++y) { 
    for (x = xmin; x <= xmax; ++x) { 
        // 检查像素是否位于三角形内
        if (pixelContainedIn2DTriangle(v0, v1, v2, x, y)) { 
            image(x,y) = triangle[i].color; 
        } 
    } 
}

我们只要在2D边界框内，判断像素属于三角形内的点，则为其上色。

Z-buffer算法：

我们的目标是产生场景的图像。我们有两种可视化程序结果的方式，一种是将渲染的图像直接显示在屏幕上，另一种是将图像保存到磁盘上，然后使用诸如Photoshop之类的程序稍后预览图像。但是在这两种情况下，我们都需要以某种方式存储正在渲染的图像，并且为此，我们在CG中使用所谓的图像或帧缓冲区。就是具有图像大小的二维颜色数组。在渲染过程开始之前，将创建帧缓冲区，并将像素全部设置为黑色。在渲染时，当对三角形进行栅格化时，如果给定像素与给定三角形重叠，则我们将该三角形的颜色存储在该像素位置的帧缓冲区中。光栅化所有三角形后，帧缓冲区将包含场景的图像。剩下要做的就是将缓冲区的内容显示到屏幕上，或将其内容保存到文件中。但是如果我们的给定像素与多个给定三角形重叠，这时候这么选择要显示的三角形呢？显然我们需要的是距离我们屏幕最近的三角形的点。我们可以采用称为Z-缓存算法的方法来获取最近的三角形的点。

Z缓冲区无非是另一个二维数组，它的维数与图像的维数相同，但是它不是颜色数组，而只是一个浮点数数组。在开始渲染图像之前，我们将该数组中的每个像素初始化为非常大的数量。当像素与三角形重叠时，我们还将读取存储在该像素位置z缓冲区中的值。该数组用于存储从相机到图像中任何像素重叠的最近三角的距离。伪代码如下：

// Z-缓存就是一个浮点的2维数组
float buffer = new float [imageWidth * imageHeight]; 
// 用非常大的数初始化该2维数组
for (uint32_t i = 0; i < imageWidth * imageHeight; ++i) 
    buffer[i] = INFINITY; 
 
for (each triangle in scene) { 
    // 投影三角形
    ... 
    // 计算2D边界盒
    ... 
    for (y = ymin; y <= ymin; ++y) { 
        for (x = xmin; x <= xmax; ++x) { 
            // 判断像素是否和三角形重叠
            float z; // 三角形上该点到摄像机距离(深度) 
            if (pixelContainedIn2DTriangle(v0, v1, v2, x, y, z)) { 
                // 如果当前的三角形的z值是最近的则更新z缓存的值，并绘制颜色
                if (z < zbuffer(x,y)) { 
                    zbuffer(x,y) = z; 
                    image(x,y) = triangle[i].color; 
                } 
            } 
        } 
    } 
}

我们介绍完光栅化框架后我们来了解具体实现的细节。

第一步我们先将三角形进行透视投影变换，这一步详细可以看https://blog.csdn.net/qq_39300235/article/details/90670282。

第二步我们在上文完成了2D边界盒的计算。

第三步我们需要了解如何判断像素位于三角形内部。

第四步了解三角形内部点的属性插值。

边缘函数

要找到像素是否与三角形重叠有许多方法，这里我们使用Juan Pineda于1988年提出，并发表在论文“多边形栅格化的并行算法”中的方法。我们将首先描述他的方法的原理。假设三角形的边缘可以看作是将2D平面（图像的平面）一分为二的线。 Pineda方法的原理是找到一个称为边缘函数的函数，这样，当我们测试该点在哪条线的哪一侧（图2中的点P）时，该函数将返回负数。在该行的左边，当正点在该行的右边时为正数；在该点正好在该行上时为零。如下图所示：

上图中我们将此方法应用于三角形的第一个边缘（由顶点v0-v1定义。请注意顺序很重要）。如果现在将相同的方法应用于其他两个边（v1-v2和v2-v0），则可以清楚地看到存在一个区域（白色三角形），其中所有点均为正。如果P实际上是像素中心的一个点，则可以使用此方法查找像素是否与三角形重叠。如果在这一点上，我们发现边缘函数为所有三个边缘返回正数，则像素包含在三角形中（或可能位于其边缘之一上），如下图所示：

要注意的是这里的顺序一定要是顺时针(如果为逆时针，则内部为负数)，Pinada使用的函数也恰好是线性的，这意味着可以递增地计算它。现在我们了解了原理，让我们找出该函数是什么。边函数定义为（对于由顶点V0和V1定义的边）：

该函数值和点P有如下关系：

E(P) > 0 ，P在右侧
E(P) = 0 ，P在线上
E(P) < 0 ，P在左侧

其实该公式正是向量 $\overrightarrow{V_{0}V_{1}}$ 和向量 $\overrightarrow{V_{0}P}$ 的叉积，叉积又可以表示平行四边形面积，叉积公式如下所示：

其中 $sin(\theta )$ 表现了正负性(我们用右手定则也可以轻松判断正负)。我们分别计算三条利用公式：

若三者都为大于0的数，则表明点P位于三角形内。代码：

bool edgeFunction(const Vec2f &a, const Vec3f &b, const Vec2f &c) 
{ 
    return ((c.x - a.x) * (b.y - a.y) - (c.y - a.y) * (b.x - a.x) >= 0); 
} 
 
bool inside = true; 
inside &= edgeFunction(V0, V1, p); 
inside &= edgeFunction(V1, V2, p); 
inside &= edgeFunction(V2, V0, p); 
 
if (inside == true) { 
    // 点P在由点V0，V1，V2构成的三角形内部
    ... 
}

重心坐标：

在进一步介绍之前，我们先了解重心坐标是什么。首先，它们是一组三个浮点数，在本课中，我们将分别表示 $\lambda _{0}$ ， $\lambda _{1}$ 和 $\lambda _{2}$ 。可以通过以下方式使用坐标定义三角形上的任何点：

通常，V0，V1和V2是三角形的顶点。这些坐标可以取任何值。而对于三角形内部（或其边缘之一）上的点。 $\lambda _{0}$ ， $\lambda _{1}$ 和 $\lambda _{2}$ 只能在[0,1]范围内，并且总和等于1。也就是说：

这是一种插值形式。有时也将它们定义为三角形顶点的权重。插值三角形的顶点以找到三角形内部的点的位置并没有太大用处。但是，该方法还可以用于在三角形的表面上插值在三角形顶点处定义的任何数量或变量。假设您在三角形的每个顶点上定义了一种颜色。假设V0为红色，V1为绿色，V2为蓝色。如下图：

您想要做的是找到如何在三角形的表面上插入这三种颜色。如果知道三角形上的点P的重心坐标，则其颜色CP（三角形顶点颜色的组合）定义为：

这是一种非常方便的技术，它将对渲染三角形有用。与三角形的顶点关联的数据称为顶点属性。这是CG中非常普遍且非常重要的技术。最常见的顶点属性是颜色，法线和纹理坐标。实际上，这意味着在定义三角形时，不仅将三角形的顶点传递给渲染器，而且将其相关的顶点属性传递给渲染器。例如，如果要渲染三角形，则可能需要颜色和法线顶点属性，这意味着每个三角形将由3个点（三角形顶点位置），3个颜色（三角形顶点的颜色）和3个法线定义（三角形顶点的法线）。法线也可以在三角形的表面内插。插值法线用于一种称为“平滑着色”的技术，该技术最早由Henri Gouraud引入。这里我们先不介绍。

我们如何找到重心坐标呢？如果点P在三角形内，那么通过查看下图可以看到，我们可以绘制三个子三角形：V0-V1-P（绿色），V1-V2-P（洋红色）和V2-V0- P（青色）。很明显，这三个子三角形的面积之和等于三角形V0-V1-V2的面积：

即：

而我们刚才的边缘函数：

正好是三角形面积的两倍。即：

那么我们可以直接用边缘函数求解重心坐标：

我们来看一下计算重心坐标的代码：

float edgeFunction(const Vec2f &a, const Vec3f &b, const Vec2f &c) 
{ 
    return (c.x - a.x) * (b.y - a.y) - (c.y - a.y) * (b.x - a.x); 
} 
 
float area = edgeFunction(v0, v1, v2); // 三角形的面积(在乘上2，下面三个都是)
float w0 = edgeFunction(v1, v2, p); // 三角形v1v2p的面积 
float w1 = edgeFunction(v2, v0, p); // 三角形v2v0p的面积 
float w2 = edgeFunction(v0, v1, p); // 三角形v0v1p的面积
 
// 判断点是否在三角形内
if (w0 >= 0 && w1 >= 0 && w2 >= 0) { 
    // 计算重心坐标
    w0 /= area; 
    w1 /= area; 
    w2 /= area; 
}

现在让我们在一个产生实际图像的程序中测试我们在本文中讲解的知识。我们假设已经投影了三角形。我们还将为三角形的每个顶点分配颜色。这是图像的形成方式。我们将遍历图像中的所有像素，并使用边缘函数方法测试它们是否与三角形重叠。对照像素的当前位置测试三角形的所有三个边缘，如果边缘函数对所有边缘返回正数，则像素与三角形重叠。然后，我们可以计算像素的重心坐标，并使用这些坐标通过对三角形每个顶点定义的颜色进行插值来对像素进行着色。帧缓冲区的结果将保存到PPM文件，该程序为：


#include <cstdio> 
#include <cstdlib> 
#include <fstream> 
 
typedef float Vec2[2]; 
typedef float Vec3[3]; 
typedef unsigned char Rgb[3]; 
 
inline 
float edgeFunction(const Vec2 &a, const Vec2 &b, const Vec2 &c) 
{ return (c[0] - a[0]) * (b[1] - a[1]) - (c[1] - a[1]) * (b[0] - a[0]); } 
 
int main(int argc, char **argv) 
{ 
    Vec2 v0 = {491.407, 411.407}; 
    Vec2 v1 = {148.593, 68.5928}; 
    Vec2 v2 = {148.593, 411.407}; 
    Vec3 c0 = {1, 0, 0}; 
    Vec3 c1 = {0, 1, 0}; 
    Vec3 c2 = {0, 0, 1}; 
 
    const uint32_t w = 512; 
    const uint32_t h = 512; 
 
    Rgb *framebuffer = new Rgb[w * h]; 
    memset(framebuffer, 0x0, w * h * 3); 
 
    float area = edgeFunction(v0, v1, v2); 
 
    for (uint32_t j = 0; j < h; ++j) { 
        for (uint32_t i = 0; i < w; ++i) { 
            Vec2 p = {i + 0.5f, j + 0.5f}; 
            float w0 = edgeFunction(v1, v2, p); 
            float w1 = edgeFunction(v2, v0, p); 
            float w2 = edgeFunction(v0, v1, p); 
            if (w0 >= 0 && w1 >= 0 && w2 >= 0) { 
                w0 /= area; 
                w1 /= area; 
                w2 /= area; 
                float r = w0 * c0[0] + w1 * c1[0] + w2 * c2[0]; 
                float g = w0 * c0[1] + w1 * c1[1] + w2 * c2[1]; 
                float b = w0 * c0[2] + w1 * c1[2] + w2 * c2[2]; 
                framebuffer[j * w + i][0] = (unsigned char)(r * 255); 
                framebuffer[j * w + i][1] = (unsigned char)(g * 255); 
                framebuffer[j * w + i][2] = (unsigned char)(b * 255); 
            } 
        } 
    } 
 
    std::ofstream ofs; 
    ofs.open("./raster2d.ppm"); 
    ofs << "P6\n" << w << " " << h << "\n255\n"; 
    ofs.write((char*)framebuffer, w * h * 3); 
    ofs.close(); 
 
    delete [] framebuffer; 
 
    return 0; 
}

程序结果为：

这和我们在OpenGL中渲染的效果几乎一摸一样。

光栅化规则：

在某些特殊情况下，一个像素可能重叠多个三角形。如图下所示：

当一个像素恰好位于两个三角形共享的边缘上时，就会发生这种情况。此类像素将通过两个三角形的覆盖率测试。如果它们是半透明的，则由于半透明对象彼此组合的方式（想象两张叠加的半透明塑料薄片），像素重叠两个三角形的地方可能会出现暗边。比不透明的纸张更不透明，并且看起来比单独的纸张更暗）。您将得到类似于下图所示的内容：

这是一条较暗的线，其中两个三角形共享一条边。该问题的解决方案是提出某种规则，以确保像素永远不会重叠两个共享边的三角形两次。我们该怎么做？大多数图形API（例如OpenGL和DirectX）都定义了一些它们称为左上角的规则。左上角的规则是，如果像素或点位于三角形内部或位于三角形的上边缘或任何被视为左边缘的边缘上，则该像素或点被视为与三角形重叠。什么是上边缘和左边缘？

上边缘是完全水平的边缘，其定义顶点在第三个边缘之上。从技术上讲，这意味着向量V [（X + 1）％3] -V [X]的y坐标等于0，并且其x坐标为正（大于0）。

左边缘本质上是上升的边缘。请记住，在我们的情况下，顶点是按顺时针顺序定义的。如果边缘的相应向量V [（X + 1）％3] -V [X]（其中X可以为0、1、2）具有y坐标为正，则认为该边缘上升。

上述图示为：

绿色标记的边是我们的上边缘和左边缘。实现代码如下：

//定义的坐标
Vec2f v0 = { ... }; 
Vec2f v1 = { ... }; 
Vec2f v2 = { ... }; 
 
//边缘函数计算
float w0 = edgeFunction(v1, v2, p); 
float w1 = edgeFunction(v2, v0, p); 
float w2 = edgeFunction(v0, v1, p); 
 
Vec2f edge0 = v2 - v1; 
Vec2f edge1 = v0 - v2; 
Vec2f edge2 = v1 - v0; 
 
bool overlaps = true; 
 
// 通过边缘函数返回值判断点是否在边上，如果点在边上，判断是否在左边缘或者上边缘。
//如果不在边上判断是否在三角形内
// 上边缘本质为向量y值为0，x值大于0，左边缘本质向量是y大于0。
overlaps &= (w0 == 0 ? ((edge0.y == 0 && edge0.x > 0) ||  edge0.y > 0) : (w0 > 0)); 
overlaps &= (w1 == 0 ? ((edge1.y == 0 && edge1.x > 0) ||  edge1.y > 0) : (w1 > 0)); 
overlaps &= (w1 == 0 ? ((edge2.y == 0 && edge2.x > 0) ||  edge2.y > 0) : (w2 > 0)); 
 
if (overlaps) { 
    //如果像素和三角形重叠，则进行之后的光栅化
    ... 
}

深度插值

我们根据前文的知识，我们再看我们的光栅化器，伪代码如下：

float *depthBuffer = new float [imageWidth * imageHeight]; 
// 初始化深度缓存
for (uint32_t y = 0; y < imageHeight; ++y) 
    for (uint32_t x = 0; x < imageWidth; ++x) 
        depthBuffer[y][x] = INFINITY; 
 
for (each triangle in scene) { 
    //投影三角形
    ... 
    // 计算三角形的边界盒
    ... 
    for (uint32_t y = bbox.min.y; y <= bbox.max.y; ++y) { 
        for (uint32_t x = bbox.min.x; x <= bbox.max.x; ++x) { 
           if (pixelOverlapsTriangle(i + 0.5, j + 0.5) { 
                // 计算三角形上点的深度值
                float z = computeDepth(...); 
                // 判断当前点的深度并记录颜色
                if (z < depthBuffer[y][x]) { 
                     // 更新深度缓存
                     depthBuffer[y][x] = z; 
                     frameBuffer[y][x] = triangleColor; 
                } 
            } 
        } 
    } 
}

我们发现还缺少深度值插值计算，我们怎么通过三角形的三个顶点插值出中间任意位置的深度值z呢？你可能会想到用上面的重心坐标向插值颜色那样插值深度值。但是这个地方直接这么使用却是行不通的。因为深度在映射之后变成了非线性，我们看下图：

我们设定屏幕到原点的距离为1，且平行于X轴。我们图中 $V_{0}V_{1}$ 上的点P投影到屏幕的 $P{}'$ 。 $\lambda_{1} =V_{0}P/V_{0}V_{1}=0.666$ ，而 $\lambda_{2} ={V}'_{0}{P}'/{V}'_{0}{V}'_{1}=0.8333$ 。两者并不相等，因为投影并不能保留距离。所以这里我们便不能直接运用求重心坐标的公式，但是方法还是一致的，只是公式有所改动：

因为我们是在投影之后进行计算，所以这里的 $\lambda$ 是投影后计算的。公式用简单的相似三角形能推出(这里就不写推倒过程了)。我们带入公式计算：

得P.z=4正好符合上图。

透视图正确的顶点属性插值

我们提到在传入顶点属性的时候可以传入颜色值和纹理坐标值，OpenGL在底层帮我们插值。颜色插值公式：

纹理坐标插值公式：

我们都可以用重心坐标去计算，但是我们的三角形通过投影后会产生形变，会出现这种情况：

在投影之后，由于不保留距离，直接获取的 $\lambda$ 会出错：

由于我们从右下角往上看，右图才是我们想要的。而左图则还是位于中心点与实际不符。找到正确的解决方案并不难。假设我们有一个三角形，在三角形的每一侧都有两个z坐标Z0和Z1，如下图所示：

如果我们连接这两个点，则可以使用线性插值对这条线上的点的z坐标进行插值。我们可以通过在三角形上分别比Z0和Z1相同的位置定义两个顶点属性C0和C1的值来执行相同的操作。从技术上讲，由于Z和C都是使用线性插值计算的，因此我们可以编写以下等式：

根据前面深度插值公式：

带入得：

我们利用这个知识点完成一个程序，代码如下：

// 编译:
// c++ -o raster3d raster3d.cpp  对于无使用投影矫正
// c++ -o raster3d raster3d.cpp -D PERSP_CORRECT  使用投影矫正

#include <cstdio> 
#include <cstdlib> 
#include <fstream> 
 
typedef float Vec2[2]; 
typedef float Vec3[3]; 
typedef unsigned char Rgb[3]; 
 
inline 
float edgeFunction(const Vec3 &a, const Vec3 &b, const Vec3 &c) 
{ return (c[0] - a[0]) * (b[1] - a[1]) - (c[1] - a[1]) * (b[0] - a[0]); } 
 
int main(int argc, char **argv) 
{ 
    Vec3 v2 = { -48, -10,  82}; 
    Vec3 v1 = {  29, -15,  44}; 
    Vec3 v0 = {  13,  34, 114}; 
    Vec3 c2 = {1, 0, 0}; 
    Vec3 c1 = {0, 1, 0}; 
    Vec3 c0 = {0, 0, 1}; 
 
    const uint32_t w = 512; 
    const uint32_t h = 512; 
 
    // 简单投影至屏幕
    v0[0] /= v0[2], v0[1] /= v0[2]; 
    v1[0] /= v1[2], v1[1] /= v1[2]; 
    v2[0] /= v2[2], v2[1] /= v2[2]; 
    // 转换至DNC坐标空间
    v0[0] = (1 + v0[0]) * 0.5 * w, v0[1] = (1 + v0[1]) * 0.5 * h; 
    v1[0] = (1 + v1[0]) * 0.5 * w, v1[1] = (1 + v1[1]) * 0.5 * h; 
    v2[0] = (1 + v2[0]) * 0.5 * w, v2[1] = (1 + v2[1]) * 0.5 * h; 
 
#ifdef PERSP_CORRECT 
    // 投影矫正
    c0[0] /= v0[2], c0[1] /= v0[2], c0[2] /= v0[2]; 
    c1[0] /= v1[2], c1[1] /= v1[2], c1[2] /= v1[2]; 
    c2[0] /= v2[2], c2[1] /= v2[2], c2[2] /= v2[2]; 
    // 与计算1/z
    v0[2] = 1 / v0[2], v1[2] = 1 / v1[2], v2[2] = 1 / v2[2]; 
#endif 
 
    Rgb *framebuffer = new Rgb[w * h]; 
    memset(framebuffer, 0x0, w * h * 3); 
 
    float area = edgeFunction(v0, v1, v2); 
 
    for (uint32_t j = 0; j < h; ++j) { 
        for (uint32_t i = 0; i < w; ++i) { 
            Vec3 p = {i + 0.5, h - j + 0.5, 0}; 
            float w0 = edgeFunction(v1, v2, p); 
            float w1 = edgeFunction(v2, v0, p); 
            float w2 = edgeFunction(v0, v1, p); 
            if (w0 >= 0 && w1 >= 0 && w2 >= 0) { 
                w0 /= area; 
                w1 /= area; 
                w2 /= area; 
                float r = w0 * c0[0] + w1 * c1[0] + w2 * c2[0]; 
                float g = w0 * c0[1] + w1 * c1[1] + w2 * c2[1]; 
                float b = w0 * c0[2] + w1 * c1[2] + w2 * c2[2]; 
#ifdef PERSP_CORRECT  //投影矫正
                float z = 1 / (w0 * v0[2] + w1 * v1[2] + w2 * v2[2]); 
                r *= z, g *= z, b *= z; 
#endif 
                framebuffer[j * w + i][0] = (unsigned char)(r * 255); 
                framebuffer[j * w + i][1] = (unsigned char)(g * 255); 
                framebuffer[j * w + i][2] = (unsigned char)(b * 255); 
            } 
        } 
    } 
 
    std::ofstream ofs; 
    ofs.open("./raster2d.ppm"); 
    ofs << "P6\n" << w << " " << h << "\n255\n"; 
    ofs.write((char*)framebuffer, w * h * 3); 
    ofs.close(); 
 
    delete [] framebuffer; 
 
    return 0; 
}

该程序运行可以对比我们是否启用投影矫正产生的结果。