递进网格算法绘制等高线

<!--<br><br>Code highlighting produced by Actipro CodeHighlighter (freeware)<br>http://www.CodeHighlighter.com/<br><br>-->/*generatescontoursusingmarchingsquares*/

/*regionsize*/

#defineX_MAX1.0
#defineY_MAX1.0
#defineX_MIN-1.0
#defineY_MIN-1.0

/*numberofcells*/

#defineN_X50
#defineN_Y50

/*contourvalue*/

#defineTHRESHOLD0.0

#include<GL/glut.h>

voiddisplay()
{
doublef(double,double);
intcell(double,double,double,double);
voidlines(int,int,int,double,double,double,double);

doubledata[N_X][N_Y];
inti,j;
intc;

glClear(GL_COLOR_BUFFER_BIT);
/*formdataarrayfromfunction*/

for(i=0;i<N_X;i++)for(j=0;j<N_Y;j++)
data[i][j]=f(X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0),Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0));

/*processeachcell*/

for(i=0;i<N_X;i++)for(j=0;j<N_Y;j++)
{
c=cell(data[i][j],data[i+1][j],data[i+1][j+1],data[i][j+1]);
lines(c,i,j,data[i][j],data[i+1][j],data[i+1][j+1],data[i][j+1]);
}
glFlush();
}

/*definefunctionf(x,y)*/

doublef(doublex,doubley)
{

doublea=0.49,b=0.5;

/*OvalsofCassini*/

return(x*x+y*y+a*a)*(x*x+y*y+a*a)-4*a*a*x*x-b*b*b*b;
}


/*definecellvertices*/

intcell(doublea,doubleb,doublec,doubled)
{
intn=0;
if(a>THRESHOLD)n+=1;
if(b>THRESHOLD)n+=8;
if(c>THRESHOLD)n+=4;
if(d>THRESHOLD)n+=2;
returnn;
}


/*drawlinesegmentsforeachcase*/

voidlines(intnum,inti,intj,doublea,doubleb,doublec,doubled)
{
voiddraw_one(int,int,int,double,double,double,double);
voiddraw_adjacent(int,int,int,double,double,double,double);
voiddraw_opposite(int,int,int,double,double,double,double);
switch(num)
{
case1:case2:case4:case7:case8:case11:case13:case14:
draw_one(num,i,j,a,b,c,d);
break;
case3:case6:case9:case12:
draw_adjacent(num,i,j,a,b,c,d);
break;
case5:case10:
draw_opposite(num,i,j,a,b,c,d);
break;
case0:case15:break;
}
}

voiddraw_one(intnum,inti,intj,doublea,doubleb,doublec,doubled)
{
doublex1,y1,x2,y2;
doubleox,oy;
doubledx,dy;
dx=(X_MAX-(X_MIN))/(N_X-1.0);
dy=(Y_MAX-(Y_MIN))/(N_Y-1.0);
ox=X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0);
oy=Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0);
switch(num)
{
case1:case14:
x1=ox;
y1=oy+dy*(THRESHOLD-a)/(d-a);
x2=ox+dx*(THRESHOLD-a)/(b-a);
y2=oy;
break;
case2:case13:
x1=ox;
y1=oy+dy*(THRESHOLD-a)/(d-a);
x2=ox+dx*(THRESHOLD-d)/(c-d);
y2=oy+dy;
break;
case4:case11:
x1=ox+dx*(THRESHOLD-d)/(c-d);
y1=oy+dy;
x2=ox+dx;
y2=oy+dy*(THRESHOLD-b)/(c-b);
break;
case7:case8:
x1=ox+dx*(THRESHOLD-a)/(b-a);
y1=oy;
x2=ox+dx;
y2=oy+dy*(THRESHOLD-b)/(c-b);
break;
}
glBegin(GL_LINES);
glVertex2d(x1,y1);
glVertex2d(x2,y2);
glEnd();
}

voiddraw_adjacent(intnum,inti,intj,doublea,doubleb,
doublec,doubled)
{
doublex1,y1,x2,y2;
doubleox,oy;
doubledx,dy;
dx=(X_MAX-(X_MIN))/(N_X-1.0);
dy=(Y_MAX-(Y_MIN))/(N_Y-1.0);
ox=X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0);
oy=Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0);
switch(num)
{
case3:case12:
x1=ox+dx*(THRESHOLD-a)/(b-a);
y1=oy;
x2=ox+dx*(THRESHOLD-d)/(c-d);
y2=oy+dy;
break;
case6:case9:
x1=ox;
y1=oy+dy*(THRESHOLD-a)/(d-a);
x2=ox+dx;
y2=oy+dy*(THRESHOLD-b)/(c-b);
break;
}
glBegin(GL_LINES);
glVertex2d(x1,y1);
glVertex2d(x2,y2);
glEnd();

}

voiddraw_opposite(intnum,inti,intj,doublea,doubleb,
doublec,doubled)
{
doublex1,y1,x2,y2,x3,y3,x4,y4;
doubleox,oy;
doubledx,dy;
dx=(X_MAX-(X_MIN))/(N_X-1.0);
dy=(Y_MAX-(Y_MIN))/(N_Y-1.0);
ox=X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0);
oy=Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0);
switch(num)
{
case5:
x1=ox;
y1=oy+dy*(THRESHOLD-a)/(d-a);
x2=ox+dx*(THRESHOLD-a)/(b-a);
y2=oy;
x3=ox+dx*(THRESHOLD-d)/(c-d);
y3=oy+dy;
x4=ox+dx;
y4=oy+dy*(THRESHOLD-b)/(c-b);
break;
case10:
x1=ox;
y1=oy+dy*(THRESHOLD-a)/(d-a);
x2=ox+dx*(THRESHOLD-d)/(c-d);
y2=oy+dy;
x3=ox+dy*(THRESHOLD-a)/(b-a);
y3=oy;
x4=ox+dx;
y4=oy+dy*(THRESHOLD-b)/(c-b);
break;
}
glBegin(GL_LINES);
glVertex2d(x1,y1);
glVertex2d(x2,y2);
glVertex2d(x3,y3);
glVertex2d(x4,y4);
glEnd();
}

voidmyReshape(intw,inth)
{
glViewport(0,0,w,h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
if(w<=h)gluOrtho2D(X_MIN,X_MAX,Y_MIN*(GLfloat)h/(GLfloat)w,Y_MAX*(GLfloat)h/(GLfloat)w);
elsegluOrtho2D(X_MIN*(GLfloat)w/(GLfloat)h,X_MAX*(GLfloat)w/(GLfloat)h,Y_MIN,Y_MAX);
glMatrixMode(GL_MODELVIEW);
}

intmain(intargc,char**argv)
{
glutInit(&argc,argv);
glutInitWindowSize(500,500);
glutCreateWindow("contourplot");
glutReshapeFunc(myReshape);
glutDisplayFunc(display);
glClearColor(0.0,0.0,0.0,1.0);
glColor3f(1.0,1.0,1.0);
glutMainLoop();
}