Bresenham算法画椭圆和斜椭圆
CG课程的第一次作业，大四才开始学CG也算是很特别【然后就迟交了一天】。
Bresenham算法用于把连续曲线投影到平面像素中，思想是只要能判断x和y哪个增量更大，就可以按x+1（或y+1）之后y（或x）是否+1来画下一个像素。判断是用x还是y的标准是斜率大于1还是小于1，在这个基础上网上能够搜到的椭圆画法做了一些优化，1.只用画一个象限内的曲线，另外三个象限直接对称过去；2.以x增量更大（斜率小于1）的情况为例，判断y是否+1的方法不是把y和y+1都算一遍，而是用y+0.5在圆内还是圆外判断；3.再进一步，每次都带入方程去算d=F(x,y+0.5)也很慢，可以考虑用算增量代替。这部分我觉得讲得比较好的是华中科技大学的这篇博客。
比较麻烦的是斜椭圆的情况，首先，按照前面的思路，要重新把公式推一遍，斜椭圆的公式是在椭圆的基础上把x和y用旋转变换替换，你博能不能插入公式啊，试试：
椭圆公式：
$$ \frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}=1 $$
旋转变换：
$$
\left[\begin{array}{cc}
\cos \beta & -\sin \beta \\
\sin \beta & \cos \beta
\end{array}\right]\left[\begin{array}{l}
x \\
y
\end{array}\right]=\left[\begin{array}{l}
x \cos \beta-y \sin \beta \\
x \sin \beta+y \cos \beta
\end{array}\right]
$$
代入得到斜椭圆公式：
$$
\frac{(x \cos \beta-y \sin \beta)^{2}}{a^{2}}+\frac{(x \sin \beta+y \cos \beta)^{2}}{b^{2}}=1
$$
在这里推荐Mathpix Snipping Tool，远离手打latex从我做起【不是】
得到曲线函数之后，按照算法核心思想，首先要判断的是切线斜率，在这里要用到隐函数求导：
$$
\frac{d y}{d x}=-\frac{F_{x}}{F_{y}}
$$
其中：
$$
F(x,y)=\frac{(x \cos \beta-y \sin \beta)^{2}}{a^{2}}+\frac{(x \sin \beta+y \cos \beta)^{2}}{b^{2}}-1
$$
为了方便编程，我去掉了上面对椭圆曲线做的三个优化，仅考虑最简单的Bresenham算法，那么公式的推导就可以到此为止了。对于优化1，在斜椭圆本来就不成立，斜椭圆不再有四个象限的对称性质，只画第一象限不能完成曲线；但是有中心对称，所以我只画了第一象限和第四象限，另外两个象限对称处理。对于优化2，最大的问题是斜椭圆在同一象限内的凹凸性不再保持不变，而F(x,y+0.5)在圆外代表我们该取y还是y+1其实是和凹凸性有关的，画个图应该很好理解。如果想要保留这个优化，就需要去计算凹凸性，也就是求隐函数二阶导，太麻烦了，我选择每次计算一下F(x,y)和F(x,y+1)。
最终代码如下，前面两个函数是助教给的生成ppm文件的代码
#define _CRT_SECURE_NO_WARNINGS #include <iostream> #include<math.h> //#include<cstdlib> using namespace std; void ppmRead(char* filename, unsigned char* data, int* w, int* h) { char header[1024]; FILE* fp = NULL; int line = 0; fp = fopen(filename, "rb"); while (line < 2) { fgets(header, 1024, fp); if (header[0] != '#') { ++line; } } sscanf(header, "%d %d\n", w, h); fgets(header, 20, fp); fread(data, (*w)*(*h) * 3, 1, fp); fclose(fp); } void ppmWrite(const char* filename, unsigned char* data, int w, int h) { FILE* fp; fp = fopen(filename, "wb"); fprintf(fp, "P6\n%d %d\n255\n", w, h); fwrite(data, w*h * 3, 1, fp); fclose(fp); } void draw_ellipse(int ex,int ey,float ra,float rb) { unsigned char final_data[400*300*3] = { 0 }; char data[400][300]; int x=0; int y=rb; float d1=rb*rb+ra*ra*(-rb+0.25); while(rb*rb*((float)x+1)<ra*ra*((float)y-0.5)) { data[x+ex][y+ey]=255; // cout<<x<<" "<<y<<endl; data[-x+ex][y+ey]=255; data[-x+ex][-y+ey]=255; data[x+ex][-y+ey]=255; if(d1<=0) { d1=d1+rb*rb*(2*(float)x+3); x=x+1; } else { d1=d1+rb*rb*(2*(float)x+3)+ra*ra*(-2*(float)y+2); x=x+1; y=y-1; } } float d2=rb*rb*((float)x+0.5)*((float)x+0.5)+ra*ra*((float)y-1)*((float)y-1)-ra*ra*rb*rb; while(y>=0) { data[x+ex][y+ey]=255; data[-x+ex][y+ey]=255; data[-x+ex][-y+ey]=255; data[x+ex][-y+ey]=255; if(d2<=0) { d2=d2+rb*rb*(2*(float)x+2)+ra*ra*(-2*(float)y+3); x=x+1; y=y-1; } else { d2=d2+ra*ra*(-2*(float)y+3); y=y-1; } } for(int i=0;i<300;i++) for(int j=0;j<400;j++) { final_data[(i*400+j)*3]=data[j][i]; final_data[(i*400+j)*3+1]=data[j][i]; final_data[(i*400+j)*3+2]=data[j][i]; } ppmWrite("test.ppm", final_data, 400, 300); } float cal_fx(float ra,float rb,float the,float x,float y) { return 2*rb*rb*(x*cos(the)-y*sin(the))*cos(the)+2*ra*ra*(x*sin(the)+y*cos(the))*sin(the); } float cal_fy(float ra,float rb,float the,float x,float y) { return -2*rb*rb*(x*cos(the)-y*sin(the))*sin(the)+2*ra*ra*(x*sin(the)+y*cos(the))*cos(the); } float cal_slope(float ra,float rb,float the,float x,float y) { return -cal_fx(ra,rb,the,x,y)/cal_fy(ra,rb,the,x,y); } float cal_f(float ra,float rb,float the,float x,float y) { return rb*rb*(x*cos(the)-y*sin(the))*(x*cos(the)-y*sin(the))+ ra*ra*(x*sin(the)+y*cos(the))*(x*sin(the)+y*cos(the))-ra*ra*rb*rb; } void draw_ellipse_bonus(int ex,int ey,float ra,float rb,float the) { unsigned char final_data[400*300*3] = { 0 }; char data[400][300]; int x=0; int y=(int)(sqrt(1.0/(sin(the)*sin(the)/ra/ra+cos(the)*cos(the)/rb/rb))+0.5); float s1=cal_slope(ra,rb,the,x,y); while(y>=0) { // system("pause"); data[x+ex][y+ey]=255; data[-x+ex][-y+ey]=255; s1=cal_slope(ra,rb,the,x,y); cout<<x<<" "<<y<<" "<<s1<<endl; if(abs(s1)<1) { x=x+1; if(abs(cal_f(ra,rb,the,x,y-1))<abs(cal_f(ra,rb,the,x,y)) &&abs(cal_f(ra,rb,the,x,y-1))<abs(cal_f(ra,rb,the,x,y+1))) { y=y-1; } if(abs(cal_f(ra,rb,the,x,y+1))<abs(cal_f(ra,rb,the,x,y)) &&abs(cal_f(ra,rb,the,x,y+1))<abs(cal_f(ra,rb,the,x,y-1))) { y=y+1; } } else { y=y-1; if(abs(cal_f(ra,rb,the,x-1,y))<abs(cal_f(ra,rb,the,x,y)) &&abs(cal_f(ra,rb,the,x-1,y))<abs(cal_f(ra,rb,the,x+1,y))) { x=x-1; } if(abs(cal_f(ra,rb,the,x+1,y))<abs(cal_f(ra,rb,the,x,y)) &&abs(cal_f(ra,rb,the,x+1,y))<abs(cal_f(ra,rb,the,x-1,y))) { x=x+1; } } } x=0; y=-(int)(sqrt(1.0/(sin(the)*sin(the)/ra/ra+cos(the)*cos(the)/rb/rb))+0.5); s1=cal_slope(ra,rb,the,x,y); while(y<=0) { // system("pause"); data[x+ex][y+ey]=255; data[-x+ex][-y+ey]=255; s1=cal_slope(ra,rb,the,x,y); cout<<x<<" "<<y<<" "<<s1<<endl; if(abs(s1)<1) { x=x+1; if(abs(cal_f(ra,rb,the,x,y-1))<abs(cal_f(ra,rb,the,x,y)) &&abs(cal_f(ra,rb,the,x,y-1))<abs(cal_f(ra,rb,the,x,y+1))) { y=y-1; } if(abs(cal_f(ra,rb,the,x,y+1))<abs(cal_f(ra,rb,the,x,y)) &&abs(cal_f(ra,rb,the,x,y+1))<abs(cal_f(ra,rb,the,x,y-1))) { y=y+1; } } else { y=y+1; if(abs(cal_f(ra,rb,the,x-1,y))<abs(cal_f(ra,rb,the,x,y)) &&abs(cal_f(ra,rb,the,x-1,y))<abs(cal_f(ra,rb,the,x+1,y))) { x=x-1; } if(abs(cal_f(ra,rb,the,x+1,y))<abs(cal_f(ra,rb,the,x,y)) &&abs(cal_f(ra,rb,the,x+1,y))<abs(cal_f(ra,rb,the,x-1,y))) { x=x+1; } } } for(int i=0;i<300;i++) for(int j=0;j<400;j++) { final_data[(i*400+j)*3]=data[j][i]; final_data[(i*400+j)*3+1]=data[j][i]; final_data[(i*400+j)*3+2]=data[j][i]; } ppmWrite("test.ppm", final_data, 400, 300); } int main() { draw_ellipse(200,150,60,70); // draw_ellipse_bonus(200,150,60,70,0.5); return 0; }
得到的ppm图片，博客园不支持上传，那就不传啦。
其实这个代码多半还是有bug的，bug在于在第一象限遍历的方式是从x=0出发（取大于0的y值），x++或者y++直到y小于等于0；在第四象限遍历的方式是从x=0出发（取小于0的y值），x++或者y--直到y大于等于0。脑跑了一下觉得好像很ok，但是总觉得并不是很放心。不过当时已经过ddl快24小时了，就先一两个小时rush一个能交的出来了。。。。
