X-Git-Url: http://git.megacz.com/?a=blobdiff_plain;ds=sidebyside;f=src%2Fedu%2Fberkeley%2Fqfat%2Fvoxel%2FMarchingCubes.java;h=46b83b6c0b80570b9b5896dc0d51a67089c7de76;hb=HEAD;hp=58434f460f81a62b3f5ab530e116227a78bd6817;hpb=a75c872663256575fa0c7330e444c0c325683003;p=anneal.git

diff --git a/src/edu/berkeley/qfat/voxel/MarchingCubes.java b/src/edu/berkeley/qfat/voxel/MarchingCubes.java
index 58434f4..46b83b6 100644
--- a/src/edu/berkeley/qfat/voxel/MarchingCubes.java
+++ b/src/edu/berkeley/qfat/voxel/MarchingCubes.java
@@ -35,18 +35,25 @@ public class MarchingCubes {
     };
 
     /** march iterates over the entire dataset, calling vMarchCube on each cube */
-    public static void march(SampledField sampledField, double threshold, int iDataSetSize, double fStepSize, Mesh mesh) {
-        int iX, iY, iZ;
+    public static void march(VoxelData voxels, Mesh mesh) { march(voxels, 0, mesh); }
+    public static void march(VoxelData voxels, double threshold, Mesh mesh) {
         int initialTriangles = mesh.numTriangles();
-        for(iX = 0; iX < iDataSetSize; iX++) {
+        double scaleX = (voxels.getMaxX() - voxels.getMinX()) / (double)voxels.getNumSamplesX();
+        double scaleY = (voxels.getMaxY() - voxels.getMinY()) / (double)voxels.getNumSamplesY();
+        double scaleZ = (voxels.getMaxZ() - voxels.getMinZ()) / (double)voxels.getNumSamplesZ();
+        for(int iX = 0; iX < voxels.getNumSamplesX(); iX++) {
             System.out.print("\r");
             for(int i=0; i<78; i++) System.out.print(' ');
             System.out.print("\r");
-            System.out.print(Math.ceil((iX/((double)iDataSetSize))*100) + "% marched, " +
+            System.out.print(Math.ceil((iX/((double)voxels.getNumSamplesX()))*100) + "% marched, " +
                              (mesh.numTriangles()-initialTriangles) + " triangles");
-            for(iY = 0; iY < iDataSetSize; iY++)
-                for(iZ = 0; iZ < iDataSetSize; iZ++)
-                    march(sampledField, mesh, threshold, iX*fStepSize, iY*fStepSize, iZ*fStepSize, fStepSize);
+            for(int iY = 0; iY < voxels.getNumSamplesY(); iY++)
+                for(int iZ = 0; iZ < voxels.getNumSamplesZ(); iZ++)
+                    /*marchCubes*/marchTetrahedrons(voxels, mesh, threshold,
+                                                    iX*scaleX + voxels.getMinX(),
+                                                    iY*scaleY + voxels.getMinY(),
+                                                    iZ*scaleZ + voxels.getMinZ(),
+                                                    scaleX, scaleY, scaleZ);
         }
         System.out.print("\r");
         for(int i=0; i<78; i++) System.out.print(' ');
@@ -55,7 +62,9 @@ public class MarchingCubes {
     }
 
     /** performs the Marching Cubes algorithm on a single cube */
-    static void march(SampledField sampledField, Mesh mesh, double threshold, double fX, double fY, double fZ, double fScale) {
+    static void marchCubes(VoxelData voxels, Mesh mesh, double threshold,
+                           double fX, double fY, double fZ,
+                           double scaleX, double scaleY, double scaleZ) {
         int iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags;
         double fOffset;
         GLvector sColor;
@@ -68,17 +77,15 @@ public class MarchingCubes {
 
         // Make a local copy of the values at the cube's corners
         for(iVertex = 0; iVertex < 8; iVertex++)
-            afCubeValue[iVertex] = sampledField.getSample(new Point(fX + a2fVertexOffset[iVertex][0]*fScale,
-                                                                    fY + a2fVertexOffset[iVertex][1]*fScale,
-                                                                    fZ + a2fVertexOffset[iVertex][2]*fScale));
+            afCubeValue[iVertex] = voxels.getSample(new Point(fX + a2fVertexOffset[iVertex][0]*scaleX,
+                                                                    fY + a2fVertexOffset[iVertex][1]*scaleY,
+                                                                    fZ + a2fVertexOffset[iVertex][2]*scaleZ));
 
         // Find which vertices are inside of the surface and which are outside
         iFlagIndex = 0;
-        for(iVertexTest = 0; iVertexTest < 8; iVertexTest++) {
-            if (afCubeValue[iVertexTest] >= threshold) {
+        for(iVertexTest = 0; iVertexTest < 8; iVertexTest++)
+            if (afCubeValue[iVertexTest] >= threshold)
                 iFlagIndex |= 1<<iVertexTest;
-            }
-        }
         
         // Find which edges are intersected by the surface
         iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];
@@ -92,15 +99,19 @@ public class MarchingCubes {
             // If there is an intersection on this edge
             if ((iEdgeFlags & (1<<iEdge))==0) continue;
             fOffset = fGetOffset(afCubeValue[ a2iEdgeConnection[iEdge][0] ], 
-                                 afCubeValue[ a2iEdgeConnection[iEdge][1] ], threshold);
+                                 afCubeValue[ a2iEdgeConnection[iEdge][1] ],
+                                 threshold,
+                                 Math.min(Math.min(scaleX, scaleY), scaleZ) * 0.1);
             
-            asEdgeVertex[iEdge].fX = fX + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][0]  +  fOffset * a2fEdgeDirection[iEdge][0]) * fScale;
-            asEdgeVertex[iEdge].fY = fY + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][1]  +  fOffset * a2fEdgeDirection[iEdge][1]) * fScale;
-            asEdgeVertex[iEdge].fZ = fZ + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][2]  +  fOffset * a2fEdgeDirection[iEdge][2]) * fScale;
+            asEdgeVertex[iEdge].fX = fX + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][0]  +  fOffset * a2fEdgeDirection[iEdge][0]) * scaleX;
+            asEdgeVertex[iEdge].fY = fY + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][1]  +  fOffset * a2fEdgeDirection[iEdge][1]) * scaleY;
+            asEdgeVertex[iEdge].fZ = fZ + (a2fVertexOffset[ a2iEdgeConnection[iEdge][0] ][2]  +  fOffset * a2fEdgeDirection[iEdge][2]) * scaleZ;
             
