X-Git-Url: http://git.megacz.com/?p=anneal.git;a=blobdiff_plain;f=src%2Fedu%2Fberkeley%2Fqfat%2Fvoxel%2FMarchingCubes.java;h=3112d8d163d46139892707a9a694e434b87c22cb;hp=cccbecf0e0eca749026ffa3116a34011d5a52fb7;hb=76e9f9a721c31fc8c2ed671ebda1640534ac2f1b;hpb=03c8f8e0d33bd7ea2940151ca474242bc9920315

diff --git a/src/edu/berkeley/qfat/voxel/MarchingCubes.java b/src/edu/berkeley/qfat/voxel/MarchingCubes.java
index cccbecf..3112d8d 100644
--- a/src/edu/berkeley/qfat/voxel/MarchingCubes.java
+++ b/src/edu/berkeley/qfat/voxel/MarchingCubes.java
@@ -23,30 +23,138 @@ import edu.berkeley.qfat.geom.*;
 //
 //////////////////////////////////////////////////////////////////////////////////////
 
-public abstract class MarchingCubes {
-
-    private int iEdgeFlags;
+public class MarchingCubes {
 
     private static class GLvector {
         public double fX;
         public double fY;
-        public double fZ;     
+        public double fZ;
+        public String toString() {
+            return "("+fX+","+fY+","+fZ+")";
+        }
     };
 
-    GLvector  sSourcePoint[] = new GLvector[3];
+    /** march iterates over the entire dataset, calling vMarchCube on each cube */
+    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();
+        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)voxels.getNumSamplesX()))*100) + "% marched, " +
+                             (mesh.numTriangles()-initialTriangles) + " triangles");
+            for(int iY = 0; iY < voxels.getNumSamplesY(); iY++)
+                for(int iZ = 0; iZ < voxels.getNumSamplesZ(); iZ++)
+                    march(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(' ');
+        System.out.print("\r");
+        System.out.println("done marching.");
+    }
+
+    /** performs the Marching Cubes algorithm on a single cube */
+    static void march(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;
+        double afCubeValue[] = new double[8];
+        GLvector asEdgeVertex[] = new GLvector[12];
+        GLvector asEdgeNorm[] = new GLvector[12];
+
+        for(int i=0; i<asEdgeVertex.length; i++) asEdgeVertex[i] = new GLvector();
+        for(int i=0; i<asEdgeNorm.length; i++) asEdgeNorm[i] = new GLvector();
 
-    public abstract double fSample(double x, double y, double z);
+        // Make a local copy of the values at the cube's corners
+        for(iVertex = 0; iVertex < 8; iVertex++)
+            afCubeValue[iVertex] = voxels.getSample(new Point(fX + a2fVertexOffset[iVertex][0]*scaleX,
+                                                                    fY + a2fVertexOffset[iVertex][1]*scaleY,
+                                                                    fZ + a2fVertexOffset[iVertex][2]*scaleZ));
 
-    //vGetNormal() finds the gradient of the scalar field at a point
-    //This gradient can be used as a very accurate vertx normal for lighting calculations
-    void vGetNormal(GLvector rfNormal, double fX, double fY, double fZ) {
-        rfNormal.fX = fSample(fX-0.01, fY, fZ) - fSample(fX+0.01, fY, fZ);
-        rfNormal.fY = fSample(fX, fY-0.01, fZ) - fSample(fX, fY+0.01, fZ);
-        rfNormal.fZ = fSample(fX, fY, fZ-0.01) - fSample(fX, fY, fZ+0.01);
+        // Find which vertices are inside of the surface and which are outside
+        iFlagIndex = 0;
+        for(iVertexTest = 0; iVertexTest < 8; iVertexTest++)
+            if (afCubeValue[iVertexTest] >= threshold)
+                iFlagIndex |= 1<<iVertexTest;
+        
+        // Find which edges are intersected by the surface
+        iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];
+        
+        // If the cube 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 < 12; iEdge++) {
+            // 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,
+                                 Math.min(Math.min(scaleX, scaleY), scaleZ) * 0.1);
+            
+            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(voxels, asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX, asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
+        }
+
+        // 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)
+                break;
+
+            Point[] points = new Point[3];
+            Vec norm = new Vec(0,0,0);
+            for(iCorner = 0; iCorner < 3; iCorner++) {
+                iVertex = a2iTriangleConnectionTable[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;
+
+            mesh.newT(points[0], points[1], points[2], norm.norm().times(-1));
+        }
+    }
+
+    /**
+     *  vGetNormal() finds the gradient of the scalar field at a point
+     *  This gradient can be used as a very accurate vertx normal for
+     *  lighting calculations
+     */
+    static void vGetNormal(VoxelData voxels, GLvector rfNormal, double fX, double fY, double fZ) {
+        rfNormal.fX =
+            voxels.getSample(new Point(fX-0.01, fY, fZ)) -
+            voxels.getSample(new Point(fX+0.01, fY, fZ));
+        rfNormal.fY =
+            voxels.getSample(new Point(fX, fY-0.01, fZ)) -
+            voxels.getSample(new Point(fX, fY+0.01, fZ));
+        rfNormal.fZ =
+            voxels.getSample(new Point(fX, fY, fZ-0.01)) -
+            voxels.getSample(new Point(fX, fY, fZ+0.01));
         vNormalizeVector(rfNormal, rfNormal);
     }
 
