LOOS  v2.3.2
ConvexHull.py
1 #!/usr/bin/env python
2 
3 import loos
4 import numpy
5 
6 try:
7  from scipy.spatial import ConvexHull as CH
8 except ImportError:
9  import sys
10  print """
11  Failure importing scipy, needed for ConvexHull calculations
12  You either need to install scipy or adjust your PYTHONPATH
13  to point to it.
14  """
15  sys.exit(1)
16 
18  def __init__(self, min_z, max_z):
19  self.min_z = float(min_z)
20  self.max_z = float(max_z)
21 
22  def __call__(self, atomicgroup):
23  ag = loos.AtomicGroup()
24  ag.periodicBox(atomicgroup.periodicBox())
25  for a in atomicgroup:
26  if self.min_z < a.coords().z() < self.max_z:
27  ag.append(a)
28  return ag
29 
30 def test_side(p, p0, p1):
31  """
32  Return val < 0 if p is to the left of the p0-p1 line segment,
33  val > 0 if p is to the right
34  val = 0 if p is on the line
35  """
36  val = (p.y()-p0.y())*(p1.x()-p0.x()) - (p.x()-p0.x())*(p1.y()-p0.y())
37  return val
38 
39 class ConvexHull:
40  def __init__(self, atomicgroup):
41  self.atoms = atomicgroup
42  self.hull = None
43  self.vertices = None
44  self.simplices = None
45 
46  def num_atoms(self):
47  return len(self.atoms)
48 
49  def update_atoms(self, atomicgroup):
50  self.atoms = atomicgroup
51 
52  def generate_hull(self):
53  points = []
54  for a in self.atoms:
55  points.append([a.coords().x(), a.coords().y()])
56 
57  points = numpy.array(points)
58  self.hull = CH(points)
59 
60  def generate_vertices(self):
61  n = self.hull.simplices
62  sorted_neighbors = []
63  sorted_neighbors.append(n[0][0])
64  sorted_neighbors.append(n[0][1])
65  prev_entered = n[0][0]
66  last_entered = n[0][1]
67  num_simplices = len(self.hull.simplices)
68  while (len(sorted_neighbors) < num_simplices):
69  # search neighbor list to find the one with the current simplex
70  # add its pair to the list
71  last_entered = sorted_neighbors[-1]
72  prev_entered = sorted_neighbors[-2]
73  for i in range(len(n)):
74  if ((n[i][0] == last_entered) and (n[i][1] != prev_entered)):
75  sorted_neighbors.append(n[i][1])
76  break
77  if ((n[i][1] == last_entered) and (n[i][0] != prev_entered)):
78  sorted_neighbors.append(n[i][0])
79  break
80  self.vertices = sorted_neighbors
81 
82 
83  def atom(self, index):
84  return self.atoms[int(index)]
85 
86  def coords(self, index):
87  return self.atoms[int(index)].coords()
88 
89  def is_inside(self, p): # p is a Gcoord
90  match = True
91  side = test_side(p, self.coords(self.vertices[0]),
92  self.coords(self.vertices[1]))
93  prev_side = side
94  for i in range(1, len(self.vertices)-1):
95  side = test_side(p, self.coords(self.vertices[i]),
96  self.coords(self.vertices[i+1]))
97  match = ( ((side >= 0) and (prev_side >=0) ) or
98  ((side <= 0) and (prev_side <=0) ) )
99  if not match: return False
100  prev_side = side
101 
102  # test the closing line segment
103  side = test_side(p, self.coords(self.vertices[-1]),
104  self.coords(self.vertices[0]))
105  match = ( ((side >= 0) and (prev_side >=0) ) or
106  ((side <= 0) and (prev_side <=0) ) )
107  return match
108 
109 
110 
111 if __name__ == '__main__':
112 
113  from numpy import random
114  ag = loos.createSystem("rhod_only.pdb")
115 
116  slicer = ZSliceSelector(-3.0,3.0)
117 
118  ag = loos.selectAtoms(ag, 'name == "CA"')
119  ag_slice = slicer(ag)
120 
121  hull = ConvexHull(ag_slice)
122  hull.generate_hull()
123  #print hull.hull.simplices
124  #print hull.hull.neighbors
125  hull.generate_vertices()
126 
127  #print hull.vertices
128  #print hull.num_atoms()
129  i = 0
130  for v in hull.vertices:
131  c = hull.coords(v)
132  print i, v, c.x(), c.y(), c.z()
133  i += 1
134 
135 
136  print hull.is_inside(loos.GCoord(0.0,0.0,0.0))
137  print hull.is_inside(loos.GCoord(20.0,0.0,0.0))
138  print hull.is_inside(loos.GCoord(0.0,20.0,0.0))
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
def coords(self, index)
Definition: ConvexHull.py:86
AtomicGroup createSystem(const std::string &filename, const std::string &filetype)
Definition: sfactories.cpp:119