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.readers; 25 26 27 import javajs.util.BS; 28 import org.jmol.jvxl.api.MeshDataServer; 29 import org.jmol.jvxl.api.VertexDataServer; 30 import org.jmol.jvxl.calc.MarchingCubes; 31 import org.jmol.jvxl.calc.MarchingSquares; 32 import org.jmol.jvxl.data.JvxlCoder; 33 import org.jmol.jvxl.data.JvxlData; 34 import org.jmol.jvxl.data.MeshData; 35 import org.jmol.jvxl.data.VolumeData; 36 import org.jmol.quantum.QuantumPlaneCalculation; 37 import org.jmol.util.BoxInfo; 38 import org.jmol.util.C; 39 import org.jmol.util.ColorEncoder; 40 import org.jmol.util.Escape; 41 42 import javajs.util.AU; 43 import javajs.util.Lst; 44 import javajs.util.SB; 45 46 import org.jmol.util.Logger; 47 48 import javajs.util.OC; 49 import javajs.util.M3; 50 import javajs.util.P3; 51 import javajs.util.P3i; 52 import javajs.util.T3; 53 import javajs.util.V3; 54 55 56 public abstract class SurfaceReader implements VertexDataServer { 57 58 /* 59 * JVXL SurfaceReader Class 60 * ---------------------- 61 * Bob Hanson, hansonr@stolaf.edu, 20 Apr 2007 62 * 63 * SurfaceReader performs four functions: 64 * 65 * 1) reading/creating volume scalar data ("voxels") 66 * 2) generating a surface (vertices and triangles) from this set 67 * based on a specific cutoff value 68 * 3) color-mapping this surface with other data 69 * 4) creating JVXL format file data for this surface 70 * 71 * In the case that the surface type does not include voxel data (EfvetReader), 72 * only steps 2 and 3 are involved, and no cutoff value is used. 73 * 74 * SurfaceReader is an ABSTRACT class, instantiated as one of the 75 * following to perform specific functions: 76 * 77 * SurfaceReader (abstract MarchingReader) 78 * | 79 * |_______VolumeDataReader (uses provided predefined data) 80 * | | 81 * | |_____IsoFxyReader (creates data as needed) 82 * | |_____IsoFxyzReader (creates data as needed) 83 * | |_____IsoMepReader (creates predefined data) 84 * | |_____IsoMOReader (creates predefined data) 85 * | |_____IsoShapeReader (creates data as needed) 86 * | |_____IsoSolventReader (creates predefined data) 87 * | |___IsoPlaneReader (predefines data) 88 * | 89 * |_______SurfaceFileReader (abstract) 90 * | 91 * |_______VolumeFileReader (abstract) 92 * | | 93 * | |______ApbsReader 94 * | |______CubeReader 95 * | |______JaguarReader 96 * | |______JvxlXmlReader 97 * | | |______JvxlReader 98 * | | 99 * | |______ElectronDensityFileReader (abstract) 100 * | |______MrcBinaryReader 101 * | |______XplorReader (version 3.1) 102 * | 103 * |_______PolygonFileReader (abstract) 104 * | 105 * |______EfvetReader 106 * |______PmeshReader 107 * 108 * The first step is to create a VolumeData structure: 109 * 110 * public final Point3f volumetricOrigin = new Point3f(); 111 * public final Vector3f[] volumetricVectors = new Vector3f[3]; 112 * public final int[] voxelCounts = new int[3]; 113 * public float[][][] voxelData; 114 * 115 * such as exists in a CUBE file. 116 * 117 * The second step is to use the Marching Cubes algorithm to 118 * create a surface set of vertices and triangles. The data structure 119 * involved for that is MeshData, containing: 120 * 121 * public int vertexCount; 122 * public Point3f[] vertices; 123 * public float[] vertexValues; 124 * public int polygonCount; 125 * public int[][] polygonIndexes; 126 * 127 * The third (optional) step is to color those vertices using 128 * a set of color index values provided by a color encoder. This 129 * data is also stored in MeshData: 130 * 131 * public short[] vertexColixes; 132 * 133 * Finally -- actually, throughout the process -- SurfaceReader 134 * creates a JvxlData structure containing the critical information 135 * that is necessary for creating Jvxl surface data files. For that, 136 * we have the JvxlData structure. 137 * 138 * Two interfaces are defined, and more should be. These include 139 * VertexDataServer and MeshDataServer. 140 * 141 * VertexDataServer 142 * ---------------- 143 * 144 * contains three methods, getSurfacePointIndex, addVertexCopy, and addTriangleCheck. 145 * 146 * These deliver MarchingCubes and MarchingSquares vertex data in 147 * return for a vertex index number that can later be used for defining 148 * a set of triangles. 149 * 150 * SurfaceReader implements this interface. 151 * 152 * 153 * MeshDataServer extends VertexDataServer 154 * --------------------------------------- 155 * 156 * contains additional methods that allow for later processing 157 * of the vertex/triangle data: 158 * 159 * public abstract void invalidateTriangles(); 160 * public abstract void fillMeshData(MeshData meshData, int mode); 161 * public abstract void notifySurfaceGenerationCompleted(); 162 * public abstract void notifySurfaceMappingCompleted(); 163 * 164 * Note that, in addition to these interfaces, some of the readers, 165 * namely IsoFxyReader, IsoMepReader,IsoMOReader, and IsoSolvenReader 166 * and (due to subclassing) IsoPlaneReader all currently require 167 * direct connections to Jmol Viewer and Atom classes. 168 * 169 * 170 * The rough outline of Jvxl files is 171 * given below: 172 * 173 174 #comments (optional) 175 info line1 176 info line2 177 -na originx originy originz [ANGSTROMS/BOHR] optional; BOHR assumed 178 n1 x y z 179 n2 x y z 180 n3 x y z 181 a1 a1.0 x y z 182 a2 a2.0 x y z 183 a3 a3.0 x y z 184 a4 a4.0 x y z 185 etc. -- na atoms 186 -ns 35 90 35 90 Jmol voxel format version 1.0 187 # more comments 188 cutoff +/-nEdges +/-nVertices [more here] 189 integer inside/outside edge data 190 ascii-encoded fractional edge data 191 ascii-encoded fractional color data 192 # optional comments 193 194 * 195 * 196 * 197 * 198 */ 199 200 protected SurfaceGenerator sg; 201 protected MeshDataServer meshDataServer; 202 203 protected Parameters params; 204 protected MeshData meshData; 205 protected JvxlData jvxlData; 206 VolumeData volumeData; 207 private String edgeData; 208 209 protected boolean haveSurfaceAtoms = false; 210 protected boolean allowSigma = false; 211 protected boolean isProgressive = false; 212 protected boolean isXLowToHigh = false; //can be overridden in some readers by --progressive 213 private float assocCutoff = 0.3f; 214 protected boolean isQuiet; 215 protected boolean isPeriodic; 216 217 boolean vertexDataOnly; 218 boolean hasColorData; 219 220 protected float dataMin = Float.MAX_VALUE; 221 protected float dataMax = -Float.MAX_VALUE; 222 protected float dataMean; 223 protected P3 xyzMin, xyzMax; 224 225 protected P3 center; 226 protected float[] anisotropy; 227 protected boolean isAnisotropic; 228 protected M3 eccentricityMatrix; 229 protected M3 eccentricityMatrixInverse; 230 protected boolean isEccentric; 231 protected float eccentricityScale; 232 protected float eccentricityRatio; 233 SurfaceReader()234 SurfaceReader() {} 235 236 /** 237 * implemented in SurfaceFileReader and 238 * 239 * @param sg 240 */ init(SurfaceGenerator sg)241 abstract void init(SurfaceGenerator sg); 242 initSR(SurfaceGenerator sg)243 void initSR(SurfaceGenerator sg) { 244 this.sg = sg; 245 params = sg.params; 246 assocCutoff = params.assocCutoff; 247 isXLowToHigh = params.isXLowToHigh; 248 center = params.center; 249 anisotropy = params.anisotropy; 250 isAnisotropic = params.isAnisotropic; 251 eccentricityMatrix = params.eccentricityMatrix; 252 eccentricityMatrixInverse = params.eccentricityMatrixInverse; 253 isEccentric = params.isEccentric; 254 eccentricityScale = params.eccentricityScale; 255 eccentricityRatio = params.eccentricityRatio; 256 marchingSquares = sg.marchingSquares; 257 meshData = sg.meshData; 258 jvxlData = sg.jvxlData; 259 setVolumeDataV(sg.volumeDataTemp); // initialize volume data to surfaceGenerator's 260 meshDataServer = sg.meshDataServer; 261 cJvxlEdgeNaN = (char) (JvxlCoder.defaultEdgeFractionBase + JvxlCoder.defaultEdgeFractionRange); 262 } 263 264 final static float ANGSTROMS_PER_BOHR = 0.5291772f; 265 final static float defaultMappedDataMin = 0f; 266 final static float defaultMappedDataMax = 1.0f; 267 final static float defaultCutoff = 0.02f; 268 269 private int edgeCount; 270 271 protected P3 volumetricOrigin; 272 protected V3[] volumetricVectors; 273 protected int[] voxelCounts; 274 protected float[][][] voxelData; 275 closeReader()276 abstract protected void closeReader(); 277 278 /** 279 * 280 * @param out 281 */ setOutputChannel(OC out)282 protected void setOutputChannel(OC out) { 283 // only for file readers 284 } 285 newVoxelDataCube()286 protected void newVoxelDataCube() { 287 volumeData.setVoxelDataAsArray(voxelData = new float[nPointsX][nPointsY][nPointsZ]); 288 } 289 setVolumeDataV(VolumeData v)290 protected void setVolumeDataV(VolumeData v) { 291 nBytes = 0; 292 volumetricOrigin = v.volumetricOrigin; 293 volumetricVectors = v.volumetricVectors; 294 voxelCounts = v.voxelCounts; 295 voxelData = v.getVoxelData(); 296 volumeData = v; 297 } 298 readVolumeParameters(boolean isMapData)299 protected abstract boolean readVolumeParameters(boolean isMapData); 300 readVolumeData(boolean isMapData)301 protected abstract boolean readVolumeData(boolean isMapData); 302 303 //////////////////////////////////////////////////////////////// 304 // CUBE/APBS/JVXL file reading stuff 305 //////////////////////////////////////////////////////////////// 306 307 protected long nBytes; 308 protected int nDataPoints; 309 protected int nPointsX, nPointsY, nPointsZ; 310 311 protected boolean isJvxl; 312 313 protected int edgeFractionBase; 314 protected int edgeFractionRange; 315 protected int colorFractionBase; 316 protected int colorFractionRange; 317 318 protected SB jvxlFileHeaderBuffer; 319 protected SB fractionData; 320 protected String jvxlEdgeDataRead = ""; 321 protected String jvxlColorDataRead = ""; 322 protected BS jvxlVoxelBitSet; 323 protected boolean jvxlDataIsColorMapped; 324 protected boolean jvxlDataIsPrecisionColor; 325 protected boolean jvxlDataIs2dContour; 326 protected boolean jvxlDataIsColorDensity; 327 protected float jvxlCutoff; 328 protected int jvxlNSurfaceInts; 329 protected char cJvxlEdgeNaN = '\0'; 330 331 protected int contourVertexCount; 332 jvxlUpdateInfo()333 void jvxlUpdateInfo() { 334 jvxlData.jvxlUpdateInfo(params.title, nBytes); 335 } 336 readAndSetVolumeParameters(boolean isMapData)337 boolean readAndSetVolumeParameters(boolean isMapData) { 338 if (!readVolumeParameters(isMapData)) 339 return false; 340 if (vertexDataOnly) 341 return true; 342 //if (volumeData.sr != null) 343 //return true; 344 return (volumeData.setUnitVectors());// || isMapData && params.thePlane != null); 345 346 // && (vertexDataOnly 347 // || isMapData && params.thePlane != null 348 // || volumeData.sr != null 349 // || volumeData.setUnitVectors())); 350 } 351 createIsosurface(boolean justForPlane)352 boolean createIsosurface(boolean justForPlane) { 353 resetIsosurface(); 354 if (params.showTiming) 355 Logger.startTimer("isosurface creation"); 356 jvxlData.cutoff = Float.NaN; 357 if (!readAndSetVolumeParameters(justForPlane)) 358 return false; 359 if (!justForPlane && !Float.isNaN(params.sigma) && !allowSigma) { 360 if (params.sigma > 0) 361 Logger 362 .error("Reader does not support SIGMA option -- using cutoff 1.6"); 363 params.cutoff = 1.6f; 364 } 365 // negative sigma just ignores the error message 366 // and means it was inserted by Jmol as a default option 367 if (params.sigma < 0) 368 params.sigma = -params.sigma; 369 nPointsX = voxelCounts[0]; 370 nPointsY = voxelCounts[1]; 371 nPointsZ = voxelCounts[2]; 372 jvxlData.isSlabbable = ((params.dataType & Parameters.IS_SLABBABLE) != 0); 373 jvxlData.insideOut = params.isInsideOut(); 374 jvxlData.isBicolorMap = params.isBicolorMap; 375 jvxlData.nPointsX = nPointsX; 376 jvxlData.nPointsY = nPointsY; 377 jvxlData.nPointsZ = nPointsZ; 378 jvxlData.jvxlVolumeDataXml = volumeData.xmlData; 379 jvxlData.voxelVolume = volumeData.voxelVolume; 380 if (justForPlane) { 381 //float[][][] voxelDataTemp = volumeData.voxelData; 382 volumeData.setMappingPlane(params.thePlane); 383 //volumeData.setDataDistanceToPlane(params.thePlane); 384 if (meshDataServer != null) 385 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); 386 params.setMapRanges(this, false); 387 generateSurfaceData(); 388 volumeData.setMappingPlane(null); 389 390 //if (volumeData != null) 391 // volumeData.voxelData = voxelDataTemp; 392 } else { 393 if (!readVolumeData(false)) 394 return false; 395 generateSurfaceData(); 396 } 397 398 if (jvxlFileHeaderBuffer == null) { 399 jvxlData.jvxlFileTitle = ""; 400 jvxlData.jvxlFileSource = null; 401 jvxlData.jvxlFileMessage = null; 402 } else { 403 String s = jvxlFileHeaderBuffer.toString(); 404 int i = s.indexOf('\n', s.indexOf('\n', s.indexOf('\n') + 1) + 1) + 1; 405 jvxlData.jvxlFileTitle = s.substring(0, i); 406 jvxlData.jvxlFileSource = params.fileName; 407 } 408 if (params.contactPair == null) 409 setBBoxAll(); 410 jvxlData.isValid = (xyzMin.x != Float.MAX_VALUE); 411 if (!params.isSilent) { 412 if (!jvxlData.isValid) 413 Logger.error("no isosurface points were found!"); 414 else 415 Logger.info("boundbox corners " + Escape.eP(xyzMin) + " " 416 + Escape.eP(xyzMax)); 417 } 418 jvxlData.boundingBox = new P3[] { xyzMin, xyzMax }; 419 jvxlData.dataMin = dataMin; 420 jvxlData.dataMax = dataMax; 421 jvxlData.cutoff = (isJvxl ? jvxlCutoff : params.cutoff); 422 jvxlData.isCutoffAbsolute = params.isCutoffAbsolute; 423 jvxlData.isModelConnected = params.isModelConnected; 424 jvxlData.pointsPerAngstrom = 1f / volumeData.volumetricVectorLengths[0]; 425 jvxlData.jvxlColorData = ""; 426 jvxlData.jvxlPlane = params.thePlane; 427 jvxlData.jvxlEdgeData = edgeData; 428 jvxlData.isBicolorMap = params.isBicolorMap; 429 jvxlData.isContoured = params.isContoured; 430 jvxlData.colorDensity = params.colorDensity; 431 jvxlData.pointSize = params.pointSize; 432 if (jvxlData.vContours != null) 433 params.nContours = jvxlData.vContours.length; 434 jvxlData.nContours = (params.contourFromZero ? params.nContours : -1 435 - params.nContours); 436 jvxlData.thisContour = params.thisContour; 437 jvxlData.nEdges = edgeCount; 438 jvxlData.edgeFractionBase = edgeFractionBase; 439 jvxlData.edgeFractionRange = edgeFractionRange; 440 jvxlData.colorFractionBase = colorFractionBase; 441 jvxlData.colorFractionRange = colorFractionRange; 442 jvxlData.jvxlDataIs2dContour = jvxlDataIs2dContour; 443 jvxlData.jvxlDataIsColorMapped = jvxlDataIsColorMapped; 444 jvxlData.jvxlDataIsColorDensity = jvxlDataIsColorDensity; 445 jvxlData.isXLowToHigh = isXLowToHigh; 446 jvxlData.vertexDataOnly = vertexDataOnly; 447 jvxlData.saveVertexCount = 0; 448 449 if (jvxlDataIsColorMapped || jvxlData.nVertexColors > 0) { 450 if (meshDataServer != null) { 451 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); 452 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_COLOR_INDEXES, 453 null); 454 } 455 jvxlData.jvxlColorData = readColorData(); 456 updateSurfaceData(); 457 if (meshDataServer != null) 458 meshDataServer.notifySurfaceMappingCompleted(); 459 } 460 if (params.showTiming) 461 Logger.checkTimer("isosurface creation", false); 462 return true; 463 } 464 resetIsosurface()465 void resetIsosurface() { 466 meshData = new MeshData(); 467 xyzMin = xyzMax = null; 468 jvxlData.isBicolorMap = params.isBicolorMap; 469 if (meshDataServer != null) 470 meshDataServer.fillMeshData(null, 0, null); 471 contourVertexCount = 0; 472 if (params.cutoff == Float.MAX_VALUE) 473 params.cutoff = defaultCutoff; 474 jvxlData.jvxlSurfaceData = ""; 475 jvxlData.jvxlEdgeData = ""; 476 jvxlData.jvxlColorData = ""; 477 //TODO: more resets of jvxlData? 478 edgeCount = 0; 479 edgeFractionBase = JvxlCoder.defaultEdgeFractionBase; 480 edgeFractionRange = JvxlCoder.defaultEdgeFractionRange; 481 colorFractionBase = JvxlCoder.defaultColorFractionBase; 482 colorFractionRange = JvxlCoder.defaultColorFractionRange; 483 params.mappedDataMin = Float.MAX_VALUE; 484 } 485 discardTempData(boolean discardAll)486 void discardTempData(boolean discardAll) { 487 discardTempDataSR(discardAll); 488 } 489 discardTempDataSR(boolean discardAll)490 protected void discardTempDataSR(boolean discardAll) { 491 if (!discardAll) 492 return; 493 voxelData = null; 494 sg.marchingSquares = marchingSquares = null; 495 marchingCubes = null; 496 } 497 initializeVolumetricData()498 void initializeVolumetricData() { 499 nPointsX = voxelCounts[0]; 500 nPointsY = voxelCounts[1]; 501 nPointsZ = voxelCounts[2]; 502 setVolumeDataV(volumeData); 503 } 504 505 // this needs to be specific for each reader readSurfaceData(boolean isMapData)506 abstract protected void readSurfaceData(boolean isMapData) throws Exception; 507 gotoAndReadVoxelData(boolean isMapData)508 protected boolean gotoAndReadVoxelData(boolean isMapData) { 509 //overloaded in jvxlReader 510 initializeVolumetricData(); 511 if (nPointsX > 0 && nPointsY > 0 && nPointsZ > 0) 512 try { 513 gotoData(params.fileIndex - 1, nPointsX * nPointsY * nPointsZ); 514 readSurfaceData(isMapData); 515 } catch (Exception e) { 516 Logger.error(e.toString()); 517 return false; 518 } 519 return true; 520 } 521 522 /** 523 * 524 * @param n 525 * @param nPoints 526 * @throws Exception 527 */ gotoData(int n, int nPoints)528 protected void gotoData(int n, int nPoints) throws Exception { 529 //only for file reader 530 } 531 readColorData()532 protected String readColorData() { 533 if (jvxlData.vertexColors == null) 534 return ""; 535 int vertexCount = jvxlData.vertexCount; 536 short[] colixes = meshData.vcs; 537 float[] vertexValues = meshData.vvs; 538 if (colixes == null || colixes.length < vertexCount) 539 meshData.vcs = colixes = new short[vertexCount]; 540 if (vertexValues == null || vertexValues.length < vertexCount) 541 meshData.vvs = vertexValues = new float[vertexCount]; 542 for (int i = 0; i < vertexCount; i++) 543 colixes[i] = C.getColix(jvxlData.vertexColors[i]); 544 return "-"; 545 } 546 547 //////////////////////////////////////////////////////////////// 548 // marching cube stuff 549 //////////////////////////////////////////////////////////////// 550 551 protected MarchingSquares marchingSquares; 552 protected MarchingCubes marchingCubes; 553 554 protected float[][] yzPlanes; 555 protected int yzCount; 556 557 protected QuantumPlaneCalculation qpc; 558 559 @Override getPlane(int x)560 public float[] getPlane(int x) { 561 return getPlaneSR(x); 562 } 563 getPlaneSR(int x)564 protected float[] getPlaneSR(int x) { 565 if (yzCount == 0) 566 initPlanes(); 567 if (qpc != null) // NCICalculation only 568 qpc.getPlane(x, yzPlanes[x % 2]); 569 return yzPlanes[x % 2]; 570 } 571 initPlanes()572 void initPlanes() { 573 yzCount = nPointsY * nPointsZ; 574 if (!isQuiet) 575 Logger.info("reading data progressively -- yzCount = " + yzCount); 576 yzPlanes = AU.newFloat2(2); 577 yzPlanes[0] = new float[yzCount]; 578 yzPlanes[1] = new float[yzCount]; 579 } 580 581 @Override getValue(int x, int y, int z, int ptyz)582 public float getValue(int x, int y, int z, int ptyz) { 583 return getValue2(x, y, z, ptyz); 584 } 585 getValue2(int x, int y, int z, int ptyz)586 protected float getValue2(int x, int y, int z, int ptyz) { 587 return (yzPlanes == null ? voxelData[x][y][z] : yzPlanes[x % 2][ptyz]); 588 } 589 generateSurfaceData()590 private void generateSurfaceData() { 591 edgeData = ""; 592 if (vertexDataOnly) { 593 try { 594 readSurfaceData(false); 595 } catch (Exception e) { 596 System.out.println(e.toString()); 597 Logger.error("Exception in SurfaceReader::readSurfaceData: " 598 + e.toString()); 599 } 600 return; 601 } 602 contourVertexCount = 0; 603 int contourType = -1; 604 marchingSquares = null; 605 606 if (params.thePlane != null || params.isContoured) { 607 marchingSquares = new MarchingSquares(this, volumeData, params.thePlane, 608 params.contoursDiscrete, params.nContours, params.thisContour, 609 params.contourFromZero); 610 contourType = marchingSquares.contourType; 611 marchingSquares.setMinMax(params.valueMappedToRed, 612 params.valueMappedToBlue); 613 } 614 params.contourType = contourType; 615 params.isXLowToHigh = isXLowToHigh; 616 marchingCubes = new MarchingCubes(this, volumeData, params, jvxlVoxelBitSet); 617 String data = marchingCubes.getEdgeData(); 618 if (params.thePlane == null) 619 edgeData = data; 620 jvxlData.setSurfaceInfoFromBitSetPts(marchingCubes.bsVoxels, 621 params.thePlane, params.mapLattice); 622 jvxlData.jvxlExcluded = params.bsExcluded; 623 if (isJvxl) 624 edgeData = jvxlEdgeDataRead; 625 postProcessVertices(); 626 } 627 postProcessVertices()628 protected void postProcessVertices() { 629 // optional 630 } 631 632 ///////////////// MarchingReader Interface Methods /////////////////// 633 634 protected final P3 ptTemp = new P3(); 635 636 @Override getSurfacePointIndexAndFraction(float cutoff, boolean isCutoffAbsolute, int x, int y, int z, P3i offset, int vA, int vB, float valueA, float valueB, T3 pointA, V3 edgeVector, boolean isContourType, float[] fReturn)637 public int getSurfacePointIndexAndFraction(float cutoff, boolean isCutoffAbsolute, 638 int x, int y, int z, P3i offset, int vA, 639 int vB, float valueA, float valueB, 640 T3 pointA, V3 edgeVector, 641 boolean isContourType, float[] fReturn) { 642 float thisValue = getSurfacePointAndFraction(cutoff, isCutoffAbsolute, valueA, 643 valueB, pointA, edgeVector, x, y, z, vA, vB, fReturn, ptTemp); 644 /* 645 * from MarchingCubes 646 * 647 * In the case of a setup for a Marching Squares calculation, 648 * we are collecting just the desired type of intersection for the 2D marching 649 * square contouring -- x, y, or z. In the case of a contoured f(x,y) surface, 650 * we take every point. 651 * 652 */ 653 654 if (marchingSquares != null && params.isContoured) 655 return marchingSquares.addContourVertex(ptTemp, cutoff); 656 int assocVertex = (assocCutoff > 0 ? (fReturn[0] < assocCutoff ? vA 657 : fReturn[0] > 1 - assocCutoff ? vB : MarchingSquares.CONTOUR_POINT) 658 : MarchingSquares.CONTOUR_POINT); 659 if (assocVertex >= 0) 660 assocVertex = marchingCubes.getLinearOffset(x, y, z, assocVertex); 661 int n = addVertexCopy(ptTemp, thisValue, assocVertex, true); 662 if (n >= 0 && params.iAddGridPoints) { 663 marchingCubes.calcVertexPoint(x, y, z, vB, ptTemp); 664 addVertexCopy(valueA < valueB ? pointA : ptTemp, Math.min(valueA, valueB), 665 MarchingSquares.EDGE_POINT, true); 666 addVertexCopy(valueA < valueB ? ptTemp : pointA, Math.max(valueA, valueB), 667 MarchingSquares.EDGE_POINT, true); 668 } 669 return n; 670 } 671 672 protected float getSurfacePointAndFraction(float cutoff, boolean isCutoffAbsolute, 673 float valueA, float valueB, T3 pointA, 674 V3 edgeVector, int x, 675 int y, int z, int vA, int vB, float[] fReturn, T3 ptReturn) { 676 // will be subclassed in many cases. 677 // JavaScript optimization: DO NOT CALL THIS method from subclassed method of the same name! 678 return getSPF(cutoff, isCutoffAbsolute, valueA, valueB, pointA, edgeVector, x, y, z, vA, vB, fReturn, ptReturn); 679 } 680 681 /** 682 * 683 * @param cutoff 684 * @param isCutoffAbsolute 685 * @param valueA 686 * @param valueB 687 * @param pointA 688 * @param edgeVector 689 * @param x TODO 690 * @param y TODO 691 * @param z TODO 692 * @param vA 693 * @param vB 694 * @param fReturn 695 * @param ptReturn 696 * @return fractional distance from A to B 697 */ 698 protected float getSPF(float cutoff, boolean isCutoffAbsolute, float valueA, 699 float valueB, T3 pointA, V3 edgeVector, int x, int y, 700 int z, int vA, int vB, float[] fReturn, T3 ptReturn) { 701 702 //JvxlReader may or may not call this 703 //IsoSolventReader overrides this for nonlinear Marching Cubes (12.1.29) 704 705 float diff = valueB - valueA; 706 float fraction = (cutoff - valueA) / diff; 707 if (isCutoffAbsolute && (fraction < 0 || fraction > 1)) 708 fraction = (-cutoff - valueA) / diff; 709 710 if (fraction < 0 || fraction > 1) { 711 //Logger.error("problem with unusual fraction=" + fraction + " cutoff=" 712 // + cutoff + " A:" + valueA + " B:" + valueB); 713 fraction = Float.NaN; 714 } 715 fReturn[0] = fraction; 716 ptReturn.scaleAdd2(fraction, edgeVector, pointA); 717 return valueA + fraction * diff; 718 } 719 720 @Override 721 public int addVertexCopy(T3 vertexXYZ, float value, int assocVertex, boolean asCopy) { 722 return addVC(vertexXYZ, value, assocVertex, asCopy); 723 } 724 725 protected int addVC(T3 vertexXYZ, float value, int assocVertex, boolean asCopy) { 726 return (Float.isNaN(value) && assocVertex != MarchingSquares.EDGE_POINT ? -1 727 : meshDataServer == null ? meshData.addVertexCopy(vertexXYZ, value, assocVertex, asCopy) : meshDataServer.addVertexCopy(vertexXYZ, value, assocVertex, asCopy)); 728 } 729 730 @Override 731 public int addTriangleCheck(int iA, int iB, int iC, int check, int iContour, 732 boolean isAbsolute, int color) { 733 if (marchingSquares != null && params.isContoured) { 734 if (color == 0) // from marching cubes 735 return marchingSquares.addTriangle(iA, iB, iC, check, iContour); 736 color = 0; // from marchingSquares 737 } 738 return (meshDataServer != null ? meshDataServer.addTriangleCheck(iA, iB, 739 iC, check, iContour, isAbsolute, color) : isAbsolute 740 && !MeshData.checkCutoff(iA, iB, iC, meshData.vvs) ? -1 741 : meshData.addTriangleCheck(iA, iB, iC, check, iContour, color)); 742 } 743 744 //////////////////////////////////////////////////////////////// 745 // color mapping methods 746 //////////////////////////////////////////////////////////////// 747 748 void colorIsosurface() { 749 if (params.isSquared && volumeData != null) 750 volumeData.filterData(true, Float.NaN); 751 /* if (params.isContoured && marchingSquares == null) { 752 // if (params.isContoured && !(jvxlDataIs2dContour || params.thePlane != null)) { 753 Logger.error("Isosurface error: Cannot contour this type of data."); 754 return; 755 } 756 */ 757 if (meshDataServer != null) { 758 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); 759 } 760 761 jvxlData.saveVertexCount = 0; 762 if (params.isContoured && marchingSquares != null) { 763 initializeMapping(); 764 params.setMapRanges(this, false); 765 marchingSquares.setMinMax(params.valueMappedToRed, 766 params.valueMappedToBlue); 767 jvxlData.saveVertexCount = marchingSquares.contourVertexCount; 768 contourVertexCount = marchingSquares 769 .generateContourData(jvxlDataIs2dContour, (params.isSquared ? 1e-8f : 1e-4f)); 770 jvxlData.contourValuesUsed = marchingSquares.contourValuesUsed; 771 minMax = marchingSquares.getMinMax(); 772 if (meshDataServer != null) 773 meshDataServer.notifySurfaceGenerationCompleted(); 774 finalizeMapping(); 775 } 776 777 applyColorScale(); 778 jvxlData.nContours = (params.contourFromZero 779 ? params.nContours : -1 - params.nContours); 780 jvxlData.thisContour = params.thisContour; 781 jvxlData.jvxlFileMessage = "mapped: min = " + params.valueMappedToRed 782 + "; max = " + params.valueMappedToBlue; 783 } 784 785 void applyColorScale() { 786 colorFractionBase = jvxlData.colorFractionBase = JvxlCoder.defaultColorFractionBase; 787 colorFractionRange = jvxlData.colorFractionRange = JvxlCoder.defaultColorFractionRange; 788 if (params.colorPhase == 0) 789 params.colorPhase = 1; 790 if (meshDataServer == null) { 791 meshData.vcs = new short[meshData.vc]; 792 } else { 793 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); 794 if (params.contactPair == null) 795 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_COLOR_INDEXES, 796 null); 797 } 798 //colorBySign is true when colorByPhase is true, but not vice-versa 799 //old: boolean saveColorData = !(params.colorByPhase && !params.isBicolorMap && !params.colorBySign); //sorry! 800 boolean saveColorData = (params.colorDensity || params.isBicolorMap 801 || params.colorBySign || !params.colorByPhase); 802 if (params.contactPair != null) 803 saveColorData = false; 804 // colors mappable always now 805 jvxlData.isJvxlPrecisionColor = true; 806 jvxlData.vertexCount = (contourVertexCount > 0 ? contourVertexCount 807 : meshData.vc); 808 jvxlData.minColorIndex = -1; 809 jvxlData.maxColorIndex = 0; 810 jvxlData.contourValues = params.contoursDiscrete; 811 jvxlData.isColorReversed = params.isColorReversed; 812 if (!params.colorDensity) 813 if (params.isBicolorMap && !params.isContoured || params.colorBySign) { 814 jvxlData.minColorIndex = C 815 .getColixTranslucent3(C.getColix(params.isColorReversed ? params.colorPos 816 : params.colorNeg), jvxlData.translucency != 0, jvxlData.translucency); 817 jvxlData.maxColorIndex = C 818 .getColixTranslucent3(C.getColix(params.isColorReversed ? params.colorNeg 819 : params.colorPos), jvxlData.translucency != 0, jvxlData.translucency); 820 } 821 jvxlData.isTruncated = (jvxlData.minColorIndex >= 0 && !params.isContoured); 822 boolean useMeshDataValues = jvxlDataIs2dContour 823 || 824 // !jvxlDataIs2dContour && (params.isContoured && jvxlData.jvxlPlane != null || 825 hasColorData || vertexDataOnly || params.colorDensity || params.isBicolorMap 826 && !params.isContoured; 827 if (!useMeshDataValues) { 828 if (haveSurfaceAtoms && meshData.vertexSource == null) 829 meshData.vertexSource = new int[meshData.vc]; 830 float min = Float.MAX_VALUE; 831 float max = -Float.MAX_VALUE; 832 float value; 833 initializeMapping(); 834 for (int i = meshData.vc; --i >= meshData.mergeVertexCount0;) { 835 /* right, so what we are doing here is setting a range within the 836 * data for which we want red-->blue, but returning the actual 837 * number so it can be encoded more precisely. This turned out to be 838 * the key to making the JVXL contours work. 839 * 840 */ 841 if (params.colorBySets) { 842 value = meshData.vertexSets[i]; 843 } else if (params.colorByPhase) { 844 value = getPhase(meshData.vs[i]); 845 //else if (jvxlDataIs2dContour) 846 //marchingSquares 847 // .getInterpolatedPixelValue(meshData.vertices[i]); 848 } else { 849 boolean needSource = haveSurfaceAtoms;//(haveSurfaceAtoms && meshData.vertexSource[i] < 0); 850 value = volumeData.lookupInterpolatedVoxelValue(meshData.vs[i], needSource); 851 //System.out.println(i + " " + meshData.vertices[i] + " " + value); 852 if (needSource) 853 meshData.vertexSource[i] = getSurfaceAtomIndex(); 854 } 855 if (value < min) 856 min = value; 857 if (value > max && value != Float.MAX_VALUE) 858 max = value; 859 meshData.vvs[i] = value; 860 } 861 if (params.rangeSelected && minMax == null) 862 minMax = new float[] { min, max }; 863 finalizeMapping(); 864 } 865 params.setMapRanges(this, true); 866 jvxlData.mappedDataMin = params.mappedDataMin; 867 jvxlData.mappedDataMax = params.mappedDataMax; 868 jvxlData.valueMappedToRed = params.valueMappedToRed; 869 jvxlData.valueMappedToBlue = params.valueMappedToBlue; 870 if (params.contactPair == null && jvxlData.vertexColors == null) 871 colorData(); 872 873 JvxlCoder.jvxlCreateColorData(jvxlData, 874 (saveColorData ? meshData.vvs : null)); 875 876 if (haveSurfaceAtoms && meshDataServer != null) 877 meshDataServer.fillMeshData(meshData, MeshData.MODE_PUT_VERTICES, null); 878 879 if (meshDataServer != null && params.colorBySets) 880 meshDataServer.fillMeshData(meshData, MeshData.MODE_PUT_SETS, null); 881 } 882 colorData()883 private void colorData() { 884 885 float[] vertexValues = meshData.vvs; 886 short[] vertexColixes = meshData.vcs; 887 meshData.pcs = null; 888 float valueBlue = jvxlData.valueMappedToBlue; 889 float valueRed = jvxlData.valueMappedToRed; 890 short minColorIndex = jvxlData.minColorIndex; 891 short maxColorIndex = jvxlData.maxColorIndex; 892 if (params.colorEncoder == null) 893 params.colorEncoder = new ColorEncoder(null, null); 894 params.colorEncoder.setRange(params.valueMappedToRed, 895 params.valueMappedToBlue, params.isColorReversed); 896 for (int i = meshData.vc; --i >= 0;) { 897 float value = vertexValues[i]; 898 if (minColorIndex >= 0) { 899 if (value <= 0) 900 vertexColixes[i] = minColorIndex; 901 else if (value > 0) 902 vertexColixes[i] = maxColorIndex; 903 } else { 904 if (value <= valueRed) 905 value = valueRed; 906 if (value >= valueBlue) 907 value = valueBlue; 908 vertexColixes[i] = params.colorEncoder.getColorIndex(value); 909 } 910 } 911 912 if ((params.nContours > 0 || jvxlData.contourValues != null) && jvxlData.contourColixes == null) { 913 int n = (jvxlData.contourValues == null ? params.nContours : jvxlData.contourValues.length); 914 short[] colors = jvxlData.contourColixes = new short[n]; 915 float[] values = jvxlData.contourValues; 916 if (values == null) 917 values = jvxlData.contourValuesUsed; 918 if (jvxlData.contourValuesUsed == null) 919 jvxlData.contourValuesUsed = (values == null ? new float[n] : values); 920 float dv = (valueBlue - valueRed) / (n + 1); 921 // n + 1 because we want n lines between n + 1 slices 922 params.colorEncoder.setRange(params.valueMappedToRed, 923 params.valueMappedToBlue, params.isColorReversed); 924 for (int i = 0; i < n; i++) { 925 float v = (values == null ? valueRed + (i + 1) * dv : values[i]); 926 jvxlData.contourValuesUsed[i] = v; 927 colors[i] = C.getColixTranslucent(params.colorEncoder.getArgb(v)); 928 } 929 //TODO -- this strips translucency 930 jvxlData.contourColors = C.getHexCodes(colors); 931 } 932 } 933 934 private final static String[] colorPhases = { "_orb", "x", "y", "z", "xy", 935 "yz", "xz", "x2-y2", "z2" }; 936 getColorPhaseIndex(String color)937 static int getColorPhaseIndex(String color) { 938 int colorPhase = -1; 939 for (int i = 0; i < colorPhases.length; i++) 940 if (color.equalsIgnoreCase(colorPhases[i])) { 941 colorPhase = i; 942 break; 943 } 944 return colorPhase; 945 } 946 getPhase(T3 pt)947 private float getPhase(T3 pt) { 948 switch (params.colorPhase) { 949 case 0: 950 case -1: 951 case 1: 952 return (pt.x > 0 ? 1 : -1); 953 case 2: 954 return (pt.y > 0 ? 1 : -1); 955 case 3: 956 return (pt.z > 0 ? 1 : -1); 957 case 4: 958 return (pt.x * pt.y > 0 ? 1 : -1); 959 case 5: 960 return (pt.y * pt.z > 0 ? 1 : -1); 961 case 6: 962 return (pt.x * pt.z > 0 ? 1 : -1); 963 case 7: 964 return (pt.x * pt.x - pt.y * pt.y > 0 ? 1 : -1); 965 case 8: 966 return (pt.z * pt.z * 2f - pt.x * pt.x - pt.y * pt.y > 0 ? 1 : -1); 967 } 968 return 1; 969 } 970 971 protected float[] minMax; 972 getMinMaxMappedValues(boolean haveData)973 public float[] getMinMaxMappedValues(boolean haveData) { 974 if (minMax != null && minMax[0] != Float.MAX_VALUE) 975 return minMax; 976 if (params.colorBySets) 977 return (minMax = new float[] { 0, Math.max(meshData.nSets - 1, 0) }); 978 float min = Float.MAX_VALUE; 979 float max = -Float.MAX_VALUE; 980 if (params.usePropertyForColorRange && params.theProperty != null) { 981 for (int i = params.theProperty.length; --i >= 0;) { 982 if (params.rangeSelected && !params.bsSelected.get(i)) 983 continue; 984 float p = params.theProperty[i]; 985 if (Float.isNaN(p)) 986 continue; 987 if (p < min) 988 min = p; 989 if (p > max) 990 max = p; 991 } 992 return (minMax = new float[] { min, max }); 993 } 994 int vertexCount = (contourVertexCount > 0 ? contourVertexCount 995 : meshData.vc); 996 T3[] vertexes = meshData.vs; 997 boolean useVertexValue = (haveData || jvxlDataIs2dContour || vertexDataOnly || params.colorDensity); 998 for (int i = meshData.mergeVertexCount0; i < vertexCount; i++) { 999 float v; 1000 if (useVertexValue) 1001 v = meshData.vvs[i]; 1002 else 1003 v = volumeData.lookupInterpolatedVoxelValue(vertexes[i], false); 1004 if (v < min) 1005 min = v; 1006 if (v > max && v != Float.MAX_VALUE) 1007 max = v; 1008 } 1009 return (minMax = new float[] { min, max }); 1010 } 1011 updateTriangles()1012 void updateTriangles() { 1013 if (meshDataServer == null) { 1014 meshData.invalidatePolygons(); 1015 } else { 1016 meshDataServer.invalidateTriangles(); 1017 } 1018 } 1019 updateSurfaceData()1020 void updateSurfaceData() { 1021 meshData.setVertexSets(true); 1022 updateTriangles(); 1023 if (params.bsExcluded[1] == null) 1024 params.bsExcluded[1] = new BS(); 1025 meshData.updateInvalidatedVertices(params.bsExcluded[1]); 1026 } 1027 1028 /** 1029 * 1030 * @param doExclude 1031 */ selectPocket(boolean doExclude)1032 public void selectPocket(boolean doExclude) { 1033 // solvent reader implements this 1034 } 1035 excludeMinimumSet()1036 void excludeMinimumSet() { 1037 if (meshDataServer != null) 1038 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); 1039 meshData.getSurfaceSet(); 1040 BS bs; 1041 for (int i = meshData.nSets; --i >= 0;) 1042 if ((bs = meshData.surfaceSet[i]) != null 1043 && bs.cardinality() < params.minSet) 1044 meshData.invalidateSurfaceSet(i); 1045 updateSurfaceData(); 1046 if (meshDataServer != null) 1047 meshDataServer.fillMeshData(meshData, MeshData.MODE_PUT_SETS, null); 1048 } 1049 excludeMaximumSet()1050 void excludeMaximumSet() { 1051 if (meshDataServer != null) 1052 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); 1053 meshData.getSurfaceSet(); 1054 BS bs; 1055 for (int i = meshData.nSets; --i >= 0;) 1056 if ((bs = meshData.surfaceSet[i]) != null 1057 && bs.cardinality() > params.maxSet) 1058 meshData.invalidateSurfaceSet(i); 1059 updateSurfaceData(); 1060 if (meshDataServer != null) 1061 meshDataServer.fillMeshData(meshData, MeshData.MODE_PUT_SETS, null); 1062 } 1063 slabIsosurface(Lst<Object[]> slabInfo)1064 public void slabIsosurface(Lst<Object[]> slabInfo) { 1065 if (meshDataServer != null) 1066 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); 1067 meshData.slabPolygonsList(slabInfo, true); 1068 if (meshDataServer != null) 1069 meshDataServer.fillMeshData(meshData, MeshData.MODE_PUT_VERTICES, null); 1070 } 1071 setVertexAnisotropy(T3 pt)1072 protected void setVertexAnisotropy(T3 pt) { 1073 pt.x *= anisotropy[0]; 1074 pt.y *= anisotropy[1]; 1075 pt.z *= anisotropy[2]; 1076 pt.add(center); 1077 } 1078 setVectorAnisotropy(T3 v)1079 protected void setVectorAnisotropy(T3 v) { 1080 haveSetAnisotropy = true; 1081 v.x *= anisotropy[0]; 1082 v.y *= anisotropy[1]; 1083 v.z *= anisotropy[2]; 1084 } 1085 1086 private boolean haveSetAnisotropy; 1087 setVolumetricAnisotropy()1088 protected void setVolumetricAnisotropy() { 1089 if (haveSetAnisotropy) 1090 return; 1091 setVolumetricOriginAnisotropy(); 1092 setVectorAnisotropy(volumetricVectors[0]); 1093 setVectorAnisotropy(volumetricVectors[1]); 1094 setVectorAnisotropy(volumetricVectors[2]); 1095 1096 } 1097 setVolumetricOriginAnisotropy()1098 protected void setVolumetricOriginAnisotropy() { 1099 volumetricOrigin.setT(center); 1100 } 1101 setBBoxAll()1102 private void setBBoxAll() { 1103 if (meshDataServer != null) 1104 meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); 1105 xyzMin = new P3(); 1106 xyzMax = new P3(); 1107 meshData.setBox(xyzMin, xyzMax); 1108 } 1109 setBBox(T3 pt, float margin)1110 protected void setBBox(T3 pt, float margin) { 1111 if (xyzMin == null) { 1112 xyzMin = P3.new3(Float.MAX_VALUE, Float.MAX_VALUE, Float.MAX_VALUE); 1113 xyzMax = P3.new3(-Float.MAX_VALUE, -Float.MAX_VALUE, -Float.MAX_VALUE); 1114 } 1115 BoxInfo.addPoint(pt, xyzMin, xyzMax, margin); 1116 } 1117 1118 /** 1119 * 1120 * @param pt 1121 * @param getSource TODO 1122 * @return value 1123 */ getValueAtPoint(T3 pt, boolean getSource)1124 public float getValueAtPoint(T3 pt, boolean getSource) { 1125 // only for readers that can support it (IsoShapeReader, AtomPropertyMapper) 1126 return 0; 1127 } 1128 initializeMapping()1129 void initializeMapping() { 1130 // initiate any iterators 1131 } 1132 finalizeMapping()1133 protected void finalizeMapping() { 1134 // release any iterators 1135 } 1136 getSurfaceAtomIndex()1137 public int getSurfaceAtomIndex() { 1138 // atomPropertyMapper 1139 return -1; 1140 } 1141 1142 } 1143 1144