IMP  2.0.0
The Integrative Modeling Platform
GeometricHash.h
Go to the documentation of this file.
1 /**
2  * \file IMP/multifit/GeometricHash.h \brief Geometric Hashing class.
3  *
4  * Copyright 2007-2013 IMP Inventors. All rights reserved.
5  *
6  */
7 
8 #ifndef IMPMULTIFIT_GEOMETRIC_HASH_H
9 #define IMPMULTIFIT_GEOMETRIC_HASH_H
10 
11 #include <vector>
12 #include <map>
13 #include <ostream>
14 #include <cmath>
15 #include <boost/array.hpp>
16 #include <IMP/multifit/multifit_config.h>
17 #include "IMP/algebra/VectorD.h"
18 
19 IMPMULTIFIT_BEGIN_NAMESPACE
20 
21 
22 /* This is the definition of Geometric Hash table.
23  * GeometricHash<T, D> stores values of type T with D-dimensional
24  * points whose coordinates are of type double inside D-dimensional
25  * buckets (cubes) */
26 template < typename T, size_t D >
27 class GeometricHash
28 {
29 public:
30  enum Distance { INF, EUCLIDEAN };
31  typedef algebra::VectorD<D> Point;
32  typedef boost::array<int, D> Bucket;
33  typedef std::pair< Point, T > ValueType;
34  typedef std::vector< ValueType > PointList;
35  typedef std::map<const Bucket, PointList> GeomMap;
36  typedef typename GeomMap::iterator iterator;
37  typedef typename GeomMap::const_iterator const_iterator;
38  typedef std::vector< const ValueType *> HashResult;
39  typedef IMP::base::Vector<T> HashResultT;
40  typedef std::vector<Bucket> BucketList;
41 
42 
43  /* Default constructor - all cubes/buckets have the same edge length */
44  GeometricHash(double radius = 1)
45  {
46  for ( size_t i=0; i<D; ++i )
47  radii_[i] = radius;
48  }
49 
50  /* This constructor allows to specify the length of cube edges per axis */
51  template < typename X >
52  GeometricHash(const X &radii)
53  {
54  for ( size_t i=0; i<D; ++i )
55  radii_[i] = radii[i];
56  }
57 
58  // For Swig
59  GeometricHash(const algebra::VectorD<D> &radii)
60  {
61  for ( size_t i=0; i<D; ++i )
62  radii_[i] = radii[i];
63  }
64 
65 
66  /* Add a point with the associated value into the hash map */
67  void add(const Point &pt, const T &data)
68  {
69  Bucket b = to_bucket(pt);
70  gmap_[b].push_back(std::make_pair(pt, data));
71  }
72 
73 
74  /* Add a point with the associated value into the hash map */
75  template < typename X >
76  void add(const X &coordinates, const T &data)
77  {
78  Point pt;
79  for ( size_t i=0; i<D; ++i )
80  pt[i] = coordinates[i];
81  Bucket b = to_bucket(pt);
82  gmap_[b].push_back(std::make_pair(pt, data));
83  }
84 
85 
86  /* Find all points that are within radius distance from center.
87  * The vector returned have elements that point to (point, value)
88  * pairs that are in the neighborhood. Distance specify which norm
89  * to use (Euclidean or l-infinity). The result is invalidated
90  * whenever add function is called. */
91  HashResult neighbors(Distance dt, const Point &center,
92  double radius) const
93  {
94  if ( dt == EUCLIDEAN )
95  return points_in_sphere(inside_sphere(center, radius));
96  else
97  return points_in_sphere(inside_sphere_inf(center, radius));
98  }
99 
100  // For Swig
101  HashResultT neighbors_data(Distance dt, const Point &center,
102  double radius) const
103  {
104  HashResult r = neighbors(dt, center, radius);
105  HashResultT tr;
106  tr.reserve(r.size());
107  for ( size_t i=0; i<r.size(); ++i )
108  tr.push_back(r[i]->second);
109  return tr;
110  }
111 
112  template < typename X >
113  HashResult neighbors(Distance dt, const X &center_coordinates,
114  double radius) const
115  {
116  Point center;
117  for ( size_t i=0; i<D; ++i )
118  center[i] = center_coordinates[i];
119  if ( dt == EUCLIDEAN )
120  return points_in_sphere(inside_sphere(center, radius));
121  else
122  return points_in_sphere(inside_sphere_inf(center, radius));
123  }
124 
125  void neighbors(Distance dt, const Point &center,
126  double radius, HashResultT &values)
127  {
128  values.clear();
129  HashResult nb = neighbors(dt, center, radius);
130  values.reserve(nb.size());
131  for ( size_t i=0; i<nb.size(); ++i )
132  values.push_back(nb[i]->second);
133  }
134 
135 
136  /* Find all points inside a cube centered in center with edges
137  * radii*2 long */
138  template < typename X >
139  HashResult neighbors(const Point &center,
140  const X &radii) const
141  {
142  return points_in_sphere(inside_cube(center, radii));
143  }
144 
145  /* Find all points inside a cube centered in center with edges
146  * 2 times bigger than the grid */
147  HashResult neighbors(const Point &center)
148  {
149  return neighbors(center, radii_);
150  }
151 
152  template < typename X >
153  HashResult neighbors(const X &center_coordinates)
154  {
155  Point center;
156  for ( size_t i=0; i<D; ++i )
157  center[i] = center_coordinates[i];
158  return neighbors(center, radii_);
159  }
160 
161  /* Return the hash map */
162  GeomMap const &Map() const
163  {
164  return gmap_;
165  }
166 
167  iterator begin()
168  {
169  return gmap_.begin();
170  }
171 
172  iterator end()
173  {
174  return gmap_.end();
175  }
176 
177  const_iterator begin() const
178  {
179  return gmap_.begin();
180  }
181 
182  const_iterator end() const
183  {
184  return gmap_.end();
185  }
186 
187  /* Return vector of buckets */
188  BucketList get_buckets() const
189  {
190  BucketList result;
191  for ( typename GeomMap::const_iterator p = gmap_.begin();
192  p != gmap_.end(); ++p )
193  result.push_back(p->first);
194  return result;
195  }
196 
197  /* Remove all data from the hash map */
198  void clear()
199  {
200  gmap_.clear();
201  }
202 
203  PointList const *bucket(const Point &pt) const
204  {
205  typename GeomMap::const_iterator p = gmap_.find(to_bucket(pt));
206  if ( p != gmap_.end() )
207  return &p->second;
208  return 0;
209  }
210 
211  double resolution(int axis = 0) const
212  {
213  return radii_[axis];
214  }
215 
216 
217 private:
218 
219  struct inside_sphere
220  {
221  inside_sphere(const Point &center, double radius)
222  : center_(center)
223  , radius_(radius)
224  , sq_radius_(radius*radius)
225  {}
226 
227  bool operator()(const Point &pt) const
228  {
229  return get_squared_distance(center_, pt) <= sq_radius_;
230  }
231 
232  double radius(int) const
233  {
234  return radius_;
235  }
236 
237  const Point &center_;
238  double radius_;
239  double sq_radius_;
240  };
241 
242  struct inside_sphere_inf
243  {
244  inside_sphere_inf(const Point &center, double radius)
245  : center_(center)
246  , radius_(radius)
247  {}
248 
249  bool operator()(const Point &pt) const
250  {
251  return inf_distance(center_, pt) <= radius_;
252  }
253 
254  double radius(int) const
255  {
256  return radius_;
257  }
258 
259  const Point &center_;
260  double radius_;
261  };
262 
263  struct inside_cube
264  {
265  template < typename X >
266  inside_cube(const Point &center, const X &radii)
267  : center_(center)
268  {
269  for ( size_t i=0; i<D; ++i )
270  radii_[i] = radii[i];
271  }
272 
273  bool operator()(const Point &pt) const
274  {
275  for ( size_t i=0; i<D; ++i )
276  if ( std::fabs(pt[i] - center_[i]) > radii_[i] )
277  return false;
278  return true;
279  }
280 
281  double radius(size_t i) const
282  {
283  return radii_[i];
284  }
285 
286  const Point &center_;
287  double radii_[D];
288  };
289 
290  template < typename Dist >
291  HashResult points_in_sphere(Dist const &dist) const
292  {
293  Point C;
294  for ( size_t i=0; i<D; ++i )
295  C[i] = dist.center_[i] - dist.radius(i);
296  Bucket l = to_bucket(C);
297  for ( size_t i=0; i<D; ++i )
298  C[i] = dist.center_[i] + dist.radius(i);
299  Bucket u = to_bucket(C);
300  HashResult result;
301  Bucket tmp;
302  points_in_sphere_rec(0, l, u, dist, tmp, result);
303  return result;
304  }
305 
306  template < typename Dist >
307  bool cube_inside_sphere(const Dist &ins, const Bucket &l) const
308  {
309  Point p = from_bucket(l);
310  return cube_inside_sphere_rec(ins, p, 0);
311  }
312 
313  template < typename Dist >
314  bool cube_inside_sphere_rec(const Dist &ins,
315  Point &p, size_t idx) const
316  {
317  if ( idx >= D )
318  return ins(p);
319  if (! cube_inside_sphere_rec(ins, p, idx + 1) )
320  return false;
321  double old = p[idx];
322  p[idx] += radii_[idx];
323  bool r = cube_inside_sphere_rec(ins, p, idx + 1);
324  p[idx] = old;
325  return r;
326  }
327 
328  template < typename Dist >
329  void points_in_sphere_rec(size_t idx, const Bucket &l, const Bucket &u,
330  const Dist &ins, Bucket &tmp,
331  HashResult &result) const
332  {
333  if ( idx >= D )
334  {
335  typename GeomMap::const_iterator p = gmap_.find(tmp);
336  if ( p != gmap_.end() )
337  {
338  const PointList &v = p->second;
339  if ( v.size() > (1<<D) && cube_inside_sphere(ins, tmp) )
340  {
341  for ( typename PointList::const_iterator q = v.begin();
342  q != v.end(); ++q )
343  result.push_back(&(*q));
344  }
345  else
346  {
347  for ( typename PointList::const_iterator q = v.begin();
348  q != v.end(); ++q )
349  if ( ins(q->first) )
350  result.push_back(&(*q));
351  }
352  }
353  return;
354  }
355  for ( int i=l[idx]; i<=u[idx]; ++i )
356  {
357  tmp[idx] = i;
358  points_in_sphere_rec(idx + 1, l, u, ins, tmp, result);
359  }
360  }
361 
362  /* Create bucket containing the given point p */
363  Bucket to_bucket(const Point &p) const
364  {
365  Bucket r;
366  for ( size_t i=0; i<D; ++i )
367  {
368  if ( p[i] < 0 )
369  r[i] = int(p[i]/radii_[i]) - 1;
370  else
371  r[i] = int(p[i]/radii_[i]);
372  }
373  return r;
374  }
375 
376  /* Find the vertex of bucket b with the smallest coordinates */
377  Point from_bucket(const Bucket &b) const
378  {
379  Point p(b.begin(), b.end());
380  for ( size_t i=0; i<D; ++i )
381  p[i] *= radii_[i];
382  return p;
383  }
384 
385  /* l-infinity distance */
386  static double inf_distance(const Point &p1, const Point &p2)
387  {
388  double r = std::fabs(p1[0] - p2[0]);
389  for ( size_t i=1; i<D; ++i )
390  {
391  double d = std::fabs(p1[i] - p2[i]);
392  if ( r < d )
393  r = d;
394  }
395  return r;
396  }
397 
398  GeomMap gmap_;
399  double radii_[D];
400 };
401 
402 IMPMULTIFIT_END_NAMESPACE
403 #endif /* IMPMULTIFIT_GEOMETRIC_HASH_H */