1
2
3
4
5
6
7
8 from __future__ import print_function
9 from rdkit import Geometry
10 from rdkit.Chem.Subshape import SubshapeObjects
11 import math
12 import numpy
13
14
16 if getattr(shapeGrid,'_indicesInSphere',None):
17 return shapeGrid._indicesInSphere
18 gridSpacing = shapeGrid.GetSpacing()
19 dX = shapeGrid.GetNumX()
20 dY = shapeGrid.GetNumY()
21 dZ = shapeGrid.GetNumZ()
22 radInGrid = int(winRad/gridSpacing)
23 indicesInSphere=[]
24 for i in range(-radInGrid,radInGrid+1):
25 for j in range(-radInGrid,radInGrid+1):
26 for k in range(-radInGrid,radInGrid+1):
27 d=int(math.sqrt(i*i + j*j + k*k ))
28 if d<=radInGrid:
29 idx = (i*dY+j)*dX+k
30 indicesInSphere.append(idx)
31 shapeGrid._indicesInSphere = indicesInSphere
32 return indicesInSphere
33
34
36 count=0
37 centroid = Geometry.Point3D(0,0,0)
38 idxI = shapeGrid.GetGridPointIndex(pt)
39 shapeGridVect = shapeGrid.GetOccupancyVect()
40
41 indicesInSphere = ComputeGridIndices(shapeGrid,winRad)
42
43 nGridPts = len(shapeGridVect)
44 for idxJ in indicesInSphere:
45 idx = idxI+idxJ;
46 if idx>=0 and idx<nGridPts:
47 wt = shapeGridVect[idx]
48 tPt = shapeGrid.GetGridPointLoc(idx)
49 centroid += tPt*wt
50 count+=wt
51 if not count:
52 raise ValueError('found no weight in sphere')
53 centroid /= count
54
55 return count,centroid
56
57
58
59
64
65
108
109
111 center = pt1+pt2
112 center /= 2.0
113 d=1e8
114 while d>shapeGrid.GetSpacing():
115 count,centroid=Geometry.ComputeGridCentroid(shapeGrid,center,winRad)
116 d = center.Distance(centroid)
117 center = centroid
118 return center
119
120
122 res = []
123 tagged = [(y,x) for x,y in enumerate(pts)]
124 while tagged:
125 head,headIdx = tagged.pop(0)
126 currSet = [head]
127 i=0
128 while i<len(tagged):
129 nbr,nbrIdx=tagged[i]
130 if head.location.Distance(nbr.location)<scale*winRad:
131 currSet.append(nbr)
132 del tagged[i]
133 else:
134 i+=1
135 pt = Geometry.Point3D(0,0,0)
136 for o in currSet:
137 pt += o.location
138 pt /= len(currSet)
139 res.append(SubshapeObjects.SkeletonPoint(location=pt))
140 return res
141
143 """ adds a set of new terminal points using a max-min algorithm
144 """
145 shapeGrid=shape.grid
146 shapeVect = shapeGrid.GetOccupancyVect()
147 nGridPts = len(shapeVect)
148
149
150 while len(pts)<targetNumber:
151 maxMin=-1
152 for i in range(nGridPts):
153 if shapeVect[i]<maxGridVal:
154 continue
155 minVal=1e8
156 posI = shapeGrid.GetGridPointLoc(i)
157 for currPt in pts:
158 dst = posI.Distance(currPt.location)
159 if dst<minVal:
160 minVal=dst
161 if minVal>maxMin:
162 maxMin=minVal
163 bestPt=posI
164 count,centroid=Geometry.ComputeGridCentroid(shapeGrid,bestPt,winRad)
165 pts.append(SubshapeObjects.SkeletonPoint(location=centroid))
166
167
169 """ find the grid point with max occupancy that is furthest from a
170 given location
171 """
172 shapeGrid=shape.grid
173 shapeVect = shapeGrid.GetOccupancyVect()
174 nGridPts = len(shapeVect)
175 dMax=-1;
176 for i in range(nGridPts):
177 if shapeVect[i]<maxGridVal:
178 continue
179 posI = shapeGrid.GetGridPointLoc(i)
180 dst = posI.Distance(loc)
181 if dst>dMax:
182 dMax=dst
183 res=posI
184
185 count,centroid=Geometry.ComputeGridCentroid(shapeGrid,res,winRad)
186 res=centroid
187 return res
188
190 """ find additional terminal points until a target number is reached
191 """
192 if len(pts)==1:
193
194
195 pt2=FindFarthestGridPoint(shape,pts[0].location,winRad,maxGridVal)
196 pts.append(SubshapeObjects.SkeletonPoint(location=pt2))
197 if len(pts)==2:
198
199 shapeGrid=shape.grid
200 pt1 = pts[0].location
201 pt2 = pts[1].location
202 center = FindGridPointBetweenPoints(pt1,pt2,shapeGrid,winRad)
203 pts.append(SubshapeObjects.SkeletonPoint(location=center))
204 if len(pts)<targetNumPts:
205 GetMoreTerminalPoints(shape,pts,winRad,maxGridVal,targetNumPts)
206
207
208
209 -def AppendSkeletonPoints(shapeGrid,termPts,winRad,stepDist,maxGridVal=3,
210 maxDistC=15.0,distTol=1.5,symFactor=1.5):
211 nTermPts = len(termPts)
212 skelPts=[]
213 shapeVect = shapeGrid.GetOccupancyVect()
214 nGridPts = len(shapeVect)
215
216 print('generate all possible')
217 for i in range(nGridPts):
218 if shapeVect[i]<maxGridVal:
219 continue
220 posI = shapeGrid.GetGridPointLoc(i)
221 ok=True
222 for pt in termPts:
223 dst = posI.Distance(pt.location)
224 if dst<stepDist:
225 ok=False
226 break
227 if ok:
228 skelPts.append(SubshapeObjects.SkeletonPoint(location=posI))
229
230 print('Compute centroids:',len(skelPts))
231 gridBoxVolume=shapeGrid.GetSpacing()**3
232 maxVol = 4.0*math.pi/3.0 * winRad**3 * maxGridVal / gridBoxVolume
233 i=0
234 while i<len(skelPts):
235 pt = skelPts[i]
236 count,centroid=Geometry.ComputeGridCentroid(shapeGrid,pt.location,winRad)
237
238 centroidPtDist=centroid.Distance(pt.location)
239 if centroidPtDist>maxDistC:
240 del skelPts[i]
241 else:
242 pt.fracVol = float(count)/maxVol
243 pt.location.x = centroid.x
244 pt.location.y = centroid.y
245 pt.location.z = centroid.z
246 i+=1
247
248 print('remove points:',len(skelPts))
249 res = termPts+skelPts
250 i=0
251 while i<len(res):
252 p=-1
253 mFrac=0.0
254 ptI = res[i]
255
256 startJ = max(i+1,nTermPts)
257 for j in range(startJ,len(res)):
258 ptJ=res[j]
259 distC = ptI.location.Distance(ptJ.location)
260 if distC<symFactor*stepDist:
261 if ptJ.fracVol>mFrac:
262 p=j
263 mFrac=ptJ.fracVol
264
265 if p>-1:
266 ptP = res.pop(p)
267 j = startJ
268 while j < len(res):
269 ptJ=res[j]
270 distC = ptI.location.Distance(ptJ.location)
271 if distC<symFactor*stepDist:
272 del res[j]
273 else:
274 j+=1
275 res.append(ptP)
276
277 i+=1
278 return res
279
280
282 shapeGridVect = shapeGrid.GetOccupancyVect()
283 nGridPts = len(shapeGridVect)
284 tmp = winRad/shapeGrid.GetSpacing()
285 radInGrid=int(tmp)
286 radInGrid2=int(tmp*tmp)
287 covMat = numpy.zeros((3,3),numpy.float64)
288
289 dX = shapeGrid.GetNumX()
290 dY = shapeGrid.GetNumY()
291 dZ = shapeGrid.GetNumZ()
292 idx = shapeGrid.GetGridPointIndex(pt.location)
293 idxZ = idx//(dX*dY)
294 rem = idx%(dX*dY)
295 idxY = rem//dX
296 idxX = rem%dX
297 totWt=0.0
298 for i in range(-radInGrid,radInGrid):
299 xi = idxX+i
300 for j in range(-radInGrid,radInGrid):
301 xj = idxY+j
302 for k in range(-radInGrid,radInGrid):
303 xk = idxZ+k
304 d2 = i*i+j*j+k*k
305 if d2>radInGrid2 and int(math.sqrt(d2))>radInGrid:
306 continue
307 gridIdx = (xk*dY+xj)*dX+xi
308 if gridIdx>=0 and gridIdx<nGridPts:
309 wtHere = shapeGridVect[gridIdx]
310 totWt += wtHere
311 ptInSphere = shapeGrid.GetGridPointLoc(gridIdx)
312 ptInSphere -= pt.location
313 covMat[0][0]+= wtHere*(ptInSphere.x*ptInSphere.x)
314 covMat[0][1]+= wtHere*(ptInSphere.x*ptInSphere.y)
315 covMat[0][2]+= wtHere*(ptInSphere.x*ptInSphere.z)
316 covMat[1][1]+= wtHere*(ptInSphere.y*ptInSphere.y)
317 covMat[1][2]+= wtHere*(ptInSphere.y*ptInSphere.z)
318 covMat[2][2]+= wtHere*(ptInSphere.z*ptInSphere.z)
319 covMat /= totWt
320 covMat[1][0] = covMat[0][1]
321 covMat[2][0] = covMat[0][2]
322 covMat[2][1] = covMat[1][2]
323
324 eVals,eVects = numpy.linalg.eigh(covMat)
325 sv = zip(eVals,numpy.transpose(eVects))
326 sv.sort(reverse=True)
327 eVals,eVects=zip(*sv)
328 pt.shapeMoments=tuple(eVals)
329 pt.shapeDirs = tuple([Geometry.Point3D(p[0],p[1],p[2]) for p in eVects])
330
331
332
333
334
335
336
337
338
339
340
341
342
352