目录
椭圆特性
- 椭圆定义
椭圆：平面内到定点F1、F2的距离之和等于常数2a（2a>|F1F2|）的动点P的轨迹。
椭圆数学表达式：
\[\tag{1} |PF1|+|PF2|=2a \]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)为椭圆上任意一点。
根据椭圆定义，有
\[\begin{aligned} & |PF1| + |PF2| = 2a \\\\ & => \sqrt{(x+c)^2+y^2}+\sqrt{(x-c)^2+y^2} = 2a \\\\ & => \sqrt{(x+c)^2+y^2} = 2a - \sqrt{(x-c)^2+y^2}, 两边平方 \\\\ & => (x+c)^2+y^2=4a^2 + (x-c)^2+y^2 - 4a\sqrt{(x-c)^2+y^2} \\\\ & => 2xc=4a^2-2xc-4a\sqrt{(x-c)^2+y^2} \\\\ & => a^2-xc=a\sqrt{(x-c)^2+y^2}, 两边平方 \\\\ & => a^4-2a^2xc+(xc)^2=a^2[(x-c)^2+y^2] \\\\ & => (a^2-c^2)x^2+a^2y^2+a^2c^2=a^4 \\\\ & => (a^2-c^2)x^2+a^2y^2=a^2(a^2-c^2) \\\\ & => x^2/a^2+y^2/(a^2-c^2)=1 \end{aligned} \]∵\(2a>2c>0\)
∴\(a^2-c^2>0\)
设\(b^2=a^2-c^2\)
则椭圆方程可写为：
a称为椭圆的半长焦距，b称为椭圆的半短焦距。
方程(2)是以椭圆的中心点（2焦点中点）为原点，建立的局部坐标系，而在世界坐标系下，原点不一定在中点。设椭圆中心点坐标\((x_c, y_c)\)，则椭圆方程：
\[\tag{3} ({x-x_c\over a})^2+({y-y_c\over b})^2=1 \]可以写成通用形式：
\[\tag{4} Ax^2+By^2+Cxy+Dx+Ey+F=0 \]其中，常数A,B,C,D,E,F可根据(3)求出。
为了便于下文讨论，我们用\(2r_x\)表示长轴，\(2r_y\)表示短轴。利用极坐标，椭圆方程可写成参数方程形式：
\[\tag{5} \begin{aligned} x=x_c+r_x\cos\theta \\\\ y=y_c+r_y\sin\theta \end{aligned} \]θ称为椭圆离心角（eccentric angle）。
椭圆离心角的几何意义，见下图：
Bresenham算法画椭圆
类似于Bresenham算法画光栅圆。
给定参数：半焦距\(r_x, r_y\), 中心点\((x_c, y_c)\)。先以中心点为原点建立局部坐标系，得到椭圆上的点(x,y)后，再将其移动到\((x_c, y_c)\)为中心点的椭圆上。如果长轴、短轴没有与x、y轴平行，就绕中心进行旋转，目前这里只考虑平移，不考虑旋转。
局部坐标系下，定义椭圆函数：
\[\tag{6} f_{ellipse}(x, y) = r_y^2x^2+r_x^2y^2-r_x^2r_y^2 \]椭圆函数可以表示 点(x,y)和椭圆位置关系：
\[\tag{7} f_{ellipse}(x, y)=\begin{cases} < 0, & (x,y)在椭圆边界内 \\\\ = 0, & (x,y)在椭圆边界上 \\\\ > 0, & (x,y)在椭圆边界外 \end{cases} \]\(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轴方向逐像素绘制。
椭圆方程：
\[f_{ellipse}(x, y)=r_y^2x^2+r_x^2y^2-r_x^2r_y^2=0 \]两边对x求导，
\[r_y^2(2x)+r_x^2(2y{dy\over dx}) - 0=0 \\\\ {dy\over dx}=-{r_y^2x \over r_x^2y} \tag{8} \]区域1、2交接处，切线斜率\({dy\over dx}=-1\)，即
\[-{r_y^2x \over r_x^2y}=-1 \\\\ => r_y^2x=r_x^2y \]区域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，有
\[\begin{aligned} p1_{k+1} &= f_{ellipse}(x_{k+1}, y_{k+1}-1/2) \\\\ & = r_y^2[(x_k+1)+1]^2 + r_x^2(y_{k+1}-1/2)^2 - r_x^2r_y^2 \end{aligned} \]做差值，可得递推公式：
\[\begin{aligned} p1_{k+1}-p1_k &= 2r_y^2(x_k+1)+r_y^2+r_x^2[(y_{k+1}-{1\over 2})^2 - (y_k-{1\over 2})^2] \\\\ &= 2r_y^2x_{k+1}+r_y^2+r_x^2[(y_{k+1}-{1\over 2})^2-(y_k-{1\over 2})^2] \end{aligned} \]而\(y_{k_1}\)值取决于\(p_k\)：
\[y_{k+1}=\begin{cases} y_k, & y_k更近<=> p1_k<0 \\\\ y_k-1, & y_k-1更近<=> p1_k>=0 \end{cases} \]因此，决策参数递推公式：
\[p_{k+1}-p_k=\begin{cases} 2x_{k+1}r_y^2+r_y^2, & p1_k<0 \\\\ 2x_{k+1}r_y^2+r_y^2-2y_{k+1}r_x^2, & p1_k >= 0 \end{cases} \]初值\(p1_0\)：起始点\((0, r_y)\)，根据(9)，可得，
\[\begin{aligned} p1_0 &= f_{ellipse}(1, r_y-{1\over 2} \\\\ & = r_y^2 + r_x^2(r_y-{1\over 2})^2-r_x^2r_y^2 \\\\ & = r_y^2 - r_x^2r_y+{1\over 4}r_x^2 \end{aligned} \]区域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)\)，那么决策参数：
\[\begin{aligned} p2_k &= f_{ellipse}(x_k+{1\over 2}, y_k-1) \\\\ & = r_y^2(x_k+{1\over 2})^2 + r_x^2(y_k-1)^2-r_x^2r_y^2 \end{aligned} \]若\(p2_k>0\)，则中点在椭圆外，下一点选择\(x_k\)点；
若\(p2_k<=0\)，则中点在椭圆内，下一点选择\(x_k+1\)点。
k取值k+1，可得
\[\begin{aligned} p2_{k+1} &= f_{ellipse}(x_{k+1}+{1\over 2}, y_{k+1}-1) \\\\ &= r_y^2(x_{k+1}+{1\over 2})^2 + r_x^2(y_{k+1}-1)^2-r_x^2r_y^2 \\\\ &= -2r_x^2(y_k-1)+r_x^2+r_y^2[(x_{k+1}+{1\over 2})^2-(x_k+{1\over 2})^2] \\\\ &= \begin{cases} -2r_x^2y_{k+1}+r_x^2, & x_k更近 <=>p2_k > 0 \\\\ -2r_x^2y_{k+1}+r_x^2+2r_y^2x_{k+1}, & x_k+1更近<=> p2_k<=0 \end{cases} \end{aligned} \]初值\(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)\)，决策参数
\[p1_{k+1}=p1_k+2r_y^2x_{k+1}-2r_x^2y_{k+1}+r_y^2 \]这里\(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)\)，决策参数
\[p2_{k+1}=p2_k+2r_y^2x_{k+1}-2r_x^2y_{k+1}+r_x^2 \]重复该步骤，直到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);
}
}
```