软光栅化渲染器学习(一)

本文为笔者入门图形学，学习光栅化渲染的记录。

SDL3框架

使用SDL3，可以轻松地创建一个窗口，作为渲染器的框架。

定义Window类如下：

class Window
{
  public:
    friend class Renderer;
    Window(const std::string &name, size_t width, size_t height)
        : width(width), height(height)
    {
        sdl_window = SDL_CreateWindow(name.c_str(), width, height, SDL_WINDOW_RESIZABLE);
        if (!sdl_window)
        {
            SDL_Log("Could not create a window: %s", SDL_GetError());
        }
        SDL_SetWindowPosition(sdl_window, SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED);
    }

    ~Window()
    {
        SDL_DestroyWindow(sdl_window);
        SDL_Quit();
    }
    void PollEvents()
    {
        SDL_Event event;
        while (SDL_PollEvent(&event))
        {
            switch (event.type)
            {
            case SDL_EVENT_QUIT:
                shouldClose = true;
                break;
            }
        }
        SDL_Delay(16);
    }
    void UpdateSurface()
    {
        SDL_UpdateWindowSurface(sdl_window);
    }
    bool ShouldClose() const { return shouldClose; }
  private:
    bool shouldClose = false;
    size_t width;
    size_t height;
    SDL_Window *sdl_window;
};

定义Renderer类：

class Renderer
{
  public:
    Renderer(Window &window)
    {
        this->window = &window;
        sdl_renderer = SDL_CreateRenderer(window.sdl_window, nullptr);
        if (!sdl_renderer)
        {
            SDL_Log("Create renderer failed: %s", SDL_GetError());
        }
        sdl_surface = SDL_CreateSurface(window.width, window.height, SDL_PIXELFORMAT_RGBA8888);
        if (!sdl_surface)
        {
            SDL_Log("Create surface failed: %s", SDL_GetError());
        }
        SDL_LockSurface(sdl_surface);
    }
    ~Renderer()
    {    
        SDL_DestroyRenderer(sdl_renderer);
    }
    void Clear(uint8_t r, uint8_t g, uint8_t b)
    {
        SDL_ClearSurface(sdl_surface, (float)r / 255, (float)g / 255, (float)b / 255, 1);
    }
    void Render()
    {
        SDL_Texture *texture = SDL_CreateTextureFromSurface(sdl_renderer, sdl_surface);
        SDL_RenderTexture(sdl_renderer, texture, NULL, NULL);
        SDL_RenderPresent(sdl_renderer);
    }
    // 设置单个像素颜色
    void SetPixel(size_t x, size_t y, uint8_t r, uint8_t g, uint8_t b)
    {
        uint32_t *pixels = static_cast<uint32_t *>(sdl_surface->pixels);
        size_t index = y * (sdl_surface->pitch / sizeof(uint32_t)) + x;

        pixels[index] = SDL_MapRGBA(
            SDL_GetPixelFormatDetails(sdl_surface->format),
            nullptr,
            r, g, b, 255);
    }
  private:
    SDL_Renderer *sdl_renderer;
    SDL_Surface *sdl_surface;
    Window *window;
};

并且在之后的渲染实现中，只会用到设置单个像素颜色的函数SetPixel。
接下来，在主函数中：

int main()
{
    Window window("CPURenderer", 1920, 1080);
    Renderer renderer(window);

    while (!window.ShouldClose())
    {
        window.PollEvents();
        window.UpdateSurface();
        renderer.Clear(200, 150, 100);
        renderer.Render();
    }

    return 0;
}

这样就实现了基本的每帧清除缓冲并渲染的功能了。

绘制直线

Bresenham 直线算法

Bresenham 直线算法用于绘制给定两端点的直线，其优势在于不用任何浮点数运算，仅通过整数的加减和比较就能完成，从而效率较高。

接下来以$k\in[0,1]$的直线为例来解释该算法原理。首先，由于x方向的差值较大，所以每次迭代时将x坐标+1，维护一个误差$e$，来控制y坐标的变化。具体而言，初始时$e$为0，每次x方向+1时，y方向实际应增加$k$，此时若绘制的点的y坐标不变，则误差$e$就不变；若绘制的点的y坐标+1,则误差$e$减1。判断y坐标的方法为：当误差$e$大于0.5时，y坐标+1；误差$e$小于0.5时，y坐标不变。

由于$k=\frac{\Delta y}{\Delta x}$，因此将上述过程中的数字都乘上$2\Delta x$，结果等效。

