光栅图形学(一):直线绘制

本来对这学期的图形学课程抱了很高的期待,可惜课程质量不尽人意。求人不如求己,于是选择自己作一点学习记录。

光栅图形学产生于GPU未产生的年代,有时要用孱弱的(也许是没有的)浮点算力解决问题,这个过程造就了各种以避开浮点运算为目的的经典算法。

直线的绘制大概是绘图最常见的过程,完成此任务的Bresenham算法可以成为中点圆、中点椭圆绘制的思路基础。

从理论上画一条线


让计算机显示一个图形,操作起来比输出到终端要麻烦许多(如果把一些封装层次非常高的库除外的话)。建模时的思路一般都是先提取问题的理论部分,考虑用什么样的思路,再设计如何实现。所以这里先不提图形API,只推导如何将直线绘制到一个理论的屏幕上。

为了简化问题,在推导时维持一个假设,姑且称为假设()(*)

给定线段的两个端点(x1,y1),  (x2,y2)(x_1,y_1),\;(x_2,y_2),假设斜率m<1m<1、前者是左端点。

其他的情况思路一致,不影响画线的思路,之后只要对公式稍加变换,就可以得到结果。

直线的描述和朴素绘制思路

数学上的直线是连续的,而光栅图形学面对的光栅扫描显示器,每次显示一个二维的像素阵列。所以将一条连续的、零宽度的直线绘制出来,就需要将其近似为一组像素点。(向量显示器就不用这样,只管给出起点终点就行了)

直线方程可以有几种形式,我看来基本上是两类:一是像一般式那样写成二(三)元一次不定方程,二是像斜截式这样写成y=f(x)y=f(x)的函数式。这也可以视作是在平面上画线的两种思路:

  • 检查整个画面范围的每一个点,判定它是否应被近似到直线上。
  • 从一个初始点开始,迭代计算出画面范围内的近似点。

也许用第一种思路绘制直线看起来有点蠢,但是它还真是暴力绘制一些奇怪曲线的可取方式[1]

以现今个人设备的算力而言,用这种方法绘制一张图像还不算过分,用来绘制动态画面绝对是不能接受的——想想每个形状都要全屏验算一遍,一个画面成千上万个三角形面片要验算成千上万次,简直是愚公移山式的任务。

那么如果给定直线

y=mx+by=mx+b

可以作这样的迭代:

  1. x=x1,  y=y1x=x_1,\;y=y_1
  2. 取一个很小的步长dxdx
  3. 计算([x+dx],f([x+dx]))([x+dx],f([x+dx]))这个点并绘制([x][x]是取整函数)
  4. 以3中的结果为基础再次迭代,直到计算结果是(x2,y2)(x_2,y_2)

无论屏幕分辨率多高,只要步长dxdx取得足够小,那么绘制的结果在xx方向上就是连续的。

这种方法无疑带来大量的浮点运算:xx方向的浮点加、yy方向上的浮点乘,还有每一步的取整。乘法已经足够耗时了,浮点乘法还要先作浮点对齐,对古老的计算机而言是不小的负担。



改进方法:数值微分法DDA

像素阵列是离散的,如果dx=1dx=1,就足够保证xx方向上连续了。那么变量xx就不需要是浮点数。

每次作迭代xn+1=xn+1x_{n+1}=x_n+1时,yy方向的增量应当是1m\frac{1}{m}

确定了yy方向的增量,每次用乘法计算yy的值就不再必要了,可以直接yn+1=yn+1my_{n+1}=y_n+\frac{1}{m}

这样一来,每次迭代的计算量减少到了两次加法。具体操作是:

  1. x=x1,  y=y1x=x_1,\;y=y_1
  2. 绘制当前点,计算下一个点(x+1,[y+1m])(x+1,[y+\frac{1}{m}])
  3. 再次迭代,直到当前点超出了给定范围

但是,DDA还是要在一个方向上,反复进行浮点运算。一方面性能开销大,另一方面如果直线很长,浮点数可能有所误差。如果想要继续改进,就势必要完全避开浮点运算。



Bresenham算法

变取整为判定

在DDA计算的过程中,yy变量保存了yy方向上线段的真实位置。这个位置本身并没有用,因为迭代的结果要么是yn+1=yny_{n+1}=y_n,要么是yn+1=yn+1y_{n+1}=y_n+1,位置只是用来判定的(取整也可以视为一种判定手段)。

