计算机图形学-抛物线的中点Bresenham算法

抛物线$f(x)=a{x}^2$的中点Bresenham算法

语言：matlab

画图：plot

1 抛物线的特征

通常定义抛物线为到一条直线（准线）和直线外一点（焦点）距离相等的点的集合。这里只讨论顶点为原点，沿纵坐标轴对称且开口向上的情况。而对于其他情况可以通过图形的平移和旋转等线性变换得到。其描述方程如下：
$$F(x,y)=y-a{x}^2(a>0)$$
与椭圆不同，抛物线是无边界的非封闭图形，若要在屏幕上绘制，必须给定坐标范围，以绘制指定抛物线的一个片段。可以在函数中设置参数$xm$，则横坐标约束其范围为$[-xm,xm]$。

抛物线关于纵坐标轴对称，故只需绘制其第一象限内的点，第二象限中的点可以通过对称得到。

为确定最大位移方向，考虑抛物线的斜率范围。在第一象限，其上一点$(x,y)$处的斜率为：
$$k(x)=\frac{\partial F(x,y)}{\partial x}=2ax\in [0,+\infty )$$

又由于斜率的变化率为：
$$\frac{dk(x)}{dx}=\frac{d(2ax)}{dx}=2a>0$$

所以在第一象限内，抛物线斜率从0开始随$x$递增至正无穷。用斜率为1的点对图形进行划分。容易解出，当斜率为1时，$x=\frac{1}{2a}$。只需沿$x$轴绘制图形，其$0<x<\frac{1}{2a}$时，最大位移方向为$x$方向；$\frac{1}{2a}<x<xm$时，最大位移方向为$y$方向。

2 算法推导过程

假定当前与抛物线距离最近者已确定为$P(x_i,y_i)$那么在抛物线前部分时，下一候选点是$P_d(x_i+1,y_i)$和$P_d(x_i+1,y_i+1)$;而在抛物线的后半部分时，下一候选点是$P_l(x_i,y_i+1)$和$P_r(x_i+1,y_i+1)$。仍然使用中点进行判别候选点。

2.1 推导前半部分的抛物线绘制公式

对于$0<x<\frac{1}{2a}$:构造判别式
$$d_{li}=F(x_i+1,y_i+0.5)=y_i+0.5-a(x_i+1{)}^2$$
若$d_{li}\ge 0$，中点在抛物线上方，应选取$P_d(x_i+1,y_i)$，反之选取$P_u(x_i+1,y_i+1)$。

误差项递推：

在$d_{li}\ge 0$时，应计算：
$$\begin{gathered} d_{l(i+1)}=F(x_i+2,y_i+0.5) \\ =y_i+0.5-a(x_i+2{)}^2 \\ =d_{li}+a[(x_i+1{)}^2-(x_i+2{)}^2] \\ =d_{li}-2a{x}_i-3a \end{gathered}$$
在$d_{li}<0$时，应计算：
$$\begin{gathered} d_{l(i+1)}=F(x_i+2,y_i+1.5) \\ =y_i+1.5-a(x_i+2{)}^2 \\ =d_{li}-2a{x}_i-3a+1 \end{gathered}$$
计算判别式的初始值。弧起点为$(0,0)$，因此第一个中点为$(1,0.5)$，对应的判别式为：
$$d_{l0}=y_0+0.5-a(x_0+1{)}^2=0.5-a$$

2.2 推导后半部分的抛物线绘制公式

对于$\frac{1}{2a}<x<xm$:构造判别式
$$d_{ri}=F(x_i+0.5,y_i+1)=y_i+1-a(x_i+0.5{)}^2$$
若$d_{ri}\ge 0$，中点在抛物线（左）上方，应选取$P_l(x_i,y_i+1)$，反之选取$P_r(x_i+1,y_i+1)$。

误差项递推：

在$d_{ri}\ge 0$时，应计算：
$$\begin{gathered} d_{r(i+1)}=F(x_i+0.5,y_i+2) \\ =y_i+2-a(x_i+0.5{)}^2 \\ =d_{ri}+1 \end{gathered}$$
在$d_{ri}<0$时，应计算：
$$\begin{gathered} d_{r(i+1)}=F(x_i+1.5,y_i+2) \\ =y_i+2-a(x_i+1.5{)}^2 \\ =d_{ri}-2a{x}_i-2a+1 \end{gathered}$$
计算判别式的初始值。弧起点横坐标为$\lceil \frac{1}{2a}\rceil$,对应的判别式为：
$$d_{r0}=y_0+1-a(x_0+0.5{)}^2=1-a\lceil \frac{1}{2a}\rceil -0.25a$$

算法代码（matlab）

function DrawBresenhamCurve(a,xm)
  div=0.5/a;
  x=0,y=0;
  d_pre=0.5-a;
  d_post=1-a*ceil(0.5/a)-0.25*a;
  while x<=xm
    plot(x,y);
    hold on;
    plot(-x,y);
    hold on;
    if x<div
      tmp=-2*a*x-3*a;
      x++;
      if d_pre<0
        y++;
        d_pre=d_pre+tmp+1;
      else
        d_pre=d_pre+tmp;
      endif
    else
      tmp=-2*a*x-2*a+1;
      y++;
      if d_post<0
        d_post=d_post+1;
      else
        d_post=d_post+tmp;
        x++;
      endif
    endif
  endwhile
end

运行情况：

测试代码：

a=0.01;
xm=100;
DrawBresenhamCurve(a,xm);
X=linspace(-xm,xm,1000);
Y=a.*X.^2.;
plot(X,Y);
grid on;

运行结果（点为算法生成的点集，曲线为所拟合的实际曲线）：