用 C 语言画直线

Milo Yip

2016 年度新知答主

画直线可能是很多人接触图形学的第一个问题，例如：

我尝试用不同技术实现画直线的方法（完整源代码在 miloyip/line ），此文简单介绍个中思路。本文的代码采用 C 语言、标准库及极简的 PNG 编码函数 svpng() ，没有使用其他 API。

1. Bresenham 算法

Bresenham直线算法 [1] 是最简单的直线光栅化（rasterization）算法。

如果像上图，直线的高度小于宽度，那么 Bresenham 直线算法会为 x 轴每个坐标填入一个像素，绘画每个像素时按斜率判断 y 是否需要调整。整个算法可以避开浮点数运算，只用整数运算实现。以下是一个简单实现：

// Modified from https://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm#C
void bresenham(int x0, int y0, int x1, int y1) {
    int dx = abs(x1 - x0), sx = x0 < x1 ? 1 : -1;
    int dy = abs(y1 - y0), sy = y0 < y1 ? 1 : -1;
    int err = (dx > dy ? dx : -dy) / 2;
    while (setpixel(x0, y0), x0 != x1 || y0 != y1) {
        int e2 = err;
        if (e2 > -dx) { err -= dy; x0 += sx; }
        if (e2 <  dy) { err += dx; y0 += sy; }
}

为了测试不同角度，我做了一个测试用例：

int main() {
    memset(img, 255, sizeof(img));
    float cx = w * 0.5f - 0.5f, cy = h * 0.5f - 0.5f;
    for (int j = 0; j < 5; j++) {
        float r1 = fminf(W, H) * (j + 0.5f) * 0.085f;
        float r2 = fminf(W, H) * (j + 1.5f) * 0.085f;
        float t = j * PI / 64.0f;
        for (int i = 1; i <= 64; i++, t += 2.0f * PI / 64.0f) {
            float ct = cosf(t), st = sinf(t);
            bresenham((int)(cx + r1 * ct), (int)(cy - r1 * st), (int)(cx + r2 * ct), (int)(cy - r2 * st));
    svpng(fopen("line_bresenham.png", "wb"), W, H, img, 0);
}

完整代码 line_bresenham.c

渲染结果及中间局部放大：

2. 采样方法

Bresenham 直线算法有 3 个问题：

不能控制直线宽度；
坐标为整数；
只能对像素写入一个颜色，只用单色会有严重的锯齿效果。

在图形学中，除了以逐个图元（如直线）方式渲染，我们还可以通过对每个像素进行采样（sampling），换句话说，我们可对整个图像逐像素询问：「这个像素的颜色是什么？」

用采样方式画直线时，我们可以用一个具有面积的形状去表示「直线」，例如是长方形。但在本文中，我们使用胶囊体（capsule）去表示直线。这样就能解决上面前两个问题：(1) 可用胶囊体半径设置直线的宽度；(2) 坐标可以为浮点数。不过，用最简单的采样方式，我们需要在每像素查询所有图元是否占有该像素，效率很低。

// 判断一点 (px, py) 是否在胶囊体（两端为(ax, ay)、(bx, by)，半径 r）之中
int capsule(float px, float py, float ax, float ay, float bx, float by, float r) {
    float pax = px - ax, pay = py - ay, bax = bx - ax, bay = by - ay;
    float h = fmaxf(fminf((pax * bax + pay * bay) / (bax * bax + bay * bay), 1.0f), 0.0f);
    float dx = pax - bax * h, dy = pay - bay * h;
    return dx * dx + dy * dy < r * r;
// 对坐标 (x, y) 进行采样
float sample(float x, float y) {
    float s = 0.0f, cx = W * 0.5f, cy = H * 0.5f;
    for (int j = 0; j < 5; j++) {
        float r1 = fmaxf(w, h) * (j + 0.5f) * 0.085f;
        float r2 = fmaxf(w, h) * (j + 1.5f) * 0.085f;
        float t = j * PI / 64.0f, r = (j + 1) * 0.5f;
        for (int i = 1; i <= 64; i++, t += 2.0f * PI / 64.0f) {
            float ct = cosf(t), st = sinf(t);
            s = fmaxf(s, capsule(x, y, cx + r1 * ct, cy - r1 * st, cx + r2 * ct, cy - r2 * st, r) ? 1.0f : 0.0f);
    return s;
int main() {
    unsigned char *p = img;
    for (int y = 0; y < H; y++)
        for (int x = 0; x < W; x++, p += 3)
            p[0] = p[1] = p[2] = (unsigned char)((1.0f - sample(x, y)) * 255);
    svpng(fopen("line_sampling.png", "wb"), W, H, img, 0);
}

完整代码 line_sampling.c

渲染结果及中间局部放大：

注意从中间往外，直线的宽度为 1, 2, 3, 4, 5。我们甚至可以渲染宽度少于一个像素的直线。

3. 超采样

所谓锯齿，其实应称作信号混叠（aliasing）。抗混叠／抗锯齿（anti-aliasing）的最简单方法是做更多采样，即超采样（supersampling）。我们可以沿用上一个例子中的采样函数 sample()，只是在每个像素内平均多个像素采样，例如以正方栅格方式做 5\times5 = 25 次采样，但计算时间也变成 25 倍。

int main() {
    unsigned char *p = img;
    int sw = 5, sh = 5;
    for (int y = 0; y < H; y++)
        for (int x = 0; x < W; x++, p += 3) {
            float sum = 0;
            for (int j = -sh / 2; j <= sh / 2; j++)
                for (int i = -sw / 2; i <= sw / 2; i++)
                    sum += sample(x + (float)i / sw, y + (float)j / sh);
            p[0] = p[1] = p[2] = (unsigned char)((1.0f - sum / (sw * sh)) * 255);
    svpng(fopen("line_supersampling.png", "wb"), W, H, img, 0);
}

完整代码 line_supersampling.c

渲染结果及中间局部放大：

这个结果可算是接近最好的了，但它的缺点就是性能极低。我们尝试另一个方法。

4. 带符号距离场

在上面的采样中，我们只用布尔值来表示一个坐标是否在形状之中。我们可以用带符号距离场（signed distance field, SDF）去表示坐标与形状的带符号距离（如正数代表坐标在形状以外，距离形状的最短距离）。[2] 采用这种方式提升矢量图形的渲染品质。

我们只需简单改变 capsule() 函数，就能把它变成一个 SDF：

int capsule(float px, float py, float ax, float ay, float bx, float by, float r) {
    float pax = px - ax, pay = py - ay, bax = bx - ax, bay = by - ay;
    float h = fmaxf(fminf((pax * bax + pay * bay) / (bax * bax + bay * bay), 1.0f), 0.0f);
    float dx = pax - bax * h, dy = pay - bay * h;
    return dx * dx + dy * dy < r * r;
float capsuleSDF(float px, float py, float ax, float ay, float bx, float by, float r) {
    float pax = px - ax, pay = py - ay, bax = bx - ax, bay = by - ay;
    float h = fmaxf(fminf((pax * bax + pay * bay) / (bax * bax + bay * bay), 1.0f), 0.0f);
    float dx = pax - bax * h, dy = pay - bay * h;
    return sqrtf(dx * dx + dy * dy) - r; // 只改变这句及返回类型
}

如果把例子中的 SDF 可视化，就会得到这种效果：

可以想像到，当一个像素的中心位于形状的边界时，SDF 应为 0，可用不透明度（opacity） 0.5 来渲染；当一个像素刚好离开形状时，SDF 为 0.5，不透明度为 0；当一个像素刚好进入形状时，SDF 为 -0.5，不透明度为 1。

因此，我们对每像素只作一次 SDF 采样，然后用 SDF 来决定不透明度：

float sample(float x, float y) {
    float s = 0.0f, cx = w * 0.5f, cy = h * 0.5f;
    for (int j = 0; j < 5; j++) {
        float r1 = fminf(W, H) * (j + 0.5f) * 0.085f;
        float r2 = fminf(W, H) * (j + 1.5f) * 0.085f;
        float t = j * PI / 64.0f, r = (j + 1) * 0.5f;
        for (int i = 1; i <= 64; i++, t += 2.0f * PI / 64.0f) {
            float ct = cosf(t), st = sinf(t);
            s = fmaxf(s, fminf(0.5f - capsuleSDF(x, y, cx + r1 * ct, cy - r1 * st, cx + r2 * ct, cy - r2 * st, r), 1.0f));
    return s;
int main() {
    unsigned char *p = img;
    for (int y = 0; y < H; y++)
        for (int x = 0; x < W; x++, p += 3)
            p[0] = p[1] = p[2] = (unsigned char)((1.0f - sample(x, y)) * 255);
    svpng(fopen("line_sdf.png", "wb"), W, H, img, 0);
}

完整代码 https:// github.com/miloyip/line /blob/master/line_sdf.c

渲染结果及中间局部放大：

这个效果几乎和超采样没有大分别，但速度快了很多。

5. 用 AABB 优化 SDF

上述的采样方法，都需要在每个像素遍历所有图元，我们可以改为另一种方式，就是对每一图元，按它的范围内采样，并把结果以 alpha 混合（alpha blending）方式写入缓冲。最简单的范围是图形的轴对齐矩形（axis-aligned bounding box, AABB）。

void alphablend(int x, int y, float alpha, float r, float g, float b) {
    unsigned char* p = img + (y * W + x) * 3;
    p[0] = (unsigned char)(p[0] * (1 - alpha) + r * alpha * 255);
    p[1] = (unsigned char)(p[1] * (1 - alpha) + g * alpha * 255);
    p[2] = (unsigned char)(p[2] * (1 - alpha) + b * alpha * 255);
void lineSDFAABB(float ax, float ay, float bx, float by, float r) {
    int x0 = (int)floorf(fminf(ax, bx) - r);
    int x1 = (int) ceilf(fmaxf(ax, bx) + r);
    int y0 = (int)floorf(fminf(ay, by) - r);
    int y1 = (int) ceilf(fmaxf(ay, by) + r);
    for (int y = y0; y <= y1; y++)
        for (int x = x0; x <= x1; x++)
            alphablend(x, y, fmaxf(fminf(0.5f - capsuleSDF(x, y, ax, ay, bx, by, r), 1.0f), 0.0f), 0.0f, 0.0f, 0.0f);
int main() {
    memset(img, 255, sizeof(img));
    float cx = w * 0.5f, cy = h * 0.5f;
    for (int j = 0; j < 5; j++) {
        float r1 = fminf(W, H) * (j + 0.5f) * 0.085f;
        float r2 = fminf(W, H) * (j + 1.5f) * 0.085f;
        float t = j * PI / 64.0f, r = (j + 1) * 0.5f;
        for (int i = 1; i <= 64; i++, t += 2.0f * PI / 64.0f) {
            float ct = cosf(t), st = sinf(t);
            lineSDFAABB(cx + r1 * ct, cy - r1 * st, cx + r2 * ct, cy - r2 * st, r);
    svpng(fopen("line_sdfaabb.png", "wb"), W, H, img, 0);
}

完整代码 line_sdfaabb.c

渲染结果及中间局部放大：

这种方式结合了光栅化的程序结构（逐一渲染图元），并在渲染每个图元中采用 SDF 采样来达到抗锯齿的效果。

视觉效果和性能对比

后三者的视觉效果几乎没有差别，但性能上就有很多差异：

$ make test
time ./line_bresenham
        0.03 real         0.02 user         0.00 sys
time ./line_sampling
        1.93 real         1.91 user         0.00 sys