-    void vNormalizeVector(GLvector rfVectorResult, GLvector rfVectorSource) {
+    static void vNormalizeVector(GLvector rfVectorResult, GLvector rfVectorSource) {
         double fOldLength;
         double fScale;
         
@@ -68,81 +176,18 @@ public abstract class MarchingCubes {
 
     // fGetOffset finds the approximate point of intersection of the surface
     // between two points with the values fValue1 and fValue2
-    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;
-        return (fValueDesired - fValue1)/fDelta;
-    }
-
-    // vMarchingCubes iterates over the entire dataset, calling vMarchCube on each cube
-    void vMarchingCubes(Mesh mesh, double targetValue, int iDataSetSize, int fStepSize) {
-        int iX, iY, iZ;
-        for(iX = 0; iX < iDataSetSize; iX++)
-            for(iY = 0; iY < iDataSetSize; iY++)
-                for(iZ = 0; iZ < iDataSetSize; iZ++)
-                    vMarchCube(mesh, targetValue, iX*fStepSize, iY*fStepSize, iZ*fStepSize, fStepSize);
-    }
-
-    //vMarchCube performs the Marching Cubes algorithm on a single cube
-    void vMarchCube(Mesh mesh, double targetValue, double fX, double fY, double fZ, double fScale) {
-        int iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags;
-        double fOffset;
-        GLvector sColor;
-        double afCubeValue[] = new double[8];
-        GLvector asEdgeVertex[] = new GLvector[12];
-        GLvector asEdgeNorm[] = new GLvector[12];
-
-        //Make a local copy of the values at the cube's corners
-        for(iVertex = 0; iVertex < 8; iVertex++)
-            afCubeValue[iVertex] = fSample(fX + a2fVertexOffset[iVertex][0]*fScale,
-                                           fY + a2fVertexOffset[iVertex][1]*fScale,
-                                           fZ + a2fVertexOffset[iVertex][2]*fScale);
-        
-        //Find which vertices are inside of the surface and which are outside
-        iFlagIndex = 0;
-        for(iVertexTest = 0; iVertexTest < 8; iVertexTest++) {
-            if(afCubeValue[iVertexTest] <= targetValue) 
-                iFlagIndex |= 1<<iVertexTest;
-        }
-        
-        //Find which edges are intersected by the surface
-        iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];
-        
-        //If the cube 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 < 12; iEdge++) {
-            //if there is an intersection on this edge
-            if ((iEdgeFlags & (1<<iEdge))!=0) {
-                fOffset = fGetOffset(afCubeValue[ a2iEdgeConnection[iEdge][0] ], 
-                                     afCubeValue[ a2iEdgeConnection[iEdge][1] ], targetValue);
-                
-                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;
-                
-                vGetNormal(asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX, asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
-            }
-        }
+        // 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;
 
-        // 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)
-                break;
-
-            for(iCorner = 0; iCorner < 3; iCorner++) {
-                iVertex = a2iTriangleConnectionTable[iFlagIndex][3*iTriangle+iCorner];
-                /*
-                vGetColor(sColor, asEdgeVertex[iVertex], asEdgeNorm[iVertex]);
-                glColor3f(sColor.fX, sColor.fY, sColor.fZ);
-                glNormal3f(asEdgeNorm[iVertex].fX,   asEdgeNorm[iVertex].fY,   asEdgeNorm[iVertex].fZ);
-                glVertex3f(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY, asEdgeVertex[iVertex].fZ);
-                */
-                //mesh.newT(
-            }
-        }
+        return (fValueDesired - fValue1)/fDelta;
     }
 
     ////////////////////////////////////////////////////////////////////////////////////////
@@ -246,7 +291,7 @@ public abstract class MarchingCubes {
     //  possible vertex states There are 12 edges.  For each entry in
     //  the table, if edge #n is intersected, then bit #n is set to 1
 
-    int aiCubeEdgeFlags[] = {
+    static final int aiCubeEdgeFlags[] = {
         0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
         0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, 
         0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
@@ -293,7 +338,7 @@ public abstract class MarchingCubes {
     //  I found this table in an example program someone wrote long
     //  ago.  It was probably generated by hand
 
-    int a2iTriangleConnectionTable[][] = {
+    static final int a2iTriangleConnectionTable[][] = {
         {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
         {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
         {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},