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一个能交的出来了。。。。