基础图形绘制-椭圆

中点画椭圆

基础

首先，我们考虑中点在 (0,0)(0,0) 的椭圆。

根据椭圆的对称性，我们只需要画出第一象限的 14\frac{1}{4} 的圆弧，即可通过对称画出整个椭圆。

image-20211013144156162

中点在 (0,0)(0,0) ，焦点在坐标轴上的椭圆的标准方程为：x2a2+y2b2=1\frac{x^2}{a^2}+\frac{y^2}{b^2}=1，隐函数为：F(x,y)=b2x2+a2y2−a2b2=0F(x,y)=b^2x^2+a^2y^2-a^2b^2=0，代入各点，有：

  1. 当 F(x,y)=b2x2+a2y2−a2b2=0F(x,y)=b^2x^2+a^2y^2-a^2b^2=0 时，点在椭圆上；
  2. 当 F(x,y)=b2x2+a2y2−a2b2>0F(x,y)=b^2x^2+a^2y^2-a^2b^2 \gt 0 时，点在椭圆外；
  3. 当 F(x,y)=b2x2+a2y2−a2b2<0F(x,y)=b^2x^2+a^2y^2-a^2b^2 \lt 0 时，点在椭圆内。

考虑画右上方的 14\frac{1}{4} 圆弧，根据法向量将椭圆弧分成A、B两部分：

image-20211013151823338

对 F(x,y)=b2x2+a2y2−a2b2=0F(x,y)=b^2x^2+a^2y^2-a^2b^2=0 分别对 xx、yy 求偏导得 (x,y)(x,y) 处的法向量方向为 (2b2x,2a2y)(2b^2x,2a^2y)，有：

  1. 当 2b2x<2a2y2b^2x\lt2a^2y 时，xx 方向法向量分量较小，水平方向较平缓，列入部分 AA；
  2. 当 2b2x>2a2y2b^2x\gt2a^2y 时，yy 方向法向量分量较小，垂直方向较平缓，列入部分 BB。

部分 A

image-20211013153438423

假设当前点在 (x,y)(x,y) ，考虑下一个点选择 (x+1,y)(x+1,y) 还是 (x+1,y−1)(x+1,y-1) 。

下一个点显然是选择 离椭圆弧更近的点 ，我们可以通过判断 (x+1,y−0.5)(x+1,y-0.5) 在 椭圆内 还是 椭圆外，来判断 (x+1,y)(x+1,y) 、 (x+1,y−1)(x+1,y-1) 哪个离圆弧更近（需要注意的是：如果 (x+1,y−0.5)(x+1,y-0.5) 恰好在椭圆上，算法选择 (x+1,y−1)(x+1,y-1) 来作为下一个点）。

也就是说，对于 (x,y)(x,y) ，我们判断

d=F(x+1,y−0.5)=b2(x+1)2+a2(y−0.5)2−a2b2=b2x2+a2y2+2b2x−a2y+0.25a2+b2−a2b2\begin{aligned} d &= F(x+1,y-0.5)\\ &= b^2(x+1)^2+a^2(y-0.5)^2-a^2b^2\\ &= b^2x^2+a^2y^2+2b^2x-a^2y+0.25a^2+b^2-a^2b^2 \end{aligned}

