1. 前沿

数学上直线是连续的,由无穷多点构成;但在计算机显示中,屏幕是离散的像素网格,只能用整数坐标的像素点去逼近连续直线。因此需要算法来确定哪些整数像素点最能代表这条直线,这就是直线绘制算法的核心作用

算法的输入是起点和终点的坐标,算法的输出是两点连成的直线经过的像素网格坐标序列

2. DDA算法

全称为数字差分分析器(Digital Differential Analyzer)算法

2.1 思路

直线的斜截式方程微分形式为

Δy=kΔx\Delta y = k \Delta x

那么有

{xi+1=xi+Δxyi+1=yi+Δy=yi+kΔx\left\{ \begin{align*} & x_{i+1} = x_i + \Delta x \\ & y_{i+1} = y_i + \Delta y = y_i + k \Delta x \end{align*} \right.

当斜率为0k<10 \leq k <1,且Δx\Delta xΔy\Delta y全为正时

{xi+1=xi+1yi+1=yi+k\left\{ \begin{align*} & x_{i+1} = x_i + 1 \\ & y_{i+1} = y_i + k \end{align*} \right.

由于像素点坐标是整数,因此yy值要四舍五入

2.4 代码参考

首先定义像素对象

class Pixel {
public:
    explicit Pixel() = default;
    Pixel(const int x_, const int y_) : x{x_}, y{y_} {}
    int x = 0;
    int y = 0;
};

简陋版代码(此代码只适合在Δx\Delta xΔy\Delta y全为正,且k0<1k \leq 0 < 1时使用)

std::vector<Pixel> DdaSimple(const int x0, const int y0, const int xn,
                             const int yn) {
    std::vector<Pixel> pixels;

    const int dx = xn - x0;
    const int dy = yn - y0;

    const double k = static_cast<double>(dy) / static_cast<double>(dx);
    for (int i = 0; i <= dx; ++i) {
        pixels.emplace_back(x0 + i, y0 + i * k);
    }
    return pixels;
}

完善版代码(考虑直线的各个方向)

std::vector<Pixel> Dda(int x0, int y0, int xn, int yn) {
    std::vector<Pixel> pixels;

    // 1. 计算x、y方向的总增量
    int dx = xn - x0;
    int dy = yn - y0;

    // 2. 确定迭代主方向
    int steps = std::max(std::abs(dx), std::abs(dy));

    // 处理特殊情况:起点=终点(无需迭代)
    if (steps == 0) {
        pixels.emplace_back(x0, y0);
        return pixels;
    }

    // 3. 计算每一步的x、y增量(浮点数,均匀步进)
    float delta_x = static_cast<float>(dx) / steps;
    float delta_y = static_cast<float>(dy) / steps;

    // 4. 初始化当前坐标(从起点开始)
    float current_x = x0;
    float current_y = y0;

    // 5. 迭代生成像素点(包含起点和终点)
    for (int i = 0; i <= steps; ++i) {
        // 四舍五入取整
        int pixel_x = static_cast<int>(std::round(current_x));
        int pixel_y = static_cast<int>(std::round(current_y));
        pixels.emplace_back(pixel_x, pixel_y);

        // 累加增量,更新当前坐标
        current_x += delta_x;
        current_y += delta_y;
    }

    return pixels;
}

3. Bresenham算法

DDA算法需要进行浮点运算,会拖慢计算速率,而此算法只进行整数运算

3.1 简单思路

假设以x轴为主方向,那么下一个像素要么是右侧的(xi+1,yi)(x_{i+1},y_{i})要么是右上的(xi+1,yi+1)(x_{i+1},y_{i+1}),判断的依据是理想直线交点离哪个更近。使用dud_u表示yi+1y_{i+1}和理想直线交点的偏差,dld_l表示yiy_{i}和理想直线交点的偏差

{du=yi+1k(xi+1)bdl=k(xi+1)+byi\left\{ \begin{align*} & d_u = y_{i+1} - k(x_i + 1) - b \\ & d_l = k(x_i + 1) +b - y_i \end{align*} \right.

计算两个值的偏差dldudl − du

dldu=2kxi2yi+2k+2b1d_l − d_u = 2k x_i − 2y_i + 2k + 2b − 1

如果偏差小于00说明交点更接近yiy_i,应该选择yiy_i;如果偏差大于00说明交点更接近yi+1y_{i + 1},应该选择yi+1y_{i + 1};如果偏差为00选择哪个都行

由于这个偏差的计算包含了浮点kk运算,为了进一步提升计算效率,接下来还需要进行转换使其只需要整数运算

上式可以写为

dldu=2ΔyΔxxi2yi+2ΔyΔx+2b1d_l − d_u = 2 \frac{\Delta y}{\Delta x} x_i − 2y_i + 2\frac{\Delta y}{\Delta x} + 2b − 1

定义决策变量pi=Δx(dldu)p_i = \Delta x(d_l − d_u),上面的公式可变换为

pi=2Δyxi2Δxyi+2Δy+Δx(2b1)p_i = 2 \Delta y x_i - 2 \Delta x y_i + 2 \Delta y + \Delta x(2b-1)

因为Δx\Delta x大于00,因此pip_i的正负号和dldud_l − d_u是一样的,可以使用pip_i的正负号作为第kk步的决策参数来确定下一个像素的坐标

pip_i的初始值p0p_0的初始值提取出来以提高计算效率

p0=2Δyx02Δxy0+2Δy+Δx(2b1)p0=2Δyx02Δx(ΔyΔxx0+b)+2Δy+Δx(2b1)p0=2ΔyΔx\begin{align*} p_0 &= 2 \Delta y x_0 - 2 \Delta x y_0 + 2 \Delta y + \Delta x(2b-1) \\ p_0 &= 2 \Delta y x_0 - 2 \Delta x (\frac{\Delta y}{\Delta x} x_0+b) + 2 \Delta y + \Delta x(2b-1) \\ p_0 &= 2 \Delta y - \Delta x \end{align*}

还需要确定怎么迭代计算pi+1p_{i+1}

pi+1=2Δyxi+12Δxynext+2Δy+Δx(2b1)pi+1=2Δyxi+2Δy2Δxynext+2Δy+Δx(2b1)pi+1=2Δyxi2Δxyi+2Δxyi2Δxynext+2Δy+Δx(2b1)+2Δypi+1=pi+2Δxyi2Δxynext+2Δypi+1=2Δy2Δx(ynextyi)\begin{align*} p_{i+1} &= 2 \Delta y x_{i+1} - 2 \Delta x y_{\text{next}} + 2 \Delta y + \Delta x(2b-1) \\ p_{i+1} &= 2 \Delta y x_i + 2 \Delta y - 2 \Delta x y_{\text{next}} + 2 \Delta y + \Delta x(2b-1) \\ p_{i+1} &= 2 \Delta y x_i - 2 \Delta x y_i + 2 \Delta x y_i - 2 \Delta x y_{\text{next}} + 2 \Delta y + \Delta x(2b-1) + 2 \Delta y \\ p_{i+1} &= p_i + 2 \Delta x y_i - 2 \Delta x y_{\text{next}} + 2 \Delta y \\ p_{i+1} &= 2 \Delta y - 2 \Delta x (y_{\text{next}} -y_i) \end{align*}

其中ynexty_{\text{next}}表示第ii步的下一个像素的yy坐标,为yiy_iyi+1y_i + 1,根据前面的内容可知它由pip_i决定

pi+1={pi+2Δy2Δx,pi0pi+2Δy,pi<0p_{i+1} = \left\{ \begin{align*} & p_i + 2 \Delta y - 2 \Delta x,p_i \geq 0\\ & p_i + 2 \Delta y, p_i < 0 \end{align*} \right.

上面的公式中,Δy\Delta yΔx\Delta x是直线端点x,yx,y轴变化量,是已知且固定不变的,因此2Δy2\Delta y2Δx2\Delta x都可以一次性计算好,后续的递推运算只需要简单的整数加减运算就可以了

3.2 进阶思路

上述方法需要分情况讨论,但是实际中需要将所有 8 个象限的情况合并成一个统一的循环逻辑

建立隐式方程

直线的隐式方程为F(x,y)=0F(x,y)=0。我们可以利用两点式方程变形得到(假设都在第一象限,且增量取绝对值 Δx,Δy):

xx0Δx=yy0Δy\frac{x-x_0}{\Delta x} = \frac{y-y_0}{\Delta y}

交叉相乘得到直线的隐式方程:

(yy0)Δx(xx0)Δy=0(y-y_0) \Delta x - (x-x_0) \Delta y = 0

我们要定义的误差项errerr,其物理意义就是当前像素点(xi,yi)(x_i,y_i)代入该方程后的残差值

erri=(yiy0)Δx(xix0)Δyerr_i = (y_i-y_0) \Delta x - (x_i-x_0) \Delta y

使用残差值的正负值可以判断点在直线的哪一侧。本质上是向量(Δx,Δy)(\Delta x,\Delta y)(xix0,yiy0)(x_i-x_0,y_i-y_0)的叉积

误差的迭代更新

观察errierr_i可知

  • 当我们在x方向前进一步,误差项减小Δy\Delta y
  • 当我们在y方向前进一步,误差项增加Δx\Delta x

初始误差的偏移

第一个点误差是00,所以将初始误差定为(x1,y1)(x_1,y_1)点的误差,即ΔxΔy\Delta x - \Delta y

但是我们不直接比较两点,而是比较中点的位置,即中点是在直线的左侧还是右侧。即判断ΔxΔy+0.5Δy\Delta x - \Delta y+0.5\Delta y是否大于零;ΔxΔy0.5Δx\Delta x - \Delta y-0.5\Delta x是否小于零;为了消除浮点值,将式子全乘以二,于是得到了2erri+Δy2err_i + \Delta y2erriΔx2err_i - \Delta x

3.3 代码参考

vector<Pixel> bresenhamLine(int x0, int y0, int xn, int yn) {
    vector<Pixel> pixels;
    // 1. 计算x/y方向总增量(带符号,体现方向)
    int dx = xn - x0;
    int dy = yn - y0;

    // 2. 确定x/y的步进方向(1=正方向,-1=反方向)
    int sx = (dx > 0) ? 1 : (dx < 0) ? -1 : 0;
    int sy = (dy > 0) ? 1 : (dy < 0) ? -1 : 0;

    // 3. 取增量绝对值(用于误差项判断,纯整数运算)
    int abs_dx = abs(dx);
    int abs_dy = abs(dy);

    // 4. 初始化核心误差项(Bresenham核心,无浮点数)
    int err = abs_dx - abs_dy;

    // 5. 迭代生成像素点(当前坐标从起点开始)
    int curr_x = x0;
    int curr_y = y0;
    while (true) {
        pixels.emplace_back(curr_x, curr_y);      // 记录当前像素(包含起点)
        if (curr_x == xn && curr_y == yn) break;  // 到达终点,退出迭代

        int err2 = 2 * err;  // 误差项*2,消去浮点数(经典优化)
        // 误差项判断:x方向步进条件
        if (err2 > -abs_dy) {
            err -= abs_dy;
            curr_x += sx;
        }
        // 误差项判断:y方向步进条件(注意是if而非else if,兼容斜率=1/-1)
        if (err2 < abs_dx) {
            err += abs_dx;
            curr_y += sy;
        }
    }
    return pixels;
}

参考文章

【图形学】直线光栅化算法(DDA算法,Bresenham算法,中点算法) - 钓一朵雪的文章 - 知乎

一文读懂Bresenham画线算法 - 一点鑫得的文章 - 知乎

Bresenham 画线算法