Main Page   Modules   Class Hierarchy   Compound List   File List   Compound Members   File Members   Related Pages   Examples  

glt/mcubes.cpp

Go to the documentation of this file.
00001 #include "mcubes.h"
00002 
00020 #include <glt/gl.h>
00021 
00022 #include <cstdlib>
00023 #include <cmath>   // Only needed for plug-in functions
00024 #include <iosfwd>
00025 using namespace std;
00026 
00027 typedef struct {
00028     float x, y, z;
00029 } Vector3;
00030 
00031 typedef struct {
00032   Vector3 p1, p2;
00033 } Pair;
00034 
00035 typedef struct {
00036     float  value;
00037     Vector3 loc;
00038 } Sample;
00039 
00040 static int x_steps = 20;
00041 static int y_steps = 20;
00042 static int z_steps = 20;
00043 
00044 static GltFunc3d func3d = NULL;
00045 
00046 static Sample **level0 = NULL;
00047 static Sample **level1 = NULL;
00048 
00049 //static FILE *f = NULL;
00050 
00051 static void write_polygon (Vector3 *poly, int size);
00052 static int  allocate_cubes (void);
00053 static void free_cubes (void);
00054 static void sample_level (Sample **level, int ix, Vector3 *vmin, Vector3 *vmax);
00055 static void tesselate_cube (int iy, int iz);
00056 static void analyze_face (Sample *cube, int a, int b, int c, int d,
00057             Pair *line, int *numline);
00058 static void interpolate (Vector3 *v, Sample *a, Sample *b);
00059 static void vect_copy (Vector3 *v1, Vector3 *v2);
00060 static int  vect_equal (Vector3 *v1, Vector3 *v2);
00061 static void vect_add (Vector3 *v1, Vector3 *v2, Vector3 *v3);
00062 static void vect_sub (Vector3 *v1, Vector3 *v2, Vector3 *v3);
00063 static void vect_scale (Vector3 *v1, Vector3 *v2, float k);
00064 
00065 
00066 int GltMarchingCubes(GltFunc3d func,
00067             float minx, float miny, float minz,
00068             float maxx, float maxy, float maxz,
00069             int xsteps, int ysteps, int zsteps)
00070 {
00071     int ix, iy, iz;
00072     Sample **temp;
00073     Vector3 vmin;
00074     Vector3 vmax;
00075 
00076     func3d = func;
00077 
00078     vmin.x = minx;
00079     vmin.y = miny;
00080     vmin.z = minz;
00081 
00082     vmax.x = maxx;
00083     vmax.y = maxy;
00084     vmax.z = maxz;
00085 
00086     x_steps = xsteps;
00087     y_steps = ysteps;
00088     z_steps = zsteps;
00089 
00090     if (!allocate_cubes())
00091         return 0;
00092 
00093     sample_level (level0, 0, &vmin, &vmax);
00094 
00095     for (ix = 1; ix <= x_steps; ix++) {
00096         sample_level (level1, ix, &vmin, &vmax);
00097 
00098         for (iy = 0; iy < y_steps; iy++)
00099             for (iz = 0; iz < z_steps; iz++)
00100                tesselate_cube (iy, iz);
00101 
00102         temp = level0;
00103         level0 = level1;
00104         level1 = temp;
00105     }
00106 
00107     free_cubes();
00108 
00109     return 1;
00110 }
00111 
00112 /* Allocate memory for the marching cubes algorithm */
00113 static int allocate_cubes()
00114 {
00115     int i;
00116 
00117     level0 = (Sample **)malloc ((y_steps+1) * sizeof(Sample *));
00118     if (level0 == NULL)
00119         return 0;
00120 
00121     for (i = 0; i < y_steps+1; i++)
00122         level0[i] = NULL;
00123 
00124     for (i = 0; i < y_steps+1; i++) {
00125         level0[i] = (Sample *)malloc ((z_steps+1) * sizeof(Sample));
00126         if (level0[i] == NULL) {
00127             free_cubes();
00128             return 0;
00129         }
00130     }
00131 
00132     level1 = (Sample **)malloc ((y_steps+1) * sizeof(Sample *));
00133     if (level1 == NULL)
00134         return 0;
00135 
00136     for (i = 0; i < y_steps+1; i++)
00137         level1[i] = NULL;
00138 
00139     for (i = 0; i < y_steps+1; i++) {
00140         level1[i] = (Sample *)malloc ((z_steps+1) * sizeof(Sample));
00141         if (level1[i] == NULL) {
00142             free_cubes();
00143             return 0;
00144         }
00145     }
00146 
00147     return 1;
00148 }
00149 
00150 
00151 static void free_cubes()
00152 {
00153     int i;
00154 
00155     if (level0 != NULL) {
00156         for (i = 0; i < y_steps+1; i++) {
00157             if (level0[i] != NULL)
00158                 free (level0[i]);
00159         }
00160 
00161         free (level0);
00162     }
00163 
00164     if (level1 != NULL) {
00165         for (i = 0; i < y_steps+1; i++) {
00166             if (level1[i] != NULL)
00167                 free (level1[i]);
00168         }
00169 
00170         free (level1);
00171     }
00172 
00173     level0 = NULL;
00174     level1 = NULL;
00175 }
00176 
00177 
00178 static void sample_level (Sample **level, int ix, Vector3 *vmin, Vector3 *vmax)
00179 {
00180     Vector3 v;
00181     int iy, iz;
00182 
00183     v.x = (float) ((float) ix/x_steps * vmin->x + (1.0 - (float) ix/x_steps)*vmax->x);
00184 
00185     for (iy = 0; iy <= y_steps; iy++) {
00186         v.y = (float) ((float) iy/y_steps * vmin->y + (1.0 - (float)iy/y_steps)*vmax->y);
00187 
00188         for (iz = 0; iz <= z_steps; iz++) {
00189             v.z =  (float) ((float) iz/z_steps * vmin->z + (1.0 - (float)iz/z_steps)*vmax->z);
00190 
00191             level[iy][iz].loc   = v;
00192             level[iy][iz].value = func3d (v.x, v.y, v.z);
00193         }
00194     }
00195 }
00196 
00197 
00198 static void tesselate_cube (int iy, int iz)
00199 {
00200     Sample cube[8];
00201     int    i, j, polysize, numline;
00202     Pair   line[12], temp;
00203     Vector3 poly[6];
00204 
00205     cube[0] = level0[iy][iz];
00206     cube[1] = level1[iy][iz];
00207     cube[2] = level0[iy+1][iz];
00208     cube[3] = level1[iy+1][iz];
00209     cube[4] = level0[iy][iz+1];
00210     cube[5] = level1[iy][iz+1];
00211     cube[6] = level0[iy+1][iz+1];
00212     cube[7] = level1[iy+1][iz+1];
00213 
00214     /* Analyze each of the 6 faces of the cube one at a time. For each
00215        face that was intersected save the intersection line */
00216 
00217     numline = 0;
00218 
00219     analyze_face (cube, 0, 1, 3, 2, line, &numline);
00220     analyze_face (cube, 4, 6, 7, 5, line, &numline);
00221     analyze_face (cube, 0, 2, 6, 4, line, &numline);
00222     analyze_face (cube, 1, 5, 7, 3, line, &numline);
00223     analyze_face (cube, 2, 3, 7, 6, line, &numline);
00224     analyze_face (cube, 0, 4, 5, 1, line, &numline);
00225 
00226     if (numline > 0) {
00227         /* Sort the line segments into polygons */
00228         polysize = 0;
00229 
00230         for (i = 0; i < numline; i++) {
00231             poly[polysize++] = line[i].p1;
00232 
00233             for (j = i+1; j < numline; j++) {
00234                 if (vect_equal (&line[j].p1, &line[i].p2)) {
00235                     temp = line[j];
00236                     line[j] = line[i+1];
00237                     line[i+1] = temp;
00238                     break;
00239                 }
00240             }
00241 
00242             if (j >= numline) {
00243                 write_polygon (poly, polysize);
00244                 polysize = 0;
00245             }
00246         }
00247     }
00248 }
00249 
00250 
00251 /* Analyze a cube face for surface intersections. Intersections are detected by
00252    looking for positive to negative field strength changes between adjacent
00253    corners. If an intersection is found then the coords of the line segment(s)
00254    are added to array line and *numline is increased */
00255 static void analyze_face (Sample *cube, int a, int b, int c, int d,
00256                           Pair *line, int *numline)
00257 {
00258     Vector3 points[4], vcenter;
00259     float center;
00260     int sign[4];
00261     int l1, l2, l3, l4;
00262     int index = 0;
00263 
00264     /* Check face edge a-b */
00265     if (cube[a].value * cube[b].value < 0.0) {
00266         sign[index] = cube[b].value > 0.0;
00267         interpolate (&points[index++], &cube[a], &cube[b]);
00268     }
00269 
00270     /* etc. */
00271     if (cube[b].value * cube[c].value < 0.0) {
00272         sign[index] = cube[c].value > 0.0;
00273         interpolate (&points[index++], &cube[b], &cube[c]);
00274     }
00275 
00276     if (cube[c].value * cube[d].value < 0.0) {
00277         sign[index] = cube[d].value > 0.0;
00278         interpolate (&points[index++], &cube[c], &cube[d]);
00279     }
00280 
00281     if (cube[d].value * cube[a].value < 0.0) {
00282         sign[index] = cube[a].value > 0.0;
00283         interpolate (&points[index++], &cube[d], &cube[a]);
00284     }
00285 
00286     if (index == 2) { /* One line segment */
00287         if (sign[0] > 0) {  /* Get the line direction right */
00288             vect_copy (&line[*numline].p1, &points[0]);
00289             vect_copy (&line[*numline].p2, &points[1]);
00290         }
00291         else {
00292             vect_copy (&line[*numline].p1, &points[1]);
00293             vect_copy (&line[*numline].p2, &points[0]);
00294         }
00295 
00296         (*numline)++;
00297     }
00298     else if (index == 4) { /* Two line segments */
00299         /* Sample the center of the face to determine which points to connect */
00300         vect_add (&vcenter, &cube[a].loc, &cube[b].loc);
00301         vect_add (&vcenter, &vcenter, &cube[c].loc);
00302         vect_add (&vcenter, &vcenter, &cube[d].loc);
00303         vect_scale (&vcenter, &vcenter, 0.25F);
00304         center = func3d (vcenter.x, vcenter.y, vcenter.z);
00305 
00306         if ((center > 0.0) != sign[0])
00307             { l1 = 0; l2 = 1; l3 = 2; l4 = 3; }
00308         else
00309             { l1 = 1; l2 = 2; l3 = 3; l4 = 0; }
00310 
00311         if (sign[l1] > 0) {
00312             vect_copy (&line[*numline].p1, &points[l1]);
00313             vect_copy (&line[*numline].p2, &points[l2]);
00314         }
00315         else {
00316             vect_copy (&line[*numline].p1, &points[l2]);
00317             vect_copy (&line[*numline].p2, &points[l1]);
00318         }
00319 
00320         (*numline)++;
00321 
00322         if (sign[l3] > 0) {
00323             vect_copy (&line[*numline].p1, &points[l3]);
00324             vect_copy (&line[*numline].p2, &points[l4]);
00325         }
00326         else {
00327             vect_copy (&line[*numline].p1, &points[l4]);
00328             vect_copy (&line[*numline].p2, &points[l3]);
00329         }
00330 
00331         (*numline)++;
00332     }
00333 }
00334 
00335 
00336 static void interpolate (Vector3 *v, Sample *a, Sample *b)
00337 {
00338     Vector3 vtemp;
00339 
00340     if (a->value < b->value) {
00341         vect_sub (&vtemp, &a->loc, &b->loc);
00342         vect_scale (&vtemp, &vtemp, a->value/(a->value - b->value));
00343         vect_sub (v, &a->loc, &vtemp);
00344     }
00345     else {
00346         vect_sub (&vtemp, &b->loc, &a->loc);
00347         vect_scale (&vtemp, &vtemp, b->value/(b->value - a->value));
00348         vect_sub (v, &b->loc, &vtemp);
00349     }
00350 }
00351 
00352 
00353 static void vect_copy (Vector3 *v1, Vector3 *v2)
00354 {
00355     v1->x = v2->x;
00356     v1->y = v2->y;
00357     v1->z = v2->z;
00358 }
00359 
00360 
00361 static int vect_equal (Vector3 *v1, Vector3 *v2)
00362 {
00363     return (v1->x == v2->x && v1->y == v2->y && v1->z == v2->z);
00364 }
00365 
00366 
00367 static void vect_add (Vector3 *v1, Vector3 *v2, Vector3 *v3)
00368 {
00369     v1->x = v2->x + v3->x;
00370     v1->y = v2->y + v3->y;
00371     v1->z = v2->z + v3->z;
00372 }
00373 
00374 
00375 static void vect_sub (Vector3 *v1, Vector3 *v2, Vector3 *v3)
00376 {
00377     v1->x = v2->x - v3->x;
00378     v1->y = v2->y - v3->y;
00379     v1->z = v2->z - v3->z;
00380 }
00381 
00382 
00383 static void vect_scale (Vector3 *v1, Vector3 *v2, float k)
00384 {
00385     v1->x = k * v2->x;
00386     v1->y = k * v2->y;
00387     v1->z = k * v2->z;
00388 }
00389 
00390 static void write_polygon (Vector3 *poly, int size)
00391 {
00392     // Use the cross-product to find a normal
00393     // Per-vertex lighting would be nicer
00394 
00395     Vector3 normal;
00396     if (size>=3)
00397     {
00398         Vector3 a,b;
00399         vect_sub(&a,poly+0,poly+1);
00400         vect_sub(&b,poly+2,poly+1);
00401         normal.x = a.y*b.z - a.z*b.y;
00402         normal.y = a.z*b.x - a.x*b.z;
00403         normal.z = a.x*b.y - a.y*b.x;
00404     }
00405     
00406     //
00407 
00408     glBegin(GL_POLYGON);
00409         glNormal3f(normal.x,normal.y,normal.z);
00410         for (int i=size-1; i>=0; i--)
00411             glVertex3f(poly[i].x,poly[i].y,poly[i].z);
00412     glEnd();
00413 }
00414 
00415 // http://atrey.karlin.mff.cuni.cz/~0rfelyus/povray/index.html
00416 
00417 template<class T> T sqr(const T &a) { return a*a; }
00418 template<class T> T cub(const T &a) { return a*a*a; }
00419 
00420 float sphere(float x, float y, float z) 
00421 { 
00422     return float( x*x + y*y + z*z - 0.95 ); 
00423 }
00424 
00425 float heart (float x, float y, float z) 
00426 {
00427     x *= 1.3f;
00428     y *= 1.3f;
00429     z *= 1.3f;
00430 
00431     return float( cub(2*sqr(x)+sqr(y)+sqr(z)-1) - (0.1*sqr(x)+sqr(y))*cub(z) );
00432 }
00433 
00434 float klein(float x,float y,float z)
00435 {
00436     x *= 5.0f;
00437     y *= 5.0f;
00438     z *= 5.0f;
00439 
00440     return float( (sqr(x)+sqr(y)+sqr(z)+2*y-1) * (sqr(sqr(x)+sqr(y)+sqr(z)-2*y-1)
00441                   -8*sqr(z)) + 16*x*z*(sqr(x)+sqr(y)+sqr(z)-2*y-1) );
00442 }
00443 
00444 float wave(float x,float y,float z)
00445 {
00446     return float( cos(15*(sqrt(sqr(x)+sqr(y))))*0.3 - z );
00447 }
00448 
00449 float sss(float x,float y,float z)
00450 {
00451     return float( sin(x*6.1) + sin(y*6.1) + sin(z*6.1) - 0.5 );
00452 }

Generated on Tue Nov 5 11:11:03 2002 for GLT by doxygen1.2.18