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一個能交的出來了。。。。