此时算法过程变为：每次x方向+1时，误差$e$增加$2\Delta y$，若绘制的点y坐标+1，则误差$e$减$2\Delta x$。判断y坐标的方法为：当误差$e$大于$\Delta x$时，y坐标+1；误差$e$小于$\Delta x$时，y坐标不变。此时我们消除了所有的浮点型变量，优化了效率。

将算法适用范围进行拓宽，对于$k\in(1,\infty)$的直线，将上述算法中的所有x与y调换即可。对于$k$为负以及两点位置关系的多种情况，可以用signx，signy规定每次迭代的x与y的走向。最终算法实现如下：

void BresenhamDrawLine(size_t x0, size_t y0, size_t x1, size_t y1, uint8_t r, uint8_t g, uint8_t b)
{
    int dx = abs((int)x1 - (int)x0);
    int dy = abs((int)y1 - (int)y0);
    int signx = (x1 > x0) ? 1 : -1;
    int signy = (y1 > y0) ? 1 : -1;
    int e = 0;
    if (dx > dy) // k < 1
    {
        for (int x = x0, y = y0; x != x1; x += signx)
        {
            SetPixel(x, y, r, g, b);
            e += 2 * dy;
            if (e >= dx)
            {
                y += signy;
                e -= 2 * dx;
            }
        }
    }
    else // k > 1
    {
        for (int x = x0, y = y0; y != y1; y += signy)
        {
            SetPixel(x, y, r, g, b);
            e += 2 * dx;
            if (e >= dy)
            {
                x += signx;
                e -= 2 * dy;
            }
        }
    }
}

Cohen-sutherland 直线裁切算法

该算法用于剔除直线在视口之外的部分。

将视口的上下左右用位掩码表示，左为0001，右为0010，下为0100，上为1000，四角的位掩码则用两者的并表示，左上为1001，右上为1010，左下为0101，右下为0110。迭代时判断，如果直线两端点对应位掩码的交不为0，说明两端点一定位于窗口同一侧，则直线必定不会经过窗口，直接舍弃，如果直线两端点的并为0，说明两端点都在窗口内，直接采纳。对于其他情况，需要进行裁切。首先找到直线两端点中位于窗口外的点，在该算法中，用两个位掩码中数字大的即可（因为窗口内的点位掩码为0，是最小）。判断该点方位，并用直线方程运算出新的点。以窗口上方的点为例，直线方程为$\frac{y-y_0}{x-x_0}=\frac{y_1-y_0}{x_1-x_0}$，令$y=y_{min}$，得到对应$x=\frac{(x_1-x_0)(y_{min}-y_0)}{y_1-y_0}+x_0$，更新该点坐标为新的$(x,y)$。

完整代码如下：

static constexpr uint8_t INSIDE = 0b0000;
static constexpr uint8_t LEFT = 0b0001;
static constexpr uint8_t RIGHT = 0b0010;
static constexpr uint8_t BOTTOM = 0b0100;
static constexpr uint8_t TOP = 0b1000;

uint8_t Renderer::ComputeOutCode(float x, float y)
{
    uint8_t code = INSIDE;

    if (x < 0)
        code |= LEFT;
    else if (x > window->width - 1)
        code |= RIGHT;
    if (y < 0)
        code |= TOP;
    else if (y > window->height - 1)
        code |= BOTTOM;

    return code;
}

// Cohen-sutherland 直线裁切算法
bool Renderer::CohenSutherlandLineClip(float &x0, float &y0, float &x1, float &y1, const float xmin, const float xmax, const float ymin, const float ymax)
{
    uint8_t outcode0 = ComputeOutCode(x0, y0);
    uint8_t outcode1 = ComputeOutCode(x1, y1);

    while (true)
    {
        if (!(outcode0 | outcode1))
        {
            return true;
        }
        else if (outcode0 & outcode1)
        {
            return false;
        }
        else
        {
            uint8_t &outcodeOut = outcode1 > outcode0 ? outcode1 : outcode0;
            float &x = outcode1 > outcode0 ? x1 : x0;
            float &y = outcode1 > outcode0 ? y1 : y0;
            if (outcodeOut & BOTTOM)
            {
                x = x0 + (x1 - x0) * (ymax - y0) / (y1 - y0);
                y = ymax;
            }
            else if (outcodeOut & TOP)
            {
                x = x0 + (x1 - x0) * (ymin - y0) / (y1 - y0);
                y = ymin;
            }
            else if (outcodeOut & RIGHT)
            {
                y = y0 + (y1 - y0) * (xmax - x0) / (x1 - x0);
                x = xmax;
            }
            else if (outcodeOut & LEFT)
            {
                y = y0 + (y1 - y0) * (xmin - x0) / (x1 - x0);
                x = xmin;
            }
            outcodeOut = ComputeOutCode(x, y);
        }
    }
    return false;
}

