/* $RCSfile$ * $Author: hansonr $ * $Date: 2007-03-30 11:40:16 -0500 (Fri, 30 Mar 2007) $ * $Revision: 7273 $ * * Copyright (C) 2007 Miguel, Bob, Jmol Development * * Contact: hansonr@stolaf.edu * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. */ package org.openscience.jvxl.simplewriter; import java.util.BitSet; import javax.vecmath.Point3i; import org.jmol.jvxl.data.JvxlCoder; import org.jmol.jvxl.data.VolumeData; import org.jmol.util.Logger; public class SimpleMarchingCubesOld { /* * SimpleMarchingCubesOld implements the original method of generating * data, using an array that is size [12] to hold each cube's edge vertex data. * and an array of size [nCubesY*nCubesZ][12] to hold a slice of cubes in memory. * Turns out it is slower by about 10% and consumes far more memory than the * new algorithm I wrote Feb 10, 2008. Saved here for posterity -- Bob Hanson * * timing: SimpleMarchingCubes with 100,100,100: * * getEdgeData: 641 ms * getEdgeData: 1625 ms * * old getEdgeData: 688 ms * old getEdgeData: 1672 ms * * An adaptation of Marching Cubes to include data slicing and the option * for progressive reading of the data. Associated SurfaceReader and VoxelData * structures are required to store the sequential values in the case of a plane * and to deliver the sequential vertex numbers in any case. * * Author: Bob Hanson, hansonr@stolaf.edu * * The "Simple" version does not create triangle data, * just the JVXL fractionData string * */ private VolumeData volumeData; private float cutoff; private boolean isCutoffAbsolute; private boolean isXLowToHigh; private StringBuffer fractionData = new StringBuffer(); private int cubeCountX, cubeCountY, cubeCountZ; private int nY, nZ; private BitSet bsVoxels = new BitSet(); public BitSet getBsVoxels() { return bsVoxels; } private int mode; private final static int MODE_CUBE = 1; private final static int MODE_BITSET = 2; private final static int MODE_GETXYZ = 3; private VoxelDataCreator vdc; public SimpleMarchingCubesOld(VoxelDataCreator vdc, VolumeData volumeData, float cutoff, boolean isCutoffAbsolute , boolean isXLowToHigh) { // when just creating a JVXL file all you really need are: // // volumeData.voxelData[x][y][z] // cutoff // this.vdc = vdc; this.volumeData = volumeData; this.cutoff = cutoff; this.isCutoffAbsolute = isCutoffAbsolute; this.isXLowToHigh = isXLowToHigh; if (vdc == null) { mode = MODE_CUBE; } else { mode = MODE_GETXYZ; } cubeCountX = volumeData.voxelCounts[0] - 1; cubeCountY = (nY = volumeData.voxelCounts[1]) - 1; cubeCountZ = (nZ = volumeData.voxelCounts[2]) - 1; yzCount = nY * nZ; setLinearOffsets(); } private final float[] vertexValues = new float[8]; private final Point3i[] vertexPoints = new Point3i[8]; { for (int i = 8; --i >= 0;) vertexPoints[i] = new Point3i(); } int edgeCount; /* Note to Jason from Bob: * * To just create a JVXL file, you need these five methods. * Their output is the fractionData string buffer and the * number of surface points * * inputs required: * * 1) volumeData.voxelData[x][y][z] * 2) cutoff * 3) values created in MarchingCubes constructor * * The first four methods are in org.jmol.jvxl.calc.MarchingCubes.java * * generateSurfaceData -- isXLowToHigh false; isContoured false * -- triangle stuff at end not needed * propagateNeighborPointIndexes -- EXACTLY as is, no changes allowed * isInside -- EXACTLY as is -- defines what "inside" means * processOneCubical -- EXACTLY as is, no changes at all * SurfaceReader.getSurfacePointIndex -- your job * -- receives the point value data and positions * -- responsible for creating the fractionData character buffer * -- just return 0 since you are not creating triangles * */ private static int[] xyPlanePts = new int[] { 0, 1, 1, 0, 0, 1, 1, 0 }; public String getEdgeData() { Logger.startTimer(); // set up the set of edge points in the YZ plane // isoPointIndixes are indices into an array of Point3f values // They will be initialized as -1 whenever a vertex is needed. // But if just creating a JVXL file, all you need to do // is set them to 0, not an index into any actual array. int[][] isoPointIndexes = new int[cubeCountY * cubeCountZ][12]; float[][] xyPlanes = (mode == MODE_GETXYZ ? new float[2][yzCount] : null); int x0, x1, xStep, ptStep, pt, ptX; if (isXLowToHigh) { x0 = 0; x1 = cubeCountX; xStep = 1; ptStep = yzCount; pt = ptX = (yzCount - 1) - nZ - 1; // we are starting at the top corner, in the next to last // cell on the next to last row of the first plane } else { x0 = cubeCountX - 1; x1 = -1; xStep = -1; ptStep = -yzCount; pt = ptX = (cubeCountX * yzCount - 1) - nZ - 1; // we are starting at the top corner, in the next to last // cell on the next to last row of the next to last plane(!) } int cellIndex0 = cubeCountY * cubeCountZ - 1; int cellIndex = cellIndex0; for (int x = x0; x != x1; x += xStep, ptX += ptStep, pt = ptX, cellIndex = cellIndex0) { if (mode == MODE_GETXYZ) { float[] plane = xyPlanes[0]; xyPlanes[0] = xyPlanes[1]; xyPlanes[1] = plane; } for (int y = cubeCountY; --y >= 0; pt--) { for (int z = cubeCountZ; --z >= 0; pt--, cellIndex--) { // set up the list of indices that need checking int[] voxelPointIndexes = propagateNeighborPointIndexes(x, y, z, pt, isoPointIndexes, cellIndex); // create the bitset mask indicating which vertices are inside. // 0xFF here means "all inside"; 0x00 means "all outside" int insideMask = 0; for (int i = 8; --i >= 0;) { // cubeVertexOffsets just gets us the specific grid point relative // to our base x,y,z cube position boolean isInside; Point3i offset = cubeVertexOffsets[i]; int pti = pt + linearOffsets[i]; switch (mode) { case MODE_GETXYZ: vertexValues[i] = getValue(i, x + offset.x, y + offset.y, z + offset.z, pti, xyPlanes[xyPlanePts[i]]); isInside = bsVoxels.get(pti); break; case MODE_BITSET: isInside = bsVoxels.get(pti); vertexValues[i] = (isInside ? 1 : 0); break; default: case MODE_CUBE: vertexValues[i] = volumeData.voxelData[x + offset.x][y + offset.y][z + offset.z]; isInside = isInside(vertexValues[i], cutoff, isCutoffAbsolute); if (isInside) bsVoxels.set(pti); } if (isInside) { insideMask |= 1 << i; } } if (insideMask == 0) { continue; } if (insideMask == 0xFF) { continue; } // This cube is straddling the cutoff. We must check all edges processOneCubical(insideMask, voxelPointIndexes, x, y, z, pt); } } } Logger.checkTimer("old getEdgeData"); return fractionData.toString(); } public static boolean isInside(float voxelValue, float max, boolean isAbsolute) { return ((max > 0 && (isAbsolute ? Math.abs(voxelValue) : voxelValue) >= max) || (max <= 0 && voxelValue <= max)); } BitSet bsValues = new BitSet(); private float getValue(@SuppressWarnings("unused") int i, int x, int y, int z, int pt, float[] tempValues) { //if (bsValues.get(pt)) //return tempValues[pt % yzCount]; bsValues.set(pt); float value = vdc.getValue(x, y, z); tempValues[pt % yzCount] = value; //System.out.println("xyz " + x + " " + y + " " + z + " v=" + value); if (isInside(value, cutoff, isCutoffAbsolute)) bsVoxels.set(pt); return value; } private final int[] nullNeighbor = { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }; private int[] propagateNeighborPointIndexes(int x, int y, int z, @SuppressWarnings("unused") int pt, int[][] isoPointIndexes, int cellIndex) { /* * * We are running through the grid points in yz planes from high x --> low x * and within those planes along strips from high y to low y * and within those strips, from high z to low z. The "leading vertex" is 0, * and the "leading edges" are {0,3,8}. * * For each such cube, edges are traversed from high to low (11-->0) * * Each edge has the potential to be "critical" and cross the surface. * Setting -1 in voxelPointIndexes indicates that this edge needs checking. * Otherwise, the crossing point for this edge is taken from the value * already determined, because it has already been determined to be critical. * * The above model, because it starts at HIGH x, requires that all x,y,z points * be in memory from the beginning. We could have instead used a progressive * streaming model, where we only pull in the slice of data that we need. In * that case, each edge corresponds to a specific pair of indices in our slice. * * Say we have a 51 x 11 x 21 block of data. This represents a 50 x 10 x 20 set * of cubes. If, instead of reading all the data, we pull in just the first two * "slices" x=0(10x20), x=1(10x20), that is just 400 points. Once a slice of * data is used, we can flush it -- it is never used again. * * When color mapping, we can do the same thing; we just have to put the verticies * into bins based on which pair of slices will be relevant, and then make sure we * process the verticies based on these bins. * * The JVXL format depends on a specific order of reading of the edge data. The * progressive model completely messes this up. The vertices will be read in the * same order around the cube, but the "leading edges" will be {0,1,9}, not {0,3,8}. * We do know which edge is which, so we could construct a progressive model from * a nonprogressive one, if necessary. * * All we are really talking about is the JVXL reader, because we can certainly * switch to progressive mode in all the other readers. * * Y * 4 --------4--------- 5 * /| /| * / | / | * / | / | * 7 8 5 | * / | / 9 * / | / | * 7 --------6--------- 6 | * | | | | * | 0 ---------0--|----- 1 X * | / | / * 11 / 10 / * | 3 | 1 * | / | / * | / | / * 3 ---------2-------- 2 * Z * * */ /* DO NOT EVER CHANGE THIS */ int[] voxelPointIndexes = isoPointIndexes[cellIndex]; boolean noYNeighbor = (y == cubeCountY - 1); int[] yNeighbor = noYNeighbor ? nullNeighbor : isoPointIndexes[cellIndex + cubeCountZ]; boolean noZNeighbor = (z == cubeCountZ - 1); int[] zNeighbor = noZNeighbor ? nullNeighbor : isoPointIndexes[cellIndex + 1]; voxelPointIndexes[0] = -1; voxelPointIndexes[2] = zNeighbor[0]; voxelPointIndexes[4] = yNeighbor[0]; voxelPointIndexes[6] = (noYNeighbor ? zNeighbor[4] : yNeighbor[2]); if (isXLowToHigh) { // reading x from low to high if (x == 0) { voxelPointIndexes[3] = -1; voxelPointIndexes[8] = -1; voxelPointIndexes[7] = yNeighbor[3]; voxelPointIndexes[11] = zNeighbor[8]; } else { voxelPointIndexes[3] = voxelPointIndexes[1]; voxelPointIndexes[7] = voxelPointIndexes[5]; voxelPointIndexes[8] = voxelPointIndexes[9]; voxelPointIndexes[11] = voxelPointIndexes[10]; } voxelPointIndexes[1] = -1; voxelPointIndexes[5] = yNeighbor[1]; voxelPointIndexes[9] = -1; voxelPointIndexes[10] = zNeighbor[9]; } else { // reading x from high to low if (x == cubeCountX - 1) { voxelPointIndexes[1] = -1; voxelPointIndexes[5] = yNeighbor[1]; voxelPointIndexes[9] = -1; voxelPointIndexes[10] = zNeighbor[9]; } else { voxelPointIndexes[1] = voxelPointIndexes[3]; voxelPointIndexes[5] = voxelPointIndexes[7]; voxelPointIndexes[9] = voxelPointIndexes[8]; voxelPointIndexes[10] = voxelPointIndexes[11]; } voxelPointIndexes[3] = -1; voxelPointIndexes[7] = yNeighbor[3]; voxelPointIndexes[8] = -1; voxelPointIndexes[11] = zNeighbor[8]; } return voxelPointIndexes; } private static final int[] Pwr2 = new int[] { 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 }; private boolean processOneCubical(int insideMask, int[] voxelPointIndexes, @SuppressWarnings("unused") int x, @SuppressWarnings("unused") int y, @SuppressWarnings("unused") int z, @SuppressWarnings("unused") int pt) { // the key to the algorithm is that we have a catalog that // maps the inside-vertex mask to an edge mask. int edgeMask = insideMaskTable[insideMask]; //for (int i =0; i < 8; i++) System.out.print("\nvpi for cell " + pt + ": vertex " + i + ": " + voxelPointIndexes[i] + " " + Integer.toBinaryString(edgeMask)); boolean isNaN = false; for (int iEdge = 12; --iEdge >= 0;) { // bit set to one means it's a relevant edge if ((edgeMask & Pwr2[iEdge]) == 0) continue; // if we have a point already, we don't need to check this edge. // for triangles, this will be an index into an array; // for just creating JVXL files, this can just be 0 if (voxelPointIndexes[iEdge] >= 0) continue; // propagated from neighbor // here's an edge that has to be checked. // get the vertex numbers 0 - 7 int vertexA = edgeVertexes[iEdge << 1]; int vertexB = edgeVertexes[(iEdge << 1) + 1]; // pick up the actual value at each vertex // this array of 8 values is updated as we go. float valueA = vertexValues[vertexA]; float valueB = vertexValues[vertexB]; // we allow for NaN values -- missing triangles if (Float.isNaN(valueA) || Float.isNaN(valueB)) isNaN = true; // the exact point position -- not important for just // creating the JVXL file. In that case, all you // need are the two values valueA and valueB and the cutoff. // from those you can define the fractional offset // here is where we get the value and assign the point for that edge // it is where the JVXL surface data line is appended voxelPointIndexes[iEdge] = edgeCount++; //System.out.println(" pt=" + pt + " edge" + iEdge + " xyz " + x + " " + y + " " + z + " vertexAB=" + vertexA + " " + vertexB + " valueAB=" + valueA + " " + valueB + " f= " + (cutoff - valueA) / (valueB - valueA)); fractionData.append(JvxlCoder.jvxlFractionAsCharacter((cutoff - valueA) / (valueB - valueA))); } return !isNaN; } final static Point3i[] cubeVertexOffsets = { new Point3i(0, 0, 0), //0 pt new Point3i(1, 0, 0), //1 pt + yz new Point3i(1, 0, 1), //2 pt + yz + 1 new Point3i(0, 0, 1), //3 pt + 1 new Point3i(0, 1, 0), //4 pt + z new Point3i(1, 1, 0), //5 pt + yz + z new Point3i(1, 1, 1), //6 pt + yz + z + 1 new Point3i(0, 1, 1) //7 pt + z + 1 }; private final int[] linearOffsets = new int[8]; int yzCount; /* set the linear offsets for unique cell ID * and for pointing into the inside/outside BitSet. * Add offset to 0: x * (nY * nZ) + y * nZ + z */ void setLinearOffsets() { linearOffsets[0] = 0; linearOffsets[1] = yzCount; linearOffsets[2] = yzCount + 1; linearOffsets[3] = 1; linearOffsets[4] = nZ; linearOffsets[5] = yzCount + nZ; linearOffsets[6] = yzCount + nZ + 1; linearOffsets[7] = nZ + 1; } public int getLinearOffset(int x, int y, int z, int offset) { return x * yzCount + y * nZ + z + linearOffsets[offset]; } /* Y * 4 --------4--------- 5 +z --------4--------- +yz+z * /| /| /| /| * / | / | / | / | * / | / | / | / | * 7 8 5 | 7 8 5 | * / | / 9 / | / 9 * / | / | / | / | * 7 --------6--------- 6 | +z+1 --------6--------- +yz+z+1| * | | | | | | | | * | 0 ---------0--|----- 1 X | 0 ---------0--|----- +yz X(outer) * | / | / | / | / * 11 / 10 / 11 / 10 / * | 3 | 1 | 3 | 1 * | / | / | / | / * | / | / | / | / * 3 ---------2-------- 2 +1 ---------2-------- +yz+1 * Z Z (inner) * * streaming data offsets * type 0: x-edges: 0 2 4 6 * type 1: y-edges: 8 9 10 11 * type 2: z-edges: 1 3 5 7 * * Data stream offsets for vertices, relative to point 0, based on reading * loops {for x {for y {for z}}} 0-->n-1 * y and z are numbers of grid points in those directions: * * 0 1 2 3 4 5 6 7 * 0 +yz +yz+1 +1 +z +yz+z +yz+z+1 +z+1 * * These are just looked up in a table. After the first set of cubes, * we are only adding points 1, 2, 5 or 6. This means that initially * we need two data slices, but after that only one (slice 1): * * base * offset 0 1 2 3 4 5 6 7 * slice[0] 0 +1 +z +z+1 * slice[1] +yz 0 +1 +z +z+1 * * slice: 0 1 1 0 0 1 1 0 * * We can request reading of two slices (2*nY*nZ data points) first, then * from then on, just nY*nZ points. "Reading" is really just being handed a * pointer into an array. Perhaps that array is already filled completely; * perhaps it is being read incrementally. * * As it is now, the JVXL data are just read into an [nX][nY][nZ] array anyway, * so we can continue to do that with NON progressive files. */ private final static byte edgeVertexes[] = { 0, 1, 1, 2, 2, 3, 3, 0, 4, 5, /*0 1 2 3 4 */ 5, 6, 6, 7, 7, 4, 0, 4, 1, 5, 2, 6, 3, 7 }; /*5 6 7 8 9 10 11 */ private final static short insideMaskTable[] = { 0x0000, 0x0109, 0x0203, 0x030A, 0x0406, 0x050F, 0x0605, 0x070C, 0x080C, 0x0905, 0x0A0F, 0x0B06, 0x0C0A, 0x0D03, 0x0E09, 0x0F00, 0x0190, 0x0099, 0x0393, 0x029A, 0x0596, 0x049F, 0x0795, 0x069C, 0x099C, 0x0895, 0x0B9F, 0x0A96, 0x0D9A, 0x0C93, 0x0F99, 0x0E90, 0x0230, 0x0339, 0x0033, 0x013A, 0x0636, 0x073F, 0x0435, 0x053C, 0x0A3C, 0x0B35, 0x083F, 0x0936, 0x0E3A, 0x0F33, 0x0C39, 0x0D30, 0x03A0, 0x02A9, 0x01A3, 0x00AA, 0x07A6, 0x06AF, 0x05A5, 0x04AC, 0x0BAC, 0x0AA5, 0x09AF, 0x08A6, 0x0FAA, 0x0EA3, 0x0DA9, 0x0CA0, 0x0460, 0x0569, 0x0663, 0x076A, 0x0066, 0x016F, 0x0265, 0x036C, 0x0C6C, 0x0D65, 0x0E6F, 0x0F66, 0x086A, 0x0963, 0x0A69, 0x0B60, 0x05F0, 0x04F9, 0x07F3, 0x06FA, 0x01F6, 0x00FF, 0x03F5, 0x02FC, 0x0DFC, 0x0CF5, 0x0FFF, 0x0EF6, 0x09FA, 0x08F3, 0x0BF9, 0x0AF0, 0x0650, 0x0759, 0x0453, 0x055A, 0x0256, 0x035F, 0x0055, 0x015C, 0x0E5C, 0x0F55, 0x0C5F, 0x0D56, 0x0A5A, 0x0B53, 0x0859, 0x0950, 0x07C0, 0x06C9, 0x05C3, 0x04CA, 0x03C6, 0x02CF, 0x01C5, 0x00CC, 0x0FCC, 0x0EC5, 0x0DCF, 0x0CC6, 0x0BCA, 0x0AC3, 0x09C9, 0x08C0, 0x08C0, 0x09C9, 0x0AC3, 0x0BCA, 0x0CC6, 0x0DCF, 0x0EC5, 0x0FCC, 0x00CC, 0x01C5, 0x02CF, 0x03C6, 0x04CA, 0x05C3, 0x06C9, 0x07C0, 0x0950, 0x0859, 0x0B53, 0x0A5A, 0x0D56, 0x0C5F, 0x0F55, 0x0E5C, 0x015C, 0x0055, 0x035F, 0x0256, 0x055A, 0x0453, 0x0759, 0x0650, 0x0AF0, 0x0BF9, 0x08F3, 0x09FA, 0x0EF6, 0x0FFF, 0x0CF5, 0x0DFC, 0x02FC, 0x03F5, 0x00FF, 0x01F6, 0x06FA, 0x07F3, 0x04F9, 0x05F0, 0x0B60, 0x0A69, 0x0963, 0x086A, 0x0F66, 0x0E6F, 0x0D65, 0x0C6C, 0x036C, 0x0265, 0x016F, 0x0066, 0x076A, 0x0663, 0x0569, 0x0460, 0x0CA0, 0x0DA9, 0x0EA3, 0x0FAA, 0x08A6, 0x09AF, 0x0AA5, 0x0BAC, 0x04AC, 0x05A5, 0x06AF, 0x07A6, 0x00AA, 0x01A3, 0x02A9, 0x03A0, 0x0D30, 0x0C39, 0x0F33, 0x0E3A, 0x0936, 0x083F, 0x0B35, 0x0A3C, 0x053C, 0x0435, 0x073F, 0x0636, 0x013A, 0x0033, 0x0339, 0x0230, 0x0E90, 0x0F99, 0x0C93, 0x0D9A, 0x0A96, 0x0B9F, 0x0895, 0x099C, 0x069C, 0x0795, 0x049F, 0x0596, 0x029A, 0x0393, 0x0099, 0x0190, 0x0F00, 0x0E09, 0x0D03, 0x0C0A, 0x0B06, 0x0A0F, 0x0905, 0x080C, 0x070C, 0x0605, 0x050F, 0x0406, 0x030A, 0x0203, 0x0109, 0x0000 }; }