From: <ha...@us...> - 2006-06-27 04:21:20
|
Revision: 5244 Author: hansonr Date: 2006-06-26 21:21:05 -0700 (Mon, 26 Jun 2006) ViewCVS: http://svn.sourceforge.net/jmol/?rev=5244&view=rev Log Message: ----------- bob200603 isosurface development -- lobes, scale, eccentricity Modified Paths: -------------- branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java Modified: branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java =================================================================== --- branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java 2006-06-26 05:27:06 UTC (rev 5243) +++ branches/bob200603/Jmol/src/org/jmol/viewer/Eval.java 2006-06-27 04:21:05 UTC (rev 5244) @@ -4160,9 +4160,9 @@ break; } if (str.equalsIgnoreCase("lobe")) { - //lobe [size] + //lobe {eccentricity} propertyName = "lobe"; - propertyValue = new Float(floatParameter(++i)); + propertyValue = getPoint4f(i + 1); break; } if (str.equalsIgnoreCase("ellipsoid")) { @@ -4207,6 +4207,11 @@ propertyName = "remappable"; break; } + if (str.equalsIgnoreCase("SCALE")) { + propertyName = "scale"; + propertyValue = new Float(floatParameter(++i)); + break; + } if (str.equalsIgnoreCase("ANISOTROPY")) { propertyName = "anisotropy"; propertyValue = getCoordinate(i + 1, false); @@ -4214,8 +4219,8 @@ break; } if (str.equalsIgnoreCase("ECCENTRICITY")) { - propertyName = "anisotropy"; - propertyValue = getCoordinate(i + 1, false); + propertyName = "eccentricity"; + propertyValue = getPoint4f(i + 1); i = pcLastExpressionInstruction; break; } Modified: branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java =================================================================== --- branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java 2006-06-26 05:27:06 UTC (rev 5243) +++ branches/bob200603/Jmol/src/org/jmol/viewer/Isosurface.java 2006-06-27 04:21:05 UTC (rev 5244) @@ -135,7 +135,10 @@ boolean isAnisotropic; boolean isEccentric; + float eccentricityScale; + float scale; Matrix3f eccentricityMatrix; + Matrix3f eccentricityMatrixInverse; boolean isSolvent; float solventRadius; @@ -213,6 +216,7 @@ colorPos = defaultColorPositive; isSphere = isOrbital = isLobe = isSolvent = isFunctionXY = false; isEccentric = isAnisotropic = false; + scale = Float.NaN; center = new Point3f(Float.MAX_VALUE, Float.MAX_VALUE, Float.MAX_VALUE); //anisotropy[0] = anisotropy[1] = anisotropy[2] = 1f; cutoff = Float.MAX_VALUE; @@ -228,6 +232,11 @@ return; } + if ("scale" == propertyName) { + scale = ((Float) value).floatValue(); + return; + } + if ("center" == propertyName) { center.set((Point3f) value); return; @@ -267,6 +276,10 @@ return; } + if ("eccentricity" == propertyName) { + setEccentricity((Point4f) value); + } + /// JVXL options if ("fileIndex" == propertyName) { @@ -329,8 +342,9 @@ boolean isCalculation = false; if ("sphere" == propertyName) { + sphere_radiusAngstroms = ((Float) value).floatValue(); isSphere = isCalculation = true; - sphere_radiusAngstroms = ((Float) value).floatValue(); + setEccentricity(new Point4f(0, 0, 1, 1)); cutoff = Float.MIN_VALUE; propertyName = "getSurface"; } @@ -364,21 +378,11 @@ } if ("lobe" == propertyName) { + setEccentricity((Point4f) value); isLobe = isCalculation = true; - float scale = lobe_sizeAngstroms = ((Float) value).floatValue(); if (cutoff == Float.MAX_VALUE) cutoff = defaultOrbitalCutoff; - if (isAnisotropic) { - anisotropy[0] *= scale; - anisotropy[1] *= scale; - anisotropy[2] *= scale; - } else { - anisotropy[2] = scale; - anisotropy[0] = anisotropy[1] = 0.5f * scale; - isAnisotropic = true; - } propertyName = "getSurface"; - } if ("solvent" == propertyName) { @@ -484,8 +488,11 @@ void setEccentricity(Point4f info) { Vector3f ecc = new Vector3f(info.x, info.y, info.z); - float c = ecc.length(); - float ab = info.w * c; + float c = (scale > 0 ? scale : info.w < 0 ? 1f : ecc.length()); + eccentricityScale = c; + float fab_c = Math.abs(info.w); + if (fab_c > 1) + eccentricityScale *= fab_c; ecc.normalize(); Vector3f z = new Vector3f(0, 0, 1); ecc.add(z); @@ -496,13 +503,14 @@ eccentricityMatrix = new Matrix3f(); eccentricityMatrix.setIdentity(); eccentricityMatrix.set(eccA); + eccentricityMatrixInverse = new Matrix3f(); + eccentricityMatrixInverse.invert(eccentricityMatrix); isEccentric = isAnisotropic = true; - anisotropy[0] = ab; - anisotropy[1] = ab; + anisotropy[0] = fab_c * c; + anisotropy[1] = fab_c * c; anisotropy[2] = c; if (center.x == Float.MAX_VALUE) center.set(0, 0, 0); - } void colorIsosurface() { @@ -574,8 +582,8 @@ jvxlReadColorData(currentMesh, "roygb"); currentMesh.colix = getDefaultColix(); currentMesh.jvxlDefinitionLine = jvxlDefinitionLine(1); - if (thePlane != null && iAddGridPoints) - addGridPoints(); + if (thePlane != null && iAddGridPoints) + addGridPointCube(); } void resetIsosurface() { @@ -596,7 +604,7 @@ edgeCount = 0; } - void addGridPoints() { + void addGridPointCube() { for (int x = 0; x < voxelCounts[0]; x += 5) for (int y = 0; y < voxelCounts[1]; y += 5) for (int z = 0; z < voxelCounts[2]; z += 5) { @@ -892,6 +900,9 @@ // update surfaceData + if (logMessages) + System.out.println(" cutoff=" + cutoff + " (" + x + " " + y + " " + + z + ") value=" + voxelValue); boolean isInside = ((cutoff > 0 && voxelValue >= cutoff) || (cutoff <= 0 && voxelValue <= cutoff)); if (inside == isInside) { dataCount++; @@ -1586,8 +1597,10 @@ float fraction = (ich - base + fracOffset) / range; if (fraction < 0f) fraction = 0f; - if (fraction >= 1f) + if (fraction > 1f) fraction = 0.999999f; + if (ich == base + range) + fraction = Float.NaN; if (logCompression) System.out.println("ffc: " + fraction + " <-- " + ich + " " + (char) ich); return fraction; @@ -1602,6 +1615,8 @@ char jvxlFractionAsCharacter(float fraction, int base, int range) { if (fraction > 0.9999f) fraction = 0.9999f; + if (Float.isNaN(fraction)) + fraction = 1.0001f; int ich = (int) (fraction * range + base); if (ich < base) ich = base; @@ -1822,7 +1837,6 @@ float vertexValue = voxelData[x + offset.x][y + offset.y][z + offset.z]; vertexValues[i] = vertexValue; - if ((cutoff > 0 && vertexValue >= cutoff) || (cutoff <= 0 && vertexValue <= cutoff)) insideMask |= 1 << i; @@ -1837,7 +1851,7 @@ continue; } ++surfaceCount; - processOneVoxel(insideMask, cutoff, voxelPointIndexes, x, y, z, + processOneCubical(insideMask, cutoff, voxelPointIndexes, x, y, z, contouring); } } @@ -1927,9 +1941,10 @@ int edgeCount; final static float assocCutoff = 0.3f; - void processOneVoxel(int insideMask, float cutoff, int[] voxelPointIndexes, - int x, int y, int z, boolean contouring) { + void processOneCubical(int insideMask, float cutoff, int[] voxelPointIndexes, + int x, int y, int z, boolean contouring) { int edgeMask = insideMaskTable[insideMask]; + boolean isNaN = false; for (int iEdge = 12; --iEdge >= 0;) { if ((edgeMask & (1 << iEdge)) == 0) continue; @@ -1940,6 +1955,8 @@ int vertexB = edgeVertexes[2 * iEdge + 1]; float valueA = vertexValues[vertexA]; float valueB = vertexValues[vertexB]; + if (Float.isNaN(valueA) || Float.isNaN(valueB)) + isNaN = true; calcVertexPoints(x, y, z, vertexA, vertexB); float fraction = calcSurfacePoint(cutoff, valueA, valueB, surfacePoints[iEdge]); @@ -1967,7 +1984,7 @@ addVertexCopy(valueA < valueB ? pointB : pointA, false, ""); } } - if (contouring) + if (contouring || isNaN) return; byte[] triangles = triangleTable[insideMask]; for (int i = triangles.length; (i -= 3) >= 0;) { @@ -1999,7 +2016,7 @@ if (Float.isNaN(fraction) || fraction < 0 || fraction > 1) { System.out.println("problem with unusual fraction=" + fraction + " cutoff=" + cutoff + " A:" + valueA + " B:" + valueB); - fraction = 0; + fraction = Float.NaN; } if (!isJvxl) fractionData += jvxlFractionAsCharacter(fraction, edgeFractionBase, @@ -2676,7 +2693,8 @@ +"\n\t" + vertexValues2d[3] + " " + vertexValues2d[2] + "\n\t" + vertexValues2d[0] + " " + vertexValues2d[1]); */ - processOnePixel(insideMask, contourCutoff, pixelPointIndexes, x, y); + processOneQuadrilateral(insideMask, contourCutoff, pixelPointIndexes, + x, y); } } @@ -2716,8 +2734,8 @@ return pixelPointIndexes; } - void processOnePixel(int insideMask, float cutoff, int[] pixelPointIndexes, - int x, int y) { + void processOneQuadrilateral(int insideMask, float cutoff, + int[] pixelPointIndexes, int x, int y) { int edgeMask = insideMaskTable2d[insideMask]; planarSquares[x * squareCountY + y].addEdgeMask(contourIndex, edgeMask, insideMask); @@ -3000,7 +3018,7 @@ void getCalcPoint(Point3f pt) { pt.sub(center); if (isEccentric) - eccentricityMatrix.transform(pt); + eccentricityMatrixInverse.transform(pt); if (isAnisotropic) { pt.x /= anisotropy[0]; pt.y /= anisotropy[1]; @@ -3054,7 +3072,14 @@ case 2: volumetricVectors[2].set(0, 0, d); volumetricOrigin.z = min - d; + if (isEccentric) { + eccentricityMatrix.transform(volumetricOrigin); + } + volumetricOrigin.add(center); } + if (isEccentric) + eccentricityMatrix.transform(volumetricVectors[index]); + unitVolumetricVectors[index].normalize(volumetricVectors[index]); } @@ -3081,16 +3106,16 @@ /////// spheres and ellipsoids ////// - int sphere_gridMax = 30; + int sphere_gridMax = 20; float sphere_ptsPerAngstrom = 10f; float sphere_radiusAngstroms; void setupSphere() { if (center.x == Float.MAX_VALUE) center.set(0, 0, 0); + float radius = sphere_radiusAngstroms * 1.1f * eccentricityScale; for (int i = 0; i < 3; i++) - setVoxelRange(i, -sphere_radiusAngstroms - 0.5f, - sphere_radiusAngstroms + 0.5f, sphere_ptsPerAngstrom, sphere_gridMax); + setVoxelRange(i, -radius, radius, sphere_ptsPerAngstrom, sphere_gridMax); jvxlFileHeader = "SPHERE \nres=" + sphere_ptsPerAngstrom + " rad=" @@ -3236,20 +3261,21 @@ //////// lobes //////// + float lobe_sizeAngstroms = 1f; + int lobe_gridMax = 21; + int lobe_ptsPerAngstrom = 10; + void setupLobe() { - psi_n = 1; - psi_l = 0; + psi_n = 2; + psi_l = 1; psi_m = 0; psi_Znuc = 6; - psi_radiusAngstroms = autoScaleOrbital(); if (center.x == Float.MAX_VALUE) center.set(0, 0, 0); - setVoxelRange(0, -psi_radiusAngstroms, - psi_radiusAngstroms, psi_ptsPerAngstrom, psi_gridMax); - setVoxelRange(1, -psi_radiusAngstroms, - psi_radiusAngstroms, psi_ptsPerAngstrom, psi_gridMax); - setVoxelRange(2, 0, - psi_radiusAngstroms, 2 * psi_ptsPerAngstrom, psi_gridMax); + float radius = lobe_sizeAngstroms / 2f * eccentricityScale; + setVoxelRange(0, -radius, radius, lobe_ptsPerAngstrom, lobe_gridMax); + setVoxelRange(1, -radius, radius, lobe_ptsPerAngstrom, lobe_gridMax); + setVoxelRange(2, 0, radius * 4.4f, lobe_ptsPerAngstrom, lobe_gridMax); jvxlFileHeader = "lobe \nn=" + psi_n + ", l=" @@ -3259,9 +3285,9 @@ + " Znuc=" + psi_Znuc + " res=" - + psi_ptsPerAngstrom + + lobe_ptsPerAngstrom + " rad=" - + psi_radiusAngstroms + + radius + (isAnisotropic ? " anisotropy=(" + anisotropy[0] + "," + anisotropy[1] + "," + anisotropy[2] + ")" : "") + "\n"; jvxlFileHeader += jvxlGetVolumeHeader(); @@ -3270,14 +3296,12 @@ calcFactors(psi_n, psi_l, psi_m); } - float lobe_sizeAngstroms; - float getLobeValue(int x, int y, int z) { voxelPtToXYZ(x, y, z, ptPsi); getCalcPoint(ptPsi); - ptPsi.z = 1.5f * lobe_sizeAngstroms - ptPsi.z; + ptPsi.z = 1.2f - ptPsi.z; float value = (float) hydrogenAtomPsiAt(ptPsi, psi_n, psi_l, psi_m); - //System.out.println(x+ " " + y + " " + z + " " + value); + //System.out.println("getlobe " + x + " " + y + " " + z + ptPsi + value); return value; } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |