从 Bresenham 画直线到画圆,其实是图形学里绕不开的一道坎。很多人第一次接触时觉得这就是“用点拼成一个圈”,但真的上手会发现,直接用 (y = \sqrt{r^2 - x^2}) 逐点算,不仅慢,而且画出来的圆坑坑洼洼,根本不能用。Bresenham 画圆算法是三条线段之后最值得吃透的扫描转换算法,也是很多嵌入式图形库和单片机显示驱动里的标准实现方式。这篇文章会从算法原理、数学推导一路讲到 MATLAB 的完整实现和调试经验,适合正在学计算机图形学、准备做课程设计,或者想在硬件上画图的开发者参考。算法本身不复杂,关键是把决策参数的更新逻辑真正理解透,而不只是背代码。
1. 为什么画圆不能“老老实实”按函数算
1.1 理想圆与像素网格之间的天然矛盾
计算机屏幕上的图像是由离散像素组成的,而圆在数学上是连续曲线。所谓画圆,本质上是找出“最接近理想圆弧”的一组像素点,这个过程叫扫描转换(Scan Conversion)。如果直接把理想圆弧覆盖到像素网格上,圆弧会穿过一系列像素的中缝,你必须决定每个位置到底点亮哪一个像素。
最朴素的想法是遍历 x 从 0 到 r,代入圆的方程算出 y 值,再取整。这个思路听起来没问题,但实际用起来有两个硬伤。第一,sqrt计算涉及浮点开方,耗时远高于整数加减运算,在 CPU 弱、没有硬件浮点单元的嵌入式设备上根本无法接受。第二,取整后像素点之间会出现明显的断点,因为圆弧在 x 轴附近变化快、在 y 轴附近变化慢,均匀步进 x 画出来的点疏密不均,还会在斜率大的区域留下空洞。
1.2 圆的八分对称性:工作量直接砍掉 87.5%
圆有一个非常好的性质:它关于坐标轴和两条对角线对称。换句话说,在圆周上随便取一个点 ((x, y)),通过对称映射可以得到另外 7 个点。这八个点合在一起,就覆盖了完整圆周。
这意味着,你只需要精确计算出第一象限 (0^\circ) 到 (45^\circ)(即从 ((0, r)) 到 ((r/\sqrt{2}, r/\sqrt{2})))这八分之一段圆弧的像素,其余七段直接通过坐标交换和取负号映射出来就行。这一下就把计算量缩小到了原来的八分之一,也让算法的主循环变得非常轻量。Bresenham 画圆算法正是基于这一点,先算八分之一圆弧,再由对称性生成完整圆。
1.3 画直线算法和画圆算法的本质区别
学过 Bresenham 画直线的人会发现,直线的决策参数很好理解:比较当前点之后“该往上还是往下”,误差项是一个关于斜率的累积偏差。画圆的时候很多人想照搬这个思路,但会遇到一个麻烦——圆弧的斜率是变化的。
在 (45^\circ) 到 (90^\circ) 这段圆弧上,x 从 0 增加到 r,每走一个像素,y 方向要不要减一,完全取决于当前位置的切线斜率。这里不能用固定斜率来决策,而是需要引入一个动态更新的误差度量。这就是 Bresenham 画圆算法比画直线多出来的核心设计:用一个关于圆心距离的符号来判断候选像素点,并让决策参数随着每一步动态递推。
2. 决策参数是怎么推出来的:核心数学原理
2.1 从圆的隐式方程出发
圆心在原点、半径为 r 的圆可以写成隐式形式:
[ f(x, y) = x^2 + y^2 - r^2 ]
这个函数的符号非常有用。对任意一点 ((x, y)):
- (f(x, y) < 0) 表示该点在圆内
- (f(x, y) = 0) 表示该点正好在圆上
- (f(x, y) > 0) 表示该点在圆外
画圆的时候,我们沿着 x 方向从 0 逐步递增,每次 x 加 1,y 的候选值只有两种可能:保持当前 y 不变,或者 y 减 1。因为在这八分之一圆弧上,x 始终增加,y 始终减少或不变化,不会出现 y 增加的情况。像素网格上的每一个格子,本质就是在判断“下一个点到底是右边那个,还是右下那个”。
2.2 中点决策:取距离圆弧更近的像素
假设当前像素点是 ((x_k, y_k)),那么下一个候选点是 E(即 ((x_k+1, y_k)))或 SE(即 ((x_k+1, y_k-1)))。为了决定选哪个,取这两个候选点的中点 (M = (x_k+1, y_k - 0.5)),计算它的隐式函数值:
[ d_k = f(x_k+1, y_k - 0.5) = (x_k+1)^2 + (y_k - 0.5)^2 - r^2 ]
如果 (d_k < 0),说明中点在圆内,圆弧从中点上方穿过,此时 E 点离圆弧更近,选 E;反之,(d_k \ge 0) 时中点落在圆外或圆上,圆弧从中点下方穿过,选 SE。这个思路就是教科书上常说的“中点画圆法”,Bresenham 画圆算法在本质上也是这个决策逻辑,只是进一步做整数化处理,避免了浮点运算。
2.3 决策参数的整数化递推:让算法只剩加减法
直接按上面的公式每步都重新计算 (d_k),虽然能工作,但每步要做平方运算,还不够“漂亮”。Bresenham 的精髓在于:根据上一步的决策结果,用一个简单的整数加减操作来更新决策参数,全程没有任何乘法和浮点数。
设当前决策参数为 (d_k),计算下一步:
- 如果选了 E(即 (x) 加 1,(y) 不变):
[ d_{k+1} = (x+2)^2 + (y - 0.5)^2 - r^2 = d_k + 2x + 3 ]
- 如果选了 SE(即 (x) 加 1,(y) 减 1):
[ d_{k+1} = (x+2)^2 + (y - 1.5)^2 - r^2 = d_k + 2x - 2y + 5 ]
初始条件:起点是 ((0, r)),因此:
[ d_0 = (0+1)^2 + (r - 0.5)^2 - r^2 = 1.25 - r ]
注意这里出现了 0.25 的小数部分。为了彻底消除浮点,可以把整个决策参数乘以 4,得到整数版的初始值和递推式。这就是市面上各种教材里那个“ (d = 3 - 2r) ”版本公式的来源。我在实际推导中发现,不同教材对 (d) 的初始化和更新式有不同的写法,有一类用 (d_0 = 1 - r)、步进分别为 (2x+3) 和 (2x-2y+5) 的变体,其实是通过不同方式消除浮点后得到的等价形式。用的时候不要混着抄,认准一个体系推导完再记忆,否则很容易在符号上翻车。
2.4 主循环的完整逻辑
整理一下,整数版 Bresenham 画圆算法主循环的流程是:
- 初始化 (x = 0),(y = r),决策参数 (d = 3 - 2r)。
- 调用对称映射,生成当前 ((x, y)) 对应的 8 个圆周像素。
- 如果 (d < 0),选择 E 点,更新 (d = d + 4x + 6),(y) 保持不变。
- 否则选择 SE 点,更新 (d = d + 4(x - y) + 10),同时 (y = y - 1)。
- (x = x + 1),检查 (x \le y) 是否成立,成立则回到步骤 2。
循环终止条件是 (x > y),因为当 (x) 超过 (y) 时,意味着我们已经走完了从 ((0, r)) 到 (45^\circ) 方向的所有像素,接下来就是对称映射的管辖范围了。
3. MATLAB 实现全流程:从函数封装到出图
3.1 函数接口设计
写 MATLAB 实现时,我建议把算法封装成一个独立函数,而不是在脚本里堆一段代码。这样既方便课程作业中使用,也能直接复用到之后做图像处理、形状识别之类的小项目里。
函数的输入参数应该包括圆心坐标(xc, yc)和半径r,返回的是圆上所有像素点的行列坐标数组。我更喜欢返回坐标数组而不是直接画图,因为这样调用方可以自己决定用plot散点图、scatter还是转成二值图像矩阵来进一步处理。
function [X, Y] = bresenhamCircle(xc, yc, r) % bresenhamCircle 使用 Bresenham 算法生成圆上像素坐标 % 输入: % xc, yc - 圆心坐标 % r - 半径(正整数) % 输出: % X, Y - 圆上像素点的横纵坐标列向量3.2 核心算法的整数实现
算法主体的实现代码如下。我用的是上面推导出的整数版决策参数体系,全程只有加减法和整数比较,这套代码可以直接移植到 C 语言或者单片机上,几乎没有改动成本。
x = 0; y = r; d = 3 - 2 * r; X = []; Y = []; while x <= y % 八对称映射:从当前八分之一圆弧点生成完整圆周的 8 个点 pts_x = [xc + x, xc - x, xc + x, xc - x, ... xc + y, xc - y, xc + y, xc - y]; pts_y = [yc + y, yc + y, yc - y, yc - y, ... yc + x, yc + x, yc - x, yc - x]; X = [X; pts_x(:)]; Y = [Y; pts_y(:)]; % 根据决策参数选择下一个点并更新参数 if d < 0 % 选择 E 点:(x+1, y) d = d + 4 * x + 6; else % 选择 SE 点:(x+1, y-1) d = d + 4 * (x - y) + 10; y = y - 1; end x = x + 1; end end有几点值得说明。第一,八对称映射的 8 个点顺序其实无所谓,因为最终都是打点,不涉及连线方向。第二,x和y的更新是有先后顺序的:先根据旧的x、y更新决策参数,再执行x = x + 1和可能的y = y - 1。如果把顺序搞反,算出来的d就对不上当前位置了,最后画出来的圆会有缺损。第三,在开始编写之前,务必确认你用的决策参数版本和初始化公式是配套的。我之前就见过有同学从两份讲义上各抄一半,导致画出来的圆在最上方出现一个明显的“台阶”。
3.3 调用与可视化:验证算法是否真的画出了圆
函数封装好之后,调用方式非常简单。下面这段代码生成一个圆心在(100, 100)、半径 50 的圆,并用散点图绘制出来:
% 生成圆上像素坐标 xc = 100; yc = 100; r = 50; [X, Y] = bresenhamCircle(xc, yc, r); % 绘制散点图 figure; plot(X, Y, 'r.', 'MarkerSize', 8); axis equal; % 关键:强制坐标轴等比例,否则圆会变形 grid on; xlim([xc - r - 10, xc + r + 10]); ylim([yc - r - 10, yc + r + 10]); title('Bresenham Circle');这里axis equal这行一定要加。MATLAB 默认会分别缩放 x 轴和 y 轴,如果不加这一句,画出来的圆会严重变形,看起来像椭圆,很多人还以为算法写错了,其实是坐标轴比例的问题。运行之后你应该能看到一个由离散点组成的、轮廓光滑的圆形,圆周上相邻点之间的间距不超过一个像素,没有明显的断口。
3.4 把算法扩展到任意位置和大小的圆
上面的函数已经支持任意圆心坐标,因为对称映射中的xc、yc会让所有坐标进行整体平移。半径r也支持较大的整数,只要不超出 MATLAB 的整数范围即可。
但有一点需要提醒:当r特别小(比如 1 或 2)时,这个算法画出的圆看起来会非常粗糙,甚至看不出是圆。这是所有像素化画圆算法的通病——在半径小于 3 个像素时,可用的像素点太少,无论怎么选点都很难形成视觉上的圆形。这种情况不要指望算法改进,而应该直接提高分辨率,或者用更大半径再缩放。
如果想把画出的圆输出为二值图像矩阵(比如为了后续做形态学处理),可以利用返回的坐标直接索引赋值:
imgSize = 201; img = zeros(imgSize, imgSize); idx = sub2ind([imgSize, imgSize], X, Y); img(idx) = 1; imshow(img, []);以圆心(100, 100)为例,只要确保圆心坐标和半径映射在图像尺寸内,就能得到一张 201×201 的圆环二值图。这是一个非常常用的小技巧,很多同学一开始不知道,都是直接plot出图,再截图转成图片,精度就丢了。
4. 多角度验证:正确性、效率与效果对比
4.1 用像素误差来验证算法的正确性
仅仅“看起来像个圆”是不够的。严谨一点,应该计算每个生成像素到理想圆心的距离,和半径r比较,看看误差有多大。Bresenham 算法保证每个像素到理想圆的误差不超过半个像素,下面这段代码可以帮助你量化这个误差:
% 计算每个像素点到圆心的实际距离 distances = sqrt((X - xc).^2 + (Y - yc).^2); errors = abs(distances - r); fprintf('最大像素误差: %.3f\n', max(errors)); fprintf('平均像素误差: %.3f\n', mean(errors));以半径 50、圆心 (100, 100) 为例,运行结果会显示最大误差在 0.5 左右,均值远小于 0.5。这个“误差不超过 0.5 像素”的结论正是算法设计的核心目标,如果实测结果大幅超过 0.5,那基本可以确定决策参数或者初始化部分写错了。
4.2 与朴素平方根方法的对比
为了直观体会 Bresenham 算法的价值,可以和“直接计算”的版本做个对比。朴素方法就是 x 从 0 到 r 遍历,每个点直接算 (y = \sqrt{r^2 - x^2}) 然后取整:
x0 = 0:r; y0 = round(sqrt(r^2 - x0.^2)); figure; plot(x0, y0, 'b.', 'MarkerSize', 8); title('Naive sqrt approach');跑一下就能看到两个明显的问题。第一,在 x 接近 r 的地方,y 的变化非常快,相邻采样点之间会直接跳跃好几个像素,画出的右半部分有明显的撕裂感。第二,自己算一遍会发现 MATLAB 里sqrt运算在循环中速度远慢于整数加减,虽然画一个圆差距不大,但如果画成千上万个圆(比如字符显示、地图渲染),性能差距会非常可观。
4.3 多个半径与不同视角的测试
我建议在验证阶段多跑几组参数,而不是只画一个圆。可以测试半径从 5 到 200 的多个圆叠加在一起,检查任意半径下有没有断点或者不对称的像素。一个常见测试是同时画同心圆,观察它们之间的间距是否均匀、每层是否有缺口:
figure; hold on; axis equal; for r = 10:10:100 [X, Y] = bresenhamCircle(150, 150, r); plot(X, Y, '.', 'MarkerSize', 6); end如果每一圈都完整且没有重叠错位,说明算法的对称映射和循环终止条件都是正确的。注意看从 (45^\circ) 附近过渡到其他象限时,像素是否平滑衔接。这里最容易出的问题就是循环条件误写成x < y导致多算一段、或者漏算对角点,最终在圆的斜对角线方向出现缺口。
4.4 算法扩展:椭圆与圆弧的切入点
掌握了 Bresenham 画圆之后,向椭圆和任意圆弧扩展会顺理成章。画椭圆的思路是把圆的决策参数变成两个方向的不等比例尺度,但推导会更繁琐,因为不再有八对称性,只有四对称性,而且切点位置会发生变化。至于画圆弧,则是在生成完整圆之后,根据起始角、终止角筛选像素点,或者直接修改主循环的遍历区间。
这一步很多人是到做绘图软件、CAD 插件或者游戏引擎里的迷你小地图时才真正需要。如果你是应付图形学课程作业,做到画圆就基本够用了,但建议把扩展思路记住,面试或者实际项目里经常会绕回这些基础算法。
5. 常见问题与排查技巧实录
5.1 圆看起来像椭圆
绝大部分情况下是显示问题,而不是算法问题。MATLABplot默认的坐标轴比例不一致,会让同样距离在屏幕上的像素长度不同,导致圆形被水平和垂直拉伸。解决方法是使用axis equal或者手动设置daspect([1 1 1])。如果已经设置了axis equal还是椭圆,检查一下圆心xc、yc和半径r是否都是整数,以及数据是否有 NaN 或者意外越界。
5.2 圆周在斜向出现缺口
这个基本上是八对称映射写错了。有些初学者会把 8 个坐标中的某个符号弄反,比如本该xc - y写成了xc + y,导致在对角线方向缺了一块或重复了一块。建议对照着圆周自己手推一遍:取 (45^\circ) 方向附近的点,比如当前(x, y) = (r/sqrt(2), r/sqrt(2))\),把 8 个映射点画在纸上,看它们是不是均匀围成一个圈。还有个隐蔽问题:当循环条件写成while x < y时,会漏掉 \(x = y\) 所对应的对角点,导致圆的 \(45^\circ\) 方向出现小缺口,应改成x <= y`。
5.3 决策参数更新后圆的上方出现“台阶”
如果圆的顶部(即 (90^\circ) 方向)出现明显的水平长线段,或者像素点排列成一个台阶形状,多半是决策参数使用了不配套的初始化和更新公式。我的排查经验是:先固定半径,用断点调试看前 5 个像素坐标,理想情况下应该是这样变化的:从 ((0, r)) 出发,x 递增,y 维持不变若干步之后再减 1,并且整个过程不会出现 y 一下子减 2 或者完全不动的情况。还可以找一个特定半径(比如 r=10)手算前几步,和代码结果逐项对照。
5.4 半径较大时 MATLAB 运行变慢
如果只是在 MATLAB 里跑,半径几千也许还撑得住,但到几万像素级别时,我上面这种X = [X; pts_x(:)]的数组动态拼接方式会带来很大的复制开销。更高效的做法是提前估算像素数量并预分配数组,或者改用循环中用索引逐步填充:
totalPoints = 8 * (floor(r / sqrt(2)) + 1); X = zeros(totalPoints, 1); Y = zeros(totalPoints, 1); idx = 1; % 在循环中逐点填入并更新 idx这个优化在小半径时感知不强,但做高性能渲染或大规模批量画圆时会拉开很大差距。毕竟 Bresenham 算法的魅力就在于极低的时间复杂度,别让 MATLAB 的数组拼接把优势浪费掉了。
5.5 在非整数网格坐标下使用
有些项目可能需要在亚像素精度下画圆,直接用 Bresenham 会很别扭。这类场景下算法本身仍然可以用,只要你把坐标先放大到整数网格再缩放回去,或者干脆改用抗锯齿画圆算法。记住 Bresenham 定位是“光栅化”,也就是把连续几何转换成离散像素,它天然工作在整数像素网格上。别试图修改它去处理浮点圆心坐标,那是舍近求远。
5.6 画出的圆整体偏移了圆心
这个问题通常不是算法造成的,而是坐标系约定不同。MATLAB 的plot默认原点在左下角,x 向右、y 向上;但如果你在读入图像并用imshow显示,图像坐标系是原点在左上角、y 向下的。同一个圆心(100, 100),在两种坐标系下的视觉效果是完全不同的,甚至可能导致圆跑到图像外。用imshow显示二值图像时,建议先确认圆心坐标和半径都落在图像矩阵有效索引范围内,或者先绘制再裁剪显示区域。
写在最后:几个实际项目中踩过的坑
我最初学这个算法时也踩了不少坑。最深刻的一点是:对称映射不仅是把 8 个点写出来那么简单,它的顺序和边界条件直接决定你是否会在斜对角方向产生一条“缝”。有一次我在移植这个算法到嵌入式 LCD 驱动时,就是因为循环条件多写了一个等号,导致圆在四个斜向全部出现微小缺口,排查了一整天才发现。用 MATLAB 验证算法时,建议把生成的点打印出来,直观对比各个象限的对称性,这是最快的排查手段。
另外一个小建议是:不要把 Bresenham 算法当成“过时的课本知识”。在很多没有 GPU、没有浮点单元的微控制器上,这套只能用整数加减法的算法到今天仍然是绘制几何图形的主力方案。真正理解了它的设计思想(用误差累积和整数化递推逼近理想几何形状),再去看反锯齿画线、椭圆生成、Bezier 曲线的离散化,都会觉得顺很多。
自己在 MATLAB 里动手把这个函数写完、跑通、画出来,比看十遍推导都管用。遇到问题也别急着搜答案,先对着推导过程手算两步,你会有完全不一样的收获。