如果能够不用取整而判定yn+1y_{n+1}的情况,那么保留这个位置就没有必要了。

现在考虑换一种判定方法:

  • 当迭代到第nn个点时,记第n+1n+1个点的真实位置为yy
  • 计算两个距离dlower=yynd_{\text{lower}}=y-y_ndupper=yn+1yd_{\text{upper}}=y_n+1-y
  • 前一个距离是真实位置到yny_n的距离,后一个是真实位置到yn+1y_n+1的距离。
  • dlower>dupperd_{\text{lower}}>d_{\text{upper}},则yyyn+1y^n+1更近,下一步应当yn+1=yn+1y^{n+1}=y^n+1,反之亦然。

只走到这一步的话,Bresenham算法也不过提供了另一种思路——直接计算dlowerd_{\text{lower}}dupperd_{\text{upper}}比进行DDA的取整操作还要麻烦。

因此,Bresenham算法的第二个任务就是用整数迭代的方式算出sgn(dlowerdupper)\mathrm{sgn}(d_{\text{lower}}-d_{\text{upper}})

辅助判定

在假设x1<x2x_1<x_2下,Δx=x2x1>0\Delta x=x_2-x_1>0成立,则sgn(Δx(dlowerdupper))=sgn(dlowerdupper)\mathrm{sgn}(\Delta x(d_{\text{lower}}-d_{\text{upper}}))=\mathrm{sgn}(d_{\text{lower}}-d_{\text{upper}})。构造一个辅助变量

pn=Δx(dlowerdupper)=2xnΔy2ynΔx+(2ΔyΔx)\begin{aligned} p_n & =\Delta x(d_{\text{lower}}-d_{\text{upper}}) \\ & = 2x_n\Delta y-2y_n\Delta x+(2\Delta y-\Delta x) \end{aligned}

再令p0=2ΔyΔxp_0=2\Delta y-\Delta x,可得{pn}\{p_n\}数列的递推关系是:

pn+1={pn+2Δy,pn>0pn+2Δy2Δx,pn0p_{n+1}=\begin{cases} p_n+2\Delta y,\quad p_n>0 \\ p_n+2\Delta y-2\Delta x,\quad p_n\geq 0 \end{cases}

大功告成,这个递推式中Δx,Δy\Delta x,\Delta y都是整数,在整个计算过程中完全不需要浮点运算。于是我们建立了在假设()(*)下的直线绘制公式。



一般化情况

在前面的推导过程都建立在假设()(*)的基础上,如果斜率m>1m>1,或者斜率是负的,递推式也应当加以变化。

基本的思路是:

  1. 将直线变换到一个坐标系xOyx'Oy',这个坐标系中,直线符合假设()(*),得到在xOyx'Oy'系下,含有x,y,mx',y',m'的计算公式。
  2. 再进行逆变换,将含有x,y,mx',y',m'的公式变换为含有x,y,mx,y,m的公式,就得到了xOyxOy系下的相应公式。

斜率为负值

对于1<m<0-1<m<0的情况,仍然保证x1<x2x_1<x_2(从左往右),考虑一个线性变换

v=[1001]v\vec{v}'=\begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}\vec{v}

这个变换就是y=yy'=-y,那么原先在xOyxOy坐标系下斜率为负值的直线,变换后在xOyx'Oy'系下斜率就是正的,可以用原来的递推式计算。

那么将计算结果变换回来,就应当左乘一个变换矩阵TT的逆矩阵T1T^{-1},而TT的逆矩阵是自身,即作变换

v=[1001]v\vec{v}=\begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}\vec{v}'

在逆变换下,Δy=Δy\Delta y=-\Delta y',从而m=mm=-m'

  • DDA算法是完全不受变换影响的。因为dx=1dx=1不变,yy方向上的迭代式变为yn+1=yn1myn+1=yn+1m-y_{n+1}=-y_n-\frac{1}{m}\Leftrightarrow y_{n+1}=y_n+\frac{1}{m},计算过程完全一样。
  • Bresenham算法则需要作一点改变:
    • 由于Δy=Δy\Delta y=-\Delta y',由pnp_n'变为pnp_n时,所有的Δy\Delta y'都变为Δy-\Delta y。也就是

      pn+1={pn2Δy,pn>0pn2Δy2Δx,pn<0p_{n+1}=\begin{cases} p_n-2\Delta y,\quad p_n>0 \\ p_n-2\Delta y-2\Delta x,\quad p_n<0 \end{cases}

    • 考虑到在m<0m<0的情况下,Δy<0\Delta y<0,这也等同于对Δy\Delta y取绝对值。

斜率大于1

这时改用变换矩阵

T=[0110]T=\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}

