Milo的编程
首发于Milo的编程
用 C 语言画直线

用 C 语言画直线

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

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

1. Bresenham 算法

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

Bresenham 直线

如果像上图,直线的高度小于宽度,那么 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 个问题:

  1. 不能控制直线宽度;
  2. 坐标为整数;
  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);
}

完整代码 github.com/miloyip/line

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

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

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
time ./line_supersampling
       47.06 real        46.85 user         0.10 sys
time ./line_sdf
        2.00 real         1.98 user         0.00 sys
time ./line_sdfaabb
        0.03 real         0.03 user         0.00 sys

分析及总结

本文介绍了 5 个直线渲染方法。但实际上,除了 Bresenham 算法,其他的方法都可用于绘画任何形状,只需要提供形状的采样函数(判断一点是否在形状之中)或 SDF 便可。这样方法也可以实际应用在 GPU 硬件上,渲染抗锯齿的矢量形状。

本文的例子只是用于教学性质,它并未有加入直线剪裁(line clipping)算法,若绘画的直线跨越画布会做成崩溃。

另外,最后的方法以 AABB 作为直线图元的范围,那么最坏情况(画布对角直线)的性能还是很差。事实上你可以把范围缩小成定向包围盒(oriented bounding box, OBB),用 Bresenham 画线方法来绘制 OBB 的外框,从而缩小了采样的范围。这个优化可留待读者自行实现。

P. S. 本文题图使用 SDF AABB 方法画线,生成自 stitchheart.c

参考

[1] Bresenham, Jack E. "Algorithm for computer control of a digital plotter." IBM Systems journal 4.1 (1965): 25-30.

[2] Green, Chris. "Improved alpha-tested magnification for vector textures and special effects." ACM SIGGRAPH 2007 courses. ACM, 2007.

编辑于 2017-10-31

文章被以下专栏收录