- Added graphics assignment 9.

graphics/ass9/source/Makefile

+CC=gcc
+
+WARNING_FLAGS=-Wall -Wextra -Werror-implicit-function-declaration -Wshadow -Wstrict-prototypes -pedantic-errors
+CFLAGS=-g -O2 -std=c99 $(WARNING_FLAGS)
+LDFLAGS=-g -lGL -lglut -lGLU
+
+.c.o:
+	$(CC) -c $(CFLAGS) $<
+
+all: main
+main: main.o marching_tetrahedra.o volume.o
+	$(CC) $(LDFLAGS) -o main main.o marching_tetrahedra.o volume.o
+
+clean:
+	rm -f *.o
+	rm -f main
+
+volume.o              : v3math.h volume.h
+volume.o              : volume.h volume.c
+v3math.o              : v3math.h
+marching_tetrahedra.o : marching_tetrahedra.h marching_tetrahedra.c
+marching_tetrahedra.o : v3math.h volume.h marching_tetrahedra.h
+main.o                : marching_tetrahedra.h v3math.h volume.h main.c

graphics/ass9/source/main.c

+/* Computer Graphics, Assignment, Volume rendering with cubes/points/isosurface
+ *
+ * Filename ........ main.c
+ * Description ..... Creates OpenGL window and draws the scene.
+ * Date ............ 29.10.2007
+ * Created by ...... Paul Melis
+ *
+ * Student name ....
+ * Student email ...
+ * Collegekaart ....
+ * Date ............
+ * Comments ........
+ *
+ * (always fill in these fields before submitting!!)
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/time.h>
+#include <GL/glut.h>
+#include <GL/gl.h>
+#include <GL/glu.h>
+#include <math.h>
+
+#include "volume.h"
+#include "marching_tetrahedra.h"
+#include "v3math.h"
+
+int             display_mode=2;            /* 0=points, 1=cubes, 2=isosurface */
+int             use_arrays=1;
+unsigned char   isovalue=128;
+unsigned char   epsilon=0;
+
+int             window;
+int             mouse_mode=0;
+int             mx, my;
+float           camDistance=100.0;
+float           camRotZ=45.0, camAzimuth=20.0;
+float           saved_camRotZ, saved_camAzimuth, saved_camDistance;
+int             regenerate_arrays=1;
+int             entered_number=0;
+int             display_fps_counter=0;
+int             backface_culling=0;
+
+struct timeval  t0, t1;
+int             frames_drawn=0;
+double          time_used=0.0;
+
+#define MAX_VERTICES_IN_ARRAY   5000000
+
+int             num_vertices_in_array=0;
+float           vertices[3*MAX_VERTICES_IN_ARRAY];
+float           normals[3*MAX_VERTICES_IN_ARRAY];
+
+void    SetupCamera(void);
+void    DrawVolumeAxes(void);
+void    DrawVolumeUsingCurrentDisplayMode(void);
+
+void ClearArrays(void)
+{
+    num_vertices_in_array = 0;
+}
+
+void AddVertexToArray(vec3 v, vec3 n)
+{
+    vertices[3*num_vertices_in_array] = v.x;
+    vertices[3*num_vertices_in_array+1] = v.y;
+    vertices[3*num_vertices_in_array+2] = v.z;
+
+    normals[3*num_vertices_in_array] = n.x;
+    normals[3*num_vertices_in_array+1] = n.y;
+    normals[3*num_vertices_in_array+2] = n.z;
+
+    num_vertices_in_array++;
+}
+
+void DrawVolumeAsPoints(void)
+{
+    int i, j, k;
+    int idx;
+
+    glBegin(GL_POINTS);
+
+    idx = 0;
+    for (k = 0; k < nz; k++)
+    {
+        for (j = 0; j < ny; j++)
+        {
+            for (i = 0; i < nx; i++)
+            {
+                if (abs(volume[idx]-isovalue) <= epsilon)
+                    glVertex3f((i+0.5)*sizex, (j+0.5)*sizey, (k+0.5)*sizez);
+
+                idx++;
+            }
+        }
+    }
+
+    glEnd();
+}
+
+void FillArrayWithPoints(void)
+{
+    int i, j, k;
+    int idx;
+
+    idx = 0;
+
+    for (k = 0; k < nz; k++)
+    {
+        for (j = 0; j < ny; j++)
+        {
+            for (i = 0; i < nx; i++)
+            {
+                if (abs(volume[idx]-isovalue) <= epsilon)
+                {
+                    AddVertexToArray(
+                        v3_create((i+0.5)*sizex, (j+0.5)*sizey, (k+0.5)*sizez),
+                        v3_create(0, 0, 0)
+                    );
+                }
+
+                idx++;
+            }
+        }
+    }
+}
+
+
+void DrawVolumeAsCubes(void)
+{
+    int i, j, k;
+    int idx;
+
+    idx = 0;
+    for (k = 0; k < nz; k++)
+    {
+        for (j = 0; j < ny; j++)
+        {
+            for (i = 0; i < nx; i++)
+            {
+                if (abs(volume[idx]-isovalue) <= epsilon)
+                {
+                    glPushMatrix();
+                    glTranslatef((i+0.5)*sizex, (j+0.5)*sizey, (k+0.5)*sizez);
+                    glScalef(sizex, sizey, sizez);
+                    glutSolidCube(1.0);
+                    glPopMatrix();
+                }
+                idx++;
+            }
+        }
+    }
+}
+
+void FillArrayWithCubes(void)
+{
+    int i, j, k;
+    int idx;
+    float   minx, miny, minz;
+    float   maxx, maxy, maxz;
+
+    idx = 0;
+
+    for (k = 0; k < nz; k++)
+    {
+        for (j = 0; j < ny; j++)
+        {
+            for (i = 0; i < nx; i++)
+            {
+                if (abs(volume[idx]-isovalue) <= epsilon)
+                {
+                    minx = i*sizex;
+                    miny = j*sizey;
+                    minz = k*sizez;
+
+                    maxx = (i+1)*sizex;
+                    maxy = (j+1)*sizey;
+                    maxz = (k+1)*sizez;
+
+                    AddVertexToArray(v3_create(minx, miny, minz), v3_create(0, 0, -1));
+                    AddVertexToArray(v3_create(maxx, miny, minz), v3_create(0, 0, -1));
+                    AddVertexToArray(v3_create(maxx, maxy, minz), v3_create(0, 0, -1));
+                    AddVertexToArray(v3_create(minx, maxy, minz), v3_create(0, 0, -1));
+
+                    AddVertexToArray(v3_create(minx, miny, maxz), v3_create(0, 0, 1));
+                    AddVertexToArray(v3_create(maxx, miny, maxz), v3_create(0, 0, 1));
+                    AddVertexToArray(v3_create(maxx, maxy, maxz), v3_create(0, 0, 1));
+                    AddVertexToArray(v3_create(minx, maxy, maxz), v3_create(0, 0, 1));
+
+
+                    AddVertexToArray(v3_create(minx, miny, minz), v3_create(0, -1, 0));
+                    AddVertexToArray(v3_create(maxx, miny, minz), v3_create(0, -1, 0));
+                    AddVertexToArray(v3_create(maxx, miny, maxz), v3_create(0, -1, 0));
+                    AddVertexToArray(v3_create(minx, miny, maxz), v3_create(0, -1, 0));
+
+                    AddVertexToArray(v3_create(minx, maxy, minz), v3_create(0, 1, 0));
+                    AddVertexToArray(v3_create(maxx, maxy, minz), v3_create(0, 1, 0));
+                    AddVertexToArray(v3_create(maxx, maxy, maxz), v3_create(0, 1, 0));
+                    AddVertexToArray(v3_create(minx, maxy, maxz), v3_create(0, 1, 0));
+
+
+                    AddVertexToArray(v3_create(minx, miny, minz), v3_create(-1, 0, 0));
+                    AddVertexToArray(v3_create(minx, maxy, minz), v3_create(-1, 0, 0));
+                    AddVertexToArray(v3_create(minx, maxy, maxz), v3_create(-1, 0, 0));
+                    AddVertexToArray(v3_create(minx, miny, maxz), v3_create(-1, 0, 0));
+
+                    AddVertexToArray(v3_create(maxx, miny, minz), v3_create(1, 0, 0));
+                    AddVertexToArray(v3_create(maxx, maxy, minz), v3_create(1, 0, 0));
+                    AddVertexToArray(v3_create(maxx, maxy, maxz), v3_create(1, 0, 0));
+                    AddVertexToArray(v3_create(maxx, miny, maxz), v3_create(1, 0, 0));
+                }
+
+                idx++;
+            }
+        }
+    }
+}
+
+void DrawVolumeAsIsosurface(void)
+{
+    /* This function can be left empty, we don't use immediate
+    rendering of the isosurface, only rendering using arrays */
+}
+
+void FillArrayWithIsosurface(void)
+{
+	int x, y, z, i, n, tri_count = 0;
+	triangle tri[12];
+	vec3 normal, *p;
+	
+	for( x = 0; x < nx-1; x++ )
+	{
+		for( y = 0; y < ny-1; y++ )
+		{
+			for( z = 0; z < nz-1; z++ )
+			{
+				n = generate_cell_triangles(tri, get_cell(x, y, z), isovalue);
+				tri_count += n;
+
+				for( i = 0; i < n; i++ )
+				{
+					p = tri[i].p;
+					normal = v3_normalize(v3_crossprod(
+						v3_subtract(p[1], p[0]),
+						v3_subtract(p[2], p[0])
+					));
+					
+					AddVertexToArray(p[0], normal);
+					AddVertexToArray(p[1], normal);
+					AddVertexToArray(p[2], normal);
+				}
+			}
+		}
+	}
+	
+	printf("generated %d triangles\n", tri_count);
+}
+
+void DrawScene(void)
+{
+    if (frames_drawn == 0)
+    {
+        /* start timing */
+        glFinish();
+        gettimeofday(&t0, NULL);
+    }
+
+    /* clear the draw buffer */
+    glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);
+
+    /* set the model-view matrix */
+    SetupCamera();
+
+    /* draw an outline of where the volume resides */
+    DrawVolumeAxes();
+
+    /* draw the volume in the chosen mode */
+    DrawVolumeUsingCurrentDisplayMode();
+
+    /* FPS timer code */
+
+    frames_drawn++;
+
+    if (frames_drawn == 30)
+    {
+        glFinish();
+        gettimeofday(&t1, NULL);
+
+        /* calculate time used for 30 frames */
+        double diff = (t1.tv_sec + t1.tv_usec/1000000.0) - (t0.tv_sec + t0.tv_usec/1000000.0);
+
+        if (display_fps_counter)
+            printf("%.1f fps (%f secs for %d frames)\n", frames_drawn/diff, diff, frames_drawn);
+
+        frames_drawn = 0;
+    }
+
+    /* finally, swap the draw buffers to make them appear on screen */
+    glutSwapBuffers();
+}
+
+/* Draw a set of colored X, Y and Z axes that match the size (and location)
+   of the volume, so we can see where output should appear */
+void
+DrawVolumeAxes(void)
+{
+    glPushAttrib(GL_LIGHTING_BIT);
+    glDisable(GL_LIGHTING);
+
+    glBegin(GL_LINES);
+
+    glColor3f(1, 0, 0);
+    glVertex3f(0, 0, 0);
+    glVertex3f(nx*sizex, 0, 0);
+
+    glColor3f(0, 1, 0);
+    glVertex3f(0, 0, 0);
+    glVertex3f(0, ny*sizey, 0);
+
+    glColor3f(0, 0, 1);
+    glVertex3f(0, 0, 0);
+    glVertex3f(0, 0, nz*sizez);
+
+    glEnd();
+
+    glPopAttrib();
+}
+
+void DrawVolumeUsingCurrentDisplayMode(void)
+{
+    if (use_arrays)
+    {
+        if (regenerate_arrays == 1)
+        {
+            printf("Filling vertex/normal arrays... ");
+            fflush(stdout);
+
+            ClearArrays();
+
+            switch (display_mode)
+            {
+            case 0:
+                FillArrayWithPoints();
+                break;
+            case 1:
+                FillArrayWithCubes();
+                break;
+            case 2:
+                FillArrayWithIsosurface();
+                break;
+            }
+
+            printf("done, %d array vertices used (1 per point, 4 per quad, 3 per triangle)\n", num_vertices_in_array);
+
+            regenerate_arrays = 0;
+        }
+
+        glEnableClientState(GL_VERTEX_ARRAY);
+        glVertexPointer(3, GL_FLOAT, 0, vertices);
+
+        if (display_mode > 0)
+        {
+            glEnableClientState(GL_NORMAL_ARRAY);
+            glNormalPointer(GL_FLOAT, 0, normals);
+        }
+    }
+
+    glColor3f(0.6, 0.0, 0.0);
+
+    if (display_mode == 0)
+        glDisable(GL_LIGHTING);
+    else
+        glEnable(GL_LIGHTING);
+
+    if (display_mode == 0)
+    {
+        /* points */
+        glPointSize(8.0);
+        if (use_arrays)
+            glDrawArrays(GL_POINTS, 0, num_vertices_in_array);
+        else
+            DrawVolumeAsPoints();
+    }
+    else if (display_mode == 1)
+    {
+        /* cubes */
+        if (use_arrays)
+            glDrawArrays(GL_QUADS, 0, num_vertices_in_array);
+        else
+        {
+            glEnable(GL_RESCALE_NORMAL);
+            DrawVolumeAsCubes();
+            glDisable(GL_RESCALE_NORMAL);
+        }
+    }
+    else if (display_mode == 2)
+    {
+        /* isosurface */
+        if (use_arrays)
+            glDrawArrays(GL_TRIANGLES, 0, num_vertices_in_array);
+        else
+            DrawVolumeAsIsosurface();
+    }
+
+    if (use_arrays)
+    {
+        glDisableClientState(GL_VERTEX_ARRAY);
+        if (display_mode > 0)
+            glDisableClientState(GL_NORMAL_ARRAY);
+    }
+}
+
+void SetupCamera(void)
+{
+    GLfloat light_position[4];
+
+    glLoadIdentity();
+
+    // Verbose, but straightforward way, of positioning the camera and lightsource.
+    // Assume the camera's final position is (cx, cy, cz).
+    // Start with c being (camDistance, 0, 0)
+    // First rotate around Y, then around Z.
+    // Now we have c at the given distance from the origin, with specified rotation angles.
+
+    float	cx, cy, cz;
+    float	t;
+    float 	beta, gamma;
+
+    // degrees -> radians
+    beta = camAzimuth / 180.0 * 3.1415926535;
+    gamma = camRotZ / 180.0 * 3.1415926535;
+
+    cx = camDistance;
+    cy = cz = 0.0;
+
+    // Rotate around Y
+    t = cx;
+    cx = cx * cos(beta) - cz * sin(beta);
+    cz = t * sin(beta) + cz * cos(beta);
+
+    // Rotate around Z
+    t = cx;
+    cx = cx * cos(gamma) + cy * sin(gamma);
+    cy = -t * sin(gamma) + cy * cos(gamma);
+
+    gluLookAt (cx, cy, cz,  0.0, 0.0, 0.0,  0.0, 0.0, 1.0);
+
+    // Finally, translate because the lookat point is at the center of the volume.
+    glTranslatef(-0.5*nx*sizex, -0.5*ny*sizey, -0.5*nz*sizez);
+
+    // Place the light source at the camera position (head light)
+    light_position[0] = cx + 0.5*nx*sizex;
+    light_position[1] = cy + 0.5*ny*sizey;
+    light_position[2] = cz + 0.5*nz*sizez;
+    light_position[3] = 1.0;
+    glLightfv(GL_LIGHT0, GL_POSITION, light_position);
+}
+
+void InitOpenGL(void)
+{
+    GLfloat light_ambient[] = { 0.4, 0.4, 0.4, 0.0 };
+    GLfloat	red[] = { 1.0, 0.0, 0.0, 1.0 };
+    GLfloat mat_no_specular[] = { 0.0, 0.0, 0.0, 1.0 };
+
+    glClearColor(0.7, 0.7, 1, 1);
+
+    glEnable(GL_COLOR_MATERIAL);
+
+    glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
+    glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat_no_specular);
+
+    glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
+    glEnable(GL_LIGHTING);
+    glEnable(GL_LIGHT0);
+    glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1);
+    glLightModeli(GL_LIGHT_MODEL_LOCAL_VIEWER, 1);
+
+    glDepthFunc(GL_LEQUAL);
+    glEnable(GL_DEPTH_TEST);
+
+    glDisable(GL_CULL_FACE);
+
+    glShadeModel(GL_FLAT);
+    //glShadeModel(GL_SMOOTH);
+}
+
+void ReSizeScene(int Width, int Height)
+{
+    if (Height==0)
+        Height=1;
+
+    glViewport(0, 0, Width, Height);
+
+    glMatrixMode(GL_PROJECTION);
+
+    glLoadIdentity();
+    gluPerspective(45.0f,(GLfloat)Width/(GLfloat)Height, 0.1f, 10000.0f);
+
+    glMatrixMode(GL_MODELVIEW);
+    glLoadIdentity ();
+}
+
+void keyPressed(unsigned char key, int x, int y)
+{
+    int	redraw=0;
+
+    /* keep gcc happy, as we don't use these parameters */
+    (void)x;
+    (void)y;
+
+    switch (key)
+    {
+        case 'q':
+            exit(-1);
+            break;
+
+        case 'b':
+            backface_culling = 1 - backface_culling;
+            if (backface_culling)
+            {
+                glEnable(GL_CULL_FACE);
+                printf("Backface culling is now enabled\n");
+            }
+            else
+            {
+                glDisable(GL_CULL_FACE);
+                printf("Backface culling is now disabled\n");
+            }
+            glutPostRedisplay();
+            break;
+
+        case 'p':
+            /* switch to points mode */
+            display_mode = 0;
+            redraw = 1;
+            break;
+        case 'c':
+            /* switch to cubes mode */
+            display_mode = 1;
+            redraw = 1;
+            break;
+        case 's':
+            /* switch to isosurface mode */
+            display_mode = 2;
+            redraw = 1;
+            break;
+
+        case '[':
+            /* decrease isovalue */
+            {
+                int newvalue = isovalue-10;
+                if (newvalue < 0)
+                    isovalue = 0;
+                else
+                    isovalue = newvalue;
+            }
+            printf("Setting isovalue to %d\n", isovalue);
+            redraw = 1;
+            break;
+        case ']':
+            /* increase isovalue */
+            {
+                int newvalue = isovalue+10;
+                if (newvalue > 255)
+                    isovalue = 255;
+                else
+                    isovalue = newvalue;
+            }
+            printf("Setting isovalue to %d\n", isovalue);
+            redraw = 1;
+            break;
+
+        /* handle isovalue/epsilon entered by number */
+        case '0':
+            entered_number = entered_number*10;
+            break;
+        case '1':
+            entered_number = entered_number*10 + 1;
+            break;
+        case '2':
+            entered_number = entered_number*10 + 2;
+            break;
+        case '3':
+            entered_number = entered_number*10 + 3;
+            break;
+        case '4':
+            entered_number = entered_number*10 + 4;
+            break;
+        case '5':
+            entered_number = entered_number*10 + 5;
+            break;
+        case '6':
+            entered_number = entered_number*10 + 6;
+            break;
+        case '7':
+            entered_number = entered_number*10 + 7;
+            break;
+        case '8':
+            entered_number = entered_number*10 + 8;
+            break;
+        case '9':
+            entered_number = entered_number*10 + 9;
+            break;
+
+        case 'i':
+            /* update isovalue with number entered */
+            isovalue = entered_number % 256;
+            entered_number = 0;
+
+            printf("Setting isovalue to %d\n", isovalue);
+            redraw = 1;
+            break;
+
+        case '{':
+            /* decrease epsilon */
+            {
+                int newvalue = epsilon - 2;
+                if (newvalue < 0)
+                    epsilon = 0;
+                else
+                    epsilon = newvalue;
+            }
+            printf("Setting epsilon to %d\n", epsilon);
+            redraw = 1;
+            break;
+        case '}':
+            /* increase epsilon */
+            {
+                int newvalue = epsilon + 2;
+                if (newvalue > 255)
+                    epsilon = 255;
+                else
+                    epsilon = newvalue;
+            }
+            printf("Setting epsilon to %d\n", epsilon);
+            redraw = 1;
+            break;
+
+        case 'e':
+            /* update epsilon with number entered */
+            epsilon = entered_number % 256;
+            entered_number = 0;
+
+            printf("Setting epsilon to %d\n", epsilon);
+            redraw = 1;
+            break;
+
+
+        case 'f':
+            /* enable/display FPS counter display */
+            display_fps_counter = 1 - display_fps_counter;
+            break;
+        case 'a':
+            use_arrays = 1 - use_arrays;
+            if (use_arrays)
+                printf("Now using arrays\n");
+            else
+                printf("Stopped using arrays\n");
+            glutPostRedisplay();
+            break;
+    }
+
+    if (redraw)
+    {
+        regenerate_arrays = 1;
+        glutPostRedisplay();
+    }
+}
+
+void specialKeyPressed(int key, int x, int y)
+{
+    /* keep gcc happy, as we don't use these parameters */
+    (void)key;
+    (void)x;
+    (void)y;
+}
+
+static void
+mouseFunc(int button, int state, int x, int y)
+{
+    // guard against both left and right buttons being pressed at the same time,
+    // by only responding when a mouse button is pressed while another one
+    // hasn't been pressed yet
+    if (state == GLUT_DOWN && mouse_mode == 0)
+    {
+        if (button == GLUT_LEFT_BUTTON)
+        {
+            mouse_mode = GLUT_LEFT_BUTTON;
+
+            saved_camRotZ = camRotZ;
+            saved_camAzimuth = camAzimuth;
+
+            mx = x;
+            my = y;
+        }
+        else if (button == GLUT_RIGHT_BUTTON)
+        {
+            mouse_mode = GLUT_RIGHT_BUTTON;
+
+            saved_camDistance = camDistance;
+
+            my = y;
+        }
+    }
+    else if (state == GLUT_UP && button == mouse_mode)
+    {
+        // pressed button released
+        mouse_mode = 0;
+    }
+
+}
+
+static void
+motionFunc(int x, int y)
+{
+    int	dx, dy;
+
+    if (mouse_mode == GLUT_LEFT_BUTTON)
+    {
+        dx = mx - x;
+        dy = my - y;
+
+        camRotZ = saved_camRotZ - dx * 0.25;
+        camAzimuth = saved_camAzimuth - dy * 0.25;
+
+        if (camAzimuth > 89.99)
+            camAzimuth = 89.99;
+        else if (camAzimuth < -89.99)
+            camAzimuth = -89.99;
+    }
+    else if (mouse_mode == GLUT_RIGHT_BUTTON)
+    {
+        dy = my - y;
+
+        camDistance = saved_camDistance - dy * sizex;
+        if (camDistance < 0.01)
+            camDistance = 0.01;
+    }
+}
+
+static void
+initialize_volume(const char *fname, unsigned char iso)
+{
+    read_volume(fname);
+
+    isovalue = iso;
+    if (epsilon > 0)
+        printf("Iso-value = %d (+/- %d)\n", isovalue, epsilon);
+    else
+        printf("Iso-value = %d\n", isovalue);
+
+    // Initial distance from camera position to center of volume
+    camDistance = 2.0 * nx * sizex;
+}
+
+int
+main(int argc, char **argv)
+{
+    glutInit(&argc, argv);
+
+    if (argc != 3)
+    {
+        printf("usage: %s datafile isovalue\n", argv[0]);
+        exit(0);
+    }
+
+    initialize_volume(argv[1], atoi(argv[2]));
+
+    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH);
+    glutInitWindowSize(640, 480);
+    glutInitWindowPosition(0, 0);
+    window = glutCreateWindow("OpenGL Framework");
+
+    glutDisplayFunc(&DrawScene);
+    glutIdleFunc(&DrawScene);
+    glutReshapeFunc(&ReSizeScene);
+    glutSpecialFunc(&specialKeyPressed);
+    glutKeyboardFunc(&keyPressed);
+    glutMouseFunc(&mouseFunc);
+    glutMotionFunc(&motionFunc);
+
+    InitOpenGL();
+
+    glutMainLoop();
+
+    return 1;
+}
+
+
+
+
+

graphics/ass9/source/marching_tetrahedra.c

+/* Computer Graphics, Assignment, Volume rendering with cubes/points/isosurface
+ *
+ * Student name ....
+ * Student email ...
+ * Collegekaart ....
+ * Date ............
+ * Comments ........
+ *
+ * (always fill in these fields before submitting!!)
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include "marching_tetrahedra.h"
+
+/* Compute a linearly interpolated position where an isosurface cuts
+   an edge between two vertices (p1 and p2), each with their own
+   scalar value (v1 and v2) */
+
+static vec3
+interpolate_points(unsigned char isovalue, vec3 p1, vec3 p2,
+	unsigned char v1, unsigned char v2)
+{
+    /* Initially, simply return the midpoint between p1 and p2.
+       So no real interpolation is done yet */
+
+    return v3_add(v3_multiply(p1, 0.5), v3_multiply(p2, 0.5));
+}
+
+/* Using the given iso-value generate triangles for the tetrahedron
+   defined by corner vertices v0, v1, v2, v3 of cell c.
+
+   Store the resulting triangles in the "triangles" array.
+
+   Return the number of triangles created (either 0, 1, or 2).
+
+   Note: the output array "triangles" should have space for at least
+         2 triangles.
+*/
+
+static int
+generate_tetrahedron_triangles(triangle *triangles, unsigned char isovalue,
+	cell c, int v0, int v1, int v2, int v3)
+{
+    unsigned char *v = c.value;
+	vec3 *p = triangles[0].p;
+	int bitsystem = 0;
+	
+	if(c.value[v0] > isovalue) bitsystem += 1;
+    if(c.value[v1] > isovalue) bitsystem += 2;
+    if(c.value[v2] > isovalue) bitsystem += 4;
+    if(c.value[v3] > isovalue) bitsystem += 8;
+	
+	switch( bitsystem )
+	{
+	case 1: case 14: // 0001
+		p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]);
+		p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]);
+		p[2] = interpolate_points(isovalue, c.p[v0], c.p[v3], v[v0], v[v3]);
+		return 1;
+	case 2: case 13: // 0010
+		p[0] = interpolate_points(isovalue, c.p[v1], c.p[v0], v[v1], v[v0]);
+		p[1] = interpolate_points(isovalue, c.p[v1], c.p[v3], v[v1], v[v3]);
+		p[2] = interpolate_points(isovalue, c.p[v1], c.p[v2], v[v1], v[v2]);
+		return 1;
+	case 4: case 11: // 0100
+		p[0] = interpolate_points(isovalue, c.p[v2], c.p[v0], v[v2], v[v0]);
+		p[1] = interpolate_points(isovalue, c.p[v2], c.p[v1], v[v2], v[v1]);
+		p[2] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]);
+		return 1;
+	case 8: case 7: // 1000
+		p[0] = interpolate_points(isovalue, c.p[v3], c.p[v0], v[v3], v[v0]);
+		p[1] = interpolate_points(isovalue, c.p[v3], c.p[v2], v[v3], v[v2]);
+		p[2] = interpolate_points(isovalue, c.p[v3], c.p[v1], v[v3], v[v1]);
+		return 1;
+	case 5: case 10: // 0101
+		p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]);
+		p[1] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]);
+		p[2] = interpolate_points(isovalue, c.p[v0], c.p[v3], v[v0], v[v3]);
+		p = triangles[1].p;
+		p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]);
+		p[1] = interpolate_points(isovalue, c.p[v1], c.p[v2], v[v1], v[v2]);
+		p[2] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]);
+		return 2;
+	case 3: case 12: // 0011
+		p[0] = interpolate_points(isovalue, c.p[v0], c.p[v3], v[v0], v[v3]);
+		p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]);
+		p[2] = interpolate_points(isovalue, c.p[v1], c.p[v3], v[v1], v[v3]);
+		p = triangles[1].p;
+		p[0] = interpolate_points(isovalue, c.p[v1], c.p[v3], v[v1], v[v3]);
+		p[2] = interpolate_points(isovalue, c.p[v1], c.p[v2], v[v1], v[v2]);
+		p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]);
+		return 2;
+	case 6: case 9: // 0110
+		p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]);
+		p[2] = interpolate_points(isovalue, c.p[v1], c.p[v3], v[v1], v[v3]);
+		p[1] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]);
+		p = triangles[1].p;
+		p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]);
+		p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]);
+		p[2] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]);
+		return 2;
+	default: // 0000
+		return 0;
+	}
+}
+
+/* Generate triangles for a single cell for the given iso-value. This function
+   should produce at most 6 * 2 triangles (for which the "triangles" array
+   should have enough space).
+
+   Use calls to generate_tetrahedron_triangles().
+
+   Return the number of triangles produced
+*/
+
+int
+generate_cell_triangles(triangle *triangles, cell c, unsigned char isovalue)
+{
+    const int points[24] = {0,1,3,7,0,2,6,7,0,1,5,7,0,2,3,7,0,4,5,7,0,4,6,7};
+	int i, tri_count = 0;
+
+	for( i = 0; i < 21; i += 4 )
+	{
+		tri_count += generate_tetrahedron_triangles(triangles + tri_count,
+			isovalue, c, points[i], points[i+1], points[i+2], points[i+3]);
+	}
+
+    return tri_count;
+}

