OpenStructure
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
spatial_organizer.hh
Go to the documentation of this file.
1 //------------------------------------------------------------------------------
2 // This file is part of the OpenStructure project <www.openstructure.org>
3 //
4 // Copyright (C) 2008-2011 by the OpenStructure authors
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License as published by the Free
8 // Software Foundation; either version 3.0 of the License, or (at your option)
9 // any later version.
10 // This library is distributed in the hope that it will be useful, but WITHOUT
11 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU Lesser General Public License
16 // along with this library; if not, write to the Free Software Foundation, Inc.,
17 // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 //------------------------------------------------------------------------------
19 #ifndef OST_SPATIAL_ORGANIZER_HI
20 #define OST_SPATIAL_ORGANIZER_HI
21 
22 #include <vector>
23 #include <map>
24 #include <cmath>
25 #include <limits>
26 
27 #include <ost/geom/geom.hh>
28 
29 namespace ost { namespace mol {
30 
32 /*
33  organizes ITEMs defined by positions
34  in a spatial hash map, into bins of
35  a defined size (as given in the ctor)
36 */
37 template <class ITEM, class VEC=geom::Vec3>
39 public:
40  typedef std::vector<ITEM> ItemList;
41 
42 private:
43  struct Index {
44  Index():
45  u(0),v(0),w(0) {}
46  Index(int uu, int vv, int ww):
47  u(uu),v(vv),w(ww) {}
48 
49  bool operator<(const Index& other) const {
50  return w!=other.w ? w<other.w : (v!=other.v ? v<other.v : u<other.u);
51  }
52 
53  int u,v,w;
54 
55  static Index Max(const Index& a, const Index& b) {
56  Index m;
57  m.u = a.u > b.u ? a.u : b.u;
58  m.v = a.v > b.v ? a.v : b.v;
59  m.w = a.w > b.w ? a.w : b.w;
60  return m;
61  }
62 
63  static Index Min(const Index& a, const Index& b) {
64  Index m;
65  m.u = a.u < b.u ? a.u : b.u;
66  m.v = a.v < b.v ? a.v : b.v;
67  m.w = a.w < b.w ? a.w : b.w;
68  return m;
69  }
70 
71  };
72 
73  struct Entry {
74  Entry(const ITEM& i, const VEC& p): item(i), pos(p) {}
75 
76  ITEM item;
77  VEC pos;
78  };
79 
80  typedef std::vector<Entry> EntryList;
81 
82  typedef std::map<Index,EntryList> ItemMap;
83 
84 public:
85 
86  SpatialOrganizer(Real delta): delta_(delta) {
87  if(delta==0.0) {
88  throw "delta cannot be zero";
89  }
90  }
91 
92  void Add(const ITEM& item, const VEC& pos) {
93  bool first = map_.empty();
94  Index indx=gen_index(pos);
95  map_[indx].push_back(Entry(item,pos));
96  if (!first) {
97  min_ = Index::Min(min_, indx);
98  max_ = Index::Max(max_, indx);
99  } else {
100  min_ = indx;
101  max_ = indx;
102  }
103  }
104  void Remove(const ITEM& item) {
105  typename ItemMap::iterator i=map_.begin();
106  for (; i!=map_.end(); ++i) {
107  for (size_t j=0; j<i->second.size(); ++j) {
108  if (i->second[j].item==item) {
109  i->second.erase(i->second.begin()+j);
110  return;
111  }
112  }
113  }
114  }
115 
116  bool HasWithin(const VEC& pos, Real dist) const {
117  Real dist2=dist*dist;
118  Index imin = Index::Max(min_, gen_index(pos-VEC(dist,dist,dist)));
119  Index imax = Index::Min(max_, gen_index(pos+VEC(dist,dist,dist)));
120  const size_t tmp = (imax.u-imin.u+1)*(imax.v-imin.v+1)*(imax.w-imin.w+1);
121  if (tmp > map_.size()) {
122  return this->has_within_all_buckets(pos, dist2);
123  }
124  for(int wc=imin.w;wc<=imax.w;++wc) {
125  for(int vc=imin.v;vc<=imax.v;++vc) {
126  for(int uc=imin.u;uc<=imax.u;++uc) {
127  typename ItemMap::const_iterator map_it = map_.find(Index(uc,vc,wc));
128  if(map_it!=map_.end()) {
129  for(typename EntryList::const_iterator entry_it = map_it->second.begin();
130  entry_it != map_it->second.end(); ++entry_it) {
131  /*
132  speed tests indicate that pre-calculating dx2 or dy2
133  and pre-checking them with an additional if gives little
134  speed improvement for very specific circumstances only,
135  but most of the time the performance is worse.
136  */
137  Real delta_x = entry_it->pos[0]-pos[0];
138  Real delta_y = entry_it->pos[1]-pos[1];
139  Real delta_z = entry_it->pos[2]-pos[2];
140  if(delta_x*delta_x+delta_y*delta_y+delta_z*delta_z<=dist2) {
141  return true;
142  }
143  }
144  }
145  }
146  }
147  }
148  return false;
149  }
150 
151  ItemList FindWithin(const VEC& pos, Real dist) const {
152  Real dist2=dist*dist;
153  Index imin = gen_index(pos-VEC(dist,dist,dist));
154  Index imax = gen_index(pos+VEC(dist,dist,dist));
155 
156  ItemList item_list;
157 
158  for(int wc=imin.w;wc<=imax.w;++wc) {
159  for(int vc=imin.v;vc<=imax.v;++vc) {
160  for(int uc=imin.u;uc<=imax.u;++uc) {
161  typename ItemMap::const_iterator map_it = map_.find(Index(uc,vc,wc));
162 
163  if(map_it!=map_.end()) {
164  for(typename EntryList::const_iterator entry_it = map_it->second.begin();
165  entry_it != map_it->second.end(); ++entry_it) {
166  /*
167  speed tests indicate that pre-calculating dx2 or dy2
168  and pre-checking them with an additional if gives little
169  speed improvement for very specific circumstances only,
170  but most of the time the performance is worse.
171  */
172  Real delta_x = entry_it->pos[0]-pos[0];
173  Real delta_y = entry_it->pos[1]-pos[1];
174  Real delta_z = entry_it->pos[2]-pos[2];
175  if(delta_x*delta_x+delta_y*delta_y+delta_z*delta_z<=dist2) {
176  item_list.push_back(entry_it->item);
177  }
178  }
179  }
180  }
181  }
182  }
183  return item_list;
184  }
185  void Clear()
186  {
187  map_.clear();
188  }
190  map_.swap(o.map_);
191  std::swap(delta_,o.delta_);
192  std::swap(min_, o.min_);
193  std::swap(max_, o.max_);
194  }
195 
196 private:
197  bool has_within_all_buckets(const VEC& pos, Real dist2) const {
198  for (typename ItemMap::const_iterator
199  i = map_.begin(), e = map_.end(); i!=e; ++i) {
200  for(typename EntryList::const_iterator j = i->second.begin();
201  j != i->second.end(); ++j) {
202  Real delta_x = j->pos[0]-pos[0];
203  Real delta_y = j->pos[1]-pos[1];
204  Real delta_z = j->pos[2]-pos[2];
205  if(delta_x*delta_x+delta_y*delta_y+delta_z*delta_z<=dist2) {
206  return true;
207  }
208  }
209  }
210  return false;
211  }
212  ItemMap map_;
213  Real delta_;
214  Index min_;
215  Index max_;
216 
217  Index gen_index(const VEC& pos) const {
218  Index nrvo(static_cast<int>(round(pos[0]/delta_)),
219  static_cast<int>(round(pos[1]/delta_)),
220  static_cast<int>(round(pos[2]/delta_)));
221  return nrvo;
222  }
223 
224  geom::Vec3 gen_middle(const Index& i) const {
225  geom::Vec3 nrvo((static_cast<Real>(i.u)+0.5)*delta_,
226  (static_cast<Real>(i.v)+0.5)*delta_,
227  (static_cast<Real>(i.w)+0.5)*delta_);
228  return nrvo;
229  }
230 
231 };
232 
233 }} // ns
234 
235 #endif
std::vector< ITEM > ItemList
bool operator<(const ResidueHandle &lhs, const ResidueHandle &rhs)
float Real
Definition: base.hh:44
void Add(const ITEM &item, const VEC &pos)
ItemList FindWithin(const VEC &pos, Real dist) const
void Swap(SpatialOrganizer &o)
Vec2 Max(const Vec2 &v1, const Vec2 &v2)
Definition: vecmat2_op.hh:150
Three dimensional vector class, using Real precision.
Definition: vec3.hh:42
#define DLLEXPORT_OST_MOL
bool HasWithin(const VEC &pos, Real dist) const
Vec2 Min(const Vec2 &v1, const Vec2 &v2)
Definition: vecmat2_op.hh:143
void Remove(const ITEM &item)