Package Pmv :: Module computeIsovalue
[hide private]
[frames] | no frames]

Source Code for Module Pmv.computeIsovalue

  1  from geomutils.surface import findComponents, meshVolume, triangleArea 
  2  try: 
  3      from UTpackages.UTisocontour import isocontour 
  4  except: 
  5      pass 
  6  try: 
  7      from UTpackages.UTblur import blur 
  8  except: 
  9      pass 
 10   
 11  try: 
 12      from mslib import MSMS 
 13  except: 
 14      pass 
 15   
 16  try: 
 17     from QSlimLib import qslimlib 
 18  except: 
 19     pass 
 20   
 21  from math import fabs 
 22  from MolKit.molecule import Atom 
 23   
 24  import Numeric 
 25  from math import sqrt 
 26   
27 -def computeIsovalue(nodes, blobbyness, densityList=3.0, gridResolution=0.5, criteria="Volume", radiiSet="united", msmsProbeRadius=1.4, computeWeightOption=0, resolutionLevel=None, gridDims=None, padding=0.0):
28 """This function is based on shapefit/surfdock.py blurSurface() by Qing Zhang. 29 30 INPUT: 31 nodes 32 blobbyness: blobbyness 33 densityList: A list of vertex densities (or just a single density value) for blur surfaces 34 e.g. [1.0, 2.0, 3.0] or 1.0 or [1.0] 35 gridResolution: Resolution of grid used for blurring 36 0.5 gives surface vertex density > 5 dots/Angstrom^2 37 0.25 gives surface vertex density > 20 dots/Angstrom^2 38 0.2 gives surface vertex density > 40 dots/Angstrom^2 39 radiiSet: Radii set for computing both MSMS and blur surfaces 40 msmsProbeRadius: Probe radius for computing MSMS 41 resolutionLevel: Level of surface resolution 42 = None: using the given blobbyness value 43 = 1 or 'very low' blobbyness = -0.1 44 = 2 or 'low' blobbyness = -0.3 45 = 3 or 'medium' blobbyness = -0.5 46 = 4 or 'high' blobbyness = -0.9 47 = 5 or 'very high' blobbyness = -3.0 48 49 padding : size of the padding around the molecule 50 51 52 """ 53 # 0. Initialization 54 55 blurSurfList = [] 56 57 isoGridResolution = 0.5 # grid resolution for finding isovalue 58 isoDensity = 5.0 # density of blur surface vertices for finding isovalue 59 msmsDensity = 20.0 # density of MSMS surface vertices for compute MSMS volume 60 #msmsFilePrefix=molName # prefix of MSMS surface filenames 61 62 63 # ... 1.1. Blur 64 print "... 1.1. Blur ......" 65 66 atoms = nodes.findType(Atom) 67 coords = atoms.coords 68 radii = atoms.radius 69 arrayf, origin, stepsize = blurCoordsRadii(coords, radii, blobbyness, 70 isoGridResolution, gridDims, padding) 71 data = isocontour.newDatasetRegFloat3D(arrayf, origin, stepsize) 72 print "##### data : ", type(data), "#######" 73 # ... 1.2. Compute area and volume of MSMS surface at msmsDensity 74 print "... 1.2. Compute MSMS volume ......" 75 msmsVolume = computeMSMSvolume(coords, radii, msmsProbeRadius, msmsDensity) 76 #msmsArea, msmsVolume = computeMSMSAreaVolumeViaCommand(molFile, radiiSet, msmsProbeRadius, dens=msmsDensity, outFilePrefix=msmsFilePrefix) 77 #os.system( "rm -rf %s.vert %s.face %s.xyzr" % ( msmsFilePrefix, msmsFilePrefix, msmsFilePrefix ) ) # remove MSMS surface files 78 # ... 1.3. Find isovalue based on either MSMS area or volume at blur isoDensity 79 print "... 1.3. Find isovalue ......" 80 res = findIsoValueMol(radii, coords, blobbyness, data, isoDensity, msmsVolume, "Volume") 81 82 target, targetValue, value, isovalue, v1, t1, n1 = res 83 isocontour.clearDataset(data) # free memory of data 84 if isovalue == 0.0: # isovalue can not be found 85 print "blurSurface(): isovalue can not be found" 86 return None 87 return isovalue
88
89 -def blurCoordsRadii(coords, radii, blobbyness=-0.1, res=0.5, 90 weights=None, dims=None, padding = 0.0):
91 """blur a set of coordinates with radii 92 """ 93 # Setup grid 94 resX = resY = resZ = res 95 if not dims: 96 minb, maxb = blur.getBoundingBox(coords, radii, blobbyness, padding) 97 Xdim = int(round( (maxb[0] - minb[0])/resX + 1)) 98 Ydim = int(round( (maxb[1] - minb[1])/resY + 1)) 99 Zdim = int(round( (maxb[2] - minb[2])/resZ + 1)) 100 else: 101 Xdim, Ydim, Zdim = dims 102 print "Xdim = %d, Ydim =%d, Zdim = %d"%(Xdim, Ydim, Zdim) 103 # Generate blur map 104 volarr, origin, span = blur.generateBlurmap( 105 coords, radii, [Xdim, Ydim, Zdim], blobbyness, weights=weights, padding=padding) 106 # Take data from blur map 107 volarr.shape = [Zdim, Ydim, Xdim] 108 volarr = Numeric.transpose(volarr).astype('f') 109 origin = Numeric.array(origin).astype('f') 110 stepsize = Numeric.array(span).astype('f') 111 arrayf = Numeric.reshape( Numeric.transpose(volarr), 112 (1, 1)+tuple(volarr.shape) ) 113 # Return data 114 return arrayf, origin, stepsize
115
116 -def computeMSMSvolume(atmCoords, atmRadii, pRadius, dens ):
117 srf = MSMS(coords=atmCoords, radii=atmRadii) 118 srf.compute(probe_radius=pRadius, density=dens) 119 vf, vi, f = srf.getTriangles() 120 121 vertices=vf[:,:3] 122 normals=vf[:,3:6] 123 triangles=f[:,:3] 124 return meshVolume(vertices, normals, triangles)
125 126 127
128 -def meshArea(verts, tri):
129 """Compute the surface area of a surface as the sum of the area of its 130 triangles. The surface is specified by vertices, indices of triangular faces 131 """ 132 areaSum = 0.0 133 for t in tri: 134 s1 = verts[t[0]] 135 s2 = verts[t[1]] 136 s3 = verts[t[2]] 137 area = triangleArea(s1,s2,s3) 138 areaSum += area 139 140 return areaSum
141
142 -def normalCorrection(norms, faces):
143 """Replace any normal, which is equal to zero, with the average 144 of its triangle neighbors' non-zero normals (triangle neighbors, as 145 the edge partners, can be derived from the faces) 146 """ 147 # Convert one-based faces temporarily to zero-based 148 newFaces = [] 149 if faces[0][0] == 1: 150 for v1, v2, v3 in faces: # vertex indices v1, v2, v3 151 newFaces.append( ( v1-1, v2-1, v3-1 ) ) 152 else: 153 newFaces = faces 154 155 # Build a neighborhood dictionary for easy neighbor search 156 neighborDict = {} 157 for v1, v2, v3 in newFaces: # vertex indices v1, v2, v3 158 # v1 159 if v1 not in neighborDict: 160 neighborDict[v1] = [] 161 if v2 not in neighborDict[v1]: 162 neighborDict[v1].append(v2) 163 if v3 not in neighborDict[v1]: 164 neighborDict[v1].append(v3) 165 # v2 166 if v2 not in neighborDict: 167 neighborDict[v2] = [] 168 if v3 not in neighborDict[v2]: 169 neighborDict[v2].append(v3) 170 if v1 not in neighborDict[v2]: 171 neighborDict[v2].append(v1) 172 # v3 173 if v3 not in neighborDict: 174 neighborDict[v3] = [] 175 if v1 not in neighborDict[v3]: 176 neighborDict[v3].append(v1) 177 if v2 not in neighborDict[v3]: 178 neighborDict[v3].append(v2) 179 180 # Find any zero-value normal and replace it 181 newNorms = [] 182 for i, eachnorm in enumerate(norms): 183 # non-zero normal 184 if not ( eachnorm[0]==0.0 and eachnorm[1]==0.0 and eachnorm[2]==0.0 ): 185 newNorms.append( [ eachnorm[0], eachnorm[1], eachnorm[2] ] ) 186 continue 187 # zero normal 188 print "normalCorrection(): zero normal at vertex %d", i 189 neighbors = neighborDict[i] 190 neighborNorms = [0.0, 0.0, 0.0] 191 Nneighbors = 0 192 for eachneighbor in neighbors: # sum up non-zero neighbors 193 eachneighborNorm = norms[eachneighbor] 194 if eachneighborNorm[0]==0.0 and eachneighborNorm[1]==0.0 and eachneighborNorm[2]==0.0: 195 continue # skip zero-normal neighbor 196 Nneighbors += 1 197 neighborNorms[0] += eachneighborNorm[0] 198 neighborNorms[1] += eachneighborNorm[1] 199 neighborNorms[2] += eachneighborNorm[2] 200 if Nneighbors == 0: # non-zero neighbors can not be found 201 print "Can not find non-zero neighbor normals for normal %d " % i 202 newNorms.append(neighborNorms) 203 continue 204 neighborNorms[0] /= Nneighbors # average 205 neighborNorms[1] /= Nneighbors 206 neighborNorms[2] /= Nneighbors 207 newNorms.append(neighborNorms) 208 209 # Return 210 return newNorms
211
212 -def findIsoValueMol(radii, coords, blobbyness, data, isoDensity, targetValue, target):
213 """Find the best isocontour value (isovalue) for a molecule 214 at one specific blobbyness by reproducing the targetValue 215 for a target (area or volume) 216 INPUT: 217 mol: molecule recognizable by MolKit 218 blobbyness: blobbyness 219 data: blurred data 220 isoDensity: density of surface vertices for finding isovalue only 221 targetValue: value of the target, to be reproduced 222 target: Area or Volume 223 OUTPUT: 224 target: (as input) 225 targetValue: (as input) 226 value: reproduced value of target 227 isovalue: the best isocontour value 228 v1d: vertices of generated blur surface 229 t1d: faces of generated blur surface 230 n1d: normals of the vertices 231 """ 232 233 234 # Initialization 235 isovalue = 1.0 # guess starting value 236 mini = 0. 237 maxi = 99. 238 cutoff = targetValue*0.01 # 99% reproduction of targetValue 239 accuracy = 0.0001 # accuracy of isovalue 240 stepsize = 0.5 # step size of increase/decrease of isovalue 241 value_adjust = 1.0 # adjustment of value 242 value = 0.0 243 244 # Optimize isovalue 245 while True: 246 print "------- isovalue = %f -------" % isovalue 247 248 # 1. Isocontour 249 print "... ... 1. Isocontour ......" 250 v1, t1, n1 = computeIsocontour(isovalue, data) 251 if len(v1)==0 or len(t1)==0: 252 value = 0.0 253 maxi = isovalue 254 isovalue = maxi - (maxi-mini)*stepsize 255 print "***** isocontoured surface has no vertices *****", len(v1), ", ", len(t1) 256 continue 257 258 # 2. Remove small components (before decimate) 259 print "... ... 2. Remove small components ......" 260 v1, t1, n1 = findComponents(v1, t1, n1, 1) 261 262 # 3. Decimate 263 print "... ... 3. Decimate ......" 264 v1d, t1d, n1d = decimate(v1, t1, n1, isoDensity) 265 if len(v1d)==0 or len(t1d)==0: 266 value = 0.0 267 maxi = isovalue 268 isovalue = maxi - (maxi-mini)*stepsize 269 continue 270 #################### Analytical normals for decimated vertices ##################### 271 normals = computeBlurCurvature(radii, coords, blobbyness, v1d) # Analytical Blur curvatures 272 n1d = normals.tolist() 273 #################################################################################### 274 275 # 4. Compute value 276 print "... ... 4. Compute value ......" 277 if target == "Area": 278 value = meshArea(v1d, t1d) 279 else: 280 n1d = normalCorrection(n1d, t1d) # replace zero normals 281 print "calling meshVolume after normalCorrection" 282 value = meshVolume(v1d, n1d, t1d) 283 284 # 5. Adjust value 285 print "... ... 5. Adjust value ......" 286 value *= value_adjust 287 288 # TEST 289 print "target=%6s, targetValue=%10.3f, value=%10.3f, isovalue=%7.3f, maxi=%7.3f, mini=%7.3f"%( 290 target, targetValue, value, isovalue, maxi, mini) 291 292 # 6. Evaluate isocontour value 293 print "... ... 6. Evaluate isovalue ......" 294 if fabs(value - targetValue) < cutoff: # Satisfy condition! 295 print "Found: target=%6s, targetValue=%10.3f, value=%10.3f, isovalue=%7.3f, maxi=%7.3f, mini=%7.3f"%( 296 target, targetValue, value, isovalue, maxi, mini) 297 break 298 299 if maxi-mini < accuracy: # can not find good isovalue 300 print "Not found: target=%6s, targetValue=%10.3f, value=%10.3f, isovalue=%7.3f, maxi=%7.3f, mini=%7.3f"%( 301 target, targetValue, value, isovalue, maxi, mini) 302 return target, targetValue, value, 0.0, v1, t1, n1 303 304 if value > targetValue: # value too big, increase isovalue 305 mini = isovalue 306 isovalue = mini + (maxi-mini)*stepsize 307 else: # value too small, decrease isovalue 308 maxi = isovalue 309 isovalue = maxi - (maxi-mini)*stepsize 310 311 # Return 312 return target, targetValue, value, isovalue, v1, t1, n1
313
314 -def computeIsocontour(isovalue, data):
315 isoc = isocontour.getContour3d(data, 0, 0, isovalue, 316 isocontour.NO_COLOR_VARIABLE) 317 318 if isoc.nvert==0 or isoc.ntri==0: 319 return [], [], [] 320 vert = Numeric.zeros((isoc.nvert,3)).astype('f') 321 norm = Numeric.zeros((isoc.nvert,3)).astype('f') 322 col = Numeric.zeros((isoc.nvert)).astype('f') 323 tri = Numeric.zeros((isoc.ntri,3)).astype('i') 324 isocontour.getContour3dData(isoc, vert, norm, col, tri, 0) 325 if len(vert) == 0 or len(tri) == 0: 326 return [], [], [] 327 return vert, tri, norm
328
329 -def decimate(vert, tri, normals, density=1.0):
330 # decimate a mesh to vertex density close to 1.0 331 if len(vert)==0 or len(tri)==0: 332 return [], [], [] 333 334 #print "before decimate: ", len(vert), " vertices, ", len(tri), " faces" 335 model = qslimlib.QSlimModel(vert, tri, bindtoface=False, 336 colors=None, norms=normals) 337 nverts = model.get_vert_count() 338 nfaces = model.get_face_count() 339 # allocate array to read decimated surface back 340 341 newverts = Numeric.zeros((nverts, 3)).astype('f') 342 newfaces = Numeric.zeros((nfaces, 3)).astype('i') 343 newnorms = Numeric.zeros((nverts, 3)).astype('f') 344 345 # print len(vert), area, len(tri), int(len(tri)*len(vert)/area) 346 area = meshArea(vert, tri) 347 tface = int(min(len(tri), int(len(tri)*area/len(vert)))*density) 348 349 model.slim_to_target(tface) 350 numFaces = model.num_valid_faces() 351 352 # print 'decimated to %d triangles from (%d)'%(numFaces, len(tri)) 353 model.outmodel(newverts, newfaces, outnorms=newnorms) 354 355 numVertices = max(newfaces.flat)+1 356 # build lookup table of vertices used in faces 357 d = {} 358 for t in newfaces[:numFaces]: 359 d[t[0]] = 1 360 d[t[1]] = 1 361 d[t[2]] = 1 362 vl = {} 363 decimVerts = Numeric.zeros((numVertices, 3)).astype('f') 364 decimFaces = Numeric.zeros((numFaces, 3)).astype('i') 365 decimNormals = Numeric.zeros((numVertices, 3)).astype('f') 366 367 nvert = 0 368 nvert = 0 369 for j, t in enumerate(newfaces[:numFaces]): 370 for i in (0,1,2): 371 n = t[i] 372 if not vl.has_key(n): 373 vl[n] = nvert 374 decimVerts[nvert] = newverts[n,:] 375 decimNormals[nvert] = newnorms[n,:] 376 nvert += 1 377 decimFaces[j] = (vl[t[0]], vl[t[1]], vl[t[2]]) 378 ## vertices=newverts[:numVertices] 379 ## for t in newfaces[:numFaces]: 380 ## for i in (0,1,2): 381 ## if not vl.has_key(t[i]): 382 ## vl[t[i]] = nvert 383 ## decimVerts.append(vertices[t[i]]) 384 ## nvert += 1 385 ## decimFaces.append( (vl[t[0]], vl[t[1]], vl[t[2]] ) ) 386 387 #print 'density after decimation', len(decimVerts)/meshArea(decimVerts, decimFaces) 388 ## norms1 = glTriangleNormals( decimVerts, decimFaces, 'PER_VERTEX') 389 #print "after decimate: ", nvert, " vertices, ", numFaces, " faces" 390 return decimVerts[:nvert], decimFaces, decimNormals[:nvert]
391
392 -def computeBlurCurvature(radii, coords, blobbiness, points):
393 # Compute Blur curvatures analytically by UTmolderivatives (since Novemeber 2005) 394 from UTpackages.UTmolderivatives import molderivatives 395 import Numeric 396 # Read input 397 npoints = len(points) 398 numberOfGaussians = len(coords) 399 gaussianCenters = Numeric.zeros((numberOfGaussians,4)).astype('d') 400 for i in range(numberOfGaussians): 401 gaussianCenters[i,0:3] = coords[i] 402 gaussianCenters[i,3] = radii[i] 403 numberOfGridDivisions = 10 # as in Readme.htm of UTmolderivatives 404 maxFunctionError = 0.001 # as in Readme.htm of UTmolderivatives 405 406 # Compute HandK, normals, k1Vector, k2Vector 407 HandK, normals, k1Vector, k2Vector = molderivatives.getGaussianCurvature(gaussianCenters, 408 numberOfGridDivisions, maxFunctionError, 409 blobbiness, points) 410 411 k1Vector = Numeric.reshape(k1Vector, (npoints, 3)) 412 k2Vector = Numeric.reshape(k2Vector, (npoints, 3)) 413 HandK = Numeric.reshape(HandK, (npoints, 2)) 414 normals = Numeric.reshape(normals, (npoints, 3)) * (-1) # change normal directions to outwards 415 return normals
416