顶点变换

模型变换

模型变换(Model Transformation)通过模型矩阵表现。模型矩阵包括SRT（缩放Scale，旋转Rotate，位移Translate）矩阵，用于将顶点从物体的局部坐标系变换到世界坐标系中。

引入齐次坐标的一个好处是支持了仿射变换，平移、旋转、缩放、剪切、反射这些变换都属于仿射变换。仿射变换即引入一个额外坐标，用3维向量和矩阵表示2维空间的变换，用4维向量和矩阵表示3维空间的变换。

SRT的顺序很重要，因为矩阵相乘没有交换律，顺序变化会导致变换结果完全不同。一般先缩放，再旋转，最后平移，这样是符合逻辑的，因为在最初的局部坐标系中，一般以物体中心作为原点，这样进行缩放和旋转不会大幅影响物体整体的位置。在具体实现时，如果将向量都认为是列向量，那么矩阵必须在向量的左边相乘，因此矩阵运算的顺序是从右往左，故应该写成位移矩阵*旋转矩阵*缩放矩阵*顶点向量的形式。

观察变换

观察变换(View Transformation)通过观察矩阵表现。观察矩阵将世界坐标系变换到观察空间，观察空间中我们约定俗成地将相机置于原点，朝向-z方向看去，并且上方向为+y方向。

世界坐标系中我们定义的相机，一般位置不在原点，或者方向不符，这时如何构造这个观察矩阵呢？

观察矩阵分为两部分，一部分为位移矩阵，将物体与相机位移到相机在原点的位置。已知相机的位置$\mathbf{P}=(t_x,t_y,t_z)$，相当于相机从原点移到世坐标系的位置，则位移矩阵反向位移即可，因此：

$$T_{view} = \begin{bmatrix} 1 & 0 & 0 & -t_x\\ 0 & 1 & 0 & -t_y\\ 0 & 0 & 1 & -t_z\\ 0 & 0 & 0 & 1 \end{bmatrix}$$

而另一部分旋转矩阵，涉及到我们如何定义相机的方向，常见有两种定义方式。第一种是用欧拉角表示，并且我们给3个欧拉角起了名字：俯仰角（Pitch）、偏航角（Yaw）、滚转角（Roll）。这组欧拉角表示从初始方向，先绕y轴旋转（Yaw），再绕x轴旋转（Pitch），最后绕z轴旋转（Roll），得到世界坐标系中的方向。这3个旋转可以组合成一个旋转矩阵，对其求逆就将方向从世界坐标系转到观察空间了。由于旋转矩阵是正交矩阵，其逆等于其转置，可以简化计算。即：

$$R_{view} = R^{-1} = (R_{roll}R_{pitch}R_{yaw})^{-1} = R_{yaw}^TR_{pitch}^TR_{roll}^T$$

$$V = R_{view}T_{view}$$

第二种用世界坐标系中相机看向的方向（或物体），以及上方向来表示。若用target表示摄像机看向的目标点，up表示世界坐标系中摄像机的上方向，则我们可以计算前、右、上3个基向量。（注意是右手坐标系）

前向向量：

$$\mathbf{F}=\text { normalize }(\text { eye }- \text { target })$$

右向量：

$$\mathbf{R}=\text{normalize}( \mathbf{F}\times \operatorname{up})$$

上向量（为确保垂直，要重新计算）：

$$\mathbf{U}=\text{normalize}(\mathbf{R}\times\mathbf{F})$$

考虑将观察空间的相机方向变换到世界坐标系中，也就是将右方向从$\begin{bmatrix}1\\0\\0\end{bmatrix}$变换到 $\mathbf{R}$，将上方向从$\begin{bmatrix}0\\1\\0\end{bmatrix}$变换到 $\mathbf{U}$，将前方向从$\begin{bmatrix}0\\0\\-1\end{bmatrix}$（因为相机朝向-z）变换到 $\mathbf{F}$。我们可以构造方程：

$$R_r \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & -1\\ \end{bmatrix} = \begin{bmatrix} \mathbf{R}_x & \mathbf{U}_x & \mathbf{F}_x\\ \mathbf{R}_y & \mathbf{U}_y & \mathbf{F}_y\\ \mathbf{R}_z & \mathbf{U}_z & \mathbf{F}_z\\ \end{bmatrix}$$

