OpenStructure
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 
56  struct Entry {
57  Entry(const ITEM& i, const VEC& p): item(i), pos(p) {}
58 
59  ITEM item;
60  VEC pos;
61  };
62 
63  typedef std::vector<Entry> EntryList;
64 
65  typedef std::map<Index,EntryList> ItemMap;
66 
67 public:
68 
69  SpatialOrganizer(Real delta): delta_(delta) {
70  if(delta==0.0) {
71  throw "delta cannot be zero";
72  }
73  }
74 
75  void Add(const ITEM& item, const VEC& pos) {
76  Index indx=gen_index(pos);
77  map_[indx].push_back(Entry(item,pos));
78  }
79 
80  void Remove(const ITEM& item) {
81  typename ItemMap::iterator i=map_.begin();
82  for (; i!=map_.end(); ++i) {
83  for (size_t j=0; j<i->second.size(); ++j) {
84  if (i->second[j].item==item) {
85  i->second.erase(i->second.begin()+j);
86  return;
87  }
88  }
89  }
90  }
91 
92  ITEM FindClosest(const VEC& pos) {
93 
94  // find closest index with non-empty itemlist
95  Real best_dist2=std::numeric_limits<Real>::max();
96  Index i0;
97  bool found_i0=false;
98  for(typename ItemMap::const_iterator map_it = map_.begin();
99  map_it!=map_.end();++map_it) {
100  Real dist2=geom::Length2(pos-gen_middle(map_it->first));
101  if(dist2<best_dist2 && !map_it->second.empty()) {
102  best_dist2=dist2;
103  i0 = map_it->first;
104  found_i0=true;
105  }
106  }
107 
108  if(!found_i0) return ITEM();
109 
110  // now find the closest item in the 3x3x3 items centered on i0
111  best_dist2=std::numeric_limits<Real>::max();
112  ITEM* best_item=NULL;
113 
114  for(int wc=i0.w-1;wc<=i0.w+1;++wc) {
115  for(int vc=i0.v-1;wc<=i0.v+1;++vc) {
116  for(int uc=i0.u-1;wc<=i0.u+1;++uc) {
117  typename ItemMap::const_iterator map_it = map_.find(Index(uc,vc,wc));
118 
119  if(map_it!=map_.end()) {
120  for(typename EntryList::iterator entry_it = map_it->second.begin();
121  entry_it != map_it->second.end(); ++entry_it) {
122  Real delta_x = entry_it->pos[0]-pos[0];
123  Real delta_y = entry_it->pos[1]-pos[1];
124  Real delta_z = entry_it->pos[2]-pos[2];
125  Real dist2=delta_x*delta_x+delta_y*delta_y+delta_z*delta_z;
126  if(dist2<best_dist2) {
127  best_dist2=dist2;
128  best_item=entry_it;
129  }
130  }
131  }
132  }
133  }
134  }
135  if(best_item) {
136  return &best_item;
137  } else {
138  return ITEM();
139  }
140  }
141 
142  ItemList FindWithin(const VEC& pos, Real dist) const {
143  Real dist2=dist*dist;
144  Index imin = gen_index(pos-VEC(dist,dist,dist));
145  Index imax = gen_index(pos+VEC(dist,dist,dist));
146 
147  ItemList item_list;
148 
149  for(int wc=imin.w;wc<=imax.w;++wc) {
150  for(int vc=imin.v;vc<=imax.v;++vc) {
151  for(int uc=imin.u;uc<=imax.u;++uc) {
152  typename ItemMap::const_iterator map_it = map_.find(Index(uc,vc,wc));
153 
154  if(map_it!=map_.end()) {
155  for(typename EntryList::const_iterator entry_it = map_it->second.begin();
156  entry_it != map_it->second.end(); ++entry_it) {
157  /*
158  speed tests indicate that pre-calculating dx2 or dy2
159  and pre-checking them with an additional if gives little
160  speed improvement for very specific circumstances only,
161  but most of the time the performance is worse.
162  */
163  Real delta_x = entry_it->pos[0]-pos[0];
164  Real delta_y = entry_it->pos[1]-pos[1];
165  Real delta_z = entry_it->pos[2]-pos[2];
166  if(delta_x*delta_x+delta_y*delta_y+delta_z*delta_z<=dist2) {
167  item_list.push_back(entry_it->item);
168  }
169  }
170  }
171  }
172  }
173  }
174  return item_list;
175  }
176  void Clear()
177  {
178  map_.clear();
179  }
181  map_.swap(o.map_);
182  std::swap(delta_,o.delta_);
183  }
184 
185 private:
186 
187  ItemMap map_;
188  Real delta_;
189 
190  Index gen_index(const VEC& pos) const {
191  Index nrvo(static_cast<int>(round(pos[0]/delta_)),
192  static_cast<int>(round(pos[1]/delta_)),
193  static_cast<int>(round(pos[2]/delta_)));
194  return nrvo;
195  }
196 
197  geom::Vec3 gen_middle(const Index& i) const {
198  geom::Vec3 nrvo((static_cast<Real>(i.u)+0.5)*delta_,
199  (static_cast<Real>(i.v)+0.5)*delta_,
200  (static_cast<Real>(i.w)+0.5)*delta_);
201  return nrvo;
202  }
203 
204 };
205 
206 }} // ns
207 
208 #endif