-            vGetNormal(sampledField, asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX, asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
+            vGetNormal(voxels, asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX, asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
         }
 
+        System.out.println();
+
         // Draw the triangles that were found.  There can be up to five per cube
         for(iTriangle = 0; iTriangle < 5; iTriangle++) {
             if(a2iTriangleConnectionTable[iFlagIndex][3*iTriangle] < 0)
@@ -115,10 +126,132 @@ public class MarchingCubes {
                 // questionable, but we do it anyways
                 norm = norm.plus(new Vec(asEdgeNorm[iVertex].fX,   asEdgeNorm[iVertex].fY,   asEdgeNorm[iVertex].fZ));
             }
-            if (points[0].equals(points[1])) continue;
-            if (points[0].equals(points[2])) continue;
-            if (points[1].equals(points[2])) continue;
-            mesh.newT(points[0], points[1], points[2], norm.norm().times(-1));
+
+            // Eliminate triangles with "length-zero" sides.
+            // Unfortunately this puts holes in the mesh.
+            if (points[0].equals(points[1]) ||
+                points[0].equals(points[2]) ||
+                points[1].equals(points[2]))
+                continue;
+
+            System.out.println("creating triangle: " + points[0] + " " + points[1] + " " + points[2]);
+            try {
+                Mesh.T t = mesh.newT(points[0], points[1], points[2], norm.norm().times(-1));
+                System.out.println("  created " + t);
+            } catch (Throwable t) { }
+        }
+    }
+
+    /**
+     *  marchTetrahedron performs the Marching Tetrahedrons algorithm
+     *  on a single tetrahedron
+     */
+    static void marchTetrahedron(VoxelData voxels, Mesh mesh, double threshold,
+                                 GLvector[] pasTetrahedronPosition,
+                                 float[] pafTetrahedronValue,
+                                 double scaleX, double scaleY, double scaleZ) {
+
+        int iEdge, iVert0, iVert1, iEdgeFlags, iTriangle, iCorner, iVertex, iFlagIndex = 0;
+        double fOffset, fInvOffset, fValue = 0.0;
+        GLvector[] asEdgeVertex = new GLvector[6];
+        GLvector[] asEdgeNorm = new GLvector[6];
+        GLvector sColor = new GLvector();
+        for(int i=0; i<asEdgeVertex.length; i++) asEdgeVertex[i] = new GLvector();
+        for(int i=0; i<asEdgeNorm.length; i++) asEdgeNorm[i] = new GLvector();
+
+        // Find which vertices are inside of the surface and which are outside
+        for(iVertex = 0; iVertex < 4; iVertex++)
+            if(pafTetrahedronValue[iVertex] <= threshold) 
+                iFlagIndex |= 1<<iVertex;
+
+        // Find which edges are intersected by the surface
+        iEdgeFlags = aiTetrahedronEdgeFlags[iFlagIndex];
+
+        // If the tetrahedron is entirely inside or outside of the surface, then there will be no intersections
+        if(iEdgeFlags == 0) return;
+
+        // Find the point of intersection of the surface with each edge
+        // Then find the normal to the surface at those points
+        for(iEdge = 0; iEdge < 6; iEdge++) {
+            // if there is an intersection on this edge
+            if( (iEdgeFlags & (1<<iEdge))==0 ) continue;
+
+            iVert0 = a2iTetrahedronEdgeConnection[iEdge][0];
+            iVert1 = a2iTetrahedronEdgeConnection[iEdge][1];
+            fOffset = fGetOffset(pafTetrahedronValue[iVert0], pafTetrahedronValue[iVert1], threshold,
+                                 Math.min(Math.min(scaleX, scaleY), scaleZ) * 0.1);
+            fInvOffset = 1.0f - fOffset;
+            
+            asEdgeVertex[iEdge].fX = fInvOffset*pasTetrahedronPosition[iVert0].fX  +  fOffset*pasTetrahedronPosition[iVert1].fX;
+            asEdgeVertex[iEdge].fY = fInvOffset*pasTetrahedronPosition[iVert0].fY  +  fOffset*pasTetrahedronPosition[iVert1].fY;
+            asEdgeVertex[iEdge].fZ = fInvOffset*pasTetrahedronPosition[iVert0].fZ  +  fOffset*pasTetrahedronPosition[iVert1].fZ;
+            
+            vGetNormal(voxels, asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX, asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
+        }
+
+        // Draw the triangles that were found.  There can be up to 2 per tetrahedron
+        for(iTriangle = 0; iTriangle < 2; iTriangle++) {
+            if(a2iTetrahedronTriangles[iFlagIndex][3*iTriangle] < 0) break;
+
+            Point[] points = new Point[3];
+            Vec norm = new Vec(0,0,0);
+            for(iCorner = 0; iCorner < 3; iCorner++) {
+                iVertex = a2iTetrahedronTriangles[iFlagIndex][3*iTriangle+iCorner];
+                points[iCorner] = new Point(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY, asEdgeVertex[iVertex].fZ);
+                // questionable, but we do it anyways
+                norm = norm.plus(new Vec(asEdgeNorm[iVertex].fX,   asEdgeNorm[iVertex].fY,   asEdgeNorm[iVertex].fZ));
+            }
+            // Eliminate triangles with "length-zero" sides.
+            // Unfortunately this puts holes in the mesh.
+            if (points[0].equals(points[1]) ||
+                points[0].equals(points[2]) ||
+                points[1].equals(points[2]))
+                continue;
+            System.out.println("creating triangle: " + points[0] + " " + points[1] + " " + points[2]);
+            Mesh.T t = mesh.newT(points[0], points[1], points[2], norm.norm().times(-1));
+            System.out.println("  created " + t);
+        }
+    }
+
+    /**
+     *  marchTetrahedrons performs the Marching Tetrahedrons algorithm on a
+     *  single cube by making six calls to marchTetrahedron
+     */
+    static void marchTetrahedrons(VoxelData voxels, Mesh mesh, double threshold,
+                                  double fX, double fY, double fZ,
+                                  double scaleX, double scaleY, double scaleZ) {
+        int iVertex, iTetrahedron, iVertexInACube;
+        GLvector[] asCubePosition = new GLvector[8];
+        float[]  afCubeValue = new float[8];
+        GLvector[] asTetrahedronPosition = new GLvector[8];
+        float[]  afTetrahedronValue = new float[4];
+
+        for(int i=0; i<asCubePosition.length; i++) asCubePosition[i] = new GLvector();
+        for(int i=0; i<asTetrahedronPosition.length; i++) asTetrahedronPosition[i] = new GLvector();
+
+        // Make a local copy of the cube's corner positions
+        for(iVertex = 0; iVertex < 8; iVertex++) {
+            asCubePosition[iVertex].fX = fX + a2fVertexOffset[iVertex][0]*scaleX;
+            asCubePosition[iVertex].fY = fY + a2fVertexOffset[iVertex][1]*scaleY;
+            asCubePosition[iVertex].fZ = fZ + a2fVertexOffset[iVertex][2]*scaleZ;
+        }
+
+        // Make a local copy of the cube's corner values
+        for(iVertex = 0; iVertex < 8; iVertex++)
+            afCubeValue[iVertex] =
+                voxels.getSample(new Point(asCubePosition[iVertex].fX,
+                                           asCubePosition[iVertex].fY,
+                                           asCubePosition[iVertex].fZ));
+
+        for(iTetrahedron = 0; iTetrahedron < 6; iTetrahedron++) {
+            for(iVertex = 0; iVertex < 4; iVertex++) {
+                iVertexInACube = a2iTetrahedronsInACube[iTetrahedron][iVertex];
+                asTetrahedronPosition[iVertex].fX = asCubePosition[iVertexInACube].fX;
+                asTetrahedronPosition[iVertex].fY = asCubePosition[iVertexInACube].fY;
+                asTetrahedronPosition[iVertex].fZ = asCubePosition[iVertexInACube].fZ;
+                afTetrahedronValue[iVertex] = afCubeValue[iVertexInACube];
+            }
+            marchTetrahedron(voxels, mesh, threshold, asTetrahedronPosition, afTetrahedronValue, scaleX, scaleY, scaleZ);
         }
     }
 
@@ -127,16 +260,16 @@ public class MarchingCubes {
      *  This gradient can be used as a very accurate vertx normal for
      *  lighting calculations
      */
-    static void vGetNormal(SampledField sampledField, GLvector rfNormal, double fX, double fY, double fZ) {
+    static void vGetNormal(VoxelData voxels, GLvector rfNormal, double fX, double fY, double fZ) {
         rfNormal.fX =
-            sampledField.getSample(new Point(fX-0.01, fY, fZ)) -
-            sampledField.getSample(new Point(fX+0.01, fY, fZ));
+            voxels.getSample(new Point(fX-0.01, fY, fZ)) -
+            voxels.getSample(new Point(fX+0.01, fY, fZ));
         rfNormal.fY =
-            sampledField.getSample(new Point(fX, fY-0.01, fZ)) -
-            sampledField.getSample(new Point(fX, fY+0.01, fZ));
+            voxels.getSample(new Point(fX, fY-0.01, fZ)) -
+            voxels.getSample(new Point(fX, fY+0.01, fZ));
         rfNormal.fZ =
-            sampledField.getSample(new Point(fX, fY, fZ-0.01)) -
-            sampledField.getSample(new Point(fX, fY, fZ+0.01));
+            voxels.getSample(new Point(fX, fY, fZ-0.01)) -
+            voxels.getSample(new Point(fX, fY, fZ+0.01));
         vNormalizeVector(rfNormal, rfNormal);
     }
 
@@ -162,9 +295,17 @@ public class MarchingCubes {
 
     // fGetOffset finds the approximate point of intersection of the surface
     // between two points with the values fValue1 and fValue2
-    static double fGetOffset(double fValue1, double fValue2, double fValueDesired) {
+    static double fGetOffset(double fValue1, double fValue2, double fValueDesired, double EPSILON) {
         double fDelta = fValue2 - fValue1;
         if(fDelta == 0.0) return 0.5;
+
+        // the following two lines are a hack; they "snap" the
+        // estimate to one grid point or the other if the distance is
+        // less than some EPSILON.  This ensures that the resulting
+        // mesh is watertight and meets the requirements of Mesh.java
+        if (Math.abs(fValueDesired-fValue1) < EPSILON) return 0;
+        if (Math.abs(fValueDesired-fValue2) < EPSILON) return 1;
+
         return (fValueDesired - fValue1)/fDelta;
     }