并解得：

$$R_r = \begin{bmatrix} \mathbf{R}_x & \mathbf{U}_x & -\mathbf{F}_x\\ \mathbf{R}_y & \mathbf{U}_y & -\mathbf{F}_y\\ \mathbf{R}_z & \mathbf{U}_z & -\mathbf{F}_z\\ \end{bmatrix}$$

或者我们可以等效认为，将相机后方向从$\begin{bmatrix}0\\0\\1\end{bmatrix}$变换到$-\mathbf{F}$，这样$\mathbf{R},\mathbf{U},-\mathbf{F}$就是世界坐标系的一组基向量，组成的矩阵$R_r$就表示该变换。

而我们需要的是从世界坐标系变换到观察空间，因此对该矩阵求逆，由于$R_r$是正交矩阵，可以求转置来求逆。

$$R_r^{-1}=R_r^T=\begin{bmatrix} \mathbf{R}_x & \mathbf{R}_y & \mathbf{R}_z\\ \mathbf{U}_x & \mathbf{U}_y & \mathbf{U}_z\\ -\mathbf{F}_x & -\mathbf{F}_y & -\mathbf{F}_z\\ \end{bmatrix}$$

最后与位移矩阵组合，不过记得要将旋转矩阵扩展为4维的形式。

$$V = R_{view}T_{view} = \begin{bmatrix}R_r^{-1}&0\\0&1\end{bmatrix}T_{view} =\begin{bmatrix} \mathbf{R}_x & \mathbf{R}_y & \mathbf{R}_z&0\\ \mathbf{U}_x & \mathbf{U}_y & \mathbf{U}_z&0\\ -\mathbf{F}_x & -\mathbf{F}_y & -\mathbf{F}_z&0\\ 0&0&0&1 \end{bmatrix}\begin{bmatrix} 1 & 0 & 0 & -t_x\\ 0 & 1 & 0 & -t_y\\ 0 & 0 & 1 & -t_z\\ 0 & 0 & 0 & 1 \end{bmatrix}\\ =\begin{bmatrix} \mathbf{R}_x & \mathbf{R}_y & \mathbf{R}_z&-t_x\mathbf{R}_x-t_y\mathbf{R}_y-t_z\mathbf{R}_z\\ \mathbf{U}_x & \mathbf{U}_y & \mathbf{U}_z&-t_x\mathbf{U}_x-t_y\mathbf{U}_y-t_z\mathbf{U}_z\\ -\mathbf{F}_x & -\mathbf{F}_y & -\mathbf{F}_z&-t_x\mathbf{F}_x-t_y\mathbf{F}_y-t_z\mathbf{F}_z\\ 0&0&0&1 \end{bmatrix}=\begin{bmatrix} \mathbf{R}_x & \mathbf{R}_y & \mathbf{R}_z&-\mathbf{P}\cdot\mathbf{R}\\ \mathbf{U}_x & \mathbf{U}_y & \mathbf{U}_z&-\mathbf{P}\cdot\mathbf{U}\\ -\mathbf{F}_x & -\mathbf{F}_y & -\mathbf{F}_z&\mathbf{P}\cdot\mathbf{F}\\ 0&0&0&1 \end{bmatrix}$$

投影变换

投影本来是指将3维空间的物体映射到2D屏幕上的过程，在经过观察变换之后，理论上可见的物体都在相机面对的方向(正对-z)，只要取该方向上的某个平面，将其后的一部分空间内的顶点映射到该平面，只保留x和y的坐标。但是由于z坐标(深度)还要参与后续的插值和深度测试，以及还有其他用处，因此我们需要同样要映射z值。OpenGL规范是以$x\in [-1,1],y\in[-1,]\,z\in[-1,1]$的归一化设备坐标(NDC)空间作为目标。

实际上，从观察空间到NDC空间经过了两步：第一步投影变换将观察空间转换为(齐次)裁剪空间，第二步透视除法将裁剪空间转换为NDC空间。

平行投影

平行投影将观察空间中指定立方体的区域直接映射到$x\in [-1,1],y\in[-1,]\,z\in[-1,1]$的范围内。在形式上，平行投影完成后已经将顶点变换至NDC空间了，但是为了与透视投影统一，之后还是要进行透视除法（$\omega=1$）。

平行投影多用于2D游戏，2.5D游戏，工程图等，因为不符合人眼透视，但保证各处物体的比例相同。

