00001 #include "mcubes.h"
00002
00020 #include <glt/gl.h>
00021
00022 #include <cstdlib>
00023 #include <cmath>
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
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
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
00215
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
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
00252
00253
00254
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
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
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) {
00287 if (sign[0] > 0) {
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) {
00299
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
00393
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
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 }