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