光栅图形学(二):中点圆和中点椭圆

用于生成圆和椭圆的中点圆(椭圆)算法实际上是是Bresenham算法的另一种实现,用代入中点的方法代替计算上下距离,更适合于圆和椭圆的方程形式。

朴素方法

圆心在原点的圆表示为

x2+y2=r2x^2+y^2=r^2

如果从xx方向逐列扫描,那么按照

y=[r2x2]y=[\sqrt{r^2-x^2}]

取整所得的近似圆,并不是屏幕上的一个封闭区域。

这种办法的不合理之处在于:每次移动只在xx方向上是等长的,弧长的变化却不相等。因此,在第一象限x>yx>y的区域内,每次xx增加1,弧长的增长却大于1,导致圆弧连不上。

直角坐标系下弧长在坐标方向上是不均匀的,极坐标系就不存在这个问题。所以换成圆的参数方程

{x=rcosθy=rsinθ\begin{cases} x=r\cos\theta \\ y=r\sin\theta \end{cases}

再对θ\theta取足够小的步长进行迭代,就能画出圆。

类似的,如果用椭圆的直角坐标方程来画,也有不连续的问题;而用椭圆的参数方程,也可以画出完整的椭圆(只要极角的步长足够小)。

这些办法都有类似暴力画直线的问题,就是计算量太大。尤其是参数方程要计算三角函数值,可能使用多项式逼近(级数展开、插值)或者Cordic算法[1]实现,算一次的计算量就很可观。



应用对称性简化计算

圆和椭圆都有一定的对称性,如果得到了一个象限的生成结果,不用再计算就可以对称到其余三个象限。

不仅如此,圆的对称性更好,只要求出半个象限就足够了,对y=xy=x作对称就可以得到另外半个象限。因此生成圆的算法也叫八分圆算法。



实现八分之一圆的 Bresenham 算法

如果能应用 Bresenham 算法的思想:在一个方向上+1+1,另一个方向上判断是不变还是±1\pm 1,就能大幅度减少计算量。

(x0,y0)=(0,r)(x_0,y_0)=(0,r),在第一象限y>xy>x的部分生成一个圆,这需要判断下一个点是(xn+1,yn)(x_n+1,y_n),还是(xn+1,yn1)(x_n+1,y_n-1)

如果还按照上一篇中直线的 Bresenham 算法,那应该计算sgn(dlowerdupper)\mathrm{sgn}(d_{\text{lower}}-d_{\text{upper}})。但是这两个距离中都含有根式,用平方之类的办法也不好消除。

那么换一种做法:取一个中点(xn+1,yn+12)(x_n+1,y_n+\frac{1}{2}),讨论这个点的位置:

  • 如果中点在圆内(靠下),那么(xn+1,yn)(x_n+1,y_n)离圆弧更近。
  • 如果中点在圆外(靠上),那么(xn+1,yn1)(x_n+1,y_n-1)离圆弧更近。
  • 如果中点在圆上,随意按一种情况处理都可以。

根据点与圆的位置关系,可以代入中点到函数

F(x,y)=x2+y2r2F(x,y)=x^2+y^2-r^2

中,根据sgn(F(xn+1,yn+12))\mathrm{sgn}(F(x_n+1,y_n+\frac{1}{2}))判断是哪一种情况。令

pn=F(xn+1,yn+12)=xn2+yn2+2xnyn54r2p_n=F(x_n+1,y_n+\frac{1}{2})=x_n^2+y_n^2+2x_n-y_n-\frac{5}{4}-r^2

pn+1p_{n+1}的可能是:

pn+1={F(xn+2,yn12),if pn<0F(xn+2,yn32),if pn>0={pn+2xn+1+1,if pn<0pn+2xn+12yn+1+1,if pn>0\begin{aligned}p_{n+1}&= \begin{cases} F(x_n+2,y_n-\frac{1}{2}),\quad\text{if } p_n<0 \\ F(x_n+2,y_n-\frac{3}{2}),\quad\text{if } p_n>0 \end{cases} \\ &=\begin{cases} p_n+2x_{n+1}+1,\quad\text{if } p_n<0 \\ p_n+2x_{n+1}-2y_{n+1}+1,\quad\text{if } p_n>0 \end{cases} \end{aligned}

再由p0=F(0,r)=54rp_0=F(0,r)=\frac{5}{4}-r,就建立了{pn}\{p_n\}的递推关系。

这里有一个细节:由于递推过程中都是整数计算,每次最少要±1\pm 1,直接取p0=1rp_0=1-r并不影响结果,还可以避免浮点运算。

上述辅助变量pn<0p_n<0,则yn+1=yny_{n+1}=y_n,反之yn+1=yn1y_{n+1}=y_n-1。从p0p_0开始逐步迭代,生成第一个区域内的像素点,同时在另外七个对称区域内作图即可。



四分之一椭圆的生成

对称性与边界

椭圆的对称性不如圆,只能用四分生成。而一个象限内,并不是在xx方向上一直递增就可以,需要划分成不均等的两个部分,分别处理。

在分界线上方迭代时,点向右移动,下方则是向下移动。那么这个分界线应当取在椭圆的切线斜率为1-1处。

将椭圆方程写成

b2x2+a2y2=a2b2b^2x^2+a^2y^2=a^2b^2

为了求切线斜率就要求导,但化成显函数求导的话,结果还是要带有根式。所以用隐函数求导法,对两边同时取全微分得到:

2b2x+2a2ydydx=02b^2x+2a^2y\cdot\frac{\mathrm{d}y}{\mathrm{d}x}=0

则令切线斜率dydx=b2xa2y=1\frac{\mathrm{d}y}{\mathrm{d}x}=-\frac{b^2x}{a^2y}=-1可以解得分界条件是

b2x=a2yb^2x=a^2y

以这条直线为界,上方的部分是不断向右迭代的过程,下方的部分是不断向下迭代的过程。


迭代计算

b2x=a2yb^2x=a^2y为界,

  • 在上方:从(0,b)(0,b)开始,有递推式

    p1k+1={p1k+2b2xk+1+b2,if p1k<0p1k+2b2xk+12a2yk+1+b2,if p1k>0p1_{k+1}=\begin{cases} p1_k+2b^2x_{k+1}+b^2,\quad\text{if }p1_k<0 \\ p1_k+2b^2x_{k+1}-2a^2y_{k+1}+b^2,\quad\text{if }p1_k>0 \end{cases}

  • 在下方:从边界点开始,有递推式

    p2k+1={p2k2a2yk+1+a2,if p2k>0p2k+2b2xk+12a2yk+1+a2,if p2k<0p2_{k+1}=\begin{cases} p2_k-2a^2y_{k+1}+a^2,\quad\text{if }p2_k>0 \\ p2_k+2b^2x_{k+1}-2a^2y_{k+1}+a^2,\quad\text{if }p2_k<0 \end{cases}

就可以生成第一象限内的四分之一椭圆。



代码实现

因为基本上和 Bresenham 算法生成直线一样,就不放代码了,可以参考我的 GitHub 仓库:

点击链接:greyishsong/basic_cg



注释


  1. 一种快速实现非线性函数计算的方法,可以只用移位和加减法计算三角函数这样的超越函数值,在硬件设计中比较常用。 ↩︎