IMP  2.2.1
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-2014 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 //! Geometric Hash table.
22 /** GeometricHash<T, D> stores values of type T with D-dimensional
23  points whose coordinates are of type double inside D-dimensional
24  buckets (cubes) */
25 template <typename T, size_t D>
27  public:
28  enum Distance {
29  INF,
30  EUCLIDEAN
31  };
32  typedef algebra::VectorD<D> Point;
33  typedef boost::array<int, D> Bucket;
34  typedef std::pair<Point, T> ValueType;
35  typedef std::vector<ValueType> PointList;
36  typedef std::map<const Bucket, PointList> GeomMap;
37  typedef typename GeomMap::iterator iterator;
38  typedef typename GeomMap::const_iterator const_iterator;
39  typedef std::vector<const ValueType *> HashResult;
41  typedef std::vector<Bucket> BucketList;
42 
43  /* Default constructor - all cubes/buckets have the same edge length */
44  GeometricHash(double radius = 1) {
45  for (size_t i = 0; i < D; ++i) radii_[i] = radius;
46  }
47 
48  /* This constructor allows to specify the length of cube edges per axis */
49  template <typename X>
50  GeometricHash(const X &radii) {
51  for (size_t i = 0; i < D; ++i) radii_[i] = radii[i];
52  }
53 
54  // For Swig
55  GeometricHash(const algebra::VectorD<D> &radii) {
56  for (size_t i = 0; i < D; ++i) radii_[i] = radii[i];
57  }
58 
59  /* Add a point with the associated value into the hash map */
60  void add(const Point &pt, const T &data) {
61  Bucket b = to_bucket(pt);
62  gmap_[b].push_back(std::make_pair(pt, data));
63  }
64 
65  /* Add a point with the associated value into the hash map */
66  template <typename X>
67  void add(const X &coordinates, const T &data) {
68  Point pt;
69  for (size_t i = 0; i < D; ++i) pt[i] = coordinates[i];
70  Bucket b = to_bucket(pt);
71  gmap_[b].push_back(std::make_pair(pt, data));
72  }
73 
74  /* Find all points that are within radius distance from center.
75  * The vector returned have elements that point to (point, value)
76  * pairs that are in the neighborhood. Distance specify which norm
77  * to use (Euclidean or l-infinity). The result is invalidated
78  * whenever add function is called. */
79  HashResult neighbors(Distance dt, const Point &center, double radius) const {
80  if (dt == EUCLIDEAN)
81  return points_in_sphere(inside_sphere(center, radius));
82  else
83  return points_in_sphere(inside_sphere_inf(center, radius));
84  }
85 
86  // For Swig
87  HashResultT neighbors_data(Distance dt, const Point &center,
88  double radius) const {
89  HashResult r = neighbors(dt, center, radius);
90  HashResultT tr;
91  tr.reserve(r.size());
92  for (size_t i = 0; i < r.size(); ++i) tr.push_back(r[i]->second);
93  return tr;
94  }
95 
96  template <typename X>
97  HashResult neighbors(Distance dt, const X &center_coordinates,
98  double radius) const {
99  Point center;
100  for (size_t i = 0; i < D; ++i) center[i] = center_coordinates[i];
101  if (dt == EUCLIDEAN)
102  return points_in_sphere(inside_sphere(center, radius));
103  else
104  return points_in_sphere(inside_sphere_inf(center, radius));
105  }
106 
107  void neighbors(Distance dt, const Point &center, double radius,
108  HashResultT &values) {
109  values.clear();
110  HashResult nb = neighbors(dt, center, radius);
111  values.reserve(nb.size());
112  for (size_t i = 0; i < nb.size(); ++i) values.push_back(nb[i]->second);
113  }
114 
115  /* Find all points inside a cube centered in center with edges
116  * radii*2 long */
117  template <typename X>
118  HashResult neighbors(const Point &center, const X &radii) const {
119  return points_in_sphere(inside_cube(center, radii));
120  }
121 
122  /* Find all points inside a cube centered in center with edges
123  * 2 times bigger than the grid */
124  HashResult neighbors(const Point &center) {
125  return neighbors(center, radii_);
126  }
127 
128  template <typename X>
129  HashResult neighbors(const X &center_coordinates) {
130  Point center;
131  for (size_t i = 0; i < D; ++i) center[i] = center_coordinates[i];
132  return neighbors(center, radii_);
133  }
134 
135  /* Return the hash map */
136  GeomMap const &Map() const { return gmap_; }
137 
138  iterator begin() { return gmap_.begin(); }
139 
140  iterator end() { return gmap_.end(); }
141 
142  const_iterator begin() const { return gmap_.begin(); }
143 
144  const_iterator end() const { return gmap_.end(); }
145 
146  /* Return vector of buckets */
147  BucketList get_buckets() const {
148  BucketList result;
149  for (typename GeomMap::const_iterator p = gmap_.begin(); p != gmap_.end();
150  ++p)
151  result.push_back(p->first);
152  return result;
153  }
154 
155  /* Remove all data from the hash map */
156  void clear() { gmap_.clear(); }
157 
158  PointList const *bucket(const Point &pt) const {
159  typename GeomMap::const_iterator p = gmap_.find(to_bucket(pt));
160  if (p != gmap_.end()) return &p->second;
161  return 0;
162  }
163 
164  double resolution(int axis = 0) const { return radii_[axis]; }
165 
166  private:
167  struct inside_sphere {
168  inside_sphere(const Point &center, double radius)
169  : center_(center), radius_(radius), sq_radius_(radius * radius) {}
170 
171  bool operator()(const Point &pt) const {
172  return get_squared_distance(center_, pt) <= sq_radius_;
173  }
174 
175  double radius(int) const { return radius_; }
176 
177  const Point &center_;
178  double radius_;
179  double sq_radius_;
180  };
181 
182  struct inside_sphere_inf {
183  inside_sphere_inf(const Point &center, double radius)
184  : center_(center), radius_(radius) {}
185 
186  bool operator()(const Point &pt) const {
187  return inf_distance(center_, pt) <= radius_;
188  }
189 
190  double radius(int) const { return radius_; }
191 
192  const Point &center_;
193  double radius_;
194  };
195 
196  struct inside_cube {
197  template <typename X>
198  inside_cube(const Point &center, const X &radii)
199  : center_(center) {
200  for (size_t i = 0; i < D; ++i) radii_[i] = radii[i];
201  }
202 
203  bool operator()(const Point &pt) const {
204  for (size_t i = 0; i < D; ++i)
205  if (std::fabs(pt[i] - center_[i]) > radii_[i]) return false;
206  return true;
207  }
208 
209  double radius(size_t i) const { return radii_[i]; }
210 
211  const Point &center_;
212  double radii_[D];
213  };
214 
215  template <typename Dist>
216  HashResult points_in_sphere(Dist const &dist) const {
217  Point C;
218  for (size_t i = 0; i < D; ++i) C[i] = dist.center_[i] - dist.radius(i);
219  Bucket l = to_bucket(C);
220  for (size_t i = 0; i < D; ++i) C[i] = dist.center_[i] + dist.radius(i);
221  Bucket u = to_bucket(C);
222  HashResult result;
223  Bucket tmp;
224  points_in_sphere_rec(0, l, u, dist, tmp, result);
225  return result;
226  }
227 
228  template <typename Dist>
229  bool cube_inside_sphere(const Dist &ins, const Bucket &l) const {
230  Point p = from_bucket(l);
231  return cube_inside_sphere_rec(ins, p, 0);
232  }
233 
234  template <typename Dist>
235  bool cube_inside_sphere_rec(const Dist &ins, Point &p, size_t idx) const {
236  if (idx >= D) return ins(p);
237  if (!cube_inside_sphere_rec(ins, p, idx + 1)) return false;
238  double old = p[idx];
239  p[idx] += radii_[idx];
240  bool r = cube_inside_sphere_rec(ins, p, idx + 1);
241  p[idx] = old;
242  return r;
243  }
244 
245  template <typename Dist>
246  void points_in_sphere_rec(size_t idx, const Bucket &l, const Bucket &u,
247  const Dist &ins, Bucket &tmp,
248  HashResult &result) const {
249  if (idx >= D) {
250  typename GeomMap::const_iterator p = gmap_.find(tmp);
251  if (p != gmap_.end()) {
252  const PointList &v = p->second;
253  if (v.size() > (1 << D) && cube_inside_sphere(ins, tmp)) {
254  for (typename PointList::const_iterator q = v.begin(); q != v.end();
255  ++q)
256  result.push_back(&(*q));
257  } else {
258  for (typename PointList::const_iterator q = v.begin(); q != v.end();
259  ++q)
260  if (ins(q->first)) result.push_back(&(*q));
261  }
262  }
263  return;
264  }
265  for (int i = l[idx]; i <= u[idx]; ++i) {
266  tmp[idx] = i;
267  points_in_sphere_rec(idx + 1, l, u, ins, tmp, result);
268  }
269  }
270 
271  /* Create bucket containing the given point p */
272  Bucket to_bucket(const Point &p) const {
273  Bucket r;
274  for (size_t i = 0; i < D; ++i) {
275  if (p[i] < 0)
276  r[i] = int(p[i] / radii_[i]) - 1;
277  else
278  r[i] = int(p[i] / radii_[i]);
279  }
280  return r;
281  }
282 
283  /* Find the vertex of bucket b with the smallest coordinates */
284  Point from_bucket(const Bucket &b) const {
285  Point p(b.begin(), b.end());
286  for (size_t i = 0; i < D; ++i) p[i] *= radii_[i];
287  return p;
288  }
289 
290  /* l-infinity distance */
291  static double inf_distance(const Point &p1, const Point &p2) {
292  double r = std::fabs(p1[0] - p2[0]);
293  for (size_t i = 1; i < D; ++i) {
294  double d = std::fabs(p1[i] - p2[i]);
295  if (r < d) r = d;
296  }
297  return r;
298  }
299 
300  GeomMap gmap_;
301  double radii_[D];
302 };
303 
304 IMPMULTIFIT_END_NAMESPACE
305 #endif /* IMPMULTIFIT_GEOMETRIC_HASH_H */
double get_squared_distance(const VectorD< D > &v1, const VectorD< D > &v2)
compute the squared distance between two vectors
Definition: VectorD.h:201
A Cartesian vector in D-dimensions.
Definition: VectorD.h:52
Simple D vector class.
Geometric Hash table.
Definition: GeometricHash.h:26