graphics/ass9/source/marching_tetrahedra.h

+#ifndef MARCHING_TETRAHEDRA_H
+#define MARCHING_TETRAHEDRA_H
+
+#include "v3math.h"
+#include "volume.h"
+
+typedef struct
+{
+    vec3    p[3];
+    vec3    n[3];
+}
+triangle;
+
+int generate_cell_triangles(triangle *triangles, cell c, unsigned char isovalue);
+
+#endif

graphics/ass9/source/v3math.h

+#ifndef V3MATH_H
+#define V3MATH_H
+
+#include <math.h>
+
+typedef struct
+{
+    float   x, y, z;
+}
+vec3;
+
+// Create a new 3-vector of floats
+static inline vec3
+v3_create(float x, float y, float z)
+{
+    vec3 res;
+    res.x = x;
+    res.y = y;
+    res.z = z;
+    return res;
+}
+
+// Return -a
+static inline vec3
+v3_negate(vec3 a)
+{
+    vec3 res;
+    res.x = - a.x;
+    res.y = - a.y;
+    res.z = - a.z;
+    return res;
+}
+
+// Return a + b
+static inline vec3
+v3_add(vec3 a,vec3 b)
+{
+    vec3 res;
+    res.x = a.x+b.x;
+    res.y = a.y+b.y;
+    res.z = a.z+b.z;
+    return res;
+}
+
+// Return a - b
+static inline vec3
+v3_subtract(vec3 a, vec3 b)
+{
+    vec3 res;
+    res.x = a.x-b.x;
+    res.y = a.y-b.y;
+    res.z = a.z-b.z;
+    return res;
+}
+
+// Return a / |a|
+static inline vec3
+v3_normalize(vec3 a)
+{
+    vec3 res;
+    double l = sqrt(a.x*a.x + a.y*a.y + a.z*a.z);
+
+    res = a;
+    res.x /= l;
+    res.y /= l;
+    res.z /= l;
+
+    return res;
+}
+
+// Return a ^ b
+static inline vec3
+v3_crossprod(vec3 a, vec3 b)
+{
+    vec3 res;
+    res.x = a.y*b.z - a.z*b.y;
+    res.y = a.z*b.x - a.x*b.z;
+    res.z = a.x*b.y - a.y*b.x;
+    return res;
+}
+
+// Return a * b
+static inline float
+v3_dotprod(vec3 a, vec3 b)
+{
+    return a.x*b.x + a.y*b.y + a.z*b.z;
+}
+
+// Return a*s
+static inline vec3
+v3_multiply(vec3 a, float s)
+{
+    vec3 res;
+    res.x = a.x*s;
+    res.y = a.y*s;
+    res.z = a.z*s;
+    return res;
+}
+
+// Return |a|
+static inline float
+v3_length(vec3 a)
+{
+    return sqrt(a.x*a.x + a.y*a.y + a.z*a.z);
+}
+
+// Return the i-th component of a, i.e. for i=0
+// this returns a.x
+static inline float
+v3_component(vec3 a, int i)
+{
+    if (i == 0)
+        return a.x;
+    else if (i == 1)
+        return a.y;
+    else
+        return a.z;
+}
+
+// Set the i-th component of a to v
+static inline vec3
+v3_set_component(vec3 a, int i, float v)
+{
+    vec3    res = a;
+
+    if (i == 0)
+        res.x = v;
+    else if (i == 1)
+        res.y = v;
+    else
+        res.z = v;
+
+    return res;
+}
+
+#endif