指定要投影的立方体区间x方向宽为$w$，y方向高为$h$，z方向$[-n,-f]$(因为相机看向-z方向)。那么我们需要完成$x:[-\frac{w}{2},\frac{w}{2}]\rarr[-1,1]$，$y:[-\frac{h}{2},\frac{h}{2}]\rarr[-1,1]$，$z:[-n, -f]\rarr[-1,1]$的映射，其中z的映射可以先将$[-n,-f]$加上$\frac{n+f}{2}$移动到$[\frac{f-n}{2},\frac{n-f}{2}]$，再进行缩放。通过简单的计算可以得出对应的矩阵为：

$$P_{ortho}=\begin{bmatrix} \frac{2}{w} & 0 & 0 & 0 \\ 0 & \frac{2}{h} & 0 & 0 \\ 0 & 0 & \frac{2}{f-n} & \frac{f+n}{f-n} \\ 0 & 0 & 0 & 1 \end{bmatrix}$$

透视投影

透视投影是符合人眼直觉，因此更加广泛使用。它将观察空间中的一个视锥(Frustum)，即一个四棱台体的空间映射到(齐次)裁剪空间内。但是这样的概念过于抽象，我们不妨一步步考虑。

首先，由于符合人眼直觉，假设最后投影的结果呈现在一张矩形的画布上，那么在相机方向的物体，光线能进入相机的部分位于一个锥体中。这个矩形画布的比例用宽高比asp描述，而后方多大角度范围的物体进入相机，则用视角(Field of View,FOV)描述。有这两个参数，我们确定了一个从相机出发，向-z方向无限延伸的四棱锥。我们将矩形画布置于该方向上的一处后，将画布后、四棱锥内的物体向前投影到画布上，这就是最基本的透视投影的效果。

接下来通过几何推导，如何实现这样的投影。

画出观察空间沿-x方向看去的样子，将画布置于距相机(原点)$n$的位置，也就是“近平面”位置。张角的大小由视角决定，以垂直FOV为例，它表示相机纵向的张角大小。
接下来，取近平面后、锥体内任意一点$M(x,y,z)$，要计算它投影在画布上的$x',y'$值。构造相似三角形，易得：

