//
//////////////////////////////////////////////////////////////////////////////////////
-public abstract class MarchingCubes {
-
- private int iEdgeFlags;
+public class MarchingCubes {
private static class GLvector {
public double fX;
public double fY;
- public double fZ;
- };
-
- GLvector sSourcePoint[] = new GLvector[3];
-
- public abstract double fSample(double x, double y, double z);
-
- //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);
- vNormalizeVector(rfNormal, rfNormal);
- }
-
- void vNormalizeVector(GLvector rfVectorResult, GLvector rfVectorSource) {
- double fOldLength;
- double fScale;
-
- fOldLength = Math.sqrt( (rfVectorSource.fX * rfVectorSource.fX) +
- (rfVectorSource.fY * rfVectorSource.fY) +
- (rfVectorSource.fZ * rfVectorSource.fZ) );
-
- if(fOldLength == 0.0) {
- rfVectorResult.fX = rfVectorSource.fX;
- rfVectorResult.fY = rfVectorSource.fY;
- rfVectorResult.fZ = rfVectorSource.fZ;
- } else {
- fScale = 1.0/fOldLength;
- rfVectorResult.fX = rfVectorSource.fX*fScale;
- rfVectorResult.fY = rfVectorSource.fY*fScale;
- rfVectorResult.fZ = rfVectorSource.fZ*fScale;
+ public double fZ;
+ public String toString() {
+ return "("+fX+","+fY+","+fZ+")";
}
- }
-
- // 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) {
- 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) {
+ /** 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;
- for(iX = 0; iX < iDataSetSize; iX++)
+ int initialTriangles = mesh.numTriangles();
+ for(iX = 0; iX < iDataSetSize; 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, " +
+ (mesh.numTriangles()-initialTriangles) + " triangles");
for(iY = 0; iY < iDataSetSize; iY++)
for(iZ = 0; iZ < iDataSetSize; iZ++)
- vMarchCube(mesh, targetValue, iX*fStepSize, iY*fStepSize, iZ*fStepSize, fStepSize);
+ march(sampledField, mesh, targetValue, 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");
+ System.out.println("done marching.");
}
- //vMarchCube performs the Marching Cubes algorithm on a single cube
- void vMarchCube(Mesh mesh, double targetValue, double fX, double fY, double fZ, double fScale) {
+ /** 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) {
int iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags;
double fOffset;
GLvector sColor;
GLvector asEdgeVertex[] = new GLvector[12];
GLvector asEdgeNorm[] = new GLvector[12];
- //Make a local copy of the values at the cube's corners
+ for(int i=0; i<asEdgeVertex.length; i++) asEdgeVertex[i] = new GLvector();
+ for(int i=0; i<asEdgeNorm.length; i++) asEdgeNorm[i] = new GLvector();
+
+ // 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
+ afCubeValue[iVertex] = sampledField.getSample(new Point(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)
+ if (afCubeValue[iVertexTest] <= targetValue) {
iFlagIndex |= 1<<iVertexTest;
+ }
}
- //Find which edges are intersected by the surface
+ // 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;
+ // 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
+ // 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);
- }
+ // 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);
+
+ 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(sampledField, 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
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];
- /*
- 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(
+ points[iCorner] = new Point(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY, asEdgeVertex[iVertex].fZ);
+
+ // questionable
+ 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);
+ }
+ }
+
+ /**
+ * 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(SampledField sampledField, 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));
+ rfNormal.fY =
+ sampledField.getSample(new Point(fX, fY-0.01, fZ)) -
+ sampledField.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));
+ vNormalizeVector(rfNormal, rfNormal);
+ }
+
+ static void vNormalizeVector(GLvector rfVectorResult, GLvector rfVectorSource) {
+ double fOldLength;
+ double fScale;
+
+ fOldLength = Math.sqrt( (rfVectorSource.fX * rfVectorSource.fX) +
+ (rfVectorSource.fY * rfVectorSource.fY) +
+ (rfVectorSource.fZ * rfVectorSource.fZ) );
+
+ if(fOldLength == 0.0) {
+ rfVectorResult.fX = rfVectorSource.fX;
+ rfVectorResult.fY = rfVectorSource.fY;
+ rfVectorResult.fZ = rfVectorSource.fZ;
+ } else {
+ fScale = 1.0/fOldLength;
+ rfVectorResult.fX = rfVectorSource.fX*fScale;
+ rfVectorResult.fY = rfVectorSource.fY*fScale;
+ rfVectorResult.fZ = rfVectorSource.fZ*fScale;
}
}
+ // 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) {
+ double fDelta = fValue2 - fValue1;
+ if(fDelta == 0.0) return 0.5;
+ return (fValueDesired - fValue1)/fDelta;
+ }
+
////////////////////////////////////////////////////////////////////////////////////////
// Tables //////////////////////////////////////////////////////////////////////////////
// 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,
// 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},
projects / anneal.git / commitdiff