graphics/ass9/source/volume.c

+/* Computer Graphics, Assignment, Volume rendering with cubes/points/isosurface
+ *
+ * Student name ....
+ * Student email ...
+ * Collegekaart ....
+ * Date ............
+ * Comments ........
+ *
+ * (always fill in these fields before submitting!!)
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "volume.h"
+
+/* The voxels of the volume dataset, stored as a one-dimensional array */
+unsigned char   *volume;
+
+/* The dimensions of the volume dataset */
+int     nx, ny, nz;
+
+/* The size of a voxel */
+float   sizex, sizey, sizez;
+
+/* Utility function to convert the index of a voxel
+   into an index in the volume array above */
+int
+voxel2idx(int i, int j, int k)
+{
+    return (k*ny + j)*nx + i;
+}
+
+/* Extract a cell from the volume, so that datapoint 0 of the
+   cell corresponds to voxel (i, j, k), datapoint 1 to voxel (i+1, j, k),
+   etc. See Figure 3 in the assignment. */
+cell
+get_cell(int i, int j, int k)
+{
+    cell c;
+
+	c.p[0] = v3_create(i, j, k);
+	c.p[1] = v3_create(i+1, j, k);
+	c.p[2] = v3_create(i, j+1, k);
+	c.p[3] = v3_create(i+1, j+1, k);
+	c.p[4] = v3_create(i, j, k+1);
+	c.p[5] = v3_create(i+1, j, k+1);
+	c.p[6] = v3_create(i, j+1, k+1);
+	c.p[7] = v3_create(i+1, j+1, k+1);
+
+	c.value[0] = volume[voxel2idx(i, j, k)];
+	c.value[1] = volume[voxel2idx(i+1, j, k)];
+	c.value[2] = volume[voxel2idx(i, j+1, k)];
+	c.value[3] = volume[voxel2idx(i+1, j+1, k)];
+	c.value[4] = volume[voxel2idx(i, j, k+1)];
+	c.value[5] = volume[voxel2idx(i+1, j, k+1)];
+	c.value[6] = volume[voxel2idx(i, j+1, k+1)];
+	c.value[7] = volume[voxel2idx(i+1, j+1, k+1)];
+	
+    return c;
+}
+
+/* Utility function to read a volume dataset from a VTK file.
+   This will store the data in the "volume" array and update the dimension
+   and size values. */
+
+void
+read_volume(const char *fname)
+{
+    FILE *f;
+    char s[256];
+    int nvoxels;
+
+    printf("Reading %s\n", fname);
+    f = fopen(fname, "rb");
+
+    if (!f)
+    {
+        fprintf(stderr, "read_volume(): Could not open file '%s' for reading!\n", fname);
+        exit(-1);
+    }
+
+    // header line
+    fgets(s, 255, f);
+
+    // comment line
+    fgets(s, 255, f);
+
+    // BINARY
+    fgets(s, 255, f);
+
+    // DATASET STRUCTURED_POINTS
+    fgets(s, 255, f);
+
+    // DIMENSIONS %d %d %d
+    fscanf(f, "%s %d %d %d\n", s, &nx, &ny, &nz);
+    printf("%d x %d x %d voxels\n", nx, ny, nz);
+
+    // ASPECT_RATIO/SPACING %f %f %f
+    fscanf(f, "%s %f %f %f\n", s, &sizex, &sizey, &sizez);
+    printf("voxel sizes: %.3f, %.3f, %.3f\n", sizex, sizey, sizez);
+
+    // ORIGIN ...
+    fgets(s, 255, f);
+
+    // POINT_DATA ...
+    fgets(s, 255, f);
+
+    // SCALARS ...
+    fgets(s, 255, f);
+
+    // LOOKUP_TABLE ...
+    fgets(s, 255, f);
+
+    // allocate memory to hold the volume data and read it from file
+    nvoxels = nx * ny * nz;
+    volume = (unsigned char*)malloc(nvoxels);
+
+    if (fread(volume, 1, nvoxels, f) < (size_t)nvoxels)
+    {
+        printf("WARNING: not all data could be read!\n");
+    }
+
+    fclose(f);
+}

graphics/ass9/source/volume.h

+#ifndef VOLUME_H
+#define VOLUME_H
+
+#include "v3math.h"
+
+/* A cell in the volume dataset, consisting of 8 neighbouring datapoints;
+   p[] contains the corner positions, value[] the associated scalar values */
+typedef struct
+{
+    vec3    p[8];
+    vec3    n[8];       // Note: we experiment with normals in the solution
+    unsigned char  value[8];
+}
+cell;
+
+/* The data points in the volume dataset, stored as a one-dimensional array */
+extern unsigned char    *volume;
+
+/* The dimensions of the volume dataset in number of voxels in each
+   dimension*/
+extern int                nx, ny, nz;
+
+/* The size of a voxel for each dimension */
+extern float              sizex, sizey, sizez;
+
+/* Utility function to convert the index of a datapoint
+   into an index in the volume array above */
+int voxel2idx(int i, int j, int k);
+
+/* Extract a cell from the volume, so that datapoint 0 of the
+   cell corresponds to voxel (i, j, k) */
+cell get_cell(int i, int j, int k);
+
+/* Utility function to read a volume dataset from a file.
+   This will store the data in the "volume" array and update the dimension
+   and size values. */
+void read_volume(const char *fname);
+
+#endif