的值，来判断下一个点的选择情况。

  1. 若 d<0d \lt 0，选择 (x+1,y)(x+1,y) 作为下一个点。
    此时，判断 (x+1,y)(x+1,y) 的下一个点的判断值为：

    d′=F((x+1)+1,(y)−0.5)=F(x+2,y−0.5)=b2(x+2)2+a2(y−0.5)2−a2b2=b2x2+a2y2+4b2x−a2y+0.25a2+4b2−a2b2=(b2x2+a2y2+2b2x−a2y+0.25a2+b2−a2b2)+2b2x+3b2=d+2b2x+3b2\begin{aligned} d' &= F((x+1)+1,(y)-0.5)\\ &= F(x+2,y-0.5)\\ &= b^2(x+2)^2+a^2(y-0.5)^2-a^2b^2\\ &= b^2x^2+a^2y^2+4b^2x-a^2y+0.25a^2+4b^2-a^2b^2\\ &= (b^2x^2+a^2y^2+2b^2x-a^2y+0.25a^2+b^2-a^2b^2)+2b^2x+3b^2\\ &= d+2b^2x+3b^2 \end{aligned}

  2. 若 d≥0d \ge 0，选择 (x+1,y−1)(x+1,y-1) 作为下一个点。
    此时，判断 (x+1,y−1)(x+1,y-1) 的下一个点的判断值为：

    d′=F((x+1)+1,(y−1)−0.5)=F(x+2,y−1.5)=b2(x+2)2+a2(y−1.5)2−a2b2=b2x2+a2y2+4b2x−3a2y+2.25a2+4b2−a2b2=(b2x2+a2y2+2b2x−a2y+0.25a2+b2−a2b2)+2b2x−2a2y+2a2+3b2=d+2b2x−2a2y+2a2+3b2\begin{aligned} d' &= F((x+1)+1,(y-1)-0.5)\\ &= F(x+2,y-1.5)\\ &= b^2(x+2)^2+a^2(y-1.5)^2-a^2b^2\\ &= b^2x^2+a^2y^2+4b^2x-3a^2y+2.25a^2+4b^2-a^2b^2\\ &= (b^2x^2+a^2y^2+2b^2x-a^2y+0.25a^2+b^2-a^2b^2)+2b^2x-2a^2y+2a^2+3b^2\\ &= d+2b^2x-2a^2y+2a^2+3b^2 \end{aligned}

这样我们就可以递推 dd 的值，不断选择点并更新判断值。

同时考虑初始点 (0,b)(0,b) ，显然有：

d0=F(0+1,b−0.5)=b2+a2(b−0.5)2−a2b2=b2+0.25a2−a2b\begin{aligned} d_0 &= F(0+1,b-0.5)\\ &= b^2+a^2(b-0.5)^2-a^2b^2\\ &= b^2+0.25a^2-a^2b \end{aligned}

这样，我们就可以画出A部分的椭圆。

部分 B

image-20211013154126540

假设当前点在 (x,y)(x,y) ，考虑下一个点选择 (x,y−1)(x,y-1) 还是 (x+1,y−1)(x+1,y-1) 。

下一个点显然是选择 离椭圆弧更近的点 ，我们可以通过判断 (x+0.5,y−1)(x+0.5,y-1) 在 椭圆外 还是 椭圆内，来判断 (x,y−1)(x,y-1) 、 (x+1,y−1)(x+1,y-1) 哪个离圆弧更近（需要注意的是：如果 (x+0.5,y−1)(x+0.5,y-1) 恰好在椭圆上，算法选择 (x,y−1)(x,y-1) 来作为下一个点）。

也就是说，对于 (x,y)(x,y) ，我们判断

d=F(x+0.5,y−1)=b2(x+0.5)2+a2(y−1)2−a2b2=b2x2+a2y2+b2x−2a2y+a2+0.25b2−a2b2\begin{aligned} d &= F(x+0.5,y-1)\\ &= b^2(x+0.5)^2+a^2(y-1)^2-a^2b^2\\ &= b^2x^2+a^2y^2+b^2x-2a^2y+a^2+0.25b^2-a^2b^2 \end{aligned}

