1
2
3
4
5 from __future__ import print_function
6
7 import math
8
9 import numpy
10
11 from rdkit import Geometry
12 from rdkit.Chem.Subshape import SubshapeObjects
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 radInGrid = int(winRad / gridSpacing)
22 indicesInSphere = []
23 for i in range(-radInGrid, radInGrid + 1):
24 for j in range(-radInGrid, radInGrid + 1):
25 for k in range(-radInGrid, radInGrid + 1):
26 d = int(math.sqrt(i * i + j * j + k * k))
27 if d <= radInGrid:
28 idx = (i * dY + j) * dX + k
29 indicesInSphere.append(idx)
30 shapeGrid._indicesInSphere = indicesInSphere
31 return indicesInSphere
32
33
35 count = 0
36 centroid = Geometry.Point3D(0, 0, 0)
37 idxI = shapeGrid.GetGridPointIndex(pt)
38 shapeGridVect = shapeGrid.GetOccupancyVect()
39
40 indicesInSphere = ComputeGridIndices(shapeGrid, winRad)
41
42 nGridPts = len(shapeGridVect)
43 for idxJ in indicesInSphere:
44 idx = idxI + idxJ
45 if idx >= 0 and idx < nGridPts:
46 wt = shapeGridVect[idx]
47 tPt = shapeGrid.GetGridPointLoc(idx)
48 centroid += tPt * wt
49 count += wt
50 if not count:
51 raise ValueError('found no weight in sphere')
52 centroid /= count
53
54 return count, centroid
55
56
61
62
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
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
190
192 """ find additional terminal points until a target number is reached
193 """
194 if len(pts) == 1:
195
196
197 pt2 = FindFarthestGridPoint(shape, pts[0].location, winRad, maxGridVal)
198 pts.append(SubshapeObjects.SkeletonPoint(location=pt2))
199 if len(pts) == 2:
200
201 shapeGrid = shape.grid
202 pt1 = pts[0].location
203 pt2 = pts[1].location
204 center = FindGridPointBetweenPoints(pt1, pt2, shapeGrid, winRad)
205 pts.append(SubshapeObjects.SkeletonPoint(location=center))
206 if len(pts) < targetNumPts:
207 GetMoreTerminalPoints(shape, pts, winRad, maxGridVal, targetNumPts)
208
209
210 -def AppendSkeletonPoints(shapeGrid, termPts, winRad, stepDist, maxGridVal=3, maxDistC=15.0,
211 distTol=1.5, symFactor=1.5, verbose=False):
212 nTermPts = len(termPts)
213 skelPts = []
214 shapeVect = shapeGrid.GetOccupancyVect()
215 nGridPts = len(shapeVect)
216
217 if verbose:
218 print('generate all possible')
219 for i in range(nGridPts):
220 if shapeVect[i] < maxGridVal:
221 continue
222 posI = shapeGrid.GetGridPointLoc(i)
223 ok = True
224 for pt in termPts:
225 dst = posI.Distance(pt.location)
226 if dst < stepDist:
227 ok = False
228 break
229 if ok:
230 skelPts.append(SubshapeObjects.SkeletonPoint(location=posI))
231
232 if verbose:
233 print('Compute centroids:', len(skelPts))
234 gridBoxVolume = shapeGrid.GetSpacing()**3
235 maxVol = 4.0 * math.pi / 3.0 * winRad**3 * maxGridVal / gridBoxVolume
236 i = 0
237 while i < len(skelPts):
238 pt = skelPts[i]
239 count, centroid = Geometry.ComputeGridCentroid(shapeGrid, pt.location, winRad)
240
241 centroidPtDist = centroid.Distance(pt.location)
242 if centroidPtDist > maxDistC:
243 del skelPts[i]
244 else:
245 pt.fracVol = float(count) / maxVol
246 pt.location.x = centroid.x
247 pt.location.y = centroid.y
248 pt.location.z = centroid.z
249 i += 1
250
251 if verbose:
252 print('remove points:', len(skelPts))
253 res = termPts + skelPts
254 i = 0
255 while i < len(res):
256 p = -1
257 mFrac = 0.0
258 ptI = res[i]
259
260 startJ = max(i + 1, nTermPts)
261 for j in range(startJ, len(res)):
262 ptJ = res[j]
263 distC = ptI.location.Distance(ptJ.location)
264 if distC < symFactor * stepDist:
265 if ptJ.fracVol > mFrac:
266 p = j
267 mFrac = ptJ.fracVol
268
269 if p > -1:
270 ptP = res.pop(p)
271 j = startJ
272 while j < len(res):
273 ptJ = res[j]
274 distC = ptI.location.Distance(ptJ.location)
275 if distC < symFactor * stepDist:
276 del res[j]
277 else:
278 j += 1
279 res.append(ptP)
280
281 i += 1
282 return res
283
284
286 shapeGridVect = shapeGrid.GetOccupancyVect()
287 nGridPts = len(shapeGridVect)
288 tmp = winRad / shapeGrid.GetSpacing()
289 radInGrid = int(tmp)
290 radInGrid2 = int(tmp * tmp)
291 covMat = numpy.zeros((3, 3), numpy.float64)
292
293 dX = shapeGrid.GetNumX()
294 dY = shapeGrid.GetNumY()
295
296 idx = shapeGrid.GetGridPointIndex(pt.location)
297 idxZ = idx // (dX * dY)
298 rem = idx % (dX * dY)
299 idxY = rem // dX
300 idxX = rem % dX
301 totWt = 0.0
302 for i in range(-radInGrid, radInGrid):
303 xi = idxX + i
304 for j in range(-radInGrid, radInGrid):
305 xj = idxY + j
306 for k in range(-radInGrid, radInGrid):
307 xk = idxZ + k
308 d2 = i * i + j * j + k * k
309 if d2 > radInGrid2 and int(math.sqrt(d2)) > radInGrid:
310 continue
311 gridIdx = (xk * dY + xj) * dX + xi
312 if gridIdx >= 0 and gridIdx < nGridPts:
313 wtHere = shapeGridVect[gridIdx]
314 totWt += wtHere
315 ptInSphere = shapeGrid.GetGridPointLoc(gridIdx)
316 ptInSphere -= pt.location
317 covMat[0][0] += wtHere * (ptInSphere.x * ptInSphere.x)
318 covMat[0][1] += wtHere * (ptInSphere.x * ptInSphere.y)
319 covMat[0][2] += wtHere * (ptInSphere.x * ptInSphere.z)
320 covMat[1][1] += wtHere * (ptInSphere.y * ptInSphere.y)
321 covMat[1][2] += wtHere * (ptInSphere.y * ptInSphere.z)
322 covMat[2][2] += wtHere * (ptInSphere.z * ptInSphere.z)
323 covMat /= totWt
324 covMat[1][0] = covMat[0][1]
325 covMat[2][0] = covMat[0][2]
326 covMat[2][1] = covMat[1][2]
327
328 eVals, eVects = numpy.linalg.eigh(covMat)
329 sv = list(zip(eVals, numpy.transpose(eVects)))
330 sv.sort(reverse=True)
331 eVals, eVects = list(zip(*sv))
332 pt.shapeMoments = tuple(eVals)
333 pt.shapeDirs = tuple([Geometry.Point3D(p[0], p[1], p[2]) for p in eVects])
334
335
336
337
338
339
340
341
342
343
344
345
355