$$\frac{x'}{n}=\frac{x}{-z}$$

即：

$$x'=-\frac{nx}{z}$$

同理可得y坐标：

$$y'=-\frac{ny}{z}$$

到此，理论上我们已经完成了投影的一大半，就是将点$M(x,y,z)$投影到近平面的$M'(-\frac{nx}{z},-\frac{ny}{z},-n)$，注意这里所有点的z值就无法区分了，都是$-n$，但我们先不管z值。只剩将$x',y'$的值映射到$[-1,1]$的区间了。由图中的几何关系得$y'\in[-n\tan\frac{FOV}{2},n\tan\frac{FOV}{2}]$，又由宽高比asp得到$x'$与$y'$的关系$\frac{x'}{y'}=\text {asp}$，故$x'\in[-\text{asp}\cdot n\tan\frac{FOV}{2},-\text{asp}\cdot n\tan\frac{FOV}{2}]$，将这两个范围映射到$[-1,1]$与平行投影类似，不再赘述。

最重要的一点是，这个投影变换要通过矩阵乘法完成，可是我们上述的变换能通过矩阵做到吗？后半部分的映射不必多说，可以通过$\begin{bmatrix} \frac{1}{\text{asp}\cdot n\tan\frac{FOV}{2}} & 0 & 0 & 0 \\ 0 & \frac{1}{n\tan\frac{FOV}{2}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}$该矩阵完成，但前半部分的变换，涉及到除以原向量的数，无法通过矩阵完成。这时就需要齐次坐标发挥用处了，引入$w$分量后，原来的目标点$M'(-\frac{nx}{z},-\frac{ny}{z},-n,1)$等效于$(nx,ny,nz,-z)$，而从$M(x,y,z,1)$到$(nx,ny,nz,-z)$的变换就可以用矩阵乘法完成，乘以$\begin{bmatrix} n & 0 & 0 & 0 \\ 0 & n & 0 & 0 \\ 0 & 0 & n & 0 \\ 0 & 0 & -1 & 0 \end{bmatrix}$该矩阵来实现。（z可以先不管）

这里解释了为什么从观察空间到NDC空间不能由透视投影一步完成，因为透视投影做的所有事情仅仅是用透视矩阵对顶点做一次乘法，而由齐次坐标转换为3维空间的点则需要后面的透视除法完成。

再来看z值，事实上z值的映射与x,y无关，如果是软光栅化渲染器，甚至可以保留观察空间中的z值，在透视除法中跳过z值，并直接用于后续的透视矫正插值和深度测试。但是这样会导致近处精度浪费，远处精度不足，从而导致Z-fighting，究其原因，是线性分布的z值不符合渲染的要求。因此我们要对z值进行非线性的映射，出于硬件优化等的考量，我们一般将z值映射到$[-1,1]$(OpenGL)或$[0,1]$(DirectX，Valkan)，这里以OpenGL为例。若将近平面映射到-1，无限远处映射到1，这是可以做到的，那么为什么我们一般还要引入“远平面”的概念呢？近远平面之间的锥体才称为视锥，因为这样的映射方便进行裁剪，将极近与极远处的物体舍弃掉，这是完全合理的，因为极远处的物体即使渲染了，也占据不了一个像素，因而无法呈现，不如直接舍弃来优化性能。

定义远平面，与相机的距离为$f$，接下来考虑$z:[-n,-f]\rarr[-1,1]$的映射。这里的映射是非线性的，因而不能仿照平行投影，同时别忘了，我们要在矩阵乘法中完成这个映射。根据矩阵已经确定的部分，列出方程：

$$\begin{bmatrix} n & 0 & 0 & 0 \\ 0 & n & 0 & 0 \\ A & B & C & D \\ 0 & 0 & -1 & 0 \end{bmatrix}\begin{bmatrix} x \\ y \\ z \\ 1 \end{bmatrix}=\begin{bmatrix} nx \\ ny \\ \text{unknown} \\ -z \end{bmatrix}\lrArr\begin{bmatrix} x' \\ y' \\ z' \\ 1 \end{bmatrix}$$

可得：

$$Ax+By+Cz+D=z'(-z)$$

我们虽不知道映射的具体过程，但是可以确定两点：

1. 顶点的x与y值对z的映射没有影响，也就是对于一个确定的$z$，任意的$x,y$都能使上述方程成立，易得：$A=B=0$。
2. 区间端点的映射关系已知，即$-n\rarr-1,-f\rarr1$，将$z$的这两个取值代入：

$$-Cn+D=-n\\ -Cf+D=f$$

解得：

$$C=\frac{n+f}{n-f}\\$$

$$D=\frac{2nf}{n-f}$$

因此完整的投影矩阵为：

$$P_{proj}=\begin{bmatrix} \frac{1}{\text{asp}\cdot n\tan\frac{FOV}{2}} & 0 & 0 & 0 \\ 0 & \frac{1}{n\tan\frac{FOV}{2}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} n & 0 & 0 & 0 \\ 0 & n & 0 & 0 \\ 0 & 0 & \frac{n+f}{n-f} & \frac{2nf}{n-f} \\ 0 & 0 & -1 & 0 \end{bmatrix}=\begin{bmatrix} \frac{1}{\text{asp}\cdot\tan\frac{FOV}{2}} & 0 & 0 & 0 \\ 0 & \frac{1}{\tan\frac{FOV}{2}} & 0 & 0 \\ 0 & 0 & \frac{n+f}{n-f} & \frac{2nf}{n-f} \\ 0 & 0 & -1 & 0 \end{bmatrix}$$

不同的标准可能导致各处有小差异，但是这样的推导过程是最重要的。

经过投影变换和透视除法这两步后，顶点的实际坐标的含义：x和y值对应投影到屏幕上的某点，z值为标准化后的位于$[0,1]$或$[-1,1]$的值，用于深度测试，而$w=-z_{view}$保留了在观察空间的真实z值，用于之后的透视矫正插值。

视口变换

视口变换将NDC空间转换为屏幕空间。这一步很简单，保留各顶点的z值和w值，进行映射$x:[-1,1]\rarr[0,\text{width}-1]$，$y:[-1,1]\rarr[0,\text{height}-1]$。需要注意的是，屏幕空间的原点一般在左上，对比NDC空间（向-z方向），需要翻转y轴。如果是CPU软光栅化，可以直接对各顶点的x,y进行运算；如果是GPU，可以写成矩阵形式。

---

到目前为止，我们已经可以用“线框模式”绘制一些基本图元了。在完成以上一系列顶点变换之后，以三角形为例，只要每两个顶点间调用上面的直线绘制算法即可。但如果想要填充三角形，必然涉及对三角形内部的逐像素操作，这就是光栅化。

光栅化

光栅化是在屏幕空间进行的。扫描线算法是更适合CPU软光栅化的算法，但是这里只提及更简单的GPU光栅化方法。

对于一个三角形，找到其3个顶点的轴对齐包围盒(AABB,Axis-Aligned Bounding Box)，也就是x,y的最值作为边界，确定一个矩形，遍历该矩形内的像素，并判断该像素是否在三角形内，若是则按照一定规则绘制该像素。

判断像素是否在三角形内可以通过几个向量的叉乘得到，将顶点分别与三个顶点组成向量，并与对应的三角形边的向量做叉乘，对应三角形边的向量必须是同一个走向首尾相连，这样就可以用叉乘判断该点是否在三边同侧（同正或同负），同侧即在三角形内。

顶点属性插值

三角形内的像素具有什么颜色，是由其三个顶点决定的，我们可以通过插值让三个顶点的属性在三角形内部平滑过渡。

对比线性插值只需一个参数$t$，就可以用$C=tA+(1-t)B$进行插值，三角形的顶点插值则需要3个参数，即重心坐标。数学定义上，一个点$P$重心坐标$\alpha,\beta,\gamma$满足$P=\alpha A+\beta B+\gamma C$，且$\alpha+\beta+\gamma=1$。而通常的计算方法是用该点分割的子三角形与整个三角形的面积比，即：

$$\alpha=\frac{S_{\triangle PBC}}{S_{\triangle ABC}},\ \beta=\frac{S_{\triangle PAC}}{S_{\triangle ABC}},\ \gamma=\frac{S_{\triangle PAB}}{S_{\triangle ABC}}$$

接着对于任意一个属性$I$，用$I_P=\alpha I_A+\beta I_B+\gamma I_C$即可得到该点属性值。

但是这样在透视情况下，并不能得到正确的视觉效果。

透视矫正插值

因为上述插值假定了三角形是“正对”着我们的，有了透视变换，屏幕空间的三角形并不都符合这点。实际的插值应该在透视投影前的世界空间内完成，但是光栅化只能在屏幕空间做，因此我们有另一种解决方法——透视矫正插值。

示意图如下：

其中$A,B,C$点位于观察空间内，由$A,B$线性插值得到C是正确插值；而$A',B',C'$是透视投影后的结果，可以认为是在NDC空间或者屏幕空间内，用$A',B'$插值得到$C'$，若使用s参数插值就是错的，$C'$的实际属性应该与$C$相同(用t参数插值)，可以明显看出两个插值的区别。

但在光栅化阶段，我们有的全部信息是$A',B',s$，以及齐次坐标的$w$保留的，观察空间的$A,B$的真实z值$Z_1,Z_2$。我们要通过这些信息推导出正确的插值参数t。

由相似三角形，易得：

$$\frac{t}{1-t}=\frac{AM}{BN}$$

$$\frac{AM}{s}=\frac{-Z_1}{d}$$

$$\frac{BN}{1-s}=\frac{-Z_2}{d}$$

因此：

$$\frac{AM}{BN}=\frac{sZ_1}{(1-s)Z_2}=\frac{t}{1-t}$$

$$\frac{1-t}{t}=\frac{1}{t}-1=\frac{(1-s)Z_2}{sZ_1}$$

$$\frac{1}{t}=\frac{sZ_1+(1-s)Z_2}{sZ_1}$$

$$t=\frac{sZ_1}{sZ_1+(1-s)Z_2}$$

解出t后，对于任意属性$I$，正确的插值如下：

$$I_t=(1-t)I_1+tI_2=\frac{(1-s)Z_2}{sZ_1+(1-s)Z_2}I_1+\frac{sZ_1}{sZ_1+(1-s)Z_2}I_2$$

当然这个式子显得有些复杂，我们可以用正确插值的z值来简化它。将$Z_1,Z_2$带入$I_1,I_2$得：

$$Z_t=\frac{Z_1Z_2}{sZ_1+(1-s)Z_2}$$

即：

$$\frac{1}{Z_t}=(1-s)\frac{1}{Z_1}+s\frac{1}{Z_2}$$

带入上述对于$I$的插值：

$$I_t=(1-s)\frac{Z_t}{Z_1}I_1+s\frac{Z_t}{Z_2}I_2=Z_t(\frac{1-s}{Z_1}I_1+\frac{s}{Z_2}I_2)$$

以上是在$yOz$平面的推导，推广至三角形内：

$$I_t=Z_t(\frac{\alpha}{Z_1}I_1+\frac{\beta}{Z_2}I_2+\frac{\gamma}{Z_3}I_3)$$

其中：

$$\frac{1}{Z_t}=\alpha\frac{1}{Z_1}+\beta\frac{1}{Z_2}+\gamma\frac{1}{Z_3}$$

这就是三维空间内的透视正确的插值方式。

虽然还没谈到纹理，在运用纹理时，对于纹理坐标的插值，透视矫正的效果很明显，如下图所示：（左边是没有用透视矫正插值，右边用了）

在实际应用中，由于已知的是并不是直接的$Z$，而是顶点的$w$（符号相反），因此我们可以改写上式：

$$I_t=w_t(\alpha\frac{I_1}{w_1}+\beta\frac{I_2}{w_2}+\gamma\frac{I_3}{w_3})$$

其中

$$\frac{1}{w_t}=\frac{\alpha}{w_1}+\frac{\beta}{w_2}+\frac{\gamma}{w_3}$$

在GPU渲染管线中，GPU在光栅化前将所有需要插值的顶点属性都除以w值，在光栅化时对这些值进行插值，得到的结果再除以$\frac{1}{w_t}$，这里的$w_t$本身也是透视矫正插值得到的。

---

结合以上，在软渲染器中绘制一个三角形的代码大致如下：

void Renderer::DrawTriangle(std::vector<Vertex> vertices, RenderMode mode)
{
    mat4f model = mat4f::GetIdentity();
    for (auto &vertex : vertices)
    {
        // MVP变换
        vertex.position = camera.GetProjectionMatrix() * camera.GetViewMatrix() * model * vertex.position;
        // 透视除法
        vertex.position.x /= vertex.position.w;
        vertex.position.y /= vertex.position.w;
        vertex.position.z /= vertex.position.w;
    }
    // 视口变换
    for (auto &vertex : vertices)
    {
        vertex.position.x = (vertex.position.x + 1.0f) * 0.5f * (window->width - 1.0f);
        vertex.position.y = window->height - (vertex.position.y + 1.0f) * 0.5f * (window->height - 1.0f);
    }
    if (mode == WIREFRAME) // 线框模式绘制
    {
        for (int i = 0; i < 3; i++)
            DrawLine(vertices[i % 3].position.x, vertices[i % 3].position.y, vertices[(i + 1) % 3].position.x, vertices[(i + 1) % 3].position.y, 0, 0, 0);
    }
    else if (mode == FILL) // 填充模式，需要光栅化
    {
        float minx = window->width - 1;
        float maxx = 0.0f;
        float miny = window->height - 1;
        float maxy = 0.0f;
        for (int i = 0; i < 3; i++)
        {
            minx = std::min(minx, vertices[i].position.x);
            miny = std::min(miny, vertices[i].position.y);
            maxx = std::max(maxx, vertices[i].position.x);
            maxy = std::max(maxy, vertices[i].position.y);
        }
        minx = std::max(0.0f, minx);
        maxx = std::min((float)(window->width - 1), maxx);
        miny = std::max(0.0f, miny);
        maxy = std::min((float)(window->height - 1), maxy); // 确定AABB包围盒
        for (int i = (int)minx; i <= (int)maxx; i++)
            for (int j = (int)miny; j <= (int)maxy; j++) // 逐像素
            {
                vec2f p(i + 0.5f, j + 0.5f);
                vec2f s0 = vertices[0].position.xy() - p;
                vec2f s1 = vertices[1].position.xy() - p;
                vec2f s2 = vertices[2].position.xy() - p;
                float Sa = vector_cross(s1, s2);
                float Sb = vector_cross(s2, s0);
                float Sc = vector_cross(s0, s1);
                if (!((Sa >= 0 && Sb >= 0 && Sc >= 0) || (Sa <= 0 && Sb <= 0 && Sc <= 0))) continue;
                float S = Sa + Sb + Sc;
                float alpha = Sa / S, beta = Sb / S, gamma = Sc / S; // 重心坐标
                float reverse_wt = alpha / vertices[0].position.w + beta / vertices[1].position.w + gamma / vertices[2].position.w;
                vec4f color = (alpha * vertices[0].color / vertices[0].position.w 
                            + beta * vertices[1].color / vertices[1].position.w 
                            + gamma * vertices[2].color / vertices[2].position.w) 
                            / reverse_wt; // 没有将属性提前除以w，在这里做
                SetPixel(i, j, color.r * 255, color.g * 255, color.b * 255);
            }
    }
}

---

参考资料：

• 单身剑法传人-软渲染
• 知乎回答:怎么理解三维投影中视锥体的近平面和远平面？