的值，来判断下一个点的选择情况。

  1. 若 d≥0d \ge 0，选择 (x,y−1)(x,y-1) 作为下一个点。
    此时，判断 (x,y−1)(x,y-1) 的下一个点的判断值为：

    d′=F((x)+0.5,(y−1)−1)=F(x+0.5,y−2)=b2(x+0.5)2+a2(y−2)2−a2b2=b2x2+a2y2+b2x−4a2y+4a2+0.25b2−a2b2=(b2x2+a2y2+b2x−2a2y+a2+0.25b2−a2b2)−2a2y+3a2=d−2a2y+3b2\begin{aligned} d' &= F((x)+0.5,(y-1)-1)\\ &= F(x+0.5,y-2)\\ &= b^2(x+0.5)^2+a^2(y-2)^2-a^2b^2\\ &= b^2x^2+a^2y^2+b^2x-4a^2y+4a^2+0.25b^2-a^2b^2\\ &= (b^2x^2+a^2y^2+b^2x-2a^2y+a^2+0.25b^2-a^2b^2)-2a^2y+3a^2\\ &= d-2a^2y+3b^2 \end{aligned}

  2. 若 d<0d \lt 0，选择 (x+1,y−1)(x+1,y-1) 作为下一个点。
    此时，判断 (x+1,y−1)(x+1,y-1) 的下一个点的判断值为：

    d′=F((x+1)+0.5,(y−1)−1)=F(x+1.5,y−2)=b2(x+1.5)2+a2(y−2)2−a2b2=b2x2+a2y2+3b2x−4a2y+4a2+2.25b2−a2b2=(b2x2+a2y2+b2x−2a2y+a2+0.25b2−a2b2)+2b2x−2a2y+3a2+2b2=d+2b2x−2a2y+3a2+2b2\begin{aligned} d' &= F((x+1)+0.5,(y-1)-1)\\ &= F(x+1.5,y-2)\\ &= b^2(x+1.5)^2+a^2(y-2)^2-a^2b^2\\ &= b^2x^2+a^2y^2+3b^2x-4a^2y+4a^2+2.25b^2-a^2b^2\\ &= (b^2x^2+a^2y^2+b^2x-2a^2y+a^2+0.25b^2-a^2b^2)+2b^2x-2a^2y+3a^2+2b^2\\ &= d+2b^2x-2a^2y+3a^2+2b^2 \end{aligned}

这样我们就可以递推 dd 的值，不断选择点并更新判断值。

同时考虑初始点，也就是 部分A 的结束点，设为 (x,y)(x,y) ，显然有：

d0=F(x+1,y−0.5)=b2(x+1)2+a2(y−0.5)2−a2b2=b2x2+a2y2+2b2x−a2y+0.25a2+b2−a2b2\begin{aligned} d_0 &= F(x+1,y-0.5)\\ &= b^2(x+1)^2+a^2(y-0.5)^2-a^2b^2\\ &= b^2x^2+a^2y^2+2b^2x-a^2y+0.25a^2+b^2-a^2b^2 \end{aligned}

这样，我们就可以画出基础的椭圆。

代码

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
static int dx[] = { 1, 1, -1, -1 };
static int dy[] = { 1, -1, 1, -1 };

static void Normalize(float& x, float& y)
{
	x = (x - (SCR_WIDTH / 2)) / (SCR_WIDTH / 2);
	y = (y - (SCR_HEIGHT / 2)) / (SCR_HEIGHT / 2);
}

static void MiddlePoint(int x0, int y0, int a, int b)
{
	int sqrA = a * a, sqrB = b * b;
	int x = 0, y = b;
	double d = sqrB + 0.25 * sqrA - sqrA * b;

	while (2 * sqrB * x < 2 * sqrA * y)
	{
		if (d < 0) d = d + 2 * sqrB * x + 3 * sqrB;
		else d = d + 2 * sqrB * x - 2 * sqrA * y + 2 * sqrA + 3 * sqrB, y--;
		x++;

		for (int i = 0; i < 4; i++)
		{
			vertices[count * 3] = x * dx[i] + x0, vertices[count * 3 + 1] = y * dy[i] + y0;
			Normalize(vertices[count * 3], vertices[count * 3 + 1]), count++;
		}
	}

	d = sqrB * (x + 0.5) * (x + 0.5) + sqrA * (y - 1) * (y - 1) - sqrA * sqrB;
	while (y >= 0)
	{
		if (d >= 0) d = d - 2 * sqrA * y + 3 * sqrB;
		else d = d + 2 * sqrB * x - 2 * sqrA * y + 3 * sqrA + 2 * sqrB, x++;

		y--;

		for (int i = 0; i < 4; i++)
		{
			vertices[count * 3] = x * dx[i] + x0, vertices[count * 3 + 1] = y * dy[i] + y0;
			Normalize(vertices[count * 3], vertices[count * 3 + 1]), count++;
		}
	}
}

