/* generates contours using marching squares */
/* region size */
#define X_MAX 1.0
#define Y_MAX 1.0
#define X_MIN -1.0
#define Y_MIN -1.0
/* number of cells */
#define N_X 50
#define N_Y 50
/* contour value */
#define THRESHOLD 0.0
#include <GL/glut.h>
void display()
{
double f(double,double);
int cell(double, double, double, double);
void lines(int, int, int, double, double, double, double);
double data[N_X][N_Y];
int i,j;
int c;
glClear(GL_COLOR_BUFFER_BIT);
/* form data array from function */
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));
/* process each cell */
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();
}
/* define function f(x,y) */
double f(double x, double y)
{
double a=0.49, b=0.5;
/* Ovals of Cassini */
return (x*x+y*y+a*a)*(x*x+y*y+a*a)-4*a*a*x*x-b*b*b*b;
}
/* define cell vertices */
int cell(double a, double b, double c , double d)
{
int n=0;
if(a>THRESHOLD) n+=1;
if(b>THRESHOLD) n+=8;
if(c>THRESHOLD) n+=4;
if(d>THRESHOLD) n+=2;
return n;
}
/* draw line segments for each case */
void lines(int num, int i, int j, double a, double b, double c, double d)
{
void draw_one(int, int, int, double, double, double, double);
void draw_adjacent(int, int, int, double, double, double, double);
void draw_opposite(int, int, int, double, double, double, double);
switch(num)
{
case 1: case 2: case 4: case 7: case 8: case 11: case 13: case 14:
draw_one(num, i,j,a,b,c,d);
break;
case 3: case 6: case 9: case 12:
draw_adjacent(num,i,j,a,b,c,d);
break;
case 5: case 10:
draw_opposite(num, i,j,a,b,c,d);
break;
case 0: case 15: break;
}
}
void draw_one(int num, int i, int j, double a, double b, double c, double d)
{
double x1, y1, x2, y2;
double ox, oy;
double dx, 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)
{
case 1: case 14:
x1=ox;
y1=oy+dy*(THRESHOLD-a)/(d-a);
x2=ox+dx*(THRESHOLD-a)/(b-a);
y2=oy;
break;
case 2: case 13:
x1=ox;
y1=oy+dy*(THRESHOLD-a)/(d-a);
x2=ox+dx*(THRESHOLD-d)/(c-d);
y2=oy+dy;
break;
case 4: case 11:
x1=ox+dx*(THRESHOLD-d)/(c-d);
y1=oy+dy;
x2=ox+dx;
y2=oy+dy*(THRESHOLD-b)/(c-b);
break;
case 7: case 8:
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();
}
void draw_adjacent(int num, int i, int j, double a, double b,
double c, double d)
{
double x1, y1, x2, y2;
double ox, oy;
double dx, 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)
{
case 3: case 12:
x1=ox+dx*(THRESHOLD-a)/(b-a);
y1=oy;
x2=ox+dx*(THRESHOLD-d)/(c-d);
y2=oy+dy;
break;
case 6: case 9:
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();
}
void draw_opposite(int num, int i, int j, double a, double b,
double c, double d)
{
double x1,y1,x2,y2,x3,y3,x4,y4;
double ox, oy;
double dx, 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;
case 10:
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();
}
void myReshape(int w, int h)
{
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);
else gluOrtho2D(X_MIN * (GLfloat) w / (GLfloat) h, X_MAX * (GLfloat) w / (GLfloat) h, Y_MIN, Y_MAX);
glMatrixMode(GL_MODELVIEW);
}
int main(int argc, char **argv)
{
glutInit(&argc, argv);
glutInitWindowSize(500, 500);
glutCreateWindow("contour plot");
glutReshapeFunc(myReshape);
glutDisplayFunc(display);
glClearColor(0.0,0.0,0.0,1.0);
glColor3f(1.0,1.0,1.0);
glutMainLoop();
}
本文转自Phinecos(洞庭散人)博客园博客，原文链接：http://www.cnblogs.com/phinecos/archive/2007/09/27/907891.html，如需转载请自行联系原作者
