8 #ifndef IMPMULTIFIT_GEOMETRIC_HASH_H
9 #define IMPMULTIFIT_GEOMETRIC_HASH_H
16 #include <IMP/multifit/multifit_config.h>
19 IMPMULTIFIT_BEGIN_NAMESPACE
25 template <
typename T,
size_t D>
33 typedef std::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;
45 for (
size_t i = 0; i < D; ++i) radii_[i] = radius;
51 for (
size_t i = 0; i < D; ++i) radii_[i] = radii[i];
56 for (
size_t i = 0; i < D; ++i) radii_[i] = radii[i];
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));
67 void add(
const X &coordinates,
const T &data) {
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));
79 HashResult neighbors(Distance dt,
const Point ¢er,
double radius)
const {
81 return points_in_sphere(inside_sphere(center, radius));
83 return points_in_sphere(inside_sphere_inf(center, radius));
88 double radius)
const {
89 HashResult r = neighbors(dt, center, radius);
92 for (
size_t i = 0; i < r.size(); ++i) tr.push_back(r[i]->second);
97 HashResult neighbors(Distance dt,
const X ¢er_coordinates,
98 double radius)
const {
100 for (
size_t i = 0; i < D; ++i) center[i] = center_coordinates[i];
102 return points_in_sphere(inside_sphere(center, radius));
104 return points_in_sphere(inside_sphere_inf(center, radius));
107 void neighbors(Distance dt,
const Point ¢er,
double radius,
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);
117 template <
typename X>
118 HashResult neighbors(
const Point ¢er,
const X &radii)
const {
119 return points_in_sphere(inside_cube(center, radii));
124 HashResult neighbors(
const Point ¢er) {
125 return neighbors(center, radii_);
128 template <
typename X>
129 HashResult neighbors(
const X ¢er_coordinates) {
131 for (
size_t i = 0; i < D; ++i) center[i] = center_coordinates[i];
132 return neighbors(center, radii_);
136 GeomMap
const &
Map()
const {
return gmap_; }
138 iterator begin() {
return gmap_.begin(); }
140 iterator end() {
return gmap_.end(); }
142 const_iterator begin()
const {
return gmap_.begin(); }
144 const_iterator end()
const {
return gmap_.end(); }
149 for (
typename GeomMap::const_iterator p = gmap_.begin(); p != gmap_.end();
151 result.push_back(p->first);
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;
164 double resolution(
int axis = 0)
const {
return radii_[axis]; }
167 struct inside_sphere {
168 inside_sphere(
const Point ¢er,
double radius)
169 : center_(center), radius_(radius), sq_radius_(radius * radius) {}
171 bool operator()(
const Point &pt)
const {
175 double radius(
int)
const {
return radius_; }
177 const Point ¢er_;
182 struct inside_sphere_inf {
183 inside_sphere_inf(
const Point ¢er,
double radius)
184 : center_(center), radius_(radius) {}
186 bool operator()(
const Point &pt)
const {
187 return inf_distance(center_, pt) <= radius_;
190 double radius(
int)
const {
return radius_; }
192 const Point ¢er_;
197 template <
typename X>
198 inside_cube(
const Point ¢er,
const X &radii)
200 for (
size_t i = 0; i < D; ++i) radii_[i] = radii[i];
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;
209 double radius(
size_t i)
const {
return radii_[i]; }
211 const Point ¢er_;
215 template <
typename Dist>
216 HashResult points_in_sphere(Dist
const &dist)
const {
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);
224 points_in_sphere_rec(0, l, u, dist, tmp, result);
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);
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;
239 p[idx] += radii_[idx];
240 bool r = cube_inside_sphere_rec(ins, p, idx + 1);
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 {
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();
256 result.push_back(&(*q));
258 for (
typename PointList::const_iterator q = v.begin(); q != v.end();
260 if (ins(q->first)) result.push_back(&(*q));
265 for (
int i = l[idx]; i <= u[idx]; ++i) {
267 points_in_sphere_rec(idx + 1, l, u, ins, tmp, result);
272 Bucket to_bucket(
const Point &p)
const {
274 for (
size_t i = 0; i < D; ++i) {
276 r[i] = int(p[i] / radii_[i]) - 1;
278 r[i] = int(p[i] / radii_[i]);
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];
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]);
304 IMPMULTIFIT_END_NAMESPACE
BucketList get_buckets() const
Return vector of buckets.
double get_squared_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the squared distance between two vectors.
GeomMap const & Map() const
Return the hash map.
A more IMP-like version of the std::vector.
A Cartesian vector in D-dimensions.
void clear()
Remove all data from the hash map.