Bresenham

Bresenham 算法是对 中点画椭圆 算法的改进。

我们发现 中点画椭圆 算法中 dd 的值涉及到过多的浮点数运算，Bresenham 将浮点数运算给取消。

观察 中点画椭圆 算法，发现最小的浮点数为 0.250.25，于是我们将所有运算扩大四倍，就可以达到消除浮点运算的效果。

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
static int dx[] = { 1, 1, -1, -1 };
static int dy[] = { 1, -1, 1, -1 };

static void Normalize(float& x, float& y)
{
	x = (x - (SCR_WIDTH / 2)) / (SCR_WIDTH / 2);
	y = (y - (SCR_HEIGHT / 2)) / (SCR_HEIGHT / 2);
}

static void Bresenham(int x0, int y0, int a, int b)
{
	int sqrA = a * a, sqrB = b * b;
	int x = 0, y = b, d = 4 * sqrB + sqrA - 4 * sqrA * b;

	while (sqrB * x < sqrA * y)
	{
		if (d < 0) d = d + 4 * (2 * sqrB * x + 3 * sqrB);
		else d = d + 4 * (2 * sqrB * x - 2 * sqrA * y + 2 * sqrA + 3 * sqrB), y--;
		x++;

		for (int i = 0; i < 4; i++)
		{
			vertices[count * 3] = x * dx[i] + x0, vertices[count * 3 + 1] = y * dy[i] + y0;
			Normalize(vertices[count * 3], vertices[count * 3 + 1]), count++;
		}
	}

	d = 4 * sqrB * x * x + 4 * sqrA * y * y + 4 * 2 * sqrB * x - sqrA * y + sqrA + 4 * sqrB - 4 * sqrA * sqrB;
	while (y >= 0)
	{
		if (d >= 0) d = d - 4 * (2 * sqrA * y + 3 * sqrB);
		else d = d + 4 * (2 * sqrB * x - 2 * sqrA * y + 3 * sqrA + 2 * sqrB), x++;

		y--;

		for (int i = 0; i < 4; i++)
		{
			vertices[count * 3] = x * dx[i] + x0, vertices[count * 3 + 1] = y * dy[i] + y0;
			Normalize(vertices[count * 3], vertices[count * 3 + 1]), count++;
		}
	}
}

OpenGL 完整代码

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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
#define GLEW_STATIC
#include <GL/glew.h>
#include <GL/GL.h>
#include <GLFW/glfw3.h>

#include <iostream>

#pragma region Setting

static GLFWwindow* window;
const unsigned int SCR_WIDTH = 800;
const unsigned int SCR_HEIGHT = 600;
const unsigned int MAX_COUNT = 800 * 600;

static void InitializeWindow()
{
	glfwInit();
	glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
	glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
	glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);

	window = glfwCreateWindow(SCR_WIDTH, SCR_HEIGHT, "Test", NULL, NULL);
	glfwMakeContextCurrent(window);
	glfwSetFramebufferSizeCallback(window, [](GLFWwindow* window, int width, int height) { glViewport(0, 0, width, height); });


	glewExperimental = GL_TRUE;
	glewInit();

}

static void ProcessInput(GLFWwindow* window)
{
	if (glfwGetKey(window, GLFW_KEY_ESCAPE) == GLFW_PRESS)
		glfwSetWindowShouldClose(window, true);
}
#pragma endregion


#pragma region InitializeVertex

static float vertices[MAX_COUNT * 3];
static unsigned int VAO, VBO;
static unsigned int count = 0;

static int dx[] = { 1, 1, -1, -1 };
static int dy[] = { 1, -1, 1, -1 };

static void Normalize(float& x, float& y)
{
	x = (x - (SCR_WIDTH / 2)) / (SCR_WIDTH / 2);
	y = (y - (SCR_HEIGHT / 2)) / (SCR_HEIGHT / 2);
}

