扫描转换椭圆与圆的绘制

扫描转换椭圆与圆的绘制

前提

我们先来设计一种数据类型

class Point{
public:
  double x,y;
  Point(){this->x = 0; this->y = 0;}
  Point(double x, double y){this->x=x;this->y=y;}
};
1
2
3
4
5
6

注意此处的类型为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$$
同样的，对于圆心在原点上的圆我们有参数方程：
{ x = R cos ⁡ θ y = R sin ⁡ θ {x=Rcosθy=Rsinθ

{x=Rcosθy=Rsinθ​

对于圆心不在原点上的圆有参数方程：
{ x = a + R cos ⁡ θ y = b + R sin ⁡ θ {x=a+Rcosθy=b+Rsinθ

{x=a+Rcosθy=b+Rsinθ​
将$\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);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

此外在接下来的计算中可能会用到三角函数的和差公式，请牢记这一点：
$$\begin{gathered} \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 \end{gathered}$$

直接计算法

直接生成有两种方式，一种是根据直角坐标系下笛卡尔方程来计算，另一种根据参数方程。

笛卡尔方程直接计算

我们扫描这样圆弧上的点。

根据方程 $x^2+y^2=R$ 可以得到：
{ x i + 1 = x i + 1 , x ∈ [ 0 , R 2 ] y i + 1 = r o u n d ( R 2 − x i + 1 2 ) {xi+1=xi+1,x∈[0,R√2]yi+1=round(√R2−x2i+1)

⎩ ⎨ ⎧​xi+1​=xi​+1,x∈[0,2 ​R​]yi+1​=round(R2−xi+12​ ​)​​
故有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);
  }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

参数方程直接计算

我们是根据离散化的参数方程来计算的

{ x = R cos ⁡ ( 2 π i n ) y = R sin ⁡ ( 2 π i n ) {x=Rcos(2πin)y=Rsin(2πin)

{x=Rcos(2πni​)y=Rsin(2πni​)​

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);
  }
  
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

DDA算法

数值微分算法(Digital Differential Analyzer, DDA)在前面绘制直线的文章中已经提到过了，现在我们来尝试使用它绘制圆。

我们考虑如下图

{ x i = R cos ⁡ α y i = R sin ⁡ α {xi=Rcosαyi=Rsinα

{xi​=Rcosαyi​=Rsinα​
我们让当前角度 $\alpha$ 沿着逆时针方向加上几个极小的角度量 $\theta$

就得到了下一个点
$$(x_{i+1},y_{i+1})=(R\cos (\alpha +\theta ),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$ 足够小的时候有
$$\begin{gathered} \cos \theta \approx 1 \\ \sin \theta \approx \theta \end{gathered}$$
故递推式可以改写为

{ x i + 1 = x i − y i θ y i + 1 = y i + x i θ {xi+1=xi−yiθyi+1=yi+xiθ

{xi+1​=xi​−yi​θyi+1​=yi​+xi​θ​

知道这一点我们可以写出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;
  }
  
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

中点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 = { d i + 2 x i + 3 , d i ≤ 0 d i + 2 ( x i − y i ) + 5 , d i > 0 d_{i+1} = {di+2xi+3,di≤0di+2(xi−yi)+5,di>0

di+1​={di​+2xi​+3,di​+2(xi​−yi​)+5,​di​≤0di​>0​

$$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++;
  }
  
}
  
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

改进

在计算中仍有小数运算，为了避免这种情况令 $d-0.25$ 代替原来的 $d$

由于为int类型，$d\le -0.25$ 可认为是 $d\le 0$

d i + 1 = { d i + 2 x i + 3 , d i ≤ 0 d i + 2 ( x i − y i ) + 5 , d i > 0 d_{i+1} = {di+2xi+3,di≤0di+2(xi−yi)+5,di>0

di+1​={di​+2xi​+3,di​+2(xi​−yi​)+5,​di​≤0di​>0​

$$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++;
  }
  
}
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

椭圆

基础知识

椭圆方程为
$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1,(a>0,b>0,a\ne b)$$

隐式方程为
$$b^2{x}^2+a^2{y}^2-a^2{b}^2=0$$
注意我们为了统一表示长轴在y上和长轴在x轴上的椭圆，并没有指明 $a>b$ 还是 $b>a$ 。

有参数方程：
{ x = a cos ⁡ θ y = b cos ⁡ θ {x=acosθy=bcosθ

{x=acosθy=bcosθ​
我们只考虑中心在原点的椭圆，如果不在原点先平移坐标系，扫描转换后再平移回去。

对于中心在原点的椭圆有如下对称性质：

• 关于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);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

直接计算法

笛卡尔方程

$$\begin{gathered} x_{i+1}=x_i+1 \\ {y}_{i+1}=round(\sqrt{\frac{a^2{b}^2-b^2{x}_{i+1}^2}{a^2}}) \end{gathered}$$

代码如下

/*
 * 利用笛卡尔方程直接绘制椭圆
 * 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);
  }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

参数方程

离散化参数方程
{ x = a cos ⁡ ( 2 π i n ) y = b sin ⁡ ( 2 π i n ) {x=acos(2πin)y=bsin(2πin)

{x=acos(2πni​)y=bsin(2πni​)​
代码如下：

/*
 * 利用参数方程直接绘制椭圆
 * 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);
  }
  
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

DDA算法

$$\begin{gathered} 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 \end{gathered}$$

当 $\theta$ 足够小的时候可认为
$$\begin{gathered} \cos \theta \approx 1 \\ \sin \theta \approx \theta \end{gathered}$$
即有
{ x i + 1 = x i − a b y i θ y i + 1 = y i + b a x i θ {xi+1=xi−abyiθyi+1=yi+baxiθ

{xi+1​=xi​−ba​yi​θyi+1​=yi​+ab​xi​θ​

/*
 * 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);
  }
  
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

中点Bresenham算法

算法原理

• 从 $(0,b)$ 到 $(a,0)$ 逆时针地确定最佳逼近于第一象限的椭圆弧像素序列
• 记切线斜率为 $k$ 。在上半部分， $-1\le k\le 0$ ，最大位移方向为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^2{b}^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^2{b}^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^2{b}^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);
  }
  
  
}
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
31
32
33
34
35
36

结束

终于写完了，累死我了。尤其是椭圆那里，程序怎么都不过，后来才发现公式推导推错了，又把推到过程重新敲了一遍。

什么你问我写的明明是伪代码，为什么还能跑。很简单，因为本来我的代码就是仿照一些接口写的，有部分甚至直接粘贴就能运行，不能运行的改个函数名加个头文件就能跑。

以下就是我在Qt中的部分测试代码，代码几乎一致，顶多加几个对象

#include "mainwindow.h"
#include "./ui_mainwindow.h"
#include 
#include 

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;
}

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

如果你发现文章中有逻辑错误，请在下方留言，与圆有关的部分代码我没有测试。