计算机图形学--中点椭圆算法原理及代码实现

椭圆的几何特性：

算法原理：

代码实现：

说明，我们这里讨论的椭圆都是对称轴平行于坐标轴的椭圆，对于其他方程较为复杂的椭圆我们不做讨论。

椭圆的几何特性：

首先我们考虑椭圆的几何特性。椭圆是抽对称图形，存在两条相互垂直的对称轴，中心位置在坐标原点的椭圆具有以下几何特性：
        1>对于一个圆心在原点的椭圆图形，椭圆上一个点，我们可以根据x=0和y=0两条对称轴，将其对称得到在圆上4个位于不同象限的点，如下图利用(xi, yi)可产生椭圆上标注的4个点。
        2>同样的，对于一个椭圆的1/4椭圆弧来讲，可以根据x=0, y=0两条对称轴，将其不断对称产生整个椭圆，如下图利用一象限椭圆弧确定的四分之椭圆可以对称得到整个椭圆。
        3>利用四分之椭圆对称得到整个椭圆的方法本质上是利用四分之椭圆上所有点对称得到整个椭圆的过程。

下面给出的椭圆绘制算法都是基于以上对椭圆几何特性的认识展开的。
说明：对于椭圆我们可以首先考虑当中心在坐标原点的情况，然后再根据圆心位置将椭圆上点的坐标进行平移，得到椭圆的实际位置。下面我们的考虑都是基于中心在坐标原点进行，实际算法中进行简单变化即可得到任意位置的椭圆。

算法原理：

中点椭圆算法的思想与中点画直线算法、中点圆形算法的思想一致：在前一已知点位置的基础上，选择下一绘制点可能像素位置中靠近椭圆的那个进行绘制。

对于对称轴平行于坐标轴的椭圆而言，假设平行于x轴方向的半轴长为a，平行于y轴方向的半轴长为b，那么我们可以获知，椭圆方程可以表示为： $F(x, y)=b^2x^2+a^2y^2-a^2b^2=0$ 。对于平面坐标系中的任意一点P(xi, yi)，我们可以记d(P)=F(xi, yi)，则存在：d(P)>0时，P点在椭圆外部；d(P)<0时，P点在椭圆内部；d(P)=0时，P点在椭圆上。基于以上认识，我们进行像素点的推算与选择。

我们考量对第一象限的四分之椭圆的顺时针绘制过程，然后对这个绘制过程中涉及到的所有应当绘制的像素点进行对称处理，得到整个椭圆。
        在对四分之椭圆进行考虑时，我们可以发现，在四分之椭圆的开始部分，x方向的变化速度大于y方向的变化速度（如上图区域1部分）；在接下来的部分，y方向的变化速度大于x方向的变化速度（如上图区域2部分）。为保证椭圆绘制像素点的密集程度相当，我们需要将这两部分分别讨论。对于第一部分，以x每次自增1为基础进行像素点的选择与绘制；对于第二部分，则以y每次自减1为基础进行像素点的选择与绘制。具体的选择与绘制过程在下面进行讨论。
a>对于x方向的变化速度大于y方向的变化速度的区域，有：
       假定当前已经确定绘制点(xpos, ypos)，那么下一个可能的绘制点为(xpos+1, ypos-1)或(xpos+1. ypos)，此时判断因子d=F(xpos+1, ypos-0.5)，则有如下结论成立：
        1>若d<0，则中点在椭圆内部，应取正右方像素，且判断因子应该变化为： $d'=F(xpos+2, ypos-0.5)\\=b^2(xpos+2)^2+a^2(ypos-0.5)^2-a^2b^2\\=b^2(xpos+1)^2+a^2(ypos-0.5)^2-a^2b^2+b^2(2xpos+3)\\=F(xpos+1, ypos-0.5)+b^2(2xpos+3)\\=d+b^2(2xpos+3)$
        即判断因子d的增量为： $b^2(2xpos+3)$
        2>若d>=0，则中点在椭圆外部（或椭圆上），应取右下方像素，且判断因子应该变化为：