static void MiddlePoint(int x0, int y0, int a, int b)
{
	int sqrA = a * a, sqrB = b * b;
	int x = 0, y = b;
	double d = sqrB + 0.25 * sqrA - sqrA * b;

	while (2 * sqrB * x < 2 * sqrA * y)
	{
		if (d < 0) d = d + 2 * sqrB * x + 3 * sqrB;
		else d = d + 2 * sqrB * x - 2 * sqrA * y + 2 * sqrA + 3 * sqrB, y--;
		x++;

		for (int i = 0; i < 4; i++)
		{
			vertices[count * 3] = x * dx[i] + x0, vertices[count * 3 + 1] = y * dy[i] + y0;
			Normalize(vertices[count * 3], vertices[count * 3 + 1]), count++;
		}
	}

	d = sqrB * x * x + sqrA * y * y + 2 * sqrB * x - sqrA * y + 0.25 * sqrA + sqrB - sqrA * sqrB;
	while (y >= 0)
	{
		if (d >= 0) d = d - 2 * sqrA * y + 3 * sqrB;
		else d = d + 2 * sqrB * x - 2 * sqrA * y + 3 * sqrA + 2 * sqrB, x++;

		y--;

		for (int i = 0; i < 4; i++)
		{
			vertices[count * 3] = x * dx[i] + x0, vertices[count * 3 + 1] = y * dy[i] + y0;
			Normalize(vertices[count * 3], vertices[count * 3 + 1]), count++;
		}
	}
}

static void Bresenham(int x0, int y0, int a, int b)
{
	int sqrA = a * a, sqrB = b * b;
	int x = 0, y = b, d = 4 * sqrB + sqrA - 4 * sqrA * b;

	while (sqrB * x < sqrA * y)
	{
		if (d < 0) d = d + 4 * (2 * sqrB * x + 3 * sqrB);
		else d = d + 4 * (2 * sqrB * x - 2 * sqrA * y + 2 * sqrA + 3 * sqrB), y--;
		x++;

		for (int i = 0; i < 4; i++)
		{
			vertices[count * 3] = x * dx[i] + x0, vertices[count * 3 + 1] = y * dy[i] + y0;
			Normalize(vertices[count * 3], vertices[count * 3 + 1]), count++;
		}
	}

	d = 4 * sqrB * x * x + 4 * sqrA * y * y + 4 * 2 * sqrB * x - sqrA * y + sqrA + 4 * sqrB - 4 * sqrA * sqrB;
	while (y >= 0)
	{
		if (d >= 0) d = d - 4 * (2 * sqrA * y + 3 * sqrB);
		else d = d + 4 * (2 * sqrB * x - 2 * sqrA * y + 3 * sqrA + 2 * sqrB), x++;

		y--;

		for (int i = 0; i < 4; i++)
		{
			vertices[count * 3] = x * dx[i] + x0, vertices[count * 3 + 1] = y * dy[i] + y0;
			Normalize(vertices[count * 3], vertices[count * 3 + 1]), count++;
		}
	}
}

static void InitializeVertex()
{
	// Generate
	glGenVertexArrays(1, &VAO);
	glGenBuffers(1, &VBO);

	// Bind
	glBindVertexArray(VAO);
	glBindBuffer(GL_ARRAY_BUFFER, VBO);

	glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);

	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(float), (void*)0);
	glEnableVertexAttribArray(0);
}

#pragma endregion


void Render()
{
	glClearColor(0.5, 0.5, 0.5, 1);
	glClear(GL_COLOR_BUFFER_BIT);

	glBindVertexArray(VAO);
	glDrawArrays(GL_POINTS, 0, count);
}

int main()
{
	int x0, y0, a, b;
	std::cin >> x0 >> y0 >> a >> b;

	InitializeWindow();

	MiddlePoint(x0, y0, a, b);
	Bresenham(x0, y0, a, b);

	InitializeVertex();

	while (!glfwWindowShouldClose(window))
	{
		ProcessInput(window);

		Render();

		glfwSwapBuffers(window);
		glfwPollEvents();
	}

	glfwTerminate();
	return 0;
}