src/edu/berkeley/qfat/voxel/MarchingCubes.java

   1 package edu.berkeley.qfat.voxel;
   2 import edu.berkeley.qfat.*;
   3 import edu.berkeley.qfat.geom.*;
   4
   5 /////////////////////////////////////////////////////////////////////////////////////
   6 // Derived from:                                                                   //
   7 // http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/marchingsource.cpp    //
   8 // (which is public domain)                                                        //
   9 /////////////////////////////////////////////////////////////////////////////////////
  10 //
  11 // Marching Cubes Example Program
  12 // by Cory Bloyd (corysama@yahoo.com)
  13 //
  14 // A simple, portable and complete implementation of the Marching Cubes
  15 // and Marching Tetrahedrons algorithms in a single source file.
  16 // There are many ways that this code could be made faster, but the
  17 // intent is for the code to be easy to understand.
  18 //
  19 // For a description of the algorithm go to
  20 // http://astronomy.swin.edu.au/pbourke/modelling/polygonise/
  21 //
  22 // This code is public domain.
  23 //
  24 //////////////////////////////////////////////////////////////////////////////////////
  25
  26 public abstract class MarchingCubes {
  27
  28     private int iEdgeFlags;
  29
  30     private static class GLvector {
  31         public double fX;
  32         public double fY;
  33         public double fZ;
  34     };
  35
  36     // These tables are used so that everything can be done in little
  37     // loops that you can look at all at once rather than in pages and
  38     // pages of unrolled code.
  39
  40     // a2fVertexOffset lists the positions, relative to vertex0, of
  41     // each of the 8 vertices of a cube
  42     static final double a2fVertexOffset[][] = {
  43         {0.0, 0.0, 0.0},{1.0, 0.0, 0.0},{1.0, 1.0, 0.0},{0.0, 1.0, 0.0},
  44         {0.0, 0.0, 1.0},{1.0, 0.0, 1.0},{1.0, 1.0, 1.0},{0.0, 1.0, 1.0}
  45     };
  46
  47     //a2iEdgeConnection lists the index of the endpoint vertices for each of the 12 edges of the cube
  48     static final int a2iEdgeConnection[][] = {
  49         {0,1}, {1,2}, {2,3}, {3,0},
  50         {4,5}, {5,6}, {6,7}, {7,4},
  51         {0,4}, {1,5}, {2,6}, {3,7}
  52     };
  53
  54     //a2fEdgeDirection lists the direction vector (vertex1-vertex0) for each edge in the cube
  55     static final double a2fEdgeDirection[][] = {
  56         {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0},
  57         {1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{-1.0, 0.0, 0.0},{0.0, -1.0, 0.0},
  58         {0.0, 0.0, 1.0},{0.0, 0.0, 1.0},{ 0.0, 0.0, 1.0},{0.0,  0.0, 1.0}
  59     };
  60
  61     // a2iTetrahedronEdgeConnection lists the index of the endpoint
  62     // vertices for each of the 6 edges of the tetrahedron
  63     static final int a2iTetrahedronEdgeConnection[][] = {
  64         {0,1},  {1,2},  {2,0},  {0,3},  {1,3},  {2,3}
  65     };
  66
  67     // a2iTetrahedronEdgeConnection lists the index of verticies from a cube
  68     // that made up each of the six tetrahedrons within the cube
  69     static final int a2iTetrahedronsInACube[][] = {
  70         {0,5,1,6},
  71         {0,1,2,6},
  72         {0,2,3,6},
  73         {0,3,7,6},
  74         {0,7,4,6},
  75         {0,4,5,6},
  76     };
  77
  78     static final double afAmbientWhite [] = {0.25, 0.25, 0.25, 1.00};
  79     static final double afAmbientRed   [] = {0.25, 0.00, 0.00, 1.00};
  80     static final double afAmbientGreen [] = {0.00, 0.25, 0.00, 1.00};
  81     static final double afAmbientBlue  [] = {0.00, 0.00, 0.25, 1.00};
  82     static final double afDiffuseWhite [] = {0.75, 0.75, 0.75, 1.00};
  83     static final double afDiffuseRed   [] = {0.75, 0.00, 0.00, 1.00};
  84     static final double afDiffuseGreen [] = {0.00, 0.75, 0.00, 1.00};
  85     static final double afDiffuseBlue  [] = {0.00, 0.00, 0.75, 1.00};
  86     static final double afSpecularWhite[] = {1.00, 1.00, 1.00, 1.00};
  87     static final double afSpecularRed  [] = {1.00, 0.25, 0.25, 1.00};
  88     static final double afSpecularGreen[] = {0.25, 1.00, 0.25, 1.00};
  89     static final double afSpecularBlue [] = {0.25, 0.25, 1.00, 1.00};
  90
  91     GLvector  sSourcePoint[] = new GLvector[3];
  92
  93     /*
  94       GLenum    ePolygonMode = GL_FILL;
  95       int     iDataSetSize = 16;
  96       double   fStepSize = 1.0/iDataSetSize;
  97       double   fTargetValue = 48.0;
  98       double   fTime = 0.0;
  99       boolean bSpin = true;
 100       boolean bMove = true;
 101       boolean bLight = true;
 102
 103
 104       void vIdle();
 105       void vDrawScene();
 106       void vResize(GLsizei, GLsizei);
 107       void vKeyboard(unsigned char cKey, int iX, int iY);
 108       void vSpecial(int iKey, int iX, int iY);
 109
 110       GLvoid vPrintHelp();
 111       GLvoid vSetTime(double fTime);
 112       double fSample1(double fX, double fY, double fZ);
 113       double fSample2(double fX, double fY, double fZ);
 114       double fSample3(double fX, double fY, double fZ);
 115       double (*fSample)(double fX, double fY, double fZ) = fSample1;
 116
 117       GLvoid vMarchingCubes();
 118       GLvoid vMarchCube1(double fX, double fY, double fZ, double fScale);
 119       GLvoid vMarchCube2(double fX, double fY, double fZ, double fScale);
 120       GLvoid (*vMarchCube)(double fX, double fY, double fZ, double fScale) = vMarchCube1;
 121
 122       void main(int argc, char **argv)
 123       {
 124       double afPropertiesAmbient [] = {0.50, 0.50, 0.50, 1.00};
 125       double afPropertiesDiffuse [] = {0.75, 0.75, 0.75, 1.00};
 126       double afPropertiesSpecular[] = {1.00, 1.00, 1.00, 1.00};
 127
 128       GLsizei iWidth = 640.0;
 129       GLsizei iHeight = 480.0;
 130
 131       glutInit(&argc, argv);
 132       glutInitWindowPosition( 0, 0);
 133       glutInitWindowSize(iWidth, iHeight);
 134       glutInitDisplayMode( GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE );
 135       glutCreateWindow( "Marching Cubes" );
 136       glutDisplayFunc( vDrawScene );
 137       glutIdleFunc( vIdle );
 138       glutReshapeFunc( vResize );
 139       glutKeyboardFunc( vKeyboard );
 140       glutSpecialFunc( vSpecial );
 141
 142       glClearColor( 0.0, 0.0, 0.0, 1.0 );
 143       glClearDepth( 1.0 );
 144
 145       glEnable(GL_DEPTH_TEST);
 146       glEnable(GL_LIGHTING);
 147       glPolygonMode(GL_FRONT_AND_BACK, ePolygonMode);
 148
 149       glLightfv( GL_LIGHT0, GL_AMBIENT,  afPropertiesAmbient);
 150       glLightfv( GL_LIGHT0, GL_DIFFUSE,  afPropertiesDiffuse);
 151       glLightfv( GL_LIGHT0, GL_SPECULAR, afPropertiesSpecular);
 152       glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, 1.0);
 153
 154       glEnable( GL_LIGHT0 );
 155
 156       glMaterialfv(GL_BACK,  GL_AMBIENT,   afAmbientGreen);
 157       glMaterialfv(GL_BACK,  GL_DIFFUSE,   afDiffuseGreen);
 158       glMaterialfv(GL_FRONT, GL_AMBIENT,   afAmbientBlue);
 159       glMaterialfv(GL_FRONT, GL_DIFFUSE,   afDiffuseBlue);
 160       glMaterialfv(GL_FRONT, GL_SPECULAR,  afSpecularWhite);
 161       glMaterialf( GL_FRONT, GL_SHININESS, 25.0);
 162
 163       vResize(iWidth, iHeight);
 164
 165       vPrintHelp();
 166       glutMainLoop();
 167       }
 168
 169       GLvoid vPrintHelp()
 170       {
 171       printf("Marching Cubes Example by Cory Bloyd (dejaspaminacan@my-deja.com)\n\n");
 172
 173       printf("+/-  increase/decrease sample density\n");
 174       printf("PageUp/PageDown  increase/decrease surface value\n");
 175       printf("s  change sample function\n");
 176       printf("c  toggle marching cubes / marching tetrahedrons\n");
 177       printf("w  wireframe on/off\n");
 178       printf("l  toggle lighting / color-by-normal\n");
 179       printf("Home  spin scene on/off\n");
 180       printf("End  source point animation on/off\n");
 181       }
 182
 183
 184       void vResize( GLsizei iWidth, GLsizei iHeight )
 185       {
 186       double fAspect, fHalfWorldSize = (1.4142135623730950488016887242097/2);
 187
 188       glViewport( 0, 0, iWidth, iHeight );
 189       glMatrixMode (GL_PROJECTION);
 190       glLoadIdentity ();
 191
 192       if(iWidth <= iHeight)
 193       {
 194       fAspect = (double)iHeight / (double)iWidth;
 195       glOrtho(-fHalfWorldSize, fHalfWorldSize, -fHalfWorldSize*fAspect,
 196       fHalfWorldSize*fAspect, -10*fHalfWorldSize, 10*fHalfWorldSize);
 197       }
 198       else
 199       {
 200       fAspect = (double)iWidth / (double)iHeight;
 201       glOrtho(-fHalfWorldSize*fAspect, fHalfWorldSize*fAspect, -fHalfWorldSize,
 202       fHalfWorldSize, -10*fHalfWorldSize, 10*fHalfWorldSize);
 203       }
 204
 205       glMatrixMode( GL_MODELVIEW );
 206       }
 207
 208       void vKeyboard(unsigned char cKey, int iX, int iY)
 209       {
 210       switch(cKey)
 211       {
 212       case 'w' :
 213       {
 214       if(ePolygonMode == GL_LINE)
 215       {
 216       ePolygonMode = GL_FILL;
 217       }
 218       else
 219       {
 220       ePolygonMode = GL_LINE;
 221       }
 222       glPolygonMode(GL_FRONT_AND_BACK, ePolygonMode);
 223       } break;
 224       case '+' :
 225       case '=' :
 226       {
 227       ++iDataSetSize;
 228       fStepSize = 1.0/iDataSetSize;
 229       } break;
 230       case '-' :
 231       {
 232       if(iDataSetSize > 1)
 233       {
 234       --iDataSetSize;
 235       fStepSize = 1.0/iDataSetSize;
 236       }
 237       } break;
 238       case 'c' :
 239       {
 240       if(vMarchCube == vMarchCube1)
 241       {
 242       vMarchCube = vMarchCube2;//Use Marching Tetrahedrons
 243       }
 244       else
 245       {
 246       vMarchCube = vMarchCube1;//Use Marching Cubes
 247       }
 248       } break;
 249       case 's' :
 250       {
 251       if(fSample == fSample1)
 252       {
 253       fSample = fSample2;
 254       }
 255       else if(fSample == fSample2)
 256       {
 257       fSample = fSample3;
 258       }
 259       else
 260       {
 261       fSample = fSample1;
 262       }
 263       } break;
 264       case 'l' :
 265       {
 266       if(bLight)
 267       {
 268       glDisable(GL_LIGHTING);//use vertex colors
 269       }
 270       else
 271       {
 272       glEnable(GL_LIGHTING);//use lit material color
 273       }
 274
 275       bLight = !bLight;
 276       };
 277       }
 278       }
 279
 280
 281       void vSpecial(int iKey, int iX, int iY)
 282       {
 283       switch(iKey)
 284       {
 285       case GLUT_KEY_PAGE_UP :
 286       {
 287       if(fTargetValue < 1000.0)
 288       {
 289       fTargetValue *= 1.1;
 290       }
 291       } break;
 292       case GLUT_KEY_PAGE_DOWN :
 293       {
 294       if(fTargetValue > 1.0)
 295       {
 296       fTargetValue /= 1.1;
 297       }
 298       } break;
 299       case GLUT_KEY_HOME :
 300       {
 301       bSpin = !bSpin;
 302       } break;
 303       case GLUT_KEY_END :
 304       {
 305       bMove = !bMove;
 306       } break;
 307       }
 308       }
 309
 310       void vIdle()
 311       {
 312       glutPostRedisplay();
 313       }
 314
 315       void vDrawScene()
 316       {
 317       static double fPitch = 0.0;
 318       static double fYaw   = 0.0;
 319       static double fTime = 0.0;
 320
 321       glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
 322
 323       glPushMatrix();
 324
 325       if(bSpin)
 326       {
 327       fPitch += 4.0;
 328       fYaw   += 2.5;
 329       }
 330       if(bMove)
 331       {
 332       fTime  += 0.025;
 333       }
 334
 335       vSetTime(fTime);
 336
 337       glTranslatef(0.0, 0.0, -1.0);
 338       glRotatef( -fPitch, 1.0, 0.0, 0.0);
 339       glRotatef(     0.0, 0.0, 1.0, 0.0);
 340       glRotatef(    fYaw, 0.0, 0.0, 1.0);
 341
 342       glPushAttrib(GL_LIGHTING_BIT);
 343       glDisable(GL_LIGHTING);
 344       glColor3f(1.0, 1.0, 1.0);
 345       glutWireCube(1.0);
 346       glPopAttrib();
 347
 348
 349       glPushMatrix();
 350       glTranslatef(-0.5, -0.5, -0.5);
 351       glBegin(GL_TRIANGLES);
 352       vMarchingCubes();
 353       glEnd();
 354       glPopMatrix();
 355
 356
 357       glPopMatrix();
 358
 359       glutSwapBuffers();
 360       }
 361
 362
 363
 364       //vGetColor generates a color from a given position and normal of a point
 365       GLvoid vGetColor(GLvector &rfColor, GLvector &rfPosition, GLvector &rfNormal)
 366       {
 367       double fX = rfNormal.fX;
 368       double fY = rfNormal.fY;
 369       double fZ = rfNormal.fZ;
 370       rfColor.fX = (fX > 0.0 ? fX : 0.0) + (fY < 0.0 ? -0.5*fY : 0.0) + (fZ < 0.0 ? -0.5*fZ : 0.0);
 371       rfColor.fY = (fY > 0.0 ? fY : 0.0) + (fZ < 0.0 ? -0.5*fZ : 0.0) + (fX < 0.0 ? -0.5*fX : 0.0);
 372       rfColor.fZ = (fZ > 0.0 ? fZ : 0.0) + (fX < 0.0 ? -0.5*fX : 0.0) + (fY < 0.0 ? -0.5*fY : 0.0);
 373       }
 374
 375
 376
 377       //Generate a sample data set.  fSample1(), fSample2() and fSample3() define three scalar fields whose
 378       // values vary by the X,Y and Z coordinates and by the fTime value set by vSetTime()
 379       GLvoid vSetTime(double fNewTime)
 380       {
 381       double fOffset;
 382       int iSourceNum;
 383
 384       for(iSourceNum = 0; iSourceNum < 3; iSourceNum++)
 385       {
 386       sSourcePoint[iSourceNum].fX = 0.5;
 387       sSourcePoint[iSourceNum].fY = 0.5;
 388       sSourcePoint[iSourceNum].fZ = 0.5;
 389       }
 390
 391       fTime = fNewTime;
 392       fOffset = 1.0 + sinf(fTime);
 393       sSourcePoint[0].fX *= fOffset;
 394       sSourcePoint[1].fY *= fOffset;
 395       sSourcePoint[2].fZ *= fOffset;
 396       }
 397
 398
 399       //fSample1 finds the distance of (fX, fY, fZ) from three moving points
 400       double fSample1(double fX, double fY, double fZ)
 401       {
 402       GLdouble fResult = 0.0;
 403       GLdouble fDx, fDy, fDz;
 404       fDx = fX - sSourcePoint[0].fX;
 405       fDy = fY - sSourcePoint[0].fY;
 406       fDz = fZ - sSourcePoint[0].fZ;
 407       fResult += 0.5/(fDx*fDx + fDy*fDy + fDz*fDz);
 408
 409       fDx = fX - sSourcePoint[1].fX;
 410       fDy = fY - sSourcePoint[1].fY;
 411       fDz = fZ - sSourcePoint[1].fZ;
 412       fResult += 1.0/(fDx*fDx + fDy*fDy + fDz*fDz);
 413
 414       fDx = fX - sSourcePoint[2].fX;
 415       fDy = fY - sSourcePoint[2].fY;
 416       fDz = fZ - sSourcePoint[2].fZ;
 417       fResult += 1.5/(fDx*fDx + fDy*fDy + fDz*fDz);
 418
 419       return fResult;
 420       }
 421
 422       //fSample2 finds the distance of (fX, fY, fZ) from three moving lines
 423       double fSample2(double fX, double fY, double fZ)
 424       {
 425       GLdouble fResult = 0.0;
 426       GLdouble fDx, fDy, fDz;
 427       fDx = fX - sSourcePoint[0].fX;
 428       fDy = fY - sSourcePoint[0].fY;
 429       fResult += 0.5/(fDx*fDx + fDy*fDy);
 430
 431       fDx = fX - sSourcePoint[1].fX;
 432       fDz = fZ - sSourcePoint[1].fZ;
 433       fResult += 0.75/(fDx*fDx + fDz*fDz);
 434
 435       fDy = fY - sSourcePoint[2].fY;
 436       fDz = fZ - sSourcePoint[2].fZ;
 437       fResult += 1.0/(fDy*fDy + fDz*fDz);
 438
 439       return fResult;
 440       }
 441
 442
 443       //fSample2 defines a height field by plugging the distance from the center into the sin and cos functions
 444       double fSample3(double fX, double fY, double fZ)
 445       {
 446       double fHeight = 20.0*(fTime + sqrt((0.5-fX)*(0.5-fX) + (0.5-fY)*(0.5-fY)));
 447       fHeight = 1.5 + 0.1*(sinf(fHeight) + cosf(fHeight));
 448       GLdouble fResult = (fHeight - fZ)*50.0;
 449
 450       return fResult;
 451       }
 452
 453
 454
 455
 456
 457
 458       //vMarchTetrahedron performs the Marching Tetrahedrons algorithm on a single tetrahedron
 459       GLvoid vMarchTetrahedron(GLvector *pasTetrahedronPosition, double *pafTetrahedronValue)
 460       {
 461       extern int aiTetrahedronEdgeFlags[16];
 462       extern int a2iTetrahedronTriangles[16][7];
 463
 464       int iEdge, iVert0, iVert1, iEdgeFlags, iTriangle, iCorner, iVertex, iFlagIndex = 0;
 465       double fOffset, fInvOffset, fValue = 0.0;
 466       GLvector asEdgeVertex[6];
 467       GLvector asEdgeNorm[6];
 468       GLvector sColor;
 469
 470       //Find which vertices are inside of the surface and which are outside
 471       for(iVertex = 0; iVertex < 4; iVertex++)
 472       {
 473       if(pafTetrahedronValue[iVertex] <= fTargetValue)
 474       iFlagIndex |= 1<<iVertex;
 475       }
 476
 477       //Find which edges are intersected by the surface
 478       iEdgeFlags = aiTetrahedronEdgeFlags[iFlagIndex];
 479
 480       //If the tetrahedron is entirely inside or outside of the surface, then there will be no intersections
 481       if(iEdgeFlags == 0)
 482       {
 483       return;
 484       }
 485       //Find the point of intersection of the surface with each edge
 486       // Then find the normal to the surface at those points
 487       for(iEdge = 0; iEdge < 6; iEdge++)
 488       {
 489       //if there is an intersection on this edge
 490       if(iEdgeFlags & (1<<iEdge))
 491       {
 492       iVert0 = a2iTetrahedronEdgeConnection[iEdge][0];
 493       iVert1 = a2iTetrahedronEdgeConnection[iEdge][1];
 494       fOffset = fGetOffset(pafTetrahedronValue[iVert0], pafTetrahedronValue[iVert1], fTargetValue);
 495       fInvOffset = 1.0 - fOffset;
 496
 497       asEdgeVertex[iEdge].fX = fInvOffset*pasTetrahedronPosition[iVert0].fX  +  fOffset*pasTetrahedronPosition[iVert1].fX;
 498       asEdgeVertex[iEdge].fY = fInvOffset*pasTetrahedronPosition[iVert0].fY  +  fOffset*pasTetrahedronPosition[iVert1].fY;
 499       asEdgeVertex[iEdge].fZ = fInvOffset*pasTetrahedronPosition[iVert0].fZ  +  fOffset*pasTetrahedronPosition[iVert1].fZ;
 500
 501       vGetNormal(asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX, asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
 502       }
 503       }
 504       //Draw the triangles that were found.  There can be up to 2 per tetrahedron
 505       for(iTriangle = 0; iTriangle < 2; iTriangle++)
 506       {
 507       if(a2iTetrahedronTriangles[iFlagIndex][3*iTriangle] < 0)
 508       break;
 509
 510       for(iCorner = 0; iCorner < 3; iCorner++)
 511       {
 512       iVertex = a2iTetrahedronTriangles[iFlagIndex][3*iTriangle+iCorner];
 513
 514       vGetColor(sColor, asEdgeVertex[iVertex], asEdgeNorm[iVertex]);
 515       glColor3f(sColor.fX, sColor.fY, sColor.fZ);
 516       glNormal3f(asEdgeNorm[iVertex].fX,   asEdgeNorm[iVertex].fY,   asEdgeNorm[iVertex].fZ);
 517       glVertex3f(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY, asEdgeVertex[iVertex].fZ);
 518       }
 519       }
 520       }
 521
 522
 523
 524       //vMarchCube2 performs the Marching Tetrahedrons algorithm on a single cube by making six calls to vMarchTetrahedron
 525       GLvoid vMarchCube2(double fX, double fY, double fZ, double fScale)
 526       {
 527       int iVertex, iTetrahedron, iVertexInACube;
 528       GLvector asCubePosition[8];
 529       double  afCubeValue[8];
 530       GLvector asTetrahedronPosition[4];
 531       double  afTetrahedronValue[4];
 532
 533       //Make a local copy of the cube's corner positions
 534       for(iVertex = 0; iVertex < 8; iVertex++)
 535       {
 536       asCubePosition[iVertex].fX = fX + a2fVertexOffset[iVertex][0]*fScale;
 537       asCubePosition[iVertex].fY = fY + a2fVertexOffset[iVertex][1]*fScale;
 538       asCubePosition[iVertex].fZ = fZ + a2fVertexOffset[iVertex][2]*fScale;
 539       }
 540
 541       //Make a local copy of the cube's corner values
 542       for(iVertex = 0; iVertex < 8; iVertex++)
 543       {
 544       afCubeValue[iVertex] = fSample(asCubePosition[iVertex].fX,
 545       asCubePosition[iVertex].fY,
 546       asCubePosition[iVertex].fZ);
 547       }
 548
 549       for(iTetrahedron = 0; iTetrahedron < 6; iTetrahedron++)
 550       {
 551       for(iVertex = 0; iVertex < 4; iVertex++)
 552       {
 553       iVertexInACube = a2iTetrahedronsInACube[iTetrahedron][iVertex];
 554       asTetrahedronPosition[iVertex].fX = asCubePosition[iVertexInACube].fX;
 555       asTetrahedronPosition[iVertex].fY = asCubePosition[iVertexInACube].fY;
 556       asTetrahedronPosition[iVertex].fZ = asCubePosition[iVertexInACube].fZ;
 557       afTetrahedronValue[iVertex] = afCubeValue[iVertexInACube];
 558       }
 559       vMarchTetrahedron(asTetrahedronPosition, afTetrahedronValue);
 560       }
 561       }
 562     */
 563
 564
 565     // vMarchingCubes iterates over the entire dataset, calling vMarchCube on each cube
 566     void vMarchingCubes(int iDataSetSize, int fStepSize)
 567     {
 568         int iX, iY, iZ;
 569         for(iX = 0; iX < iDataSetSize; iX++)
 570             for(iY = 0; iY < iDataSetSize; iY++)
 571                 for(iZ = 0; iZ < iDataSetSize; iZ++)
 572                     {
 573                         vMarchCube(iX*fStepSize, iY*fStepSize, iZ*fStepSize, fStepSize);
 574                     }
 575     }
 576
 577     public abstract double fSample(double x, double y, double z);
 578     public abstract double fTargetValue();
 579
 580     //vGetNormal() finds the gradient of the scalar field at a point
 581     //This gradient can be used as a very accurate vertx normal for lighting calculations
 582     void vGetNormal(GLvector rfNormal, double fX, double fY, double fZ) {
 583         rfNormal.fX = fSample(fX-0.01, fY, fZ) - fSample(fX+0.01, fY, fZ);
 584         rfNormal.fY = fSample(fX, fY-0.01, fZ) - fSample(fX, fY+0.01, fZ);
 585         rfNormal.fZ = fSample(fX, fY, fZ-0.01) - fSample(fX, fY, fZ+0.01);
 586         vNormalizeVector(rfNormal, rfNormal);
 587     }
 588
 589     void vNormalizeVector(GLvector rfVectorResult, GLvector rfVectorSource) {
 590         double fOldLength;
 591         double fScale;
 592
 593         fOldLength = Math.sqrt( (rfVectorSource.fX * rfVectorSource.fX) +
 594                                 (rfVectorSource.fY * rfVectorSource.fY) +
 595                                 (rfVectorSource.fZ * rfVectorSource.fZ) );
 596
 597         if(fOldLength == 0.0) {
 598             rfVectorResult.fX = rfVectorSource.fX;
 599             rfVectorResult.fY = rfVectorSource.fY;
 600             rfVectorResult.fZ = rfVectorSource.fZ;
 601         } else {
 602             fScale = 1.0/fOldLength;
 603             rfVectorResult.fX = rfVectorSource.fX*fScale;
 604             rfVectorResult.fY = rfVectorSource.fY*fScale;
 605             rfVectorResult.fZ = rfVectorSource.fZ*fScale;
 606         }
 607     }
 608
 609     // fGetOffset finds the approximate point of intersection of the surface
 610     // between two points with the values fValue1 and fValue2
 611     double fGetOffset(double fValue1, double fValue2, double fValueDesired) {
 612         double fDelta = fValue2 - fValue1;
 613         if(fDelta == 0.0) return 0.5;
 614         return (fValueDesired - fValue1)/fDelta;
 615     }
 616
 617     //vMarchCube performs the Marching Cubes algorithm on a single cube
 618     void vMarchCube(double fX, double fY, double fZ, double fScale) {
 619         int iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags;
 620         double fOffset;
 621         GLvector sColor;
 622         double afCubeValue[] = new double[8];
 623         GLvector asEdgeVertex[] = new GLvector[12];
 624         GLvector asEdgeNorm[] = new GLvector[12];
 625
 626         //Make a local copy of the values at the cube's corners
 627         for(iVertex = 0; iVertex < 8; iVertex++)
 628             afCubeValue[iVertex] = fSample(fX + a2fVertexOffset[iVertex][0]*fScale,
 629                                            fY + a2fVertexOffset[iVertex][1]*fScale,
 630                                            fZ + a2fVertexOffset[iVertex][2]*fScale);
 631
 632         //Find which vertices are inside of the surface and which are outside
 633         iFlagIndex = 0;
 634         for(iVertexTest = 0; iVertexTest < 8; iVertexTest++) {
 635             if(afCubeValue[iVertexTest] <= fTargetValue())
 636                 iFlagIndex |= 1<<iVertexTest;
 637         }
 638
 639         //Find which edges are intersected by the surface
 640         iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];
 641
 642         //If the cube is entirely inside or outside of the surface, then there will be no intersections
 643         if(iEdgeFlags == 0) return;
 644
 645         //Find the point of intersection of the surface with each edge
 646         //Then find the normal to the surface at those points
 647         for(iEdge = 0; iEdge < 12; iEdge++) {
 648             //if there is an intersection on this edge
 649             if ((iEdgeFlags & (1<<iEdge))!=0) {
 650                 fOffset = fGetOffset(afCubeValue[ a2iEdgeConnection[iEdge][0] ],
 651                                      afCubeValue[ a2iEdgeConnection[iEdge][1] ], fTargetValue());
 652
 653                 asEdgeVertex[iEdge].fX = fX + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][0]  +  fOffset * a2fEdgeDirection[iEdge][0]) * fScale;
 654                 asEdgeVertex[iEdge].fY = fY + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][1]  +  fOffset * a2fEdgeDirection[iEdge][1]) * fScale;
 655                 asEdgeVertex[iEdge].fZ = fZ + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][2]  +  fOffset * a2fEdgeDirection[iEdge][2]) * fScale;
 656
 657                 vGetNormal(asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX, asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
 658             }
 659         }
 660
 661         // Draw the triangles that were found.  There can be up to five per cube
 662         for(iTriangle = 0; iTriangle < 5; iTriangle++) {
 663             if(a2iTriangleConnectionTable[iFlagIndex][3*iTriangle] < 0)
 664                 break;
 665
 666             for(iCorner = 0; iCorner < 3; iCorner++) {
 667                 iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner];
 668                 /*
 669                 vGetColor(sColor, asEdgeVertex[iVertex], asEdgeNorm[iVertex]);
 670                 glColor3f(sColor.fX, sColor.fY, sColor.fZ);
 671                 glNormal3f(asEdgeNorm[iVertex].fX,   asEdgeNorm[iVertex].fY,   asEdgeNorm[iVertex].fZ);
 672                 glVertex3f(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY, asEdgeVertex[iVertex].fZ);
 673                 */
 674                 //mesh.newT(
 675             }
 676         }
 677     }
 678
 679
 680     // For any edge, if one vertex is inside of the surface and the
 681     //  other is outside of the surface then the edge intersects the
 682     //  surface For each of the 4 vertices of the tetrahedron can be
 683     //  two possible states : either inside or outside of the surface
 684     //  For any tetrahedron the are 2^4=16 possible sets of vertex
 685     //  states This table lists the edges intersected by the surface
 686     //  for all 16 possible vertex states There are 6 edges.  For each
 687     //  entry in the table, if edge #n is intersected, then bit #n is
 688     //  set to 1
 689
 690     static int aiTetrahedronEdgeFlags[] = {
 691         0x00, 0x0d, 0x13, 0x1e, 0x26, 0x2b, 0x35, 0x38, 0x38, 0x35, 0x2b, 0x26, 0x1e, 0x13, 0x0d, 0x00,
 692     };
 693
 694
 695     // For each of the possible vertex states listed in
 696     // aiTetrahedronEdgeFlags there is a specific triangulation of the
 697     // edge intersection points.  a2iTetrahedronTriangles lists all of
 698     // them in the form of 0-2 edge triples with the list terminated
 699     // by the invalid value -1.
 700     //
 701     // I generated this table by hand
 702
 703     static int a2iTetrahedronTriangles[][] = {
 704         {-1, -1, -1, -1, -1, -1, -1},
 705         { 0,  3,  2, -1, -1, -1, -1},
 706         { 0,  1,  4, -1, -1, -1, -1},
 707         { 1,  4,  2,  2,  4,  3, -1},
 708
 709         { 1,  2,  5, -1, -1, -1, -1},
 710         { 0,  3,  5,  0,  5,  1, -1},
 711         { 0,  2,  5,  0,  5,  4, -1},
 712         { 5,  4,  3, -1, -1, -1, -1},
 713
 714         { 3,  4,  5, -1, -1, -1, -1},
 715         { 4,  5,  0,  5,  2,  0, -1},
 716         { 1,  5,  0,  5,  3,  0, -1},
 717         { 5,  2,  1, -1, -1, -1, -1},
 718
 719         { 3,  4,  2,  2,  4,  1, -1},
 720         { 4,  1,  0, -1, -1, -1, -1},
 721         { 2,  3,  0, -1, -1, -1, -1},
 722         {-1, -1, -1, -1, -1, -1, -1},
 723     };
 724
 725     // For any edge, if one vertex is inside of the surface and the
 726     //  other is outside of the surface then the edge intersects the
 727     //  surface For each of the 8 vertices of the cube can be two
 728     //  possible states : either inside or outside of the surface For
 729     //  any cube the are 2^8=256 possible sets of vertex states This
 730     //  table lists the edges intersected by the surface for all 256
 731     //  possible vertex states There are 12 edges.  For each entry in
 732     //  the table, if edge #n is intersected, then bit #n is set to 1
 733
 734     int aiCubeEdgeFlags[] = {
 735         0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
 736         0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
 737         0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
 738         0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
 739         0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
 740         0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
 741         0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
 742         0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
 743         0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
 744         0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
 745         0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
 746         0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460,
 747         0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
 748         0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230,
 749         0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190,
 750         0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000
 751     };
 752
 753     //  For each of the possible vertex states listed in
 754     //  aiCubeEdgeFlags there is a specific triangulation of the edge
 755     //  intersection points.  a2iTriangleConnectionTable lists all of
 756     //  them in the form of 0-5 edge triples with the list terminated
 757     //  by the invalid value -1.  For example:
 758     //  a2iTriangleConnectionTable[3] list the 2 triangles formed when
 759     //  corner[0] and corner[1] are inside of the surface, but the
 760     //  rest of the cube is not.
 761     //
 762     //  I found this table in an example program someone wrote long
 763     //  ago.  It was probably generated by hand
 764
 765     int a2iTriangleConnectionTable[][] = {
 766         {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 767         {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 768         {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 769         {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 770         {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 771         {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 772         {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 773         {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
 774         {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 775         {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 776         {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 777         {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
 778         {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 779         {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
 780         {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
 781         {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 782         {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 783         {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 784         {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 785         {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
 786         {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 787         {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
 788         {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
 789         {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
 790         {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 791         {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
 792         {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
 793         {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
 794         {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
 795         {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
 796         {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
 797         {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
 798         {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 799         {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 800         {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 801         {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
 802         {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 803         {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
 804         {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
 805         {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
 806         {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 807         {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
 808         {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
 809         {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
 810         {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
 811         {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
 812         {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
 813         {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
 814         {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 815         {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
 816         {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
 817         {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 818         {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
 819         {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
 820         {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
 821         {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
 822         {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
 823         {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
 824         {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
 825         {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
 826         {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
 827         {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
 828         {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
 829         {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 830         {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 831         {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 832         {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 833         {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
 834         {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 835         {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
 836         {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
 837         {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
 838         {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 839         {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
 840         {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
 841         {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
 842         {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
 843         {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
 844         {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
 845         {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
 846         {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 847         {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
 848         {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
 849         {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
 850         {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
 851         {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
 852         {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
 853         {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
 854         {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
 855         {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
 856         {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
 857         {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
 858         {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
 859         {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
 860         {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
 861         {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
 862         {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 863         {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
 864         {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
 865         {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
 866         {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
 867         {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
 868         {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 869         {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
 870         {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
 871         {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
 872         {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
 873         {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
 874         {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
 875         {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
 876         {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
 877         {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 878         {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
 879         {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
 880         {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
 881         {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
 882         {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
 883         {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
 884         {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
 885         {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 886         {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
 887         {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
 888         {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
 889         {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
 890         {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
 891         {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 892         {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
 893         {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 894         {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 895         {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 896         {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 897         {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
 898         {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 899         {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
 900         {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
 901         {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
 902         {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 903         {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
 904         {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
 905         {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
 906         {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
 907         {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
 908         {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
 909         {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
 910         {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 911         {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
 912         {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
 913         {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
 914         {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
 915         {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
 916         {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
 917         {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
 918         {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
 919         {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 920         {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
 921         {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
 922         {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
 923         {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
 924         {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
 925         {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 926         {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 927         {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
 928         {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
 929         {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
 930         {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
 931         {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
 932         {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
 933         {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
 934         {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
 935         {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
 936         {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
 937         {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
 938         {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
 939         {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
 940         {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
 941         {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
 942         {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
 943         {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
 944         {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
 945         {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
 946         {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
 947         {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
 948         {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
 949         {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
 950         {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
 951         {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
 952         {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
 953         {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 954         {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
 955         {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
 956         {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 957         {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 958         {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 959         {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
 960         {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
 961         {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
 962         {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
 963         {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
 964         {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
 965         {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
 966         {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
 967         {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
 968         {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
 969         {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
 970         {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 971         {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
 972         {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
 973         {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 974         {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
 975         {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
 976         {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
 977         {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
 978         {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
 979         {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
 980         {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
 981         {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 982         {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
 983         {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
 984         {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
 985         {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
 986         {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
 987         {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 988         {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
 989         {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
 990         {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
 991         {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
 992         {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
 993         {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
 994         {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
 995         {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
 996         {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
 997         {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
 998         {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
 999         {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
1000         {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
1001         {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1002         {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
1003         {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
1004         {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1005         {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1006         {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1007         {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
1008         {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
1009         {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1010         {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
1011         {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
1012         {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1013         {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1014         {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
1015         {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1016         {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
1017         {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1018         {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1019         {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1020         {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
1021         {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
1022     };
1023
1024 }