二维变换概念应用于对象几何描述并改变其位置、方向或者大小的变换叫做几何变换，有时候也被叫做建模变换。而本文仅讨论平面中的几何变换，即二维变换。矩阵表示和齐次坐标对于普通的2x2矩阵，我们总是要将平移项与其它变换对应的矩阵写成不同规格，为了统一形式且方便运算，我们需要将2x2的矩阵扩展到3x3。此时，二为坐标必须用三元向量来表示。标准实现技术是将二维坐标$(x,y)$扩充到三维$(x_h,y_h,h)$，这称为齐次坐标，这个过程就被叫做齐次化。对于每一个维度有$$ x = \frac{x_h}{h},y=\frac{y_h}{h} $$其中非零值$h$被称为齐次参数。显然对于齐次参数，可以有无数个非零值，同样也意味着有无数个等价的齐次表达式。既然如此，为了方便计算，不妨令$h=1$。平移变换对于平移变换，我们可有参数方程$$ \begin{cases} x^{'}=x+\delta x \\ y^{'} = y+\delta y \end{cases} $$我们将方程组转换为矩阵的形式$$ \begin{bmatrix} x'& y'& 1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0\\ 0& 1& 0\\ \delta x& \delta y & 1 \end{bmatrix} $$旋转变换对于一个复杂的图形的旋转，我们可以看做是多个点在同时旋转。因此只需要研究出一个点的旋转变换方法就可以了。我们不妨令点为平面任意一点，其绕原点进行旋转变换，该点的运动轨迹一定是一个以该点到原点连线为半径的圆弧。那么问题就简单了，对于旋转任意角度，我们只需要用圆的参数方程就能搞定，$$ \begin{cases} x = \cos{\theta} \\ y = \sin{\theta} \end{cases} $$对于我们假设旋转了$\alpha$的弧度（逆时针为正方向），则有$$ \begin{cases} x' = \cos{(\theta+\alpha)}=\cos{\theta}\cos{\alpha}-\sin{\theta}\sin{\alpha}=x\cos{\alpha}-y\sin{\alpha} \\ y' = \sin{(\theta+\alpha)}=\sin{\theta}\cos{\alpha}+\cos{\theta}\sin{\alpha}=x\sin{\alpha}+y\cos{\alpha} \end{cases} $$很显然对应矩阵形式为$$ \begin{bmatrix} x'&y'&1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} \cos{\alpha}& \sin{\alpha}& 0 \\ -\sin{\alpha}& \cos{\alpha}& 0 \\ 0& 0& 1 \end{bmatrix} $$那么问题来了，如果所绕的旋转点不是原点怎么办呢？前面已经讲过平移变换了，只需要平移坐标系原点至该点（移轴），再进行旋转，最后平移回去就行了。不妨令被围绕点坐标为$P(x_1,y_1)$，则该流程的矩阵运算如下：$$ \begin{bmatrix} x'& y'& 1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0\\ 0& 1& 0\\ -x_1& -y_1 & 1 \end{bmatrix} \cdot \begin{bmatrix} \cos{\alpha}& \sin{\alpha}& 0 \\ -\sin{\alpha}& \cos{\alpha}& 0 \\ 0& 0& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0\\ 0& 1& 0\\ x_1& y_1 & 1 \end{bmatrix} $$如果你对平移坐标系难以理解，不妨将坐标轴与平面想象成两个分离的东西。以坐标轴往右移动为例，则对于平面上的点来讲，就相当于坐标轴不动点向左平移。（平移坐标系只是移动的轴，不带平面上其他点，否则你会得到相反的结果！）缩放变换首先给出参数方程$$ \begin{cases} x' = s_xx \\ y' = s_yy \end{cases} $$对应矩阵运算为$$ \begin{bmatrix} x'&y'&1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} s_x& 0& 0 \\ 0& s_y& 0 \\ 0& 0& 1 \end{bmatrix} $$对称变换关于对称变换可以是轴对称或者点对称。我们先来看关于y轴对称：只需要纵坐标不变，横坐标取相反数即可。矩阵运算如下$$ \begin{bmatrix} x'&y'&1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} -1& 0& 0 \\ 0& 1& 0 \\ 0& 0& 1 \end{bmatrix} $$同理可得关于y轴对称矩阵运算$$ \begin{bmatrix} x'&y'&1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0 \\ 0& -1& 0 \\ 0& 0& 1 \end{bmatrix} $$那么对于原点对称就有$$ \begin{bmatrix} x'&y'&1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} -1& 0& 0 \\ 0& -1& 0 \\ 0& 0& 1 \end{bmatrix} $$如果是对于平面中任意一点$P(x_1,y_1)$对称，则平移坐标系原点至该点，然后进行关于原点的对称，再平移回去就行。$$ \begin{bmatrix} x'& y'& 1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0\\ 0& 1& 0\\ -x_1& -y_1 & 1 \end{bmatrix} \cdot \begin{bmatrix} -1& 0& 0 \\ 0& -1& 0 \\ 0& 0& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0\\ 0& 1& 0\\ x_1& y_1 & 1 \end{bmatrix} $$如果是关于平面内任意一条直线对称，那比较麻烦了，你需要先平移坐标系到直线的一点，并旋转坐标系使x轴与直线重合，然后进行关于x轴的对称变换，再旋转回去，然后再平移回去。特别地，如果直线斜率不存在，平移后直线与y轴重合，那么直接进行关于y轴对称然后再平移回去即可。注意该过程的顺序，因为矩阵并不满足交换律！这里实际上是有两个大坑的，我们不妨假设直线经过的两点分别为$$ A(x_1,y_1),B(x_2,y_2) $$不失一般性，我们令$$ x_1<=x_2 $$我们在上面已经说过直线斜率不存在的情况了，这里仅讨论直线斜率存在的情况。直线斜率为$$ k = \frac{\delta x}{\delta y}=\frac{y_2-y_1}{x_2-x_1} $$我们想要得到直线与x轴的夹角，可以利用反正切函数$$ \theta = \arctan{k} $$但是事实就是如此吗？注意反正切函数的取值范围为$(-\frac{\pi}{2},\frac{\pi}{2})$，而直线倾斜角范围是$[0,\pi]$，即使挖去了$\frac{\pi}{2}$这个点（因为我们不讨论这种情况），范围仍不一致！我们记倾斜角为$\alpha$，则有$$ \alpha= \begin{cases} \theta,& 0\le \theta<\frac{\pi}{2} \\ \pi+\theta, &\theta<0 \end{cases} $$现在开始进入第二个坑了，我们进行旋转倾斜角的时候，是顺时针还是逆时针？我们不妨站在坐标系的角度来看，逆时针旋转倾斜角的角度，就相当于顺时针旋转图形这个角度。也就是说，我们在进行第一次旋转变换时输入的角度参数应该为$-\alpha$综上所述，矩阵运算为$$ \begin{bmatrix} x'& y'& 1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0\\ 0& 1& 0\\ -x_1& -y_1 & 1 \end{bmatrix} \cdot \begin{bmatrix} \cos{\alpha}& -\sin{\alpha}& 0 \\ \sin{\alpha}& \cos{\alpha}& 0 \\ 0& 0& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0 \\ 0& -1& 0 \\ 0& 0& 1 \end{bmatrix} \cdot \begin{bmatrix} \cos{\alpha}& \sin{\alpha}& 0 \\ -\sin{\alpha}& \cos{\alpha}& 0 \\ 0& 0& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0\\ 0& 1& 0\\ x_1& y_1 & 1 \end{bmatrix} $$错切变换错切变换实际上是物体在投影平面上非垂直投影的结果。一般为水平错切和垂直错切，也可以同时对两个方向进行错切。我们先来看看沿y轴方向的错切显然有方程$$ \begin{cases} x' = x \\ y' = y+x\tan{\theta} \end{cases} $$有矩阵运算$$ \begin{bmatrix} x'&y'&1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& \tan{\theta}& 0 \\ 0& 1& 0 \\ 0& 0& 1 \end{bmatrix} $$同理可得沿着x方向的错切矩阵运算$$ \begin{bmatrix} x'&y'&1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& 0& 0 \\ \tan{\theta}& 1& 0 \\ 0& 0& 1 \end{bmatrix} $$同时沿着x方向和y方向的错切$$ \begin{bmatrix} x'&y'&1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} 1& \tan{\alpha}& 0 \\ \tan{\theta}& 1& 0 \\ 0& 0& 1 \end{bmatrix} \space \alpha,\theta \in (-\frac{\pi}{2},\frac{\pi}{2}) $$刚体变换如果一个矩阵仅包含平移参数和旋转参数，那么它就是一个刚体变换矩阵$$ \begin{bmatrix} x'& y'& 1 \end{bmatrix} = \begin{bmatrix} x& y& 1 \end{bmatrix} \cdot \begin{bmatrix} r_{xx}& r_{yx}& 0\\ r_{xy}& r_{yy}& 0\\ tr_{x}& tr_{y}& 1 \end{bmatrix} $$对于刚体变换左上角的2x2矩阵满足正交矩阵的特性。也就是说对于两个列向量$\begin{bmatrix}r_{xx}& r_{xy} \end{bmatrix}^{T}$和$\begin{bmatrix}r_{yx}& r_{yy} \end{bmatrix}^{T}$（或者两个行向量）形成单位向量的正交组，这样的向量也称为正交向量组。对于每个向量具有单位长度，并且数量积为0：$$ r_{xx}^2+r_{xy}^2=r_{yx}^2+r_{yy}^2=1 \\ r_{xx}r_{yx}+r_{xy}r_{yy} = 0 $$如果这些向量通过旋转矩阵进行变换，则可得出x方向和y方向的单位向量$$ \begin{bmatrix} r_{xx}& r_{xy}& 1 \end{bmatrix} \cdot \begin{bmatrix} r_{xx}& r_{yx}& 0\\ r_{xy}& r_{yy}& 0\\ 0& 0& 1 \end{bmatrix} =\begin{bmatrix} 1& 0& 1 \end{bmatrix} $$$$ \begin{bmatrix} r_{yx}& r_{yy}& 1 \end{bmatrix} \cdot \begin{bmatrix} r_{xx}& r_{yx}& 0\\ r_{xy}& r_{yy}& 0\\ 0& 0& 1 \end{bmatrix} =\begin{bmatrix} 0& 1& 1 \end{bmatrix} $$代码对于二维变换的代码相对来讲较麻烦一些，我不太喜欢OpenGL所带的矩阵运算操作，因此我自己写了一套简单的矩阵模版类。由于代码相对较长，我这里仅放出这几个变换的代码片段。class Polygon { public: Transform::Matrix<GLdouble> *vex; Polygon(Pts vex) { this->vex = new Transform::Matrix<GLdouble>(vex.size(), 3, {{}}, 0); for (int i = 0; i < vex.size(); i++) { (*this->vex)[i][0] = vex[i].x; (*this->vex)[i][1] = vex[i].y; (*this->vex)[i][2] = 1; } } ~Polygon() { delete vex; } Polygon *Translate(const Point &end) { vector<vector<GLdouble>> temp{{1, 0, 0}, {0, 1, 0}, {end.x, end.y, 1}}; Transform::Matrix<GLdouble> t(3, 3, temp); (*vex) = (*vex) * t; return this; } Polygon *Draw() { // cout << "Start" << endl; glBegin(GL_POLYGON); for (GLint i = 0; i < vex->CountRow(); i++) { // cout << i << ": " << (*vex)[i][0] << " " << (*vex)[i][1] << endl; glVertex2d((*vex)[i][0], (*vex)[i][1]); } glEnd(); return this; } Polygon *Rotate(Point p, GLdouble radian) { Transform::Matrix<GLdouble> t1(3, 3, {{1, 0, 0}, {0, 1, 0}, {-p.x, -p.y, 1}}); Transform::Matrix<GLdouble> t2(3, 3, {{cos(radian), sin(radian), 0}, {-sin(radian), cos(radian), 0}, {0, 0, 1}}); Transform::Matrix<GLdouble> t3(3, 3, {{1, 0, 0}, {0, 1, 0}, {p.x, p.y, 1}}); (*vex) = (*vex) * t1 * t2 * t3; return this; } Polygon *Scale(GLdouble x, GLdouble y) { if (x <= 0 || y <= 0) { cout << "WARNING: Invalid parameters" << endl; return this; } Transform::Matrix<GLdouble> t(3, 3, { {x, 0, 0}, {0, y, 0}, {0, 0, 1}, }); (*vex) = (*vex) * t; return this; } Polygon *PointReflect(const Point &p) { Transform::Matrix<GLdouble> t1(3, 3, {{1, 0, 0}, {0, 1, 0}, {-p.x, -p.y, 1}}); Transform::Matrix<GLdouble> t2(3, 3, { {-1, 0, 0}, {0, -1, 0}, {0, 0, 1}, }); Transform::Matrix<GLdouble> t3(3, 3, {{1, 0, 0}, {0, 1, 0}, {p.x, p.y, 1}}); *vex = (*vex) * t1 * t2 * t3; return this; } Polygon *LineReflect(const Point &start, const Point &end, bool showAxis = false) { GLdouble x1 = start.x, x2 = end.x, y1 = start.y, y2 = end.y; if (x1 > x2) swap(x1, x2), swap(y1, y2); if (showAxis) { glBegin(GL_LINES); glVertex2d(x1, y1); glVertex2d(x2, y2); glEnd(); } this->Translate({-x1, -y1}); GLdouble dy = y2 - y1; GLdouble dx = x2 - x1; if (dx == 0) { Transform::Matrix<GLdouble> t(3, 3, {{-1, 0, 0}, {0, 1, 0}, {0, 0, 1}}); *vex = (*vex) * t; this->Translate({x1, y1}); return this; } GLdouble rad = atan(dy / dx); if (rad < 0) rad = 4 * atan(1) + rad; this->Rotate({0, 0}, -rad); Transform::Matrix<GLdouble> t1(3, 3, {{1, 0, 0}, {0, -1, 0}, {0, 0, 1}}); *vex = (*vex) * t1; this->Rotate({0, 0}, rad); this->Translate({x1, y1}); return this; } Polygon *Shear(const Point &sh) { Transform::Matrix<GLdouble> t1(3, 3, {{1, sh.y, 0}, {sh.x, 1, 0}, {0, 0, 1}}); *vex = (*vex) * t1; return this; } };完整代码请参考我的Github仓库: 2D变换-Github
扫描转换椭圆与圆的绘制前提我们先来设计一种数据类型class Point{ public: double x,y; Point(){this->x = 0; this->y = 0;} Point(double x, double y){this->x=x;this->y=y;} };注意此处的类型为double这一点很重要，因为在某些迭代过程中设置为int会导致坐标数值不变，然后进入死循环或者绘制不出来。如果需要int类型，那么无需担心，double到int之间会有一个隐式转换，直接将Point类型中属性值赋给对应的变量即可。（如果你想知道哪些情况下使用int会无法绘制请参考椭圆部分的DDA算法）圆基础知识在直角坐标系中，圆心在原点上的圆有如下方程：$$ x^2+y^2 = R^2 $$其中R为圆的半径。若圆心不在原点上，则有更一般形式的方程。假设圆心为(a,b)则圆的方程为：$$ (x-a)^2+(y-b)^2 = R^2 $$同样的，对于圆心在原点上的圆我们有参数方程：$$ \begin{cases} x=R\cos \theta \\ y = R\sin \theta \end{cases} $$对于圆心不在原点上的圆有参数方程：$$ \begin{cases} x=a+R\cos \theta \\ y = b+R\sin \theta \end{cases} $$将$\theta$离散化可得$$ \theta = 2\pi \frac{i}{n},i\in[0,n] $$然而我们在计算的时候并不喜欢圆心不在原点上的圆，所以在进行扫描转换之前，我们先平移坐标系原点至圆心，当扫描转换完成后再平移回去。故而在接下来的讨论中，我们只使用圆心在原点上的方程。此外，圆还有如下几种对称性质：关于 $y=\pm x$ 对称关于x轴对称关于y轴对称这几条对称轴将圆分割成了八等份，也就是说，给出圆周上任意一个点，就能通过数学计算得到其余的7个点。这对于圆的绘制就简单了许多，我们只需要计算出八分之一的圆弧，就能得到整个圆。在开始之前我们不妨约定这样一个函数（C++风格伪代码，实际代码根据你自己的API来调用相应的函数）/* 以原点为圆心，通过八分之一圆上的一点绘制其余七个对称点 * p为绘制的点 * center为圆心 * color为绘制颜色 */ void draw_circle_point(const Point& p, const Point& center, Color color){ int x = p.x; int y = p.y; int ox = center.x, oy = center.y; setPixel(x+ox,y+oy,color); //绘制像素 setPixel(y+ox,x+oy,color); setPixel(-y+ox,x+oy,color); setPixel(-x+ox,y+oy,color); setPixel(-x+ox,-y+oy,color); setPixel(-y+ox,-x+oy,color); setPixel(y+ox,-x+oy,color); setPixel(x+ox,-y+oy,color); }此外在接下来的计算中可能会用到三角函数的和差公式，请牢记这一点：$$ \sin{(\alpha \pm \beta)}=\sin{\alpha}\cos{\beta}\pm \cos{\alpha}\sin{\beta} \\ \cos{(\alpha \pm \beta)} = \cos{\alpha}\cos{\beta} \mp \sin{\alpha}\sin{\beta} $$直接计算法直接生成有两种方式，一种是根据直角坐标系下笛卡尔方程来计算，另一种根据参数方程。笛卡尔方程直接计算我们扫描这样圆弧上的点。根据方程 $x^2+y^2=R$ 可以得到：$$ \begin{cases} x_{i+1} = x_i+1,x\in [0,\frac{R}{\sqrt{2}}] \\ \\ y_{i+1} = round(\sqrt{R^2-x_{i+1}^2}) & \end{cases} $$故有C++风格伪代码const double SQRT_2 = sqrt(2); /* * 利用笛卡尔方程直接绘制圆 * circle_center为圆心 * R为圆半径 * color为绘制的颜色 */ void draw_circle_common(const Point& circle_center, double R, Color color){ double xr = R/SQRT_2; Point p; for(p.x=0;p.x<=xr;p.x++){ p.y=sqrt(R*R-p.x*p.x); draw_circle_point(p, circle_center, color); } }参数方程直接计算我们是根据离散化的参数方程来计算的$$ \begin{cases} x=R\cos (2\pi \frac{i}{n}) \\ y = R\sin (2\pi \frac{i}{n}) \end{cases} $$C++风格伪代码如下const double PI = 4.0*atan(1.0); //圆周率 /* * 利用参数方程直接绘制圆 * circle_center为圆心 * R为圆半径 * color为绘制的颜色 * n为离散化程度，n越大越精细但越耗时 */ void draw_circle_param(const Point& circle_center, double R, Color color,int n){ int k; Point p; for(k=0;k<=n;k++){ p.x = R*cos(2*PI*k/n); p.y = R*sin(2*PI*k/n); draw_circle_point(p, circle_center, color); } } DDA算法数值微分算法(Digital Differential Analyzer, DDA)在前面绘制直线的文章中已经提到过了，现在我们来尝试使用它绘制圆。我们考虑如下图$$ \begin{cases} x_i=R\cos \alpha \\ y_i = R\sin \alpha \end{cases} $$我们让当前角度 $\alpha$ 沿着逆时针方向加上几个极小的角度量 $\theta$就得到了下一个点$$ (x_{i+1},y_{i+1}) = (R\cos{(\alpha+\theta)},\space R\sin{(\alpha + \theta)}) $$所以，对于x有如下递推式$$ x_{i+1} = R\cos{(\alpha+\theta)} = R\cos\alpha\cos\theta - R\sin\alpha\sin\theta = x_i\cos\theta-y_i\sin\theta $$同样对于y也有如下递推式$$ y_{i+1} = R\cos{(\alpha + \theta)} = y_i\cos\theta+x_i\sin\theta $$又因为当 $\theta$ 足够小的时候有$$ \cos\theta \approx 1\\ \sin\theta \approx \theta $$故递推式可以改写为$$ \begin{cases} x_{i+1} = x_i - y_i\theta \\ y_{i+1} = y_i+x_i\theta \end{cases} $$知道这一点我们可以写出C++风格伪代码：/* * DDA画圆 * circle_center为圆心 * R为圆半径 * color为绘制的颜色 */ void draw_circle_dda(const Point& circle_center, double R, Color color){ Point p; p.x = R; p.y = 0; int x_old; double delta=1.0/R, theta=0; for(theta=0; theta<=PI/4; theta+=delta){ draw_circle_point(p, circle_center, color); x_old = p.x; p.x -= p.y*delta; p.y += x_old*delta; } }中点Bresenham算法基本思路圆的隐式方程如下$$ x^2 + y^2 - R^2 = 0 $$由此构造二元函数$$ F(x,y) = x^2 + y^2 -R^2 $$然后根据这个函数来判别点与圆的位置关系：若 $F(x,y)=0$ ，则点在圆上若 $F(x,y)>0$ ，则点在圆外若 $F(x,y)<0$ ， 则点在圆内我们考虑如下圆弧对于点$(x_i,y_i)$，有下一可选择两点中点坐标 $(x_M,y_M)$有：当 $F(x_M,y_M)<0$ 时，取点 $(x_i+1,y_i)$当 $F(x_M,y_M)>0$ 时，取点 $(x_i+1,y_i-1)$当 $F(x_M,y_M)=0$ 时，约定取点 $(x_i+1,y_i)$构造判别式$$ d_i = F(x_M,y_M) = (x_i+1)^2 + (y_i-0.5)^2 - R^2 $$当 $d_i \le 0$的时候，下一中点坐标为$$ (x_i + 2, y_i-0.5) $$得到递推式$$ d_{i+1} = F(x_i+2,y_i-0.5)=d_i+2x_i+3 $$同理对于 $d_i>0$ 的时候有下一中点$$ (x_i+1,y_i-1) $$递推式如下：$$ d_{i+1} = F(x_i+1,y_i-1)=d_i+2(x_i-y_i)+5 $$我们还得求出第一次迭代时候的初始值 $d_0$$$ d_0 = F(1,R-0.5) = 1.25-R $$综上所述$$ d_{i+1} = \begin{cases} d_i+2x_i+3, &d_i\le 0 \\ d_i+2(x_i-y_i) + 5, &d_i>0 \end{cases} $$$$ d_0 = 1.25-R $$代码如下（当然还是伪代码）：/* * 中点Bresenham画圆 * circle_center为圆心 * R为圆半径 * color为绘制的颜色 */ void draw_circle_midb(const Point& circle_center, double R, Color color){ Point p; p.x = 0, p.y = R; double d = 1.25-R; while(p.x<=p.y){ draw_circle_point(p, circle_center, color); if(d<=0) d+=2*p.x+3; else{ d += 2*(p.x-p.y)+5; p.y--; } p.x++; } } 改进在计算中仍有小数运算，为了避免这种情况令 $d-0.25$ 代替原来的 $d$由于为int类型，$d\le -0.25$ 可认为是 $d\le 0$$$ d_{i+1} = \begin{cases} d_i+2x_i+3, &d_i\le 0 \\ d_i+2(x_i-y_i) + 5, &d_i>0 \end{cases} $$$$ d_0 = 1-R $$代码就变成了（其实改动不是很大）/* * 改进中点Bresenham画圆 * circle_center为圆心 * R为圆半径 * color为绘制的颜色 */ void draw_circle_midb(const Point& circle_center, double R, Color color){ Point p; p.x = 0, p.y = R; int d = 1-R; while(p.x<=p.y){ draw_circle_point(p, circle_center, color); if(d<=0) d+=2*p.x+3; else{ d += 2*(p.x-p.y)+5; p.y--; } p.x++; } }椭圆基础知识椭圆方程为$$ \frac{x^2}{a^2}+\frac{y^2}{b^2} = 1, (a>0,b>0,a\ne b) $$隐式方程为$$ b^2x^2 + a^2y^2 -a^2b^2 = 0 $$注意我们为了统一表示长轴在y上和长轴在x轴上的椭圆，并没有指明 $a>b$ 还是 $b>a$ 。有参数方程：$$ \begin{cases} x = a\cos\theta \\ y = b\cos\theta \end{cases} $$我们只考虑中心在原点的椭圆，如果不在原点先平移坐标系，扫描转换后再平移回去。对于中心在原点的椭圆有如下对称性质：关于x轴对称关于y轴对称故我们用四分法来绘制椭圆。也就是说我们只需有四分之一的弧线，即可绘制出整个椭圆。不妨约定一个函数/* 以原点为中心，通过四分之一弧线上的一点绘制其余四个对称点 * p为绘制的点 * center为椭圆中心 * color为绘制颜色 */ void draw_oval_point(const Point& p, const Point& center, Color color){ int x = p.x; int y = p.y; int ox = center.x, oy = center.y; setPixel(x+ox,y+oy,color); //绘制像素 setPixel(-x+ox,y+oy,color); setPixel(-x+ox,-y+oy,color); setPixel(x+ox,-y+oy,color); }直接计算法笛卡尔方程$$ x_{i+1} = x_i+1 \\ y_{i+1} = round(\sqrt{\frac{a^2b^2-b^2x_{i+1}^2}{a^2}}) $$代码如下/* * 利用笛卡尔方程直接绘制椭圆 * center为中心点 * a为标准方程x分数线下的参数（此处不一定为长轴长度一般，也可能是短轴长一半） * b为标准方程下y分数线下的参数（同理不一定为短轴长一半） * color 为颜色 */ void draw_oval_common(const Point& center, double a, double b, Color color){ double aa = a*a; double bb = b*b; double temp = aa*bb; Point p; for(p.x=0;p.x<=a;p.x++){ p.y = sqrt((temp-bb*p.x*p.x)/aa); draw_oval_point(p,center,color); } }参数方程离散化参数方程$$ \begin{cases} x=a\cos (2\pi \frac{i}{n}) \\ y = b\sin (2\pi \frac{i}{n}) \end{cases} $$代码如下：/* * 利用参数方程直接绘制椭圆 * center为中心点 * a为标准方程x分数线下的参数（此处不一定为长轴长度一般，也可能是短轴长一半） * b为标准方程下y分数线下的参数（同理不一定为短轴长一半） * color 为颜色 * n为分割参数，越大越精细但越慢 */ void draw_oval_param(const Point& center, double a, double b, Color color,int n){ int k; Point p; for(k=0;k<=n;k++){ p.x = a*cos(2*PI*k/n); p.y = b*sin(2*PI*k/n); draw_oval_point(p, center, color); } } DDA算法$$ x_{i+1} = a\cos{(\alpha + \theta)} =a\cos\alpha\cos\theta - a\sin\alpha\sin\theta = x_i\cos\theta-\frac{a}{b}y_i\sin\theta \\ y_{i+1} = b\cos{(\alpha + \theta)} = y_i\cos\theta+\frac{b}{a}x_i\sin\theta $$当 $\theta$ 足够小的时候可认为$$ \cos\theta \approx 1\\ \sin\theta \approx \theta $$即有$$ \begin{cases} x_{i+1} = x_i - \frac{a}{b} y_i\theta \\ y_{i+1} = y_i+\frac{b}{a} x_i\theta \end{cases} $$/* * DDA绘制椭圆 * center为中心点 * a为标准方程x分数线下的参数（此处不一定为长轴长度一般，也可能是短轴长一半） * b为标准方程下y分数线下的参数（同理不一定为短轴长一半） * color 为颜色 */ void draw_circle_dda(const Point& center, double a, double b, Color color){ Point p; p.x = a; p.y = 0; double d = 1/a; for(double t=0;t<=PI/2;t+=d){ draw_oval_point(p,center,color); double old_x = p.x; p.x -= p.y/b; p.y += old_x*b/(a*a); } }中点Bresenham算法算法原理从 $(0,b)$ 到 $(a,0)$ 逆时针地确定最佳逼近于第一象限的椭圆弧像素序列记切线斜率为 $k$ 。在上半部分， $-1\le k \le0$ ，最大位移方向为x；在下半部分， $k\le-1$ ，最大位移方向为y。分别对这两部分判断中点与椭圆的关系我们先来看第一象限的上半部分与圆同理，判别式为：$$ d_i = F(x_M,y_M) = b^2(x_i+1)+a^2(y_i-0.5)^2 - a^2b^2 $$若 $d_i\le 0$ ，取 $(x_i+1,y_i)$若 $d_i>0$， 取 $(x_i+1,y_i-1)$当 $d_i \le 0$ 的时候，下一中点为 $(x_i+2,y_i-0.5)$ ，有递推式：$$ d_{i+1} = F(x_i+2,y_i-0.5) = d_i + b^2(2x_i + 3) $$当 $d_i>0$ 的时候，下一中点为 $(x_i-2,y_i-1.5)$ ，有递推式：$$ d_{i+1} = d_i + b^2(2x_i+3)+a^2(-2y_i+2) $$初始值为$$ d_0 = b^2+a^2(-b+0.25) $$然后再来看看下半部分同样有判别式$$ d_i = b^2(x_i+0.5)^2+a^2(y_i-1)^2-a^2b^2 $$当$d_i\le 0$ 时，取点 $(x_i+1,y_i-1)$ 且有下一中点 $(x_i+1.5, y_i-2)$$$ d_{i+1} = d_i+b^2(2x_i+2)+a^2(-2y_i+3) $$当 $d_i>0$ 时，取 $(x_i,y_i-1)$ ，有下一中点 $(x_i+0.5,y_i-2)$$$ d_{i+1} = d_i + a^2(-2y_i+3) $$但初始值如何计算？既然无法直接得到，我们干脆利用上半部分最后一个点来计算下半部分初始值！$$ d_0 = b^2(x+0.5)^2+a^2(y-1)^2-a^2b^2 $$代码void draw_oval_midb(const Point& center, double a, double b, Color color){ Point p; double d1 = b*b+a*a*(-b+0.25); p.x = 0; p.y = b; draw_oval_point(p,center,color); while(b*b*(p.x+1)<a*a*(p.y-0.5)){ if(d1<0){ d1+=b*b*(2*p.x+3); p.x++; }else{ d1+=b*b*(2*p.x+3)+a*a*(-2*p.y+2); p.x++;p.y--; } draw_oval_point(p,center,color); } //下半部分 double d2 = sqrt(b*(p.x+0.5))+a*(p.y-1)-a*b; while(p.y>0){ if(d2<0){ d2 += b*b*(2*p.x+2)+a*a*(-2*p.y+3); p.x++;p.y--; }else{ d2+=a*a*(-2*p.y+3); p.y--; } draw_oval_point(p,center,color); } }结束终于写完了，累死我了。尤其是椭圆那里，程序怎么都不过，后来才发现公式推导推错了，又把推到过程重新敲了一遍。什么你问我写的明明是伪代码，为什么还能跑。很简单，因为本来我的代码就是仿照一些接口写的，有部分甚至直接粘贴就能运行，不能运行的改个函数名加个头文件就能跑。以下就是我在Qt中的部分测试代码，代码几乎一致，顶多加几个对象#include "mainwindow.h" #include "./ui_mainwindow.h" #include <qmath.h> #include <QPainter> const double PI = 4*atan(1.0); MainWindow::MainWindow(QWidget *parent) : QMainWindow(parent) , ui(new Ui::MainWindow) { ui->setupUi(this); } void MainWindow::draw_oval_point(const Point &p, const Point &center){ int x = p.x; int y = p.y; int ox = center.x; int oy = center.y; this->image->setPixel(x+ox,y+oy,qRgb(255,0,0)); this->image->setPixel(-x+ox,y+oy,qRgb(255,0,0)); this->image->setPixel(x+ox,-y+oy,qRgb(255,0,0)); this->image->setPixel(-x+ox,-y+oy,qRgb(255,0,0)); } void MainWindow::draw_oval_dda(const Point& center,double a,double b){ Point p; double d1 = b*b+a*a*(-b+0.25); p.x = 0; p.y = b; draw_oval_point(p,center); while(b*b*(p.x+1)<a*a*(p.y-0.5)){ if(d1<0){ d1+=b*b*(2*p.x+3); p.x++; }else{ d1+=b*b*(2*p.x+3)+a*a*(-2*p.y+2); p.x++;p.y--; } draw_oval_point(p,center); } //下半部分 double d2 = sqrt(b*(p.x+0.5))+a*(p.y-1)-a*b; while(p.y>0){ if(d2<0){ d2 += b*b*(2*p.x+2)+a*a*(-2*p.y+3); p.x++;p.y--; }else{ d2+=a*a*(-2*p.y+3); p.y--; } draw_oval_point(p,center); } } void MainWindow::paintEvent(QPaintEvent *event){ QPainter painter(this); this->image = new QImage(1000, 600, QImage::Format_RGB32); this->image->fill(Qt::white); Point center(400,200); this->draw_oval_dda(center,300,100); painter.drawImage(QPoint(0, 0), *image); } MainWindow::~MainWindow() { delete image; delete ui; } 如果你发现文章中有逻辑错误，请在下方留言，与圆有关的部分代码我没有测试。
数值微分DDA算法与Bresenham算法画直线前提最近在学计算机图形学画线算法，主要学了两种数值微分算法(Digital Differential Analyzer, DDA)中点Bresenham算法直线我们既然要画线就必然要研究一些直线的基本性质。在中学就学过直线的斜截式方程为（我们先不考虑斜率不存在的情况，对与这种情况我们非常容易处理，只需单独拿出来判断下就行）$$ y=kx+b $$其中k为直线斜率，b为直线截距。在计算机中，我们如果要画一条线段，只需要取直线上的部分的点即可。我们不妨考虑线段的起始点和终点，其中起点坐标为：$$ (x_0,y_0) $$终点坐标为：$$ (x_{end},y_{end}) $$我们很轻松能得到斜率和截距与两点的关系式：$$ k = \frac{y_{end}-y_{0}}{x_{end}-x_{0s}} $$$$ b = y_0 - kx_0 $$不妨分别给定x方向和y方向的增量$\Delta{x}$和$\Delta{y}$$$ \Delta{y}=k\Delta{x} $$变形可得$$ k= \frac{\Delta{y}}{\Delta{x}} $$在光栅算法中，我们规定x和y方向上变化步长（变化量的绝对值）大的那一方向为最大位移方向。（沿正还是负根据实际情况判断）DDA算法原理取无穷小量$\varepsilon >0$，由于直线的一阶导数是连续的，而且对于$\Delta{x}$和$\Delta{y}$是成比例的，因此可通过在当前位置$(x_i,y_i)$，分别加上$\varepsilon x$和$\varepsilon y$来求出下一点的坐标$(x_{i+1},y_{i+1})$，即有方程$$ \begin{cases} x_{i+1} = x_{i}+\varepsilon \Delta x \\ y_{i+1} = y_{i}+\varepsilon \Delta y \\ \end{cases} $$对于计算机来讲这些点都是离散的，也就是说我们不可能取真正意义上的无穷小量。对于屏幕来讲，最小的单位就是一个像素，也就是说对于x和y方向上我们只能移动正整数个像素或者不位移，而最小的正整数为1，我们要绘制下一个点的坐标必然要在某个方向上移动1。通俗一点来讲，我们取x、y其中一个方向移动一个像素，以x为例：我们在x轴上移动一个像素大小，那么根据直线方程$x_i+1$代入得到的y值很可能是个小数，我们需要通过向上取整或者向下取整来找一个与这个小数近似的整数来作为近似的y坐标。接下来就很简单了，我们仅需要计算一下恰好比这个小数大的整数和恰好比这个小数小的整数离该小数的距离即可（四舍五入，哪个整数离得近就选哪个整数）。那么问题来了，我们究竟是选择x方向移动1，还是y方向移动1？回到我们的实际问题上来：我们仅想要一些近似于实际直线上的点，越近似越好！如果x变化比y快，也就是说x为最大位移方向，当我们选择移动y一个单位，x变化就可能非常大，带来的误差就可能大。y为最大位移方向，选择移动x同理。那么结论就有了：为降低误差，选最大位移方向为移动方向，移动1个单位的问题就转变成了找最大位移方向问题。当x为最大位移方向的时候$$ max(|\Delta x|,|\Delta y|) = |\Delta x| $$也就是说$$ 0 \le |k|\le 1 \\ \varepsilon |\Delta x| = 1 $$我们可以得出$$ \varepsilon |\Delta y| =\frac{|\Delta y|}{|\Delta x|} $$那么代入方程组可得$$ \begin{cases} x_{i+1} = x_{i}+\varepsilon |\Delta x| = x_{i}+\frac{\Delta x}{|\Delta x|}=x_i \pm 1\\ y_{i+1} = y_{i} + \varepsilon |\Delta y| = y_{i}+\frac{\Delta y}{|\Delta x|} = y_{i} \pm k \end{cases} $$同理我们得出y为最大方向的时候方程组如下$$ \begin{cases} x_{i+1} = x_{i} \pm \frac{1}{k} \\ y_{i+1} = y_{i} \pm 1 \end{cases} $$对得到的$(x_{i+1},y_{i+1})$的坐标四舍五入即为我们要绘制的近似点的坐标。上述两个方程组正负号处取正号还是负号与$\Delta x$和$\Delta y$有关系，例如对于初始点x坐标比终点x坐标小，斜率在$(0,1]$上的直线，最大位移方向为x方向，从起点到终点绘制x是不断增大的，所以$x_i \pm 1$应该为$x_i + 1$；结合斜率可知y也是增大的，则$y_i \pm k$应该变为$y_i + k$。不妨将起点和终点坐标交换一下，此时$\Delta x$和$\Delta y$都为负，那么就变成了$x_i - 1$和$y_i - k$。代码C++风格的伪代码如下：/// @brief DDA画直线 /// @param start 起始点 /// @param end 终点 void dda_line(Point& start,Point& end) { int dx = end.x-start.x; int dy = end.y - start.y; int steps = abs(dx)>abs(dy)?abs(dx):abs(dy); //迭代次数 double x_incre = dx/(double)steps; //x每次变化量 double y_incre = dy/(double)steps; //y每次变化量 double x=start.x,y=start.y; //绘制像素,round为四舍五入函数 set_pixel(round(x),round(y),rgb(0,0,255)); for(int i=0;i<steps;i++){ x+=x_incre; y+=y_incre; set_pixel(round(x),round(y),rgb(0,0,255)); } }特点DDA算法简单直观，容易实现，但是计算机并不太擅长计算浮点数，每一步都得四舍五入取整消耗资源相对较大，不利于硬件实现。中点Bresenham算法原理对于直线$y=kx+b$移项可得隐式方程$$ y-kx-b=0 $$不妨令二元函数$F(x,y)$有如下形式$$ F(x,y)=y-kx-b $$对于点$(x_n,y_n)$，当$F(x_n,y_n)=0$时，该点在直线上；当$F(x_n,y_n)<0$时，点在直线下方；当$F(x_n,y_n)>0$时，点在直线上方。我们分斜率区间为$(-\infty ,-1],(-1,0],(0,1](1,+\infty )$四种情况讨论事实上当斜率$k =\pm 1$的时候x与y方向移动速度是一样的，我们选择x为最大移动方向和选择y为最大移动方向是一样的，因此上方给出的这四个区间有些可以在这个地方取等号。另外，当$x=0$时，无论0被划分为$(-1,0]$还是$[0,1]$其最大位移方向都是x。在接下来的分类讨论中可能会出现如$k\le-1$和$-1\le k \le0$两种情况都能取到-1之类的描述，对于0和1同理，这并不妨碍结论的严谨性，我们知道这一点即可。在只考虑x沿正方向增长的条件下（要求该方向上起始点横坐标比终点横坐标小，如果大的话我们只调换下这两个点坐标让小的点作为起始点，大的点作为终点，再按照下方流程绘制就行了），我们来看$0\le k\le1$的情况：此时x为最大位移方向，不妨设$(x_i,y_i)$为绘制过程中的某一点，则下一点有$(x_{i}+1,y_i)$和$(x_i+1,y_i + 1)$两种选择。我们根据中点公式$$ x_{mid} = \frac{x_{a}+x_{b}}{2} \\ y_{mid} = \frac{y_a+y_b}{2} $$得到这两个可选择点的中点坐标为$$ (x_i+1,y_i+0.5) $$将其带回二元函数$F(x,y)$得到判别式$$ d_i = F(x_i+1,y_i+0.5)=y_i + 0.5 - k(x_i + 1)-b=y_i-kx_i-b-0.5-k $$当$d_i<0$的时候说明中点在直线下方，$(x_i + 1,y_i+1)$离直线较近，我们应该选它；当$d_i>0$当时候，说明中点在直线上方，此时我们应该选择$(x_i+1,y_i)$然后我们再来看看绘制下一个点此时我们如果是从$(x_i+1,y_i+1)$出发，则下一步可选择的两个点分别为$(x_i+2,y_i+1)$和$(x_i+2,y_i+2)$，然后求出中点$$ (x_i+2,y_i+1.5) $$代入公式可得$$ d_{i+1} = F(x_i+2,y_i+1.5)=y_i-kx_i-b-2k-1.5=(y_i-kx_i-b+0.5-k)+1-k=d_i+1-k $$即$$ d_{i+1}=d_i+1-k $$如果从$(x_i+1,y_i)$出发，则可选择两个点分别为$(x_i+2,y_i+1)$和$(x_i+2,y_i)$中点坐标为$$ (x_i+2,y_i+0.5) $$那么判别式为$$ d_{i+1} = F(x_i+2,y_i+0.5)=y_i+0.5-k(x_i+2)-b=d_i-k $$由此可以得出$d_i$与$d_{i+1}$之间的关系递推式$$ d_{i+1}=\begin{cases} d_i-k,& d_i\ge 0 \\ d_i-k+1,&d_i < 0 \end{cases} $$对于递推公式的初始值$d_0$，将$(x_0+1,y_0+0.5)$代入可得$$ d_0 = F(x_0+1,y_0+0.5)=(y_0-kx_0-b)+0.5-k=0+0.5-k=0.5-k $$同理，我们可知$k\ge 1$的情况：最大位移方向为y方向$$ d_i = F(x_i+0.5,y_i+1)=y_i+1-k(x_i+0.5)-b \\ d_0 = F(x_0+0.5,y_0+1) = (y_i-kx_0-b)-0.5k+1=1-0.5k $$有递推式$$ d_{i+1} = \begin {cases} d_i-k+1,& d_i\ge 0 \\ d_i+1, & d_i<0 \end{cases} $$再来看$-1\le k\le 0$的情况：最大位移方向为x方向，但y为负增长$$ d_i = F(x_i+1,y_i-0.5) = y_i-0.5-k(x_i+1)-b \\ d_0 = F(x_0+1,y_0-0.5) = (y_0-kx_0-b)-k-0.5 $$递推式为$$ d_{i+1} = \begin{cases} d_i -k, & d_i<0 \\ d_i -k -1, & d_i\ge0 \end{cases} $$最后看看$k\le -1$的情况此时y为最大位移方向，y为负增长$$ d_i = F(x_i+0.5,y_i-1) = y_i-1-k(x_i+0.5)-b \\ d_0 = F(x_0+0.5,y_0-1) = (y_0-kx_0-b)-0.5k-1 = -1-0.5k $$依然有递推式$$ d_{i+1} = \begin{cases} d_i - k - 1, & d_i<0 \\ d_i - 1, & d_i\ge 0 \end{cases} $$优化在DDA中我们就提到过计算机不喜欢计算浮点数和乘除法，这会降低执行效率不利于硬件实现，但是我们在Bresenham算法中看到依然会有这些情况存在，所以需要对算法进行一定的优化。替换原有的d利用$2\Delta xd$替换原有的$d$当$0\le k\le1$时$$ d_0 = \Delta x - 2\Delta y \\ d_{i+1} = \begin{cases} d_i+2\Delta x - 2\Delta y, & d_i < 0 \\ d_i-2\Delta y , & d_i\ge 0 \end{cases} $$当$k\ge 1$时$$ d_0 = 2\Delta x -\Delta y \\ d_{i+1} = \begin{cases} d_i + 2\Delta x, & d_i<0 \\ d_i + 2\Delta x - 2\Delta y, & d_i\ge 0 \end{cases} $$当$ -1\le k \le 0$时$$ d_0 = -2\Delta x - 2\Delta y \\ d_{i+1} = \begin{cases} d_i - 2\Delta y, & d_i<0 \\ d_i - 2\Delta y - 2\Delta x, & d_i \ge 0 \end{cases} $$当$k\le -1$时$$ d_0 = -2\Delta x - \Delta y \\ d_{i+1} = \begin{cases} d_i - 2\Delta x - 2\Delta y,& d_i<0 \\ d_i - 2\Delta x, &d_i \ge 0 \end{cases} $$这里注意一点，我们递推式里的$d_i$并没有写成$2\Delta x d$的形式，这是因为我们在$d_0$的时候已经乘进去，这里的$d_i$已经是变化后的了代码C++风格伪代码/// @brief Breseham中点画线算法 /// @param start 起始点 /// @param end 终止点 void mid_bresenham_line(Point &start, Point &end) { // 如果起点横坐标大于终点，则交换两点（保证x正增长） if (start.x > end.x) { int temp = start.x; start.x = end.x; end.x = temp; temp = start.y; start.y = end.y; end.y = temp; } int x = start.x; int y = start.y; // 绘制图形 set_pixel(x, y, rgb(0, 0, 255)); int dx = end.x - start.x; int dy = end.y = start.y; int d; // 判别式 if (dy <= dx && dy >= 0) { x += 1; d = dx - 2 * dy; if (d < 0) { y += 1; set_pixel(x, y, rgb(0, 0, 255)); d = d + 2 * dx - 2 * dy; } else { set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dy; } } else if (dy >= dx) { y += 1; d = 2 * dx - dy; if (d < 0) { x += 1; set_pixel(x, y, rgb(0, 0, 255)); d = d + 2 * dx; } else { set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dy - 2 * dx; } } else if (dy >= -dx && dy < 0) { x += 1; d = -2 * dx - 2 * dy; if (d < 0) { set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dy; } else { y -= 1; set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dy - 2 * dx; } } else { y -= 1; d = -2 * dx - dy; if (d < 0) { x += 1; set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dx - 2 * dy; } else { set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dx; } } while (x < end.x) { if (dy <= dx && dy >= 0) { x += 1; if (d < 0) { y += 1; set_pixel(x, y, rgb(0, 0, 255)); d = d + 2 * dx - 2 * dy; } else { set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dy; } } else if (dy >= dx) { y += 1; if (d < 0) { x += 1; set_pixel(x, y, rgb(0, 0, 255)); d = d + 2 * dx; } else { set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dy - 2 * dx; } } else if (dy >= -dx && dy < 0) { x += 1; if (d < 0) { set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dy; } else { y -= 1; set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dy - 2 * dx; } } else { y -= 1; if (d < 0) { x += 1; set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dx - 2 * dy; } else { set_pixel(x, y, rgb(0, 0, 255)); d = d - 2 * dx; } } } }改进的Bresenham算法原理在中点Bresenham中，我们是通过判断中点与直线位置关系来进行选择点位的。现在我们换一种思路，仅仅使用增量整数计算。当斜率$0\le k\le 1$时算法步骤大致描述如下：输入直线的两端点计算初始值 $\Delta x$ 、 $\Delta y$ 、 $d=0$ 、 $x=x_0$ 、 $y = y_0$绘制点 $(x,y)$$d$ 更新为 $d+k$，判断d的符号。若 $d>0.5$ ，则 $(x,y)$ 更新为 $(x+1,y+1)$，同时将 $d$ 更新为 $d-1$ ;否则 $(x,y)$ 更新为 $(x+1,y)$当直线没有画完时，重复3、4两步我们每次都要判断 $d>0.5$ ，为了方便计算，不妨令 $e = d-0.5$对于该种情况下的关于 $x$ 和 $y$ 的递推公式如下$$ \begin{cases} x_{i+1} = x_i+1 \\ y_{i+1} = \begin{cases} y_i+1,& e>0,\\ y_i,& e \le 0 \end{cases} \end{cases} $$对于e有如下操作：$e_0 = -0.5$每走一步在原有基础上加 $k$ ，即 $e_{i+1} = e_i+k$若 $e_{i+1}>0$ ， $e_{i+1} = e_{i+1}-1$类似于中点Bresenham算法，我们在这里也要优化掉小数和除法运算：即用$2\Delta x e$来替换$e$所以对于e就变成了：$e_0 = -\Delta x$每走一步 $e_{i+1}=e_i + 2\Delta y$若 $e_{i+1} > 0,e_{i+1}=e_{i+1}-2\Delta x$此时算法步骤为：输入直线两端点计算初始值 $\Delta x$ 、 $\Delta y$ 、 $e = -\Delta x$ 、 $x=x_0$和$y = y_0$绘制点 $(x,y)$$e$更新至$e+2\Delta y$ ，判断e符号。若 $e>0$ ，则 $(x,y)$ 更新为 $(x+1,y+1)$ ，同时将 $e$ 更新为 $e-2\Delta x$ ，否则为 $(x+1,y)$重复步骤3、4，直至画完直线注意上述只给出了斜率k大于等于0小于等于1点步骤，其他范围同理可得代码为了简洁，值给出 $0\le k \le 1$ 的流程（CPP风格伪代码）/// @brief Bresenham画直线（此函数仅为斜率0<=k<=1） /// @param start 起始点 /// @param end 终点 void bresenham_line(Point& start, Point& end){ int dx = fabs(end.x-start.x),dy=(end.y-start.y); int double_dy = 2*dy; int double_dy_mins_dx = 2*(dy-dx); int double_dx = 2*dx; int e = -dx; int x,y; //比较大小，判断哪个是起点（确保起点横坐标小于终点） if(start.x>end.x){ x = end.x; y = end.y; end.x = start.x; }else{ x = start.x; y = start.y; } set_pixel(x,y,rgb(0,0,255)); while(x<end.x){ ++x; e+=double_dy; if(e>0){ x++,y++; e -= double_dx; } else{ ++x; } set_pixel(x,y,rgb(0,0,255)); } }结束代码虽然是伪代码，主要是因为我当前环境下没有set_pixel这个函数的实现，事实上这个API应该是由系统提供的，在Windows下好像是有（不确定），但当前环境为MacOS所以就没尝试，直接以伪代码的形式给出了。（如果Win下有这个函数的话应该导入相应的库然后以正确的函数名代替set_pixel就能运行）另外在写改进的Bresenham算法那一部分的介绍时，发现再写斜率为其它三种情况的推导和代码有些太过于麻烦了，所以仅给出了一种情况的，其它情况同理可得。另外，对于终点Bresenham算法的代码部分用了多层if嵌套，这实际上是个不好的习惯，应该将其逻辑拆分成好几个函数逐层细化，不过此处主要介绍原理，就没考虑这方面问题。当时写代码时为了不让IDE弹出烦人的报错，我自己定义了点的类型Pointtypedef struct Point{ int x; int y; } Point;最后，DDA和Bresenham不仅仅适用于画直线，曲线也可以，之后会写一篇画椭圆或圆的文章。