Commit 4d11ffe5 authored by Taddeüs Kroes's avatar Taddeüs Kroes

- Worked on graohics ass9.

parent 3909dfd7
...@@ -229,18 +229,22 @@ void FillArrayWithIsosurface(void) ...@@ -229,18 +229,22 @@ void FillArrayWithIsosurface(void)
triangle tri[12]; triangle tri[12];
vec3 normal, *p; vec3 normal, *p;
// Loop through cells
for( x = 0; x < nx-1; x++ ) for( x = 0; x < nx-1; x++ )
{ {
for( y = 0; y < ny-1; y++ ) for( y = 0; y < ny-1; y++ )
{ {
for( z = 0; z < nz-1; z++ ) for( z = 0; z < nz-1; z++ )
{ {
// Calculate triangles in cell
n = generate_cell_triangles(tri, get_cell(x, y, z), isovalue); n = generate_cell_triangles(tri, get_cell(x, y, z), isovalue);
tri_count += n; tri_count += n;
for( i = 0; i < n; i++ ) for( i = 0; i < n; i++ )
{ {
p = tri[i].p; p = tri[i].p;
// Calculate the triangle normal and set it to every corner
normal = v3_normalize(v3_crossprod( normal = v3_normalize(v3_crossprod(
v3_subtract(p[1], p[0]), v3_subtract(p[1], p[0]),
v3_subtract(p[2], p[0]) v3_subtract(p[2], p[0])
......
...@@ -12,6 +12,7 @@ ...@@ -12,6 +12,7 @@
#include <stdlib.h> #include <stdlib.h>
#include <stdio.h> #include <stdio.h>
#include <assert.h> #include <assert.h>
#include <math.h>
#include "marching_tetrahedra.h" #include "marching_tetrahedra.h"
/* Compute a linearly interpolated position where an isosurface cuts /* Compute a linearly interpolated position where an isosurface cuts
...@@ -22,10 +23,14 @@ static vec3 ...@@ -22,10 +23,14 @@ static vec3
interpolate_points(unsigned char isovalue, vec3 p1, vec3 p2, interpolate_points(unsigned char isovalue, vec3 p1, vec3 p2,
unsigned char v1, unsigned char v2) unsigned char v1, unsigned char v2)
{ {
/* Initially, simply return the midpoint between p1 and p2. //float diff1 = fabsf((float)(isovalue - v1)),
So no real interpolation is done yet */ // diff2 = fabsf((float)(isovalue - v2)),
// total = diff1 + diff2;
return v3_add(v3_multiply(p1, 0.5), v3_multiply(p2, 0.5)); //vec3 v = v3_add(v3_multiply(p1, diff1/total), v3_multiply(p2, diff2/total));
vec3 v = v3_add(v3_multiply(p1, 0.5), v3_multiply(p2, 0.5));
return v3_create(v.x * sizex, v.y * sizey, v.z * sizez);
} }
/* Using the given iso-value generate triangles for the tetrahedron /* Using the given iso-value generate triangles for the tetrahedron
...@@ -39,71 +44,94 @@ interpolate_points(unsigned char isovalue, vec3 p1, vec3 p2, ...@@ -39,71 +44,94 @@ interpolate_points(unsigned char isovalue, vec3 p1, vec3 p2,
2 triangles. 2 triangles.
*/ */
// Some abbreviating definitions that are used in the function below
#define ADD_VEC(i,a,b) (triangles[tri_count].p[(i)] = \
interpolate_points(isovalue, c.p[(a)], c.p[(b)], v[(a)], v[(b)]))
#define ADD_TRI(a,b,c,d,e,f) \
ADD_VEC(0,(a),(b)); ADD_VEC(1,(c),(d)); ADD_VEC(2,(e),(f)); tri_count++
static int static int
generate_tetrahedron_triangles(triangle *triangles, unsigned char isovalue, generate_tetrahedron_triangles(triangle *triangles, unsigned char isovalue,
cell c, int v0, int v1, int v2, int v3) cell c, int v0, int v1, int v2, int v3)
{ {
unsigned char *v = c.value; unsigned char *v = c.value;
vec3 *p = triangles[0].p; //vec3 *p = triangles[0].p;
int bitsystem = 0; int tri_count = 0;
if(c.value[v0] > isovalue) bitsystem += 1; // Use a 4-bit bitmask to determine the order of "black/white" points of
if(c.value[v1] > isovalue) bitsystem += 2; // the tetrahedon. Generate 0, 1 or 2 triangle(s) according to the case.
if(c.value[v2] > isovalue) bitsystem += 4; switch( (v[v0] <= isovalue ? 1 : 0) | (v[v1] <= isovalue ? 1<<1 : 0)
if(c.value[v3] > isovalue) bitsystem += 8; | (v[v2] <= isovalue ? 1<<2 : 0) | (v[v3] <= isovalue ? 1<<3 : 0) )
switch( bitsystem )
{ {
case 1: case 14: // 0001 case 0x1: case 0xE: // 0001 or 1110
p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]); ADD_TRI(v0, v1, v0, v2, v0, v3);
p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]); break;
p[2] = interpolate_points(isovalue, c.p[v0], c.p[v3], v[v0], v[v3]); //p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]);
return 1; //p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]);
case 2: case 13: // 0010 //p[2] = interpolate_points(isovalue, c.p[v0], c.p[v3], v[v0], v[v3]);
p[0] = interpolate_points(isovalue, c.p[v1], c.p[v0], v[v1], v[v0]); //return 1;
p[1] = interpolate_points(isovalue, c.p[v1], c.p[v3], v[v1], v[v3]); case 0x2: case 0xD: // 0010 or 1101
p[2] = interpolate_points(isovalue, c.p[v1], c.p[v2], v[v1], v[v2]); ADD_TRI(v1, v0, v1, v3, v1, v2);
return 1; break;
case 4: case 11: // 0100 //p[0] = interpolate_points(isovalue, c.p[v1], c.p[v0], v[v1], v[v0]);
p[0] = interpolate_points(isovalue, c.p[v2], c.p[v0], v[v2], v[v0]); //p[1] = interpolate_points(isovalue, c.p[v1], c.p[v3], v[v1], v[v3]);
p[1] = interpolate_points(isovalue, c.p[v2], c.p[v1], v[v2], v[v1]); //p[2] = 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 1;
return 1; case 0x4: case 0xB: // 0100 or 1011
case 8: case 7: // 1000 ADD_TRI(v2, v0, v2, v1, v2, v3);
p[0] = interpolate_points(isovalue, c.p[v3], c.p[v0], v[v3], v[v0]); break;
p[1] = interpolate_points(isovalue, c.p[v3], c.p[v2], v[v3], v[v2]); //p[0] = interpolate_points(isovalue, c.p[v2], c.p[v0], v[v2], v[v0]);
p[2] = interpolate_points(isovalue, c.p[v3], c.p[v1], v[v3], v[v1]); //p[1] = interpolate_points(isovalue, c.p[v2], c.p[v1], v[v2], v[v1]);
return 1; //p[2] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]);
case 5: case 10: // 0101 //return 1;
p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]); case 0x8: case 0x7: // 1000 or 0111
p[1] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]); ADD_TRI(v3, v0, v3, v2, v3, v1);
p[2] = interpolate_points(isovalue, c.p[v0], c.p[v3], v[v0], v[v3]); break;
p = triangles[1].p; //p[0] = interpolate_points(isovalue, c.p[v3], c.p[v0], v[v3], v[v0]);
p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]); //p[1] = interpolate_points(isovalue, c.p[v3], c.p[v2], v[v3], v[v2]);
p[1] = interpolate_points(isovalue, c.p[v1], c.p[v2], v[v1], v[v2]); //p[2] = interpolate_points(isovalue, c.p[v3], c.p[v1], v[v3], v[v1]);
p[2] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]); //return 1;
return 2; case 0x3: case 0xC: // 0011 or 1100
case 3: case 12: // 0011 ADD_TRI(v0, v3, v0, v2, v1, v3);
p[0] = interpolate_points(isovalue, c.p[v0], c.p[v3], v[v0], v[v3]); ADD_TRI(v1, v3, v0, v2, v1, v2);
p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]); break;
p[2] = interpolate_points(isovalue, c.p[v1], c.p[v3], v[v1], v[v3]); //p[0] = interpolate_points(isovalue, c.p[v0], c.p[v3], v[v0], v[v3]);
p = triangles[1].p; //p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]);
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[v3], v[v1], v[v3]);
p[2] = interpolate_points(isovalue, c.p[v1], c.p[v2], v[v1], v[v2]); //p = triangles[1].p;
p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]); //p[0] = interpolate_points(isovalue, c.p[v1], c.p[v3], v[v1], v[v3]);
return 2; //p[2] = interpolate_points(isovalue, c.p[v1], c.p[v2], v[v1], v[v2]);
case 6: case 9: // 0110 //p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]);
p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]); //return 2;
p[2] = interpolate_points(isovalue, c.p[v1], c.p[v3], v[v1], v[v3]); case 0x5: case 0xA: // 0101 or 1010
p[1] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]); ADD_TRI(v0, v1, v2, v3, v0, v3);
p = triangles[1].p; ADD_TRI(v0, v1, v1, v2, v2, v3);
p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]); break;
p[1] = interpolate_points(isovalue, c.p[v0], c.p[v2], v[v0], v[v2]); //p[0] = interpolate_points(isovalue, c.p[v0], c.p[v1], v[v0], v[v1]);
p[2] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]); //p[1] = interpolate_points(isovalue, c.p[v2], c.p[v3], v[v2], v[v3]);
return 2; //p[2] = interpolate_points(isovalue, c.p[v0], c.p[v3], v[v0], v[v3]);
default: // 0000 //p = triangles[1].p;
return 0; //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 0x6: case 0x9: // 0110 or 1001
ADD_TRI(v0, v1, v2, v3, v1, v3);
ADD_TRI(v0, v1, v0, v2, v2, v3);
break;
//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 or 1111
break;
} }
return tri_count;
} }
/* Generate triangles for a single cell for the given iso-value. This function /* Generate triangles for a single cell for the given iso-value. This function
...@@ -118,9 +146,11 @@ generate_tetrahedron_triangles(triangle *triangles, unsigned char isovalue, ...@@ -118,9 +146,11 @@ generate_tetrahedron_triangles(triangle *triangles, unsigned char isovalue,
int int
generate_cell_triangles(triangle *triangles, cell c, unsigned char isovalue) generate_cell_triangles(triangle *triangles, cell c, unsigned char isovalue)
{ {
// Indices to cell corners for each tetrahedon (6 * 4 indices)
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}; 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; int i, tri_count = 0;
// Generate triangles for each tetrahedon
for( i = 0; i < 21; i += 4 ) for( i = 0; i < 21; i += 4 )
{ {
tri_count += generate_tetrahedron_triangles(triangles + tri_count, tri_count += generate_tetrahedron_triangles(triangles + tri_count,
......
...@@ -38,7 +38,7 @@ cell ...@@ -38,7 +38,7 @@ cell
get_cell(int i, int j, int k) get_cell(int i, int j, int k)
{ {
cell c; cell c;
c.p[0] = v3_create(i, j, k); c.p[0] = v3_create(i, j, k);
c.p[1] = v3_create(i+1, j, k); c.p[1] = v3_create(i+1, j, k);
c.p[2] = v3_create(i, j+1, k); c.p[2] = v3_create(i, j+1, k);
...@@ -47,7 +47,7 @@ get_cell(int i, int j, int k) ...@@ -47,7 +47,7 @@ get_cell(int i, int j, int k)
c.p[5] = v3_create(i+1, 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[6] = v3_create(i, j+1, k+1);
c.p[7] = v3_create(i+1, 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[0] = volume[voxel2idx(i, j, k)];
c.value[1] = volume[voxel2idx(i+1, j, k)]; c.value[1] = volume[voxel2idx(i+1, j, k)];
c.value[2] = volume[voxel2idx(i, j+1, k)]; c.value[2] = volume[voxel2idx(i, j+1, k)];
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment