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