};
/** 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) {
- 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, targetValue, 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(' ');
}
/** 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 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;
// 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] <= targetValue) {
+ 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 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,
+ 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)
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));
}
- 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;
+
+ 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);
}
}
* 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);
}
// 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;
}