圆的两种生成算法(角度微分法、Bresenham算法)
1.角度微分法的原理
圆的角度微分法是用圆的内接正多边形来逼近该圆。
若我们设圆的参数方程为:
{ x = R ⋅ cos φ y = R ⋅ sin φ \begin{cases} x=R·\text{cos} \varphi \\ y=R·\text{sin} \varphi \\ \end{cases} {x=R⋅cosφy=R⋅sinφ
其中 φ \varphi φ ( 0 ° ≤ φ ≤ 360 ° 0 \degree \le \varphi \le 360 \degree 0°≤φ≤360°)为旋转角, 现把该圆 n n n 等分,用 n n n 条直线来逐次逼近该圆,则各直线的端点坐标计算如下:
设旋转角
φ
\varphi
φ 的起始角、终止角分别为
α
,
β
\alpha, \beta
α,β,且满足关系式
0
°
≤
α
≤
β
≤
360
°
0 \degree \le \alpha \le \beta \le 360 \degree
0°≤α≤β≤360°。
则对圆进行
n
n
n 等分时,其对应的角度增量为:
Δ
φ
=
2
π
n
\varDelta \varphi = \frac{2\pi}{n}
Δφ=n2π。
于是,第
i
i
i 个点所对应的角度与端点坐标为(其中
i
=
0
,
1
,
2
,
…
,
n
i = 0 , 1 , 2 , \dots , n
i=0,1,2,…,n):
{ φ i = α + i Δ φ x i = R cos φ i y i = R sin φ i \begin{cases} \varphi_i=\alpha + i \varDelta \varphi \\ x_i = R\text{cos} \varphi_i \\ y_i=R\text{sin} \varphi_i \\ \end{cases} ⎩ ⎨ ⎧φi=α+iΔφxi=Rcosφiyi=Rsinφi
上式即为生成圆的基本公式之一。我们不难看出,
n
n
n 越大,精度越高,但计算也越 复杂。下一步是如何确定
R
R
R 和
n
n
n 的关系?
可以这样认为:如果圆的内接多边形的边线与理想圆弧之间的距离在一个刻度单位之内时,其近似显示效果是能接受的。如下图所示:
即:如果线段
a
b
ab
ab 的距离在一个刻度单位内都是能被接受的,即要求
a
b
<
1
ab < 1
ab<1。
那
a
b
ab
ab 的距离等于多少呢?显然
a
b
=
R
−
R
cos
Δ
φ
<
1
ab = R − R\text{cos} \varDelta \varphi < 1
ab=R−RcosΔφ<1。
对于
cos
Δ
φ
\text{cos} \varDelta \varphi
cosΔφ,可按泰勒级数展开取前两项,得到
cos
Δ
φ
=
1
−
Δ
φ
2
2
\text{cos} \varDelta \varphi=1-\frac{{\varDelta \varphi}^2}{2}
cosΔφ=1−2Δφ2。
将其带入
a
b
ab
ab 中即得到
R
−
R
(
1
−
Δ
φ
2
2
)
<
1
R-R \left(1-\frac{{\varDelta \varphi}^2}{2} \right)<1
R−R(1−2Δφ2)<1,解之可得
Δ
φ
<
2
R
\varDelta \varphi<\sqrt\frac{2}{R}
Δφ<R2。
将
Δ
φ
=
2
π
n
\varDelta \varphi = \frac{2 \pi}{n}
Δφ=n2π 带入上式可得
n
>
π
2
R
n>\pi\sqrt{2R}
n>π2R。
用经验公式取代得到
n
=
0.1
R
+
15
n = 0.1R + 15
n=0.1R+15。
这就是圆的半径与其内接正多边形边数之间的关系 。
2.角度微分法的实现(基于matlab)
下面给出利用matlab实现该算法的完整代码:
function main %测试主函数
clc
clear;
DifferentiationOfCircle(100,'r') %调用角度微分法
function DifferentiationOfCircle(R,color) %圆的角度微分算法
n = 0.1*R+15; %计算出需要微分的次数
delta = 2*pi/n; %从而得到每次需要增加的角度
alpha = delta; %初始化alpha
x0 = R,y0 = 0; %初始化x0,y0
hold on
while(n>0)
x = R*cos(alpha);
y = R*sin(alpha);
plot([x0,x],[y0,y]); %绘制点(x0,y0)和(x,y)之间的线段
alpha = alpha + delta; %更新alpha
x0 = x , y0 = y; %更新x0、y0
n = n-1; %n减1
end
grid minor
hold off
在matlab中运行上述代码,实验结果如下:
可以看出,角度微分法在绘圆时具有明显的分段现象。当然,我们可以通过加大 n n n 的值来适当减轻分段效果。不过这样的代价是会使得程序需花费更多的时间来进行微分。 所以这才有了 n = 0.1 R + 15 n = 0.1R + 15 n=0.1R+15 这一公式的出现。
3.Bresenham 算法的原理
圆可被定义为到给定中心位置 O O O(圆心)距离为 r r r 的点集。圆心位于原点的圆有四条对称轴 x = 0 , y = 0 , x = y x=0,y=0, x=y x=0,y=0,x=y 和 x = − y x=-y x=−y。若已知圆弧上一点 ( x , y ) (x,y) (x,y),可以得到其关于四条对称轴的其它 7 个点,这种性质称为八分对称性(如下图所示)。因此,只要扫描转换八分之一圆弧,就可以求出整个圆弧的象素集。
所以接下来我们的重心就放在:如何将某段圆弧尽可能描绘清楚?
对于上图中,点 ( y , x ) (y,x) (y,x) 所在那一段圆弧而言,其特征是 x x x 坐标变化率大于 y y y 坐标变化 率,且随着 x x x 的增加, y y y 在减小。因此逼近该圆弧的绘点步进规律是:每推算一次,其 x x x 坐标都加 1,而 y y y 坐标根据相应误差决定是否减 1。我们易知该圆中最上方的那个点的坐标为 ( 0 , R ) (0,R) (0,R),若将其设为 ( x 0 , y 0 ) (x_0,y_0) (x0,y0),则接下来的一点的坐标只有以下两种可能:
- 情况一: ( x 0 + 1 , y 0 ) (x_0+1,y_0) (x0+1,y0)
- 情况二: ( x 0 + 1 , y 0 − 1 ) (x_0+1,y_0-1) (x0+1,y0−1)
现在我们的任务就是确定到底取谁。
为便于分析,设圆的方程为: d ( x , y ) = x 2 + y 2 − R 2 d(x,y)=x^2+y^2-R^2 d(x,y)=x2+y2−R2,则对给定的任意点 ( x 0 , y 0 ) (x_0,y_0) (x0,y0), d ( x 0 , y 0 ) d(x_0,y_0) d(x0,y0) 的正负就能反应该点位于圆内还是圆外:
- 若 d ( x 0 , y 0 ) < 0 d(x_0,y_0)<0 d(x0,y0)<0,则表示该点位于圆内;
- 若 d ( x 0 , y 0 ) = 0 d(x_0,y_0)=0 d(x0,y0)=0,则表示该点位于圆上;
- 若 d ( x 0 , y 0 ) > 0 d(x_0,y_0)>0 d(x0,y0)>0,则表示该点位于圆外;
在确定下一个点的位置坐标时,我们实际上并不关心它到底在圆内还是圆外,而是关心谁距离圆的那一段圆弧更近。因此,我们可以将这两个点带入圆的定义式 中,即得到 d ( x 0 + 1 , y 0 ) , d ( x 0 + 1 , y 0 − 1 ) d(x_0+1,y_0), d(x_0+1,y_0-1) d(x0+1,y0),d(x0+1,y0−1),通过加绝对值的方式来获得这两点距离圆弧的绝对距离,最后取其中的较小值作为下一个点的位置即可。
因此,如果我们设点 U i ( x i − 1 + 1 , y i − 1 ) , D i ( x i − 1 + 1 , y i − 1 − 1 ) U_i(x_{i-1}+1,y_{i-1}), D_i(x_{i-1}+1,y_{i-1}-1) Ui(xi−1+1,yi−1),Di(xi−1+1,yi−1−1),则有:
{ d ( U i ) = ( x i − 1 + 1 ) 2 + y i − 1 2 − R 2 d ( D i ) = ( x i − 1 + 1 ) 2 + ( y i − 1 − 1 ) 2 − R 2 \begin{cases} d(U_i)=(x_{i-1}+1)^2+y_{i-1}^2-R^2 \\ d(D_i)=(x_{i-1}+1)^2+(y_{i-1}-1)^2-R^2 \\ \end{cases} {d(Ui)=(xi−1+1)2+yi−12−R2d(Di)=(xi−1+1)2+(yi−1−1)2−R2
则可以将圆的 Bresenham 算法的偏差判别式定义为:
d i = ∣ d ( U i ) ∣ − ∣ d ( D i ) ∣ d_i=\mid d(U_i)\mid - \mid d(D_i) \mid di=∣d(Ui)∣−∣d(Di)∣
其原则是:
- 当 d i < 0 d_i<0 di<0 时,选 U i U_i Ui作为下一个绘制的像素点
- 当 d i ≥ 0 d_i \ge 0 di≥0 时,选 D i D_i Di作为下一个绘制的像素点
由于存在约束: d ( U i ) ≥ 0 , d ( D i ) ≤ 0 d(U_i) \ge 0, d(D_i) \le 0 d(Ui)≥0,d(Di)≤0 ,故可将上面偏差判别式的绝对值去掉,得到:
d i = d ( U i ) + d ( D i ) d_i=d(U_i) + d(D_i) di=d(Ui)+d(Di)
根据这个判别式,我们可通过一个循环来将八分之一圆弧绘制出。这之后,就能同时通过坐标变换将整个圆画出(圆的八分对称性)。但是,这个式子里的运算偏多,我们需要对其进行优化,以使整个绘制过程尽量以迭代运算为主,从而减少复杂的乘幂运算次数。
首先是将 U i , D i U_i, D_i Ui,Di 的坐标带入判别式得到:
d i = 2 ( x i − 1 + 1 ) 2 + y i − 1 2 + ( y i − 1 − 1 ) 2 − 2 R 2 d_i = 2(x_{i-1}+1)^2+y_{i-1}^2+(y_{i-1}-1)^2-2R^2 di=2(xi−1+1)2+yi−12+(yi−1−1)2−2R2
然后增加下标值可得到:
d i + 1 = 2 ( x i + 1 ) 2 + y i 2 + ( y i − 1 ) 2 − 2 R 2 d_{i+1} = 2(x_i+1)^2+y_i^2+(y_i-1)^2-2R^2 di+1=2(xi+1)2+yi2+(yi−1)2−2R2
将其作差,得到:
Δ d = 2 ( ( x i + 1 ) 2 − ( x i − 1 + 1 ) 2 ) + ( y i 2 − y i − 1 2 ) + ( ( y i − 1 ) 2 − ( y i − 1 − 1 ) 2 ) \varDelta d = 2\left((x_i+1)^2-(x_{i-1}+1)^2\right)+\left( y_i^2-y_{i-1}^2 \right)+\left( (y_i-1)^2-(y_{i-1}-1)^2 \right) Δd=2((xi+1)2−(xi−1+1)2)+(yi2−yi−12)+((yi−1)2−(yi−1−1)2)
在当前选择不同的点时,该差值的结果分别如下:
- 若当前选择的下一个点是 U i U_i Ui,此时 x i = x i − 1 + 1 , y i = y i − 1 x_i = x_{i-1}+1, y_i=y_{i-1} xi=xi−1+1,yi=yi−1,带入上式得到 Δ d = 4 x i − 1 + 6 \varDelta d=4x_{i-1}+6 Δd=4xi−1+6
- 若当前选择的下一个点是 D i D_i Di,此时 x i = x i − 1 + 1 , y i = y i − 1 − 1 x_i = x_{i-1}+1, y_i=y_{i-1}-1 xi=xi−1+1,yi=yi−1−1,带入上式得到 Δ d = 4 ( x i − 1 − y i − 1 ) + 10 \varDelta d=4(x_{i-1}-y_{i-1})+10 Δd=4(xi−1−yi−1)+10
通过这样的方式,可将原来对 d i d_i di 的算术求值转换为迭代求值(即通过当前的选择算出 Δ d \varDelta d Δd 进而直接得到下一步的选择),从而大大减少了程序的运算时间。
4.Bresenham 算法的实现(基于matlab)
下面给出利用matlab实现该算法的完整代码:
function main %测试函数
clc
clear;
Bresenham(100) %调用Bresenham算法
function Plot(x,y) %此函数会将某个点集在对应的8个区域中所构成的弧线绘制出,从而构成一个圆
plot(x,y,'r') %用红色的线绘制第一象限的圆弧
plot(y,x,'r')
plot(-x,y,'b') %用蓝色的线绘制第二象限的圆弧
plot(-y,x,'b')
plot(-x,-y','y') %用黄色的线绘制第三象限的圆弧
plot(-y,-x,'y')
plot(x,-y,'g') %用绿色的线绘制第四象限的圆弧
plot(y,-x','g')
function Bresenham(R) %Bresenham算法
d = 0 , n = 1; %初始情况下d=0
x(1) = 0 , y(1) = R;
while(x < y)
if(d < 0)
d = d + 4*x(n) + 6;
y(n+1) = y(n);
else
d = d + 4*(x(n)-y(n)) + 10;
y(n+1) = y(n) - 1;
end
x(n+1) = x(n) + 1;
n = n + 1;
end
hold on;
Plot(x,y);
grid minor;
hold off;
在matlab中运行上述代码,实验结果如下:
可以看出,Bresenham 算法在绘圆时不会出现明显的分段现象,但是却存在“抖动”弧。出现这一现象的原因是由于 plot 函数是对像素点进行连结绘制,这个函数并不具有平滑拟合的效果。但如果仅对此算法在描绘圆弧的点上进行评估的话,该算法的效果明显是优于微分算法的。