ContourLine

Posted 2020-08-03 GIS开发

#define MULTI_PLOT true //Determine whether or not to plot multiple iterations.

#define X_MAX 1.0   // Define extent of reference plane, used in call to gluOrtho2D(...)
#define Y_MAX 1.0
#define X_MIN -1.0
#define Y_MIN -1.0
#define N_X 640      // Number of cells in x and y dimensions
#define N_Y 480
#define OX (-1 + i * dx)
#define OY (1 - j * dy)

#include <stdio.h>
#include <GL/glut.h>
#include <iostream>>

double h = 0.0;    // Height of reference plane in which contour is to be plotted --  f(x,y)=h
double dx = 2 / (double)(N_X - 1);
double dy = 2 / (double)(N_Y - 1);

double data[N_X][N_Y];   // Array to hold function values (a,b,c,d) at corners of each cell
double f(double x, double y);    // Function whose contour is to be plotted
int cell(double a, double b, double c, double d);  // Helper ftn to determine cell type from f corner values
void lines(int, int, int, double, double, double, double);  // Helper ftn to draw correct line in a cell
void display();

void multidisplay(){
    glClear(GL_COLOR_BUFFER_BIT);

    h = 0;
    for (; h < 1; h += .01){
        display();
    }
}
void display()   // Display callback function - called when window appears; drawing done here
{
    if (!MULTI_PLOT) glClear(GL_COLOR_BUFFER_BIT);

    // form data array from function
    for (int i = 0; i < N_X; i++)
        for (int j = 0; j < N_Y; j++)
            data[i][j] = f(OX, OY);

    // double loop to process each cell (loop through each cell i,j; determine cell type, draw line in cell)
    int k = 0;
    for (int i = 0; i < N_X - 1; i++){
        for (int j = 0; j < N_Y - 1; j++){
            int c = cell(data[i][j], data[i + 1][j], data[i][j + 1], data[i + 1][j + 1]);  // return cell type c (0-15)
            lines(c, i, j, data[i][j], data[i + 1][j], data[i][j + 1], data[i + 1][j + 1]);  // draw correct line
            k++;
        }
    }


    glFlush();
}

double f(double x, double y)   // Function to be plotted, example: Ovals of Cassini
{
    double a = 0.48, b = 0.5;
    return (x*x + y*y + a*a)*(x*x + y*y + a*a) - 4 * a*a*x*x - b*b*b*b;
}

void init()
{
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluOrtho2D(X_MIN, X_MAX, Y_MIN, Y_MAX);
    glMatrixMode(GL_MODELVIEW);
}

int main(int argc, char **argv)
{
    glutInit(&argc, argv);
    glutInitWindowPosition(50, 50);
    glutInitWindowSize(500, 500);
    glutCreateWindow("Oval of Cassini Contour plot");
    init();
    glutDisplayFunc(MULTI_PLOT ? multidisplay : display);           // White drawing color
    glClearColor(0.0, 0.0, 0.0, 1.0);   // Black background
    glColor3f(1.0, 1.0, 1.0);
    glutMainLoop();
}

int cell(double a, double b, double c, double d)  // Determine and return cell type
{
    int n = 0;
    if (a > h) n += 1;
    if (b > h) n += 8;
    if (c > h) n += 4;
    if (d > h) n += 2;
    return n;
}

void draw_one(int n, int i, int j, double a, double b, double c, double d){
    double x1 = 0, x2 = 0, y1 = 0, y2 = 0;
    switch (n){
    case 1:
    case 14:	x1 = OX;
        y1 = OY - (dy * ((h - a) / (c - a)));
        x2 = OX + (dx * ((h - a) / (b - a)));
        y2 = OY;
        //printf("1: (%f, %f)\n", (h - a) / (d - a), (h - a) / (b - a));
        break;
    case 4:
    case 11: 	x1 = OX;
        y1 = OY - (dy * ((h - a) / (c - a)));
        x2 = OX + (dx * ((h - c) / (d - c)));
        y2 = OY - dy;
        //printf("2: (%f, %f)\n", (h - c)/(b - c), (h - c)/(a - c));
        break;
    case 2:
    case 13: 	x1 = OX + dx;
        y1 = OY - (dy * ((h - b) / (d - b)));
        x2 = OX + (dx * ((h - c) / (d - c)));
        y2 = OY - dy;
       // printf("3: (%f, %f)\n", (h - d) / (a - d), (h - d) / (c - d));
        break;
    case 7:
    case 8: 	x1 = OX + dx;
        y1 = OY - (dy * ((h - b) / (d - b)));
        x2 = OX + (dx * ((h - a) / (b - a)));
        y2 = OY;
       // printf("4: (%f, %f)\n", (h - b) / (c - b), (h - b) / (d - b));
        break;
    }
    glBegin(GL_LINES);
    glVertex2d(x1, y1);
    glVertex2d(x2, y2);
    glEnd();
}

void draw_opposite(int n, int i, int j, double a, double b, double c, double d){
    double x1 = 0, x2 = 0, y1 = 0, y2 = 0;
    switch (n){
    case 5:
    case 10:	y1 = OY;
        x1 = OX + (dx * ((h - a) / (b - a)));
        y2 = OY - dy;
        x2 = OX + (dx * ((h - c) / (d - c)));
        break;
    case 6:
    case 9:		x2 = OX;
        y1 = OY - (dy * ((h - b) / (d - b)));
        x1 = OX + dx;
        y2 = OY - (dy * ((h - a) / (c - a)));
        break;
    }
    glBegin(GL_LINES);
    glVertex2d(x1, y1);
    glVertex2d(x2, y2);
    glEnd();

}

void lines(int n, int i, int j, double a, double b, double c, double d) // Draw correct line
{
    // Your code goes here
    // n is cell type (0-15)
    // i,j specifes which cell
    // a,b,c,d are function values at cell corners (from data array)
    switch (n){
    case 1:
    case 2:
    case 4:
    case 7:
    case 8:
    case 11:
    case 13:
    case 14:
        draw_one(n, i, j, a, b, c, d);
        break;
    case 5:
    case 6:
    case 9:
    case 10:	draw_opposite(n, i, j, a, b, c, d);
        break;
    case 0:
    case 15:	break;
    default:
        //exit(2);
        break;
    }
}