MarchingCubes.java: snap points to grid if they are close; this eliminates many mesh...

[anneal.git] / src / edu / berkeley / qfat / voxel / MarchingCubes.java

diff --git a/src/edu/berkeley/qfat/voxel/MarchingCubes.java b/src/edu/berkeley/qfat/voxel/MarchingCubes.java
index bda9531..a365585 100644 (file)
--- a/src/edu/berkeley/qfat/voxel/MarchingCubes.java
+++ b/src/edu/berkeley/qfat/voxel/MarchingCubes.java
@@ -35,7 +35,7 @@ public class MarchingCubes {
     };
 
     /** march iterates over the entire dataset, calling vMarchCube on each cube */
     };
 
     /** march iterates over the entire dataset, calling vMarchCube on each cube */

-    public static void march(SampledField sampledField, double targetValue, int iDataSetSize, double fStepSize, Mesh mesh) {
+    public static void march(SampledField sampledField, double threshold, int iDataSetSize, double fStepSize, Mesh mesh) {

         int iX, iY, iZ;
         int initialTriangles = mesh.numTriangles();
         for(iX = 0; iX < iDataSetSize; iX++) {
         int iX, iY, iZ;
         int initialTriangles = mesh.numTriangles();
         for(iX = 0; iX < iDataSetSize; iX++) {

@@ -46,7 +46,7 @@ public class MarchingCubes {
                              (mesh.numTriangles()-initialTriangles) + " triangles");
             for(iY = 0; iY < iDataSetSize; iY++)
                 for(iZ = 0; iZ < iDataSetSize; iZ++)
                              (mesh.numTriangles()-initialTriangles) + " triangles");
             for(iY = 0; iY < iDataSetSize; iY++)
                 for(iZ = 0; iZ < iDataSetSize; iZ++)

-                    march(sampledField, mesh, targetValue, iX*fStepSize, iY*fStepSize, iZ*fStepSize, fStepSize);
+                    march(sampledField, mesh, threshold, iX*fStepSize, iY*fStepSize, iZ*fStepSize, fStepSize);

         }
         System.out.print("\r");
         for(int i=0; i<78; i++) System.out.print(' ');
         }
         System.out.print("\r");
         for(int i=0; i<78; i++) System.out.print(' ');

@@ -55,7 +55,7 @@ public class MarchingCubes {
     }
 
     /** performs the Marching Cubes algorithm on a single cube */
     }
 
     /** performs the Marching Cubes algorithm on a single cube */

-    static void march(SampledField sampledField, Mesh mesh, double targetValue, double fX, double fY, double fZ, double fScale) {
+    static void march(SampledField sampledField, Mesh mesh, double threshold, double fX, double fY, double fZ, double fScale) {

         int iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags;
         double fOffset;
         GLvector sColor;
         int iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags;
         double fOffset;
         GLvector sColor;

@@ -74,11 +74,9 @@ public class MarchingCubes {
 
         // Find which vertices are inside of the surface and which are outside
         iFlagIndex = 0;
 
         // Find which vertices are inside of the surface and which are outside
         iFlagIndex = 0;

-        for(iVertexTest = 0; iVertexTest < 8; iVertexTest++) {
-            if (afCubeValue[iVertexTest] <= targetValue) {
+        for(iVertexTest = 0; iVertexTest < 8; iVertexTest++)
+            if (afCubeValue[iVertexTest] >= threshold)

                 iFlagIndex |= 1<<iVertexTest;
                 iFlagIndex |= 1<<iVertexTest;

-            }
-        }

         
         // Find which edges are intersected by the surface
         iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];
         
         // Find which edges are intersected by the surface
         iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];

@@ -92,7 +90,9 @@ public class MarchingCubes {
             // If there is an intersection on this edge
             if ((iEdgeFlags & (1<<iEdge))==0) continue;
             fOffset = fGetOffset(afCubeValue[ a2iEdgeConnection[iEdge][0] ], 
             // If there is an intersection on this edge
             if ((iEdgeFlags & (1<<iEdge))==0) continue;
             fOffset = fGetOffset(afCubeValue[ a2iEdgeConnection[iEdge][0] ], 

-                                 afCubeValue[ a2iEdgeConnection[iEdge][1] ], targetValue);
+                                 afCubeValue[ a2iEdgeConnection[iEdge][1] ],
+                                 threshold,
+                                 fScale * 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].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;

@@ -112,13 +112,18 @@ public class MarchingCubes {
                 iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner];
                 points[iCorner] = new Point(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY, asEdgeVertex[iVertex].fZ);
 
                 iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner];
                 points[iCorner] = new Point(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY, asEdgeVertex[iVertex].fZ);
 

-                // questionable
+                // questionable, but we do it anyways

                 norm = norm.plus(new Vec(asEdgeNorm[iVertex].fX,   asEdgeNorm[iVertex].fY,   asEdgeNorm[iVertex].fZ));
             }
                 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(), 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;
+
+            mesh.newT(points[0], points[1], points[2], norm.norm().times(-1));

         }
     }
 
         }
     }
 

@@ -162,9 +167,17 @@ public class MarchingCubes {
 
     // fGetOffset finds the approximate point of intersection of the surface
     // between two points with the values fValue1 and fValue2
 
     // 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;
         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;
     }
 
         return (fValueDesired - fValue1)/fDelta;
     }