$d'=F(xpos+2, ypos-1.5)\\=b^2(xpos+2)^2+a^2(ypos-1.5)^2-a^2b^2\\=b^2(xpos+1)^2+a^2(ypos-0.5)^2-a^2b^2+b^2(2xpos+3)+a^2(-2ypos+2)\\=F(xpos+1, ypos-0.5)+b^2(2xpos+3)+a^2(-2ypos+2)\\=d+b^2(2xpos+3)+a^2(-2ypos+2)$
        及判断因子d的增量为： $b^2(2xpos+3)+a^2(-2ypos+2)$
b>对于y方向的变化速度大于x方向的变化速度的区域，有：
        假定当前已经确定绘制点(xpos, ypos)，那么下一个可能的绘制点为(xpos, ypos-1)或(xpos+1. ypos-1)，此时判断因子d=F(xpos+0.5, ypos-1)，则有如下结论成立：
        1>若d<0，则中点在椭圆内部，应取右下方像素，且判断因子应该变化为： $d'=F(xpos+1.5, ypos-2)\\=b^2(xpos+1.5)+a^2(ypos-2)-a^2b^2\\=b^2(xpos+0.5)^2+a^2(ypos-1)^2-a^2b^2+b^2(2xpos+2)+a^2(-2ypos+3)\\=F(xpos+0.5, ypos-1)+b^2(2xpos+2)+a^2(-2ypos+3)\\=d+b^2(2xpos+2)+a^2(-2ypos+3)$
        即判断因子d的增量为： $b^2(2xpos+2)+a^2(-2ypos+3)$
        2>若d>=0，则中点在椭圆外部（或椭圆上），应取正下方像素，且判断因子应该变化为： $d'=F(xpos+0.5, ypos-2)\\=b^2(xpos+0.5)^2+a^2(ypos-2)^2-a^2b^2\\=b^2(xpos+0.5)^2+a^2(ypos-1)^2-a^2b^2+a^2(-2ypos+3)\\=F(xpos+0.5, ypos-1)+a^2(-2ypos+3)\\=d+a^2(-2ypos+3)$
         即判断因子d的增量为： $a^2(-2ypos+3)$
c>易知在斜率>-1的位置存在对于x方向的变化速度大于y方向的变化速度，即满足条件a；在斜率<-1的位置存在对于x方向的变化速度大于y方向的变化速度，即满足条件b。对于椭圆上任意一点，其斜率为 $\frac{dy}{dx}=-\frac{2b^2x}{2a^2y}=-\frac{b^2x}{a^2y}$ 。那么当 $-\frac{b^2x}{a^2y}>-1$ 即 $a^2y>b^2x$ 时，条件a成立，进行a条件下的像素点推算与绘制；否则条件b成立，进行对应推算与绘制。
d>易知，初始情况下，判断因子 $d=F(1, b-0.5)=b^2+a^2(b-0.5)^2-a^2b^2=b^2+a^2(0.25-b)$

代码实现：

//计算
void MidPointEllipse(const vector<pair<int, int>>& center, int a, int b, unsigned long color)
{//以下是椭圆的绘制int xpos = 0, ypos = b;int a2 = a * a, b2 = b * b;int d = b2 + a2 * (0.25 - b);while (a2 * ypos > b2 * xpos){quarterDraw(xpos, ypos, center, color);if (d < 0){d = d + b2 * ((xpos << 1) + 3);xpos += 1;}else{d = d + b2 * ((xpos << 1) + 3) + a2 * (-(ypos << 1) + 2);xpos += 1, ypos -= 1;}}d = b2 * (xpos + 0.5) * (xpos + 0.5) + a2 * (ypos - 1) * (ypos - 1) - a2 * b2;while (ypos > 0){quarterDraw(xpos, ypos, center, color);if (d < 0){d = d + b2 * ((xpos << 1) + 2) + a2 * (-(ypos << 1) + 3);xpos += 1, ypos -= 1;}else{d = d + a2 * (-(ypos << 1) + 3);ypos -= 1;}}
}
//画点
void quarterDraw(const int& xpos, const int& ypos, const vector<pair<int,int>>& center, unsigned long& color)
{CClientDC dc(mView); dc.SetPixel(center.first+ xpos, center.second+ ypos, color);dc.SetPixel(center.first+ xpos, center.second- ypos, color);                                                                    /**/dc.SetPixel(center.first- xpos, center.second+ ypos, color);dc.SetPixel(center.first- xpos, center.second- ypos, color);
}