From: <ha...@us...> - 2011-01-02 21:02:28
|
Revision: 14929 http://jmol.svn.sourceforge.net/jmol/?rev=14929&view=rev Author: hansonr Date: 2011-01-02 21:02:20 +0000 (Sun, 02 Jan 2011) Log Message: ----------- 12.1.29 VERSION (really!) # new feature: isosurface SLAB within x.y {point} # -- smooth circular slabbing # -- negative x.y --> "not within" (creates a circular hole at this position) # code: isosurface MOLECULAR/SOLVENT x.y completed; "FULL" option removed # code: (new concept) nonlinear Marching Cubes for IsoSolventReader # bug fix: isosurface molecular/solvent improvements # bug fix: logFile/logCommands/logGestures should not be saved in state # bug fix: log file cannot be turned off by setting value to "" # bug fix: user-defined function calls should be able to look like Jmol commands, without () Modified Paths: -------------- trunk/Jmol/src/org/jmol/jvxl/data/JvxlData.java trunk/Jmol/src/org/jmol/jvxl/data/MeshData.java trunk/Jmol/src/org/jmol/jvxl/readers/IsoSolventReader.java trunk/Jmol/src/org/jmol/jvxl/readers/SurfaceReader.java trunk/Jmol/src/org/jmol/script/ScriptEvaluator.java trunk/Jmol/src/org/jmol/shapespecial/Draw.java trunk/Jmol/src/org/jmol/shapesurface/Isosurface.java trunk/Jmol/src/org/jmol/util/MeshSurface.java trunk/Jmol/src/org/jmol/viewer/Jmol.properties Modified: trunk/Jmol/src/org/jmol/jvxl/data/JvxlData.java =================================================================== --- trunk/Jmol/src/org/jmol/jvxl/data/JvxlData.java 2011-01-02 19:19:17 UTC (rev 14928) +++ trunk/Jmol/src/org/jmol/jvxl/data/JvxlData.java 2011-01-02 21:02:20 UTC (rev 14929) @@ -168,14 +168,16 @@ this.nBytes = nBytes; } - public void updateSurfaceData(float[] vertexValues, int vertexCount, int vertexIncrement, char isNaN) { - if (jvxlEdgeData.length() == 0) - return; - char[] chars = jvxlEdgeData.toCharArray(); - for (int i = 0, ipt = 0; i < vertexCount; i+= vertexIncrement, ipt++) + public static String updateSurfaceData(String edgeData, float[] vertexValues, + int vertexCount, int vertexIncrement, + char isNaN) { + if (edgeData.length() == 0) + return ""; + char[] chars = edgeData.toCharArray(); + for (int i = 0, ipt = 0; i < vertexCount; i += vertexIncrement, ipt++) if (Float.isNaN(vertexValues[i])) - chars[ipt] = isNaN; - jvxlEdgeData = String.copyValueOf(chars); + chars[ipt] = isNaN; + return String.copyValueOf(chars); } } Modified: trunk/Jmol/src/org/jmol/jvxl/data/MeshData.java =================================================================== --- trunk/Jmol/src/org/jmol/jvxl/data/MeshData.java 2011-01-02 19:19:17 UTC (rev 14928) +++ trunk/Jmol/src/org/jmol/jvxl/data/MeshData.java 2011-01-02 21:02:20 UTC (rev 14929) @@ -201,6 +201,32 @@ return surfaceSet; } + private class SSet { + BitSet bs; + int n; + + protected SSet (BitSet bs) { + this.bs = bs; + n = bs.cardinality(); + } + } + + protected class SortSet implements Comparator<SSet> { + + public int compare(SSet o1, SSet o2) { + return (o1.n > o2.n ? -1 : o1.n < o2.n ? 1 : 0); + } + } + + private void sortSurfaceSets() { + SSet[] sets = new SSet[nSets]; + for (int i = 0; i < nSets; i++) + sets[i] = new SSet(surfaceSet[i]); + Arrays.sort(sets, new SortSet()); + for (int i = 0; i < nSets; i++) + surfaceSet[i] = sets[i].bs; + } + public void setVertexSets(boolean onlyIfNull) { if (surfaceSet == null) return; @@ -228,32 +254,6 @@ vertexSets[j] = i; } - private class SSet { - BitSet bs; - int n; - - protected SSet (BitSet bs) { - this.bs = bs; - n = bs.cardinality(); - } -} - - protected class SortSet implements Comparator<SSet> { - - public int compare(SSet o1, SSet o2) { - return (o1.n > o2.n ? -1 : o1.n < o2.n ? 1 : 0); - } - } - - private void sortSurfaceSets() { - SSet[] sets = new SSet[nSets]; - for (int i = 0; i < nSets; i++) - sets[i] = new SSet(surfaceSet[i]); - Arrays.sort(sets, new SortSet()); - for (int i = 0; i < nSets; i++) - surfaceSet[i] = sets[i].bs; - } - private int findSet(int vertex) { for (int i = 0; i < nSets; i++) if (surfaceSet[i] != null && surfaceSet[i].get(vertex)) Modified: trunk/Jmol/src/org/jmol/jvxl/readers/IsoSolventReader.java =================================================================== --- trunk/Jmol/src/org/jmol/jvxl/readers/IsoSolventReader.java 2011-01-02 19:19:17 UTC (rev 14928) +++ trunk/Jmol/src/org/jmol/jvxl/readers/IsoSolventReader.java 2011-01-02 21:02:20 UTC (rev 14929) @@ -40,6 +40,7 @@ import org.jmol.util.BitSetUtil; import org.jmol.util.Logger; import org.jmol.util.Measure; +import org.jmol.util.MeshSurface; import org.jmol.api.AtomIndexIterator; import org.jmol.jvxl.data.MeshData; @@ -134,6 +135,7 @@ private float cavityRadius; private float envelopeRadius; + private Point3f[] dots; private boolean doCalculateTroughs; private boolean isCavity, isPocket; @@ -176,7 +178,7 @@ doUseIterator = doCalculateTroughs; getAtoms(Float.NaN, false, true); if (isCavity || isPocket) - meshData.dots = meshDataServer.calculateGeodesicSurface(bsMySelected, + dots = meshDataServer.calculateGeodesicSurface(bsMySelected, envelopeRadius); setHeader("solvent/molecular surface", params.calculationType); @@ -328,30 +330,15 @@ int vB = marchingCubes.getLinearOffset(x, y, z, vB0); isSurfacePoint = (bsSurfaceVoxels != null && (bsSurfaceVoxels.get(vA) || bsSurfaceVoxels .get(vB))); - + if (testLinear || voxelRadii == null || voxelRadii[vA] == 0 || voxelRadii[vA] != voxelRadii[vB]) return super.getSurfacePointAndFraction(cutoff, isCutoffAbsolute, valueA, valueB, pointA, edgeVector, x, y, z, vA, vB, fReturn, ptReturn); - float d = edgeVector.length(); - double r = voxelRadii[vA]; - double ra = Math.abs(r + valueA) / d; - double rb = Math.abs(r + valueB) / d; - double ra2 = ra * ra; - double q = ra2 - rb * rb + 1; - r /= d; - double p = 4 * (r * r - ra2); - double factor = (ra < rb ? 1 : -1); - float diff = valueB - valueA; - - float fraction = (float) (((q) + factor * Math.sqrt(q * q + p)) / 2); - if (fraction < 0 || fraction > 1) { - Logger.error("problem with unusual fraction=" + fraction + " cutoff=" - + cutoff + " A:" + valueA + " B:" + valueB); - fraction = Float.NaN; - } - fReturn[0] = fraction; + float fraction = fReturn[0] = MeshSurface.getSphericalInterpolationFraction( + voxelRadii[vA], valueA, valueB, edgeVector.length()); ptReturn.scaleAdd(fraction, edgeVector, pointA); + float diff = valueB - valueA; /* //for testing -- linear with ttest.xyz: // float fractionLinear = (cutoff - valueA) / diff; @@ -377,14 +364,12 @@ @Override public void selectPocket(boolean doExclude) { - if (meshDataServer == null) - return; //can't do this without help! - meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); + if (meshDataServer != null) + meshDataServer.fillMeshData(meshData, MeshData.MODE_GET_VERTICES, null); //mark VERTICES for proximity to surface Point3f[] v = meshData.vertices; int nVertices = meshData.vertexCount; float[] vv = meshData.vertexValues; - Point3f[] dots = meshData.dots; int nDots = dots.length; for (int i = 0; i < nVertices; i++) { for (int j = 0; j < nDots; j++) { @@ -415,8 +400,10 @@ updateSurfaceData(); if (!doExclude) meshData.surfaceSet = null; - meshDataServer.fillMeshData(meshData, MeshData.MODE_PUT_SETS, null); - meshData = new MeshData(); + if (meshDataServer != null) { + meshDataServer.fillMeshData(meshData, MeshData.MODE_PUT_SETS, null); + meshData = new MeshData(); + } } @Override @@ -447,7 +434,7 @@ //3) rerun the calculation to mark a solvent around these! BitSet bs = new BitSet(nPointsX * nPointsY * nPointsZ); int i = 0; - int nDots = meshData.dots.length; + int nDots = dots.length; int n = 0; float d; float r2 = envelopeRadius;// - cavityRadius; @@ -458,7 +445,7 @@ if ((d = voxelData[x][y][z]) < Float.MAX_VALUE && d >= cavityRadius) { volumeData.voxelPtToXYZ(x, y, z, ptXyzTemp); for (int j = 0; j < nDots; j++) { - if (meshData.dots[j].distance(ptXyzTemp) < r2) + if (dots[j].distance(ptXyzTemp) < r2) continue out; } bs.set(i); Modified: trunk/Jmol/src/org/jmol/jvxl/readers/SurfaceReader.java =================================================================== --- trunk/Jmol/src/org/jmol/jvxl/readers/SurfaceReader.java 2011-01-02 19:19:17 UTC (rev 14928) +++ trunk/Jmol/src/org/jmol/jvxl/readers/SurfaceReader.java 2011-01-02 21:02:20 UTC (rev 14929) @@ -886,7 +886,7 @@ void updateSurfaceData() { meshData.setVertexSets(true); updateTriangles(); - jvxlData.updateSurfaceData(meshData.vertexValues, + jvxlData.jvxlEdgeData = edgeData = JvxlData.updateSurfaceData(edgeData, meshData.vertexValues, meshData.vertexCount, meshData.vertexIncrement, cJvxlEdgeNaN); } Modified: trunk/Jmol/src/org/jmol/script/ScriptEvaluator.java =================================================================== --- trunk/Jmol/src/org/jmol/script/ScriptEvaluator.java 2011-01-02 19:19:17 UTC (rev 14928) +++ trunk/Jmol/src/org/jmol/script/ScriptEvaluator.java 2011-01-02 21:02:20 UTC (rev 14929) @@ -14780,10 +14780,15 @@ sb.append(" ").append(getToken(i).value).append(" "); int tok = tokAt(i + 1); Point4f plane = null; + float d = 0; + Point3f pt = null; switch (tok) { case Token.within: i++; - data = getPointArray(++i, 4); + if (isFloatParameter(++i)) + data = new Object[] { Float.valueOf(d = floatParameter(i)), pt = centerParameter(++i) }; + else + data = getPointArray(i, 4); break; case Token.boundbox: data = BoxInfo.getCriticalPoints(viewer.getBoundBoxVertices(), null); @@ -14831,10 +14836,12 @@ data = plane; } if (sb != null) { - if (plane == null) + if (plane != null) + sb.append(Escape.escape(plane)); + else if (pt != null) + sb.append("within ").append(d).append(" ").append(Escape.escape(pt)); + else sb.append("within ").append(Escape.escape(data)); - else - sb.append(Escape.escape(plane)); } return data; } Modified: trunk/Jmol/src/org/jmol/shapespecial/Draw.java =================================================================== --- trunk/Jmol/src/org/jmol/shapespecial/Draw.java 2011-01-02 19:19:17 UTC (rev 14928) +++ trunk/Jmol/src/org/jmol/shapespecial/Draw.java 2011-01-02 21:02:20 UTC (rev 14929) @@ -631,7 +631,7 @@ private void setSlabData() { if (plane != null) { - slabData.getIntersection(plane, null, false, true); + slabData.getIntersection(plane, null, 0, null, false, true); polygon = new ArrayList<Object>(); polygon.add(slabData.vertices); polygon.add(slabData.polygonIndexes); Modified: trunk/Jmol/src/org/jmol/shapesurface/Isosurface.java =================================================================== --- trunk/Jmol/src/org/jmol/shapesurface/Isosurface.java 2011-01-02 19:19:17 UTC (rev 14928) +++ trunk/Jmol/src/org/jmol/shapesurface/Isosurface.java 2011-01-02 21:02:20 UTC (rev 14929) @@ -521,7 +521,7 @@ if (mesh == null) return false; data[3] = Integer.valueOf(mesh.modelIndex); - return mesh.getIntersection((Point4f) data[1], (List<Point3f[]>) data[2], false, false); + return mesh.getIntersection((Point4f) data[1], null, 0, (List<Point3f[]>) data[2], false, false); } if (property == "getBoundingBox") { String id = (String) data[0]; Modified: trunk/Jmol/src/org/jmol/util/MeshSurface.java =================================================================== --- trunk/Jmol/src/org/jmol/util/MeshSurface.java 2011-01-02 19:19:17 UTC (rev 14928) +++ trunk/Jmol/src/org/jmol/util/MeshSurface.java 2011-01-02 21:02:20 UTC (rev 14929) @@ -28,7 +28,6 @@ public BitSet bsFaces; public Point3f ptOffset; public float scale3d; - public Point3f[] dots; public float[] vertexValues; public BitSet[] surfaceSet; public int[] vertexSets; @@ -155,32 +154,55 @@ public void slabPolygons(Object slabbingObject, boolean andCap) { if (slabbingObject instanceof Point4f) { - getIntersection((Point4f) slabbingObject, null, andCap, false); + getIntersection((Point4f) slabbingObject, null, 0, null, andCap, false); return; } if (slabbingObject instanceof Point3f[]) { Point4f[] faces = BoxInfo.getFacesFromCriticalPoints((Point3f[]) slabbingObject); for (int i = 0; i < faces.length; i++) - getIntersection(faces[i], null, andCap, false); + getIntersection(faces[i], null, 0, null, andCap, false); return; } + if (slabbingObject instanceof Object[]) { + Object[] o = (Object[]) slabbingObject; + float distance = ((Float) o[0]).floatValue(); + Point3f center = (Point3f) o[1]; + getIntersection(null, center, distance, null, andCap, false); + } } - public boolean getIntersection(Point4f plane, List<Point3f[]> vData, boolean andCap, boolean doClean) { + public boolean getIntersection(Point4f plane, Point3f ptCenter, + float distance, List<Point3f[]> vData, + boolean andCap, boolean doClean) { boolean isSlab = (vData == null); + if (ptCenter != null) + andCap = false; // can only cap faces Point3f[] pts; int iD, iE; - + double absD = Math.abs(distance); List<int[]> iPts = (andCap ? new ArrayList<int[]>() : null); for (int i = polygonIndexes.length; --i >= 0;) { if (!setABC(i)) continue; int check1 = polygonIndexes[i][3]; - int check2 = (checkCount == 2 ? polygonIndexes[i][4]: 0); - Point3f vA, vB, vC; - float d1 = Measure.distanceToPlane(plane, vA = vertices[iA]); - float d2 = Measure.distanceToPlane(plane, vB = vertices[iB]); - float d3 = Measure.distanceToPlane(plane, vC = vertices[iC]); + int check2 = (checkCount == 2 ? polygonIndexes[i][4] : 0); + Point3f vA = vertices[iA]; + Point3f vB = vertices[iB]; + Point3f vC = vertices[iC]; + float d1, d2, d3; + if (ptCenter == null) { + d1 = Measure.distanceToPlane(plane, vA); + d2 = Measure.distanceToPlane(plane, vB); + d3 = Measure.distanceToPlane(plane, vC); + } else if (distance < 0) { + d1 = -vA.distance(ptCenter) - distance; + d2 = -vB.distance(ptCenter) - distance; + d3 = -vC.distance(ptCenter) - distance; + } else { + d1 = vA.distance(ptCenter) - distance; + d2 = vB.distance(ptCenter) - distance; + d3 = vC.distance(ptCenter) - distance; + } int test1 = (d1 < 0 ? 1 : 0) + (d2 < 0 ? 2 : 0) + (d3 < 0 ? 4 : 0); int test2 = (d1 >= 0 ? 1 : 0) + (d2 >= 0 ? 2 : 0) + (d3 >= 0 ? 4 : 0); pts = null; @@ -192,20 +214,32 @@ case 1: case 6: // BC on same side - pts = new Point3f[] { interpolatePoint(vA, vB, -d1, d2), - interpolatePoint(vA, vC, -d1, d3)}; + if (ptCenter == null) + pts = new Point3f[] { interpolatePoint(vA, vB, -d1, d2), + interpolatePoint(vA, vC, -d1, d3) }; + else + pts = new Point3f[] { interpolateSphere(vA, vB, -d1, -d2, absD), + interpolateSphere(vA, vC, -d1, -d3, absD) }; break; case 2: case 5: //AC on same side - pts = new Point3f[] { interpolatePoint(vB, vA, -d2, d1), - interpolatePoint(vB, vC, -d2, d3)}; + if (ptCenter == null) + pts = new Point3f[] { interpolatePoint(vB, vA, -d2, d1), + interpolatePoint(vB, vC, -d2, d3) }; + else + pts = new Point3f[] { interpolateSphere(vB, vA, -d2, -d1, absD), + interpolateSphere(vB, vC, -d2, -d3, absD) }; break; case 3: case 4: //AB on same side need A-C, B-C - pts = new Point3f[] { interpolatePoint(vC, vA, -d3, d1), - interpolatePoint(vC, vB, -d3, d2)}; + if (ptCenter == null) + pts = new Point3f[] { interpolatePoint(vC, vA, -d3, d1), + interpolatePoint(vC, vB, -d3, d2) }; + else + pts = new Point3f[] { interpolateSphere(vC, vA, -d3, -d1, absD), + interpolateSphere(vC, vB, -d3, -d2, absD) }; break; } if (isSlab) { @@ -222,35 +256,35 @@ break; case 1: // BC on side to keep - iD = addVertexCopy(pts[1], vertexValues[iA]); //AC - iE = addVertexCopy(pts[0], vertexValues[iA]); //AB + iD = addVertexCopy(pts[1], vertexValues[iA]); //AC + iE = addVertexCopy(pts[0], vertexValues[iA]); //AB addTriangleCheck(iE, iB, iC, check1 & 3, check2, 0); addTriangleCheck(iD, iE, iC, check1 & 4 | 1, check2, 0); break; case 2: // AC on side to keep - iD = addVertexCopy(pts[0], vertexValues[iB]); //AB - iE = addVertexCopy(pts[1], vertexValues[iB]); //BC + iD = addVertexCopy(pts[0], vertexValues[iB]); //AB + iE = addVertexCopy(pts[1], vertexValues[iB]); //BC addTriangleCheck(iA, iD, iC, check1 & 5, check2, 0); addTriangleCheck(iD, iE, iC, check1 & 2 | 1, check2, 0); break; case 3: //AB on side to toss - iD = addVertexCopy(pts[0], vertexValues[iA]); //AC - iE = addVertexCopy(pts[1], vertexValues[iB]); //BC + iD = addVertexCopy(pts[0], vertexValues[iA]); //AC + iE = addVertexCopy(pts[1], vertexValues[iB]); //BC addTriangleCheck(iD, iE, iC, check1 & 6 | 1, check2, 0); break; case 4: //AB on side to keep - iD = addVertexCopy(pts[1], vertexValues[iC]); //BC - iE = addVertexCopy(pts[0], vertexValues[iC]); //AC + iD = addVertexCopy(pts[1], vertexValues[iC]); //BC + iE = addVertexCopy(pts[0], vertexValues[iC]); //AC addTriangleCheck(iA, iB, iE, check1 & 5, check2, 0); addTriangleCheck(iE, iB, iD, check1 & 2 | 4, check2, 0); break; case 5: //AC on side to toss - iD = addVertexCopy(pts[1], vertexValues[iC]); //BC - iE = addVertexCopy(pts[0], vertexValues[iA]); //AB + iD = addVertexCopy(pts[1], vertexValues[iC]); //BC + iE = addVertexCopy(pts[0], vertexValues[iA]); //AB addTriangleCheck(iE, iB, iD, check1 & 3 | 4, check2, 0); break; case 6: @@ -262,7 +296,7 @@ } polygonIndexes[i] = null; if (andCap && iD > 0) - iPts.add(new int[] {iD, iE}); + iPts.add(new int[] { iD, iE }); } else if (pts != null) { vData.add(pts); } @@ -322,15 +356,36 @@ return false; } + private Point3f interpolateSphere(Point3f v1, Point3f v2, float d1, float d2, + double absD) { + return interpolateFraction(v1, v2, getSphericalInterpolationFraction(absD, d1, + d2, v1.distance(v2))); + } + private static Point3f interpolatePoint(Point3f v1, Point3f v2, float d1, float d2) { - float f = d1 / (d1 + d2); + return interpolateFraction(v1, v2, d1 / (d1 + d2)); + } + + private static Point3f interpolateFraction(Point3f v1, Point3f v2, float f) { if (f < 0.0001) f = 0; else if (f > 0.9999) f = 1; return new Point3f(v1.x + (v2.x - v1.x) * f, v1.y + (v2.y - v1.y) * f, - v1.z + (v2.z - v1.z) * f); + v1.z + (v2.z - v1.z) * f); } + + public static float getSphericalInterpolationFraction(double r, double valueA, + double valueB, double d) { + double ra = Math.abs(r + valueA) / d; + double rb = Math.abs(r + valueB) / d; + r /= d; + double ra2 = ra * ra; + double q = ra2 - rb * rb + 1; + double p = 4 * (r * r - ra2); + double factor = (ra < rb ? 1 : -1); + return (float) (((q) + factor * Math.sqrt(q * q + p)) / 2); + } } Modified: trunk/Jmol/src/org/jmol/viewer/Jmol.properties =================================================================== --- trunk/Jmol/src/org/jmol/viewer/Jmol.properties 2011-01-02 19:19:17 UTC (rev 14928) +++ trunk/Jmol/src/org/jmol/viewer/Jmol.properties 2011-01-02 21:02:20 UTC (rev 14929) @@ -1,12 +1,11 @@ # Developers: to add a description of changes you have made, # add it on a line starting with # below the "version=..." line -version=12.1.29_dev +version=12.1.29 -# TODO: fix jvxl write of new molecular to exclude unreal surface fragments # new feature: isosurface SLAB within x.y {point} # -- smooth circular slabbing -# -- negative x.y --> "not within" (creates a hole at this position) +# -- negative x.y --> "not within" (creates a circular hole at this position) # code: isosurface MOLECULAR/SOLVENT x.y completed; "FULL" option removed # code: (new concept) nonlinear Marching Cubes for IsoSolventReader # bug fix: isosurface molecular/solvent improvements This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |