椭圆特性
- 椭圆定义
椭圆：平面内到定点F1、F2的距离之和等于常数2a（2a>|F1F2|）的动点P的轨迹。
椭圆数学表达式：
F1、F2称为椭圆的2个焦点，两焦点之间距离2c（|F1F2|=2c）称为焦距。
椭圆与两焦点连线相交所得弦为长轴，长2a，a也称为半长轴的长；椭圆与两焦点垂线相交所得弦为短轴，长2b，b也称为半短轴的长。
平面上动点P轨迹形成椭圆示意图：
- 椭圆方程
以焦点F1、F2中点为原点，以F1->F2方向为x轴正方向，建立直角坐标系（局部坐标系）。记F2(c,0)，则F1(-c, 0)。设P(x, y)为椭圆上任意一点。
根据椭圆定义，有
∵\(2a>2c>0\)
∴\(a^2-c^2>0\)
设\(b^2=a^2-c^2\)
则椭圆方程可写为：
a称为椭圆的半长焦距，b称为椭圆的半短焦距。
方程(2)是以椭圆的中心点（2焦点中点）为原点，建立的局部坐标系，而在世界坐标系下，原点不一定在中点。设椭圆中心点坐标\((x_c, y_c)\)，则椭圆方程：
可以写成通用形式：
其中，常数A,B,C,D,E,F可根据(3)求出。
为了便于下文讨论，我们用\(2r_x\)表示长轴，\(2r_y\)表示短轴。利用极坐标，椭圆方程可写成参数方程形式：
θ称为椭圆离心角（eccentric angle）。
椭圆离心角的几何意义，见下图：
Bresenham算法画椭圆
类似于Bresenham算法画光栅圆。
给定参数：半焦距\(r_x, r_y\), 中心点\((x_c, y_c)\)。先以中心点为原点建立局部坐标系，得到椭圆上的点(x,y)后，再将其移动到\((x_c, y_c)\)为中心点的椭圆上。如果长轴、短轴没有与x、y轴平行，就绕中心进行旋转，目前这里只考虑平移，不考虑旋转。
局部坐标系下，定义椭圆函数：
椭圆函数可以表示 点(x,y)和椭圆位置关系：
\(f_{ellipse}(x, y)\)可用来作为算法决策参数。
- 以x轴 or y轴逐像素计算？
本质是看x轴、y轴，哪个方向变化大，就按那个轴方向逐个像素计算点位置。如，画直线时，取决于斜率|m|与1的关系；画圆时，x、y轴方向变化相同，都可以。而画椭圆时，如何判断？
椭圆的切线斜率可以代表沿着曲线的x、y轴方向变化大小。
下图是\(r_y<r_x\)时，根据切削斜率对椭圆第一象限进行划分的示意图：
在区域1，斜率0~-1，所以\(|{dy\over dx}|<1\)，可按+x轴方向逐像素绘制；
在区域2，斜率-1~-∞，所以\(|{dy\over dx}|>1\)，可按-y轴方向逐像素绘制。
椭圆方程：
两边对x求导，
区域1、2交接处，切线斜率\({dy\over dx}=-1\)，即
区域1：\({dy\over dx}=-{r_y^2x \over r_x^2y}>-1\)，即\(r_y^2x < r_x^2y\)
区域2：\(r_y^2>=r_x^2y\)
区域1
- 计算下一个点位置
因为在第一象限，y随着x增加而减少。当进行到第k步，像素点位置\((x_k, y_k)\)，区域1中，下一个像素点位置只能是\((x_k+1, y_k)\) or \((x_k+1, y_k-1)\)，可以用2备选点的中点与椭圆位置关系，判断哪个点距离椭圆更近，从而选择近点。
因此，区域1决策参数：
由(7)可知，如果\(p1_k<0\)，那么中点位于椭圆内，\(y_k\)点更接近椭圆；否则，\(y_k-1\)点更接近。
取k为k+1，有
做差值，可得递推公式：
而\(y_{k_1}\)值取决于\(p_k\)：
因此，决策参数递推公式：
初值\(p1_0\)：起始点\((0, r_y)\)，根据(9)，可得，
区域2
从区域1与2的边界（切线斜率=-1）点开始，按-y方向逐像素取样。
假设进行到第k步，对应像素点\((x_k, y_k)\)，那么第k+1步可选点：\((x_k, y_k-1)\) or \((x_k+1, y_k-1)\)。
2备选的中点\(f_{ellipse}(x_k+{1\over 2}, y_k-1)\)，那么决策参数：
若\(p2_k>0\)，则中点在椭圆外，下一点选择\(x_k\)点；
若\(p2_k<=0\)，则中点在椭圆内，下一点选择\(x_k+1\)点。
k取值k+1，可得
初值\(p2_0\)：区域2的起始点\((x_0, y_0)\)就是区域1的结束点。可以在区域1画完后，记录结束点作为区域2起始点。
- 对称性求其他3象限椭圆
求出椭圆第一象限部分一点(x, y)后，可用对称性得到另外3个象限的点：(-x, y) (y, -x) (-x, -y)。
算法步骤
- 输入\(r_x, r_y\)、椭圆中心\((x_c, y_c)\)，得到椭圆（中心在原点）上的第一个点 \((x_0, y_0)=(0, r_y)\)；
- 计算区域1决策参数初值，
- 在区域1中，从k=0开始，对每个\(x_k\)位置的像素点进行计算：
如果\(p1_k<0\)，则椭圆下一个点为\((x_k+1, y_k)\)，决策参数
如果\(p1_k>=0\)，则下一个\((x_k+1, y_k-1)\)，决策参数
这里\(x_{k+1}=x_k+1, y_{k+1}=y_k-1\)
重复该步骤，直到\(2r_y^2x>=2r_x^2y\)
- 记录区域1得到的最后一个点\((x_0, y_0)\)，用例计算区域2决策参数初值，
- 在区域2中，从k=0开始，对每个\(y_k\)位置的像素点进行计算：
如果\(p2_k>0\)，则椭圆下一点\((x_k, y_k-1)\)，决策参数
如果\(p2_k<=0\)，则椭圆下一点\((x_k+1, y_k-1)\)，决策参数
重复该步骤，直到y=0（终点为\((r_x, 0)\)）
6. 通过对称性获得其他3个象限的对称点；
7. 将得到的椭圆点位置(x, y)由原点为椭圆中心平移到\((x_c, y_c)\)，得到最终像素点位置，有\(x=x+x_c, y=y+y_c\)。
算法程序
void setPixel(int x, int y)
{
glBegin(GL_POINTS);
glVertex2i(x, y);
glEnd();
}
// 浮点数向上取整
inline int round(const float a)
{
return (int)(a + 0.5);
}
// 绘制椭圆上的点(x, y)及其对称点, 椭圆中心(xc, yc)
// (x, y)是原点在椭圆中心的局部坐标系下的坐标
void ellipsePlotPoints(int xc, int yc, int x, int y)
{
setPixel(xc + x, yc + y);
// 对称点
setPixel(xc - x, yc + y);
setPixel(xc + x, yc - y);
setPixel(xc - x, yc - y);
}
// 绘制中心点在(xc, yc)的椭圆, x轴半焦距rx, y轴半焦距ry
// 使用Bresenham算法思想, 根据2备选像素点位置的中点与椭圆位置关系,
// 判断并选择距离椭圆更近的点作为下一个点
void ellipseMidpoint(int xc, int yc, int rx, int ry)
{
int rx2 = rx * rx;
int ry2 = ry * ry;
int p;
int x = 0, y = ry; // (区域1)初始点
int px = 0, py = 2 * rx2 * y; // 像素计算终止条件项
// 绘制椭圆起始点
ellipsePlotPoints(xc, yc, x, y);
// 区域1
p = round(ry2 - rx2 * ry + 0.25 * rx2);
while (px < py) {
x++;
px += 2 * ry2; // 通过累加来计算 2ry^2 x[k+1], 而不是每次都用乘法
if (p < 0) {
p += ry2 + px;
}
else {
y--;
py -= 2 * rx2;
p += ry2 + px - py;
}
ellipsePlotPoints(xc, yc, x, y);
}
// 区域2
p = round(ry2 * (x + 0.5) * (x + 0.5) + rx2 * (y - 1) * (y - 1) - rx2 *
ry2);
while (y > 0) {
y--;
py -= 2 * rx2;
if (p > 0) {
p += rx2 - py;
}
else {
x++;
px += 2 * ry2;
p += rx2 - py + px;
}
ellipsePlotPoints(xc, yc, x, y);
}
}
```