xx轴和yy轴互换,从而在xOyx'Oy'系中,m<1m'<1。这个逆变换也满足T1=TT^{-1}=T,那么公式的逆变换结果,就是原公式中所有的xxyy互换。





程序实现

DDA

DDA的思路非常直接,实现起来也容易。首先在输入点对时作好判定,处理出两个方向的增量:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
void DDALine::inputPoints()
{
std::scanf("%d%d%d%d", &x1, &y1, &x2, &y2);
float m = ((float)y2-y1) / (x2-x1);
if (std::fabs(m) < 1)
{
if (x1 > x2)
{
std::swap(x1, x2);
std::swap(y1, y2);
}
dx = 1.0;
dy = m;
}
else
{
if (y1 > y2)
{
std::swap(x1, x2);
std::swap(y1, y2);
}
dx = 1.0 / m;
dy = 1.0;
}
x = (float)x1;
y = (float)y1;
}

这里保证是从左往右(或者从下往上)画线,即有一个方向的增量是1。

然后就是不断迭代更新点的位置:

1
2
3
4
5
6
7
8
void DDALine::updateStatus()
{
if ((int)std::round(x) == x2 && (int)std::round(y) == y2)
return;
x += dx;
y += dy;
std::printf("(%d, %d)\n", (int)std::round(x), (int)std::round(y));
}

注意,std::round这个函数至少要C++ 11才支持,否则需要自己实现。

Bresenham

输入时需要注意一些标记的处理,顺便把dx(Δx\Delta x)和dy(Δy\Delta y)直接倍增,省得迭代的时候再乘。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
void BresenhamLine::inputPoints()
{
std::scanf("%d%d%d%d", &x1, &y1, &x2, &y2);
float m = ((float)y2-y1) / (x2-x1);
if (std::fabs(m) < 1)
{
if (x1 > x2)
std::swap(x1, x2), std::swap(y1, y2);
dx = x2 - x1;
dy = std::abs(y2 - y1); //这里就是对不同情况的公式作了变换之后的结果
p = dy*2 - dx;
dy *= 2;
dx *= 2;
slope = 1; //用来标示斜率的绝对值<1
}
else
{
if (y1 > y2)
std::swap(x1, x2), std::swap(y1, y2);
dx = std::abs(x2 - x1);
dy = y2 - y1;
p = dx*2 - dy;
dx *= 2;
dy *= 2;
slope = 0;
}
x = x1;
y = y1;
sgn = (m > 0.0f)? 1 : -1;
}

在更新状态时,按斜率绝对值不同分类:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
void BresenhamLine::updateStatus()
{
if (x == x2 && y == y2)
return;
if (slope)
{
++x;
if (p < 0)
p = p + dy;
else
{
y += sgn;
p = p + dy - dx;
}
}
else //这就是对斜率>1的情况通过逆变换得到的公式
{
++y;
if (p < 0)
p = p + dx;
else
{
x += sgn;
p = p + dx - dy;
}
}
std::printf("(%d, %d)\n", x, y);
}

渲染循环

OpenGL是在一个渲染循环中进行绘制,对帧的刷新是不会停止的。为了在绘制完一条线段后保持画面,可以引入几个状态,将计算新的点视为更新状态的操作,在计算完所有点之后进入等待状态。(所以上面的方法都叫做updateStatus

完整的代码可以在我的GitHub找到:找到项目greyishsong/basic_cg3_line目录。这些代码中加入了更新状态和进行判定的部分。

项目链接:greyishsong/basic_cg





注释

[1] 比如说我见过一个叫做GrafEq的二维函数图像绘制软件,可以绘制像sinx+cosy+sinxcosy+ysinx>0\sin x+\cos y+\sin\lfloor x\cos y+y\sin x\rfloor>0这样的奇怪不等式(有些甚至能难倒MATLAB、Mathmatica这样的软件),似乎就是用这样的方法来完成。