1 /* $RCSfile$ 2 * $Author: hansonr $ 3 * $Date: 2007-03-30 11:40:16 -0500 (Fri, 30 Mar 2007) $ 4 * $Revision: 7273 $ 5 * 6 * Copyright (C) 2007 Miguel, Bob, Jmol Development 7 * 8 * Contact: hansonr@stolaf.edu 9 * 10 * This library is free software; you can redistribute it and/or 11 * modify it under the terms of the GNU Lesser General Public 12 * License as published by the Free Software Foundation; either 13 * version 2.1 of the License, or (at your option) any later version. 14 * 15 * This library is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 * Lesser General License for more details. 19 * 20 * You should have received a copy of the GNU Lesser General Public 21 * License along with this library; if not, write to the Free Software 22 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 23 */ 24 package org.jmol.jvxl.calc; 25 26 27 28 import javajs.util.BS; 29 import org.jmol.jvxl.api.VertexDataServer; 30 import org.jmol.jvxl.data.JvxlCoder; 31 import org.jmol.jvxl.data.VolumeData; 32 import org.jmol.jvxl.readers.Parameters; 33 34 import javajs.util.SB; 35 import javajs.util.P3; 36 import javajs.util.P3i; 37 import javajs.util.P4; 38 import org.jmol.util.TriangleData; 39 import javajs.util.V3; 40 41 public class MarchingCubes extends TriangleData { 42 43 /* 44 * An adaptation of Marching Cubes that includes data slicing. 45 * Associated SurfaceReader and VoxelData structures are required 46 * to store the sequential values in the case of a plane 47 * and to deliver the sequential vertex numbers in any case. 48 * 49 * Author: Bob Hanson, hansonr@stolaf.edu 50 * 51 * inputs: surfaceReader: interface to other methods 52 * volumeData, containing all information relating to origin, axes, etc. 53 * 54 * MODE_PLANE volume data will be read progressively, by yz planes, starting from x = 0. 55 * This is a very efficient way to read the data -- only two active planes, 56 * which are swapped. 57 * MODE_CUBE volumeData.voxelData possibly with 3D voxel data 58 * MODE_JVXL bsVoxels -- optional alternative, JVXL-type BitSet indicating 59 * which vertices are inside and which outside 60 * in the order for(x 0 to nX){for(y 0 to nY){for(z 0 to nZ){}}} 61 * 62 * surfaceReader.getSurfacePointIndex gets the actual Point3f vertex points and 63 * returns the fraction information. It is exported here, because the JVXL reader may use this 64 * as an opportunity to deliver the fractional information, thus providing the 65 * actual position of the vertex point on the isosurface. 66 * 67 * outputs: Point3f positions of isosurface intersections with grid via SurfaceReader.addVertex() 68 * Point4f triangle data via surfaceReader.addTriangleCheck() 69 * bsVoxels -- returned same (JVXL file data) or filled (otherwise) 70 * edgeData -- encoded fraction data as a string 71 * 72 * bsExcludedVertices -- an option to exclude vertices based on having NaN values 73 * bsExcludedTriangles -- an option to exclude triangles based on position in space 74 * 75 */ 76 77 protected VertexDataServer surfaceReader; 78 protected VolumeData volumeData; 79 protected int contourType; 80 protected boolean isContoured; 81 protected float cutoff; 82 protected boolean isCutoffAbsolute; 83 protected boolean isSquared; 84 protected boolean isXLowToHigh; 85 86 protected int cubeCountX, cubeCountY, cubeCountZ; 87 protected int nY, nZ; 88 protected int yzCount; 89 90 protected boolean colorDensity; 91 protected boolean integrateSquared = true; 92 public BS bsVoxels; 93 protected BS bsExcludedVertices; 94 protected BS bsExcludedTriangles; 95 protected BS bsExcludedPlanes; 96 97 protected SB edgeData = new SB(); 98 99 private boolean excludePartialCubes = true; // original way 100 MarchingCubes()101 public MarchingCubes() { 102 // as triangleServer 103 } 104 MarchingCubes(VertexDataServer surfaceReader, VolumeData volumeData, Parameters params, BS bsVoxels)105 public MarchingCubes(VertexDataServer surfaceReader, VolumeData volumeData, 106 Parameters params, BS bsVoxels) { 107 108 // If just creating a JVXL file, see org.openscience.jmol.jvxl.simplewriter.SimpleMarchingCubes.java 109 // 110 111 // setting this false could upset reading 112 // older Jmol version files -- will need a flag IN the file for this if we do it 113 excludePartialCubes = true; 114 115 this.surfaceReader = surfaceReader; 116 this.bsVoxels = bsVoxels; 117 BS[] bsExcluded = params.bsExcluded; 118 bsExcludedVertices = (bsExcluded[0] == null ? bsExcluded[0] = new BS() : bsExcluded[0]); 119 bsExcludedPlanes = (bsExcluded[2] == null ? bsExcluded[2] = new BS() : bsExcluded[2]); 120 bsExcludedTriangles = (bsExcluded[3] == null ? bsExcluded[3] = new BS() : bsExcluded[3]); 121 // TODO -- need not be setting up planes for simple value-readers 122 123 mode = (volumeData.getVoxelData() != null || volumeData.mappingPlane != null ? MODE_CUBE 124 : bsVoxels != null ? MODE_JVXL : MODE_PLANES); 125 setParameters(volumeData, params); 126 } 127 setParameters(VolumeData volumeData, Parameters params)128 protected void setParameters(VolumeData volumeData, Parameters params) { 129 this.volumeData = volumeData; 130 colorDensity = params.colorDensity; 131 isContoured = params.thePlane == null && params.isContoured && !colorDensity; 132 cutoff = params.cutoff; 133 isCutoffAbsolute = params.isCutoffAbsolute; 134 contourType = params.contourType; 135 isSquared = params.isSquared; 136 isXLowToHigh = params.isXLowToHigh; 137 138 cubeCountX = volumeData.voxelCounts[0] - 1; 139 cubeCountY = volumeData.voxelCounts[1] - 1; 140 cubeCountZ = volumeData.voxelCounts[2] - 1; 141 volumeData.getYzCount(); 142 if (params.mapLattice != null) { 143 // extend plane out to full lattice, if specified 144 cubeCountX *= Math.abs(params.mapLattice.x); 145 cubeCountY *= Math.abs(params.mapLattice.y); 146 cubeCountZ *= Math.abs(params.mapLattice.z); 147 } 148 nY = cubeCountY + 1; 149 nZ = cubeCountZ + 1; 150 yzCount = nY * nZ; 151 if (bsVoxels == null) 152 bsVoxels = new BS(); 153 edgeVertexPointers = (isXLowToHigh ? edgeVertexPointersLowToHigh : edgeVertexPointersHighToLow); 154 edgeVertexPlanes = (isXLowToHigh ? edgeVertexPlanesLowToHigh : edgeVertexPlanesHighToLow); 155 isoPointIndexPlanes = new int[2][yzCount][3]; 156 yzPlanes = (mode == MODE_PLANES ? new float[2][yzCount] : null); 157 setLinearOffsets(); 158 calcVoxelVertexVectors(); 159 } 160 161 protected int mode; 162 protected final static int MODE_CUBE = 1; 163 protected final static int MODE_JVXL = 2; 164 protected final static int MODE_PLANES = 3; 165 166 protected final float[] vertexValues = new float[8]; 167 168 protected int edgeCount; 169 170 protected final V3[] voxelVertexVectors = new V3[8]; 171 protected final V3[] edgeVectors = new V3[12]; 172 { 173 for (int i = 12; --i >= 0;) 174 edgeVectors[i] = new V3(); 175 } 176 calcVoxelVertexVectors()177 protected void calcVoxelVertexVectors() { 178 for (int i = 8; --i >= 0;) 179 volumeData.transform(cubeVertexVectors[i], 180 voxelVertexVectors[i] = new V3()); 181 for (int i = 12; --i >= 0;) 182 edgeVectors[i].sub2(voxelVertexVectors[edgeVertexes[i + i + 1]], 183 voxelVertexVectors[edgeVertexes[i + i]]); 184 } 185 186 /* see also org.openscience.jmol.jvxl.simplewriter.SimpleMarchingCubes.java 187 * for a streamlined version of this method that does simple writing 188 * of JVXL files from cube data. 189 * 190 */ 191 192 protected static int[] yzPlanePts = new int[] { 193 0, 1, 1, 0, 194 0, 1, 1, 0 195 }; 196 protected final int[] edgePointIndexes = new int[12]; 197 protected int[][][] isoPointIndexPlanes; 198 protected float[][] yzPlanes; 199 resetIndexPlane(int[][] plane)200 protected int[][] resetIndexPlane(int[][] plane) { 201 for (int i = 0; i < yzCount; i++) 202 for (int j = 0; j < 3; j++) 203 plane[i][j] = Integer.MIN_VALUE; 204 return plane; 205 } 206 207 private P4 mappingPlane; 208 private boolean allInside; 209 private boolean isInside; 210 private P3i offset; 211 private float[][][] voxelData; 212 getEdgeData()213 public String getEdgeData() { 214 215 if (cubeCountX < 0 || cubeCountY < 0 || cubeCountZ < 0) 216 return ""; 217 mappingPlane = volumeData.mappingPlane; 218 219 // Logger.startTimer(); 220 221 /* 222 * The (new, Jmol 11.7.26) Marching Cubes code creates the 223 * isoPointIndexes[2][nY * nZ][3] array that holds two slices of edge data. 224 * Each edge is assigned a specific vertex, such that each vertex may have 225 * up to 3 associated edges. 226 * 227 * Feb 10, 2009 -- Bob Hanson 228 */ 229 230 //int insideCount = 0, outsideCount = 0, surfaceCount = 0; 231 edgeCount = 0; 232 233 int x0, x1, xStep, ptStep, pt, ptX; 234 if (isXLowToHigh) { 235 // used for progressive plane readers 236 // we are starting at the top corner, in the next to last 237 // cell on the next to last row of the first plane 238 x0 = 0; 239 x1 = cubeCountX + (colorDensity ? 1 : 0); 240 if (colorDensity) { 241 x1 = cubeCountX + 1; 242 ptX = yzCount - 1; 243 } else { 244 x1 = cubeCountX; 245 ptX = (yzCount - 1) - nZ - 1; 246 } 247 xStep = 1; 248 ptStep = yzCount; 249 } else { 250 // NOTE: isXLowToHigh FALSE is not compatible with progressive plane reading. 251 // One should never need to do that anyway, because one can always 252 // define the cube axes such that they go "from low to high in X" 253 // we are starting at the top corner, in the next to last 254 // cell on the next to last row of the next to last plane(!) 255 if (colorDensity) { 256 x0 = cubeCountX; 257 ptX = (cubeCountX + 1) * yzCount - 1; 258 } else { 259 x0 = cubeCountX - 1; 260 ptX = (cubeCountX * yzCount - 1) - nZ - 1; 261 } 262 x1 = -1; 263 xStep = -1; 264 ptStep = -yzCount; 265 } 266 pt = ptX; 267 resetIndexPlane(isoPointIndexPlanes[1]); 268 voxelData = null; 269 int y1 = cubeCountY + (colorDensity ? 1 : 0); 270 int z1 = cubeCountZ + (colorDensity ? 1 : 0); 271 switch (mode) { 272 case MODE_PLANES: 273 getPlane(x0, false); 274 break; 275 case MODE_CUBE: 276 voxelData = volumeData.getVoxelData(); 277 break; 278 } 279 allInside = (colorDensity && (cutoff == 0 280 || mode == MODE_JVXL && bsVoxels.nextSetBit(0) < 0)); 281 boolean colorDensityAll = (colorDensity && cutoff == 0); 282 float v = 0; 283 for (int x = x0; x != x1; x += xStep, ptX += ptStep, pt = ptX) { 284 285 // we swap planes of grid data when 286 // obtaining the grid data point by point 287 288 if (mode == MODE_PLANES) { 289 // for a progressive reader, we read the next two planes 290 // for x = 0, 2, 4, 6... 291 // say we have 50 points; then we have 49 cubes 292 // we definely what planes 0-49, but we will be processing those 293 // [0 1] [1 2] [2 3] ... [48 49] 294 // so our last x is 48, and there is no problem with (x + xStep) being == x1. 295 // fixed in Jmol 12.1.51 -- was "x + xStep != x1" 296 if (x + xStep <= x1) 297 getPlane(x + xStep, true); 298 } 299 300 // now scan the plane of cubicals 301 302 if (bsExcludedPlanes.get(x) && bsExcludedPlanes.get(x + xStep)) 303 continue; 304 305 if (colorDensity) { 306 // MODE_JVXL not allowed here 307 for (int y = y1; --y >= 0;) 308 for (int z = z1; --z >= 0; pt--) { 309 // 0 cutoff read as "show grid points only" 310 v = getValue(x, y, z, pt, 0); 311 if (colorDensityAll || isInside) { 312 // xyz corner is inside, so add this point 313 // TODO - we are missing the outside edges of these 314 // in the Y and Z direction 315 addVertex(x, y, z, pt, v); 316 } 317 } 318 continue; 319 } 320 321 // we swap the edge vertex index planes 322 323 int[][] indexPlane = isoPointIndexPlanes[0]; 324 isoPointIndexPlanes[0] = isoPointIndexPlanes[1]; 325 isoPointIndexPlanes[1] = resetIndexPlane(indexPlane); 326 327 boolean noValues = true; 328 for (int y = y1; --y >= 0; pt--) { 329 for (int z = z1; --z >= 0; pt--) { 330 int insideMask = 0; 331 // create the bitset mask indicating which vertices are inside. 332 // 0xFF here means "all inside"; 0x00 means "all outside" 333 for (int i = 8; --i >= 0;) { 334 v = getValue(x, y, z, pt, i); 335 if (isInside) 336 insideMask |= Pwr2[i]; 337 } 338 // last value for i is 0, the x,y,z point itself 339 // so we need only check it for noValues -- on THIS plane 340 if (noValues && !Float.isNaN(v)) 341 noValues = false; 342 if (insideMask == 0) { 343 //++outsideCount; 344 continue; 345 } 346 if (insideMask == 0xFF) { 347 //++insideCount; 348 continue; 349 } 350 //++surfaceCount; 351 352 // This cube is straddling the cutoff. We must check all edges 353 // Note that we do not process it if it has an NaN values 354 if (processOneCubical(insideMask, x, y, z, pt) && !isContoured 355 && !colorDensity) { 356 processTriangles(insideMask); 357 } 358 } 359 } 360 if (noValues) { 361 bsExcludedPlanes.set(x); 362 } 363 } 364 365 return edgeData.toString(); 366 } 367 getValue(int x, int y, int z, int pt, int i)368 private float getValue(int x, int y, int z, int pt, int i) { 369 float v; 370 371 // cubeVertexOffsets just gets us 372 // the specific grid point relative 373 // to our base x,y,z cube position 374 375 offset = cubeVertexOffsets[i]; 376 int pti = pt + linearOffsets[i]; 377 switch (mode) { 378 case MODE_PLANES: 379 v = vertexValues[i] = getValueArray(x + offset.x, y + offset.y, z 380 + offset.z, pti, yzPlanes[yzPlanePts[i]]); 381 isInside = (allInside || bsVoxels.get(pti)); 382 break; 383 case MODE_JVXL: 384 isInside = (allInside || bsVoxels.get(pti)); 385 v = vertexValues[i] = (bsExcludedVertices.get(pti) ? Float.NaN 386 : isInside ? 1 : 0); 387 break; 388 default: 389 case MODE_CUBE: 390 //if (i == 0 && y == 0 && z == 0) 391 //dumpPlane(x, null); 392 if (mappingPlane == null) { 393 v = vertexValues[i] = voxelData[x + offset.x][y 394 + offset.y][z + offset.z]; 395 } else { 396 volumeData.voxelPtToXYZ(x + offset.x, y + offset.y, z 397 + offset.z, pt0); 398 v = vertexValues[i] = volumeData.distanceToMappingPlane(pt0); 399 } 400 if (isSquared) 401 vertexValues[i] *= vertexValues[i]; 402 isInside = (allInside ? true : isInside(vertexValues[i], cutoff, isCutoffAbsolute)); 403 if (isInside) 404 bsVoxels.set(pti); 405 } 406 return v; 407 } 408 getPlane(int i, boolean andSwap)409 private void getPlane(int i, boolean andSwap) { 410 if (i < 0 || i > cubeCountX) 411 return; 412 /*float[] p = */surfaceReader.getPlane(i); 413 //dumpPlane(i, p); 414 if (andSwap) { 415 float[] plane = yzPlanes[0]; 416 yzPlanes[0] = yzPlanes[1]; 417 yzPlanes[1] = plane; 418 } 419 } 420 // private void dumpPlane(int n, float[] plane) { 421 // float test = 0; 422 // if (plane == null) 423 // for (int y = 0; y <= cubeCountY; y++) 424 // for (int z = 0; z <= cubeCountZ; z++) 425 // test += volumeData.voxelData[n][y][z]; 426 // else 427 // for (int i = 0; i < yzCount; i++) 428 // test += plane[i]; 429 // System.out.println(isXLowToHigh + " test " + n + "=" + test); 430 // } 431 processTriangles(int insideMask)432 protected void processTriangles(int insideMask) { 433 434 // the inside mask serves to define the triangles necessary 435 // if just creating JVXL files, this step is unnecessary 436 437 byte[] triangles = triangleTable2[insideMask]; 438 for (int i = triangles.length; (i -= 4) >= 0;) 439 addTriangle(triangles[i], triangles[i + 1], triangles[i + 2], 440 triangles[i + 3]); 441 } 442 addVertex(int x, int y, int z, int pti, float value)443 protected void addVertex(int x, int y, int z, int pti, float value) { 444 volumeData.voxelPtToXYZ(x, y, z, pt0); 445 if (surfaceReader.addVertexCopy(pt0, value, -4, true) < 0) 446 bsExcludedVertices.set(pti); 447 } 448 449 protected int nTriangles; addTriangle(int ia, int ib, int ic, int edgeType)450 protected void addTriangle(int ia, int ib, int ic, int edgeType) { 451 if (!bsExcludedTriangles.get(nTriangles) && 452 surfaceReader.addTriangleCheck(edgePointIndexes[ia], 453 edgePointIndexes[ib], edgePointIndexes[ic], 454 edgeType, 0, isCutoffAbsolute, 0) < 0) { 455 bsExcludedTriangles.set(nTriangles); 456 } 457 nTriangles++; 458 } 459 460 protected BS bsValues = new BS(); 461 getValueArray(int x, int y, int z, int pt, float[] tempValues)462 protected float getValueArray(int x, int y, int z, int pt, float[] tempValues) { 463 int ptyz = pt % yzCount; 464 //if (bsValues.get(pt)) 465 //return tempValues[ptyz]; 466 bsValues.set(pt); 467 float value = surfaceReader.getValue(x, y, z, ptyz); 468 if (isSquared) 469 value *= value; 470 tempValues[ptyz] = value; 471 if (isInside(value, cutoff, isCutoffAbsolute)) 472 bsVoxels.set(pt); 473 return value; 474 } 475 isInside(float voxelValue, float max, boolean isAbsolute)476 public static boolean isInside(float voxelValue, float max, boolean isAbsolute) { 477 return ((max > 0 && (isAbsolute ? Math.abs(voxelValue) : voxelValue) >= max) || (max <= 0 && voxelValue <= max)); 478 } 479 480 protected final P3 pt0 = new P3(); 481 protected final P3 pointA = new P3(); 482 483 protected final static int[] edgeVertexPointersLowToHigh = new int[] { 484 1, 1, 2, 0, 485 5, 5, 6, 4, 486 0, 1, 2, 3 487 }; 488 489 protected final static int[] edgeVertexPointersHighToLow = new int[] { 490 0, 1, 3, 0, 491 4, 5, 7, 4, 492 0, 1, 2, 3 493 }; 494 495 protected int[] edgeVertexPointers; 496 497 protected final static int[] edgeVertexPlanesLowToHigh = new int[] { 498 1, 1, 1, 0, 499 1, 1, 1, 0, 500 0, 1, 1, 0 501 }; // from high to low, only edges 3, 7, 8, and 11 are from plane 0 502 503 protected final static int[] edgeVertexPlanesHighToLow = new int[] { 504 1, 0, 1, 1, 505 1, 0, 1, 1, 506 1, 0, 0, 1 507 }; //from high to low, only edges 1, 5, 9, and 10 are from plane 0 508 509 protected int[] edgeVertexPlanes; 510 processOneCubical(int insideMask, int x, int y, int z, int pt)511 protected boolean processOneCubical(int insideMask, int x, int y, int z, int pt) { 512 513 /* 514 * The key to the algorithm is that we have a catalog that 515 * maps the inside-vertex mask to an edge mask, and then 516 * each edge is associated with a specific vertex. 517 * 518 * Each cube vertex may be associated with from 0 to 3 edges, 519 * depending upon where it lies in the overall cube of data. 520 * 521 * When scanning X from low to high, the "leading vertex" is 522 * vertex 1 and edgeVertexPlanes[1]. Edges 0, 1, and 9 are 523 * associated with vertex 1, and others are associated similarly. 524 * 525 * When scanning X from high to low, the "leading vertex" is 526 * vertex 0 and edgeVertexPlanes[1]. Edges 0, 3, and 8 are 527 * associated with vertex 0, and others are associated similarly. 528 * 529 * edgePointIndexes[iEdge] tracks the vertex index for this 530 * specific cubical so that triangles can be created properly. 531 * 532 * 533 * Y 534 * 4 --------4--------- 5 535 * /| /| 536 * / | / | 537 * / | / | 538 * 7 8 5 | 539 * / | / 9 540 * / | / | 541 * 7 --------6--------- 6 | 542 * | | | | 543 * | 0 ---------0--|----- 1 X 544 * | / | / 545 * 11 / 10 / 546 * | 3 | 1 547 * | / | / 548 * | / | / 549 * 3 ---------2-------- 2 550 * Z 551 * / / 552 * edgeVertexPlanes[0] [1] (scanning x low to high) 553 * edgeVertexPlanes[1] [0] (scanning x high to low) 554 * 555 */ 556 557 558 559 int edgeMask = insideMaskTable[insideMask]; 560 boolean isNaN = false; 561 for (int iEdge = 12; --iEdge >= 0;) { 562 563 // bit set to one means it's a relevant edge 564 565 int xEdge = Pwr2[iEdge]; 566 if ((edgeMask & xEdge) == 0) 567 continue; 568 569 // if we have a point already, we don't need to check this edge. 570 // for triangles, this will be an index into an array; 571 // for just creating JVXL files, this can just be 0 572 573 int iPlane = edgeVertexPlanes[iEdge]; 574 int iPt = (pt + linearOffsets[edgeVertexPointers[iEdge]]) % yzCount; 575 int iType = edgeTypeTable[iEdge]; 576 int index = edgePointIndexes[iEdge] = isoPointIndexPlanes[iPlane][iPt][iType]; 577 if (index != Integer.MIN_VALUE) { 578 if (index == -1) 579 isNaN = excludePartialCubes; // -- problem with older Jmol files? 580 // this says, "If any point on the cube is NaN, then 581 //don't process the cube. 582 continue; // propagated from neighbor 583 } 584 // here's an edge that has to be checked. 585 586 // get the vertex numbers 0 - 7 587 588 int vertexA = edgeVertexes[iEdge << 1]; 589 int vertexB = edgeVertexes[(iEdge << 1) + 1]; 590 591 // pick up the actual value at each vertex 592 // this array of 8 values is updated as we go. 593 594 float valueA = vertexValues[vertexA]; 595 float valueB = vertexValues[vertexB]; 596 597 // we allow for NaN values -- missing triangles 598 599 600 // the exact point position -- not important for just 601 // creating the JVXL file. In that case, all you 602 // need are the two values valueA and valueB and the cutoff. 603 // from those you can define the fractional offset 604 605 // here is where we get the value and assign the point for that edge 606 // it is where the JVXL surface data line is appended 607 608 calcVertexPoint(x, y, z, vertexA, pointA); 609 610 edgeCount++; 611 612 int i = edgePointIndexes[iEdge] = isoPointIndexPlanes[iPlane][iPt][iType] = surfaceReader 613 .getSurfacePointIndexAndFraction(cutoff, isCutoffAbsolute, x, y, z, 614 cubeVertexOffsets[vertexA], vertexA, vertexB, valueA, valueB, 615 pointA, edgeVectors[iEdge], iType == contourType, fReturn); 616 617 addEdgeData(i < 0 ? Float.NaN : fReturn[0]); 618 619 // If the fraction returns NaN, this is because one of the end points 620 // is NaN; if the point index returns -1, then the point has been 621 // excluded by the meshDataServer (Isosurface) for some other reason, 622 // for example because it is outside the limits of the box. 623 624 if (Float.isNaN(fReturn[0]) || i < 0) 625 isNaN = excludePartialCubes; 626 } 627 return !isNaN; 628 } 629 630 protected void addEdgeData(float f) { 631 char ch = JvxlCoder.jvxlFractionAsCharacter(f); 632 edgeData.appendC(ch); 633 } 634 635 protected float[] fReturn = new float[1]; 636 637 public void calcVertexPoint(int x, int y, int z, int vertex, P3 pt) { 638 volumeData.voxelPtToXYZ(x, y, z, pt0); 639 pt.add2(pt0, voxelVertexVectors[vertex]); 640 } 641 642 protected final static V3[] cubeVertexVectors = { 643 V3.new3(0, 0, 0), 644 V3.new3(1, 0, 0), 645 V3.new3(1, 0, 1), 646 V3.new3(0, 0, 1), 647 V3.new3(0, 1, 0), 648 V3.new3(1, 1, 0), 649 V3.new3(1, 1, 1), 650 V3.new3(0, 1, 1) }; 651 652 653 /* Y 654 * 4 --------4--------- 5 +z --------4--------- +yz+z 655 * /| /| /| /| 656 * / | / | / | / | 657 * / | / | / | / | 658 * 7 8 5 | 7 8 5 | 659 * / | / 9 / | / 9 660 * / | / | / | / | 661 * 7 --------6--------- 6 | +z+1 --------6--------- +yz+z+1| 662 * | | | | | | | | 663 * | 0 ---------0--|----- 1 X | 0 ---------0--|----- +yz X(outer) 664 * | / | / | / | / 665 * 11 / 10 / 11 / 10 / 666 * | 3 | 1 | 3 | 1 667 * | / | / | / | / 668 * | / | / | / | / 669 * 3 ---------2-------- 2 +1 ---------2-------- +yz+1 670 * Z Z (inner) 671 * 672 * streaming data offsets 673 * type 0: x-edges: 0 2 4 6 674 * type 1: y-edges: 8 9 10 11 675 * type 2: z-edges: 1 3 5 7 676 * 677 * Data stream offsets for vertices, relative to point 0, based on reading 678 * loops {for x {for y {for z}}} 0-->n-1 679 * y and z are numbers of grid points in those directions: 680 * 681 * 0 1 2 3 4 5 6 7 682 * 0 +yz +yz+1 +1 +z +yz+z +yz+z+1 +z+1 683 * 684 * These are just looked up in a table. After the first set of cubes, 685 * we are only adding points 1, 2, 5 or 6. This means that initially 686 * we need two data slices, but after that only one (slice 1): 687 * 688 * base 689 * offset 0 1 2 3 4 5 6 7 690 * slice[0] 0 +1 +z +z+1 691 * slice[1] +yz 0 +1 +z +z+1 692 * 693 * slice: 0 1 1 0 0 1 1 0 694 * 695 * We can request reading of two slices (2*nY*nZ data points) first, then 696 * from then on, just nY*nZ points. "Reading" is really just being handed a 697 * pointer into an array. Perhaps that array is already filled completely; 698 * perhaps it is being read incrementally. 699 * 700 * As it is now, the JVXL data are read into a BitSet 701 * so we can continue to do that with NON progressive files. 702 * 703 * 704 */ 705 706 protected final static int edgeTypeTable[] = { 707 0, 2, 0, 2, 708 0, 2, 0, 2, 709 1, 1, 1, 1 }; 710 // 0=along X, 1=along Y, 2=along Z 711 712 protected final int[] linearOffsets = new int[8]; 713 714 /* 715 * set the linear offsets for generating a unique cell ID, 716 * for pointing into the inside/outside BitSet, 717 * and for finding the associated vertex for an edge. 718 * 719 */ 720 721 protected void setLinearOffsets() { 722 linearOffsets[0] = 0; 723 linearOffsets[1] = yzCount; 724 linearOffsets[2] = yzCount + 1; 725 linearOffsets[3] = 1; 726 linearOffsets[4] = nZ; 727 linearOffsets[5] = yzCount + nZ; 728 linearOffsets[6] = yzCount + nZ + 1; 729 linearOffsets[7] = nZ + 1; 730 } 731 732 public int getLinearOffset(int x, int y, int z, int offset) { 733 return x * yzCount + y * nZ + z + linearOffsets[offset]; 734 } 735 736 protected final static short insideMaskTable[] = { 0x0000, 0x0109, 0x0203, 737 0x030A, 0x0406, 0x050F, 0x0605, 0x070C, 0x080C, 0x0905, 0x0A0F, 0x0B06, 738 0x0C0A, 0x0D03, 0x0E09, 0x0F00, 0x0190, 0x0099, 0x0393, 0x029A, 0x0596, 739 0x049F, 0x0795, 0x069C, 0x099C, 0x0895, 0x0B9F, 0x0A96, 0x0D9A, 0x0C93, 740 0x0F99, 0x0E90, 0x0230, 0x0339, 0x0033, 0x013A, 0x0636, 0x073F, 0x0435, 741 0x053C, 0x0A3C, 0x0B35, 0x083F, 0x0936, 0x0E3A, 0x0F33, 0x0C39, 0x0D30, 742 0x03A0, 0x02A9, 0x01A3, 0x00AA, 0x07A6, 0x06AF, 0x05A5, 0x04AC, 0x0BAC, 743 0x0AA5, 0x09AF, 0x08A6, 0x0FAA, 0x0EA3, 0x0DA9, 0x0CA0, 0x0460, 0x0569, 744 0x0663, 0x076A, 0x0066, 0x016F, 0x0265, 0x036C, 0x0C6C, 0x0D65, 0x0E6F, 745 0x0F66, 0x086A, 0x0963, 0x0A69, 0x0B60, 0x05F0, 0x04F9, 0x07F3, 0x06FA, 746 0x01F6, 0x00FF, 0x03F5, 0x02FC, 0x0DFC, 0x0CF5, 0x0FFF, 0x0EF6, 0x09FA, 747 0x08F3, 0x0BF9, 0x0AF0, 0x0650, 0x0759, 0x0453, 0x055A, 0x0256, 0x035F, 748 0x0055, 0x015C, 0x0E5C, 0x0F55, 0x0C5F, 0x0D56, 0x0A5A, 0x0B53, 0x0859, 749 0x0950, 0x07C0, 0x06C9, 0x05C3, 0x04CA, 0x03C6, 0x02CF, 0x01C5, 0x00CC, 750 0x0FCC, 0x0EC5, 0x0DCF, 0x0CC6, 0x0BCA, 0x0AC3, 0x09C9, 0x08C0, 0x08C0, 751 0x09C9, 0x0AC3, 0x0BCA, 0x0CC6, 0x0DCF, 0x0EC5, 0x0FCC, 0x00CC, 0x01C5, 752 0x02CF, 0x03C6, 0x04CA, 0x05C3, 0x06C9, 0x07C0, 0x0950, 0x0859, 0x0B53, 753 0x0A5A, 0x0D56, 0x0C5F, 0x0F55, 0x0E5C, 0x015C, 0x0055, 0x035F, 0x0256, 754 0x055A, 0x0453, 0x0759, 0x0650, 0x0AF0, 0x0BF9, 0x08F3, 0x09FA, 0x0EF6, 755 0x0FFF, 0x0CF5, 0x0DFC, 0x02FC, 0x03F5, 0x00FF, 0x01F6, 0x06FA, 0x07F3, 756 0x04F9, 0x05F0, 0x0B60, 0x0A69, 0x0963, 0x086A, 0x0F66, 0x0E6F, 0x0D65, 757 0x0C6C, 0x036C, 0x0265, 0x016F, 0x0066, 0x076A, 0x0663, 0x0569, 0x0460, 758 0x0CA0, 0x0DA9, 0x0EA3, 0x0FAA, 0x08A6, 0x09AF, 0x0AA5, 0x0BAC, 0x04AC, 759 0x05A5, 0x06AF, 0x07A6, 0x00AA, 0x01A3, 0x02A9, 0x03A0, 0x0D30, 0x0C39, 760 0x0F33, 0x0E3A, 0x0936, 0x083F, 0x0B35, 0x0A3C, 0x053C, 0x0435, 0x073F, 761 0x0636, 0x013A, 0x0033, 0x0339, 0x0230, 0x0E90, 0x0F99, 0x0C93, 0x0D9A, 762 0x0A96, 0x0B9F, 0x0895, 0x099C, 0x069C, 0x0795, 0x049F, 0x0596, 0x029A, 763 0x0393, 0x0099, 0x0190, 0x0F00, 0x0E09, 0x0D03, 0x0C0A, 0x0B06, 0x0A0F, 764 0x0905, 0x080C, 0x070C, 0x0605, 0x050F, 0x0406, 0x030A, 0x0203, 0x0109, 765 0x0000 }; 766 767 768 } 769