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