IMP logo
IMP Reference Guide  develop.29b24c6ad6,2021/05/15
The Integrative Modeling Platform
grid_storages.h
Go to the documentation of this file.
1 /**
2  * \file IMP/algebra/grid_storages.h
3  * \brief A class to represent a voxel grid.
4  *
5  * Copyright 2007-2021 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPALGEBRA_GRID_STORAGES_H
10 #define IMPALGEBRA_GRID_STORAGES_H
11 
12 #include <IMP/algebra/algebra_config.h>
13 #include <IMP/bracket_macros.h>
14 #include <IMP/types.h>
15 #include "internal/grid_apply.h"
16 #include "grid_indexes.h"
17 #include "grid_ranges.h"
18 #include "Vector3D.h"
19 #include "BoundingBoxD.h"
20 #include <boost/iterator/transform_iterator.hpp>
21 #include <boost/unordered_map.hpp>
22 #include <IMP/Vector.h>
23 
24 #include <limits>
25 
26 IMPALGEBRA_BEGIN_NAMESPACE
27 
28 /** Store a grid as a densely packed set of voxels.
29  The mapping from GridIndex3D to index in the data array is:
30  \code
31  i[2]*BoundedGridRangeD::get_number_of_voxels(0)
32  *BoundedGridRange3D::get_number_of_voxels(1)
33  + i[1]*BoundedGridRange3D::get_number_of_voxels(0)+i[0]
34  \endcode
35  \see Grid3D
36 */
37 template <int D, class VT>
39  typedef boost::scoped_array<VT> Data;
40  Data data_;
41  unsigned int extent_;
42  VT default_;
43 
44  template <class I>
45  unsigned int index(const I &i) const {
46  unsigned int ii = 0;
47  for (int d = D - 1; d >= 0; --d) {
48  unsigned int cur = i[d];
49  for (int ld = d - 1; ld >= 0; --ld) {
51  }
52  ii += cur;
53  }
55  if (D == 3) {
56  unsigned int check =
60  IMP_UNUSED(check);
61  IMP_USAGE_CHECK(check == ii, "Wrong value returned");
62  }
63  }
64  IMP_INTERNAL_CHECK(ii < get_extent(), "Invalid grid index " << i);
65  return ii;
66  }
67 #ifndef SWIG
68  struct NonDefault {
69  VT default_;
70  NonDefault(const VT &def) : default_(def) {}
71  template <class P>
72  bool operator()(const P &def) const {
73  return def.second != default_;
74  }
75  };
76 #endif
77  unsigned int get_extent() const { return extent_; }
78  void copy_from(const DenseGridStorageD &o) {
79  default_ = o.default_;
80  extent_ = o.extent_;
81  data_.reset(new VT[extent_]);
82  std::copy(o.data_.get(), o.data_.get() + o.extent_, data_.get());
84  }
85  void set_number_of_voxels(Ints dims) {
86  extent_ = 1;
87  for (unsigned int i = 0; i < dims.size(); ++i) {
88  extent_ *= dims[i];
89  }
90  data_.reset(new VT[extent_]);
91  std::fill(data_.get(), data_.get() + get_extent(), default_);
92  // BoundedGridRangeD<D>::set_number_of_voxels(dims);
93  }
94 
95  public:
97  typedef VT Value;
98  DenseGridStorageD(const Ints &counts, const VT &default_value = VT())
99  : BoundedGridRangeD<D>(counts), default_(default_value) {
100  set_number_of_voxels(counts);
101  }
102  IMP_BRACKET(VT, GridIndexD<D>, true, return data_[index(i)]);
103  /** \name Direct indexing
104  One can directly access a particular voxel based on its index
105  in the array of all voxels. This can be faster than using a
106  GridIndexD.
107  @{
108  */
109  IMP_BRACKET(VT, unsigned int, i < extent_, return data_[i]);
110 /** @}
111  */
112 #ifndef IMP_DOXYGEN
113  DenseGridStorageD(const VT &default_value = VT())
114  : extent_(0), default_(default_value) {}
115  static bool get_is_dense() { return true; }
116 #endif
118 #ifndef SWIG
119  const VT *get_raw_data() const { return data_.get(); }
120  VT *get_raw_data() { return data_.get(); }
121 #endif
122 
123 #ifndef IMP_DOXYGEN
124  GridIndexD<D> add_voxel(const ExtendedGridIndexD<D> &, const VT &) {
125  IMP_FAILURE("Cannot add voxels to dense grid");
126  }
127 #endif
128 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
129  VT &get_voxel_always(const ExtendedGridIndexD<D> &i) {
130  GridIndexD<D> gi(i.begin(), i.end());
131  return operator[](gi);
132  }
133  const VT &get_value_always(const ExtendedGridIndexD<D> &i) const {
134  GridIndexD<D> gi(i.begin(), i.end());
135  return operator[](gi);
136  }
137 #endif
138 
139 /** \name All voxel iterators
140  The value type is VT.
141  @{
142 */
143 #ifndef SWIG
144  typedef VT *AllVoxelIterator;
145  typedef const VT *AllVoxelConstIterator;
146  AllVoxelIterator all_voxels_begin() { return data_.get(); }
147  AllVoxelIterator all_voxels_end() { return data_.get() + get_extent(); }
148  AllVoxelConstIterator all_voxels_begin() const { return data_.get(); }
149  AllVoxelConstIterator all_voxels_end() const {
150  return data_.get() + get_extent();
151  }
152 #endif
153  Vector<VT> get_all_voxels() const {
154  return Vector<VT>(all_voxels_begin(), all_voxels_end());
155  }
156  /** @} */
157  template <class Functor, class Grid>
158  Functor apply(const Grid &g, const Functor &fi) const {
159  Functor f = fi;
160  typename Grid::ExtendedIndex lb(typename Grid::ExtendedIndex::Filled(),
161  g.get_dimension(), 0);
162  typename Grid::ExtendedIndex ub(BoundedGridRangeD<D>::get_end_index());
163  typename Grid::Vector corner = g.get_bounding_box().get_corner(0);
164  typename Grid::Vector cell = g.get_unit_cell();
165  typename Grid::Index cur;
166  typename Grid::Vector center;
167  internal::GridApplier<Functor, Grid, D - 1>::apply(g, lb, ub, corner, cell,
168  cur, center, f);
169  return f;
170  }
171 };
172 
173 /** Store a grid as a sparse set of voxels (only the voxels which have
174  been added are actually stored). The
175  get_has_index() functions allow one to tell if a voxel has been added.
176  \unstable{SparseGridStorageD}
177 
178  Base should be one of BoundedGridRangeD or UnboundedGridRangeD.
179  \see Grid3D
180 */
181 template <int D, class VT, class Base,
182  class Map = typename boost::unordered_map<GridIndexD<D>, VT> >
183 class SparseGridStorageD : public Base {
184  typedef Map Data;
185  struct GetIndex {
186  typedef GridIndexD<D> result_type;
187  typedef typename Data::const_iterator::value_type argument_type;
188  template <class T>
189  GridIndexD<D> operator()(const T &t) const {
190  return t.first;
191  }
192  };
193  struct ItHelper {
194  const SparseGridStorageD<D, VT, Base> *stor_;
195  ItHelper(const SparseGridStorageD<D, VT, Base> *stor) : stor_(stor) {}
196  bool get_is_good(const ExtendedGridIndexD<D> &ei) {
197  /*std::cout << "Checking " << ei << " getting "
198  << stor_->get_has_index(ei) << std::endl;*/
199  return stor_->get_has_index(ei);
200  }
201  typedef GridIndexD<D> ReturnType;
202  ReturnType get_return(const ExtendedGridIndexD<D> &ei) const {
203  return stor_->get_index(ei);
204  }
205  ItHelper() : stor_(nullptr) {}
206  };
207 
208  Data data_;
209  VT default_;
210 
211  public:
212  typedef VT Value;
213  SparseGridStorageD(const Ints &counts, const VT &default_value)
214  : Base(counts), default_(default_value) {}
215  IMP_SHOWABLE_INLINE(SparseGridStorage3D, out << "Sparse grid with "
216  << data_.size() << " cells set");
217  //! Add a voxel to the storage, this voxel will now have a GridIndex3D
219  IMP_USAGE_CHECK(Base::get_has_index(i), "Out of grid domain " << i);
220  GridIndexD<D> ret(i.begin(), i.end());
221  data_[ret] = gi;
222  return ret;
223  }
224  void remove_voxel(const GridIndexD<D> &gi) { data_.erase(gi); }
225 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
226  SparseGridStorageD(const VT &def) : default_(def) {}
227  static bool get_is_dense() { return false; }
228  using Base::get_number_of_voxels;
229 #endif
230  unsigned int get_number_of_voxels() const { return data_.size(); }
231  //! Return true if the voxel has been added
232  bool get_has_index(const ExtendedGridIndexD<D> &i) const {
233  return data_.find(GridIndexD<D>(i.begin(), i.end())) != data_.end();
234  }
235  //! requires get_has_index(i) is true.
237  IMP_USAGE_CHECK(get_has_index(i), "Index is not a valid "
238  << "voxel " << i);
239  return GridIndexD<D>(i.begin(), i.end());
240  }
241 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
242  VT &get_voxel_always(const ExtendedGridIndexD<D> &i) {
243  GridIndexD<D> gi(i.begin(), i.end());
244  typename Map::iterator it = data_.find(gi);
245  if (it == data_.end()) {
246  return data_.insert(std::make_pair(gi, default_)).first->second;
247  } else {
248  return it->second;
249  }
250  }
251  const VT &get_value_always(const ExtendedGridIndexD<D> &i) const {
252  GridIndexD<D> gi(i.begin(), i.end());
253  typename Map::const_iterator it = data_.find(gi);
254  if (it == data_.end()) {
255  return default_;
256  } else {
257  return it->second;
258  }
259  }
260 #endif
261  /** \name Operator []
262  Operator[] isn't very useful at the moment as it can only
263  be used with a cell which has already been set. This
264  behavior/the existence of these functions is likely to change.
265  @{
266  */
267  IMP_BRACKET(VT, GridIndexD<D>, true,
268  IMP_USAGE_CHECK(data_.find(i) != data_.end(), "Invalid index "
269  << i);
270  return data_.find(i)->second);
271 /** @} */
272 
273 /** \name Iterators through set cells
274  Iterate through the voxels which have been set. The value
275  type is a pair of GridIndex3D and VT.
276  @{
277 */
278 #ifndef SWIG
279  typedef typename Data::const_iterator AllConstIterator;
280  AllConstIterator all_begin() const { return data_.begin(); }
281  AllConstIterator all_end() const { return data_.end(); }
282  typedef boost::transform_iterator<GetIndex, AllConstIterator>
283  AllIndexIterator;
284  AllIndexIterator all_indexes_begin() const {
285  return boost::make_transform_iterator(all_begin(), GetIndex());
286  }
287  AllIndexIterator all_indexes_end() const {
288  return boost::make_transform_iterator(all_end(), GetIndex());
289  }
290 #endif
291 /** @} */
292 
293 /** \name Index Iterators
294 
295  Iterate through a range of actual indexes. The value
296  type for the iterator is an GridIndex3D.
297 
298  The range is defined by a pair of indexes. It includes
299  all indexes in the axis aligned box defined by lb
300  as the lower corner and the second as the ub. That is, if
301  lb is \f$(l_x, l_y, l_z)\f$ and ub is \f$(u_x, u_y, u_z)\f$,
302  then the range includes all
303  indexes \f$(i_x, i_y, i_z)\f$ such that \f$l_x \leq i_x \leq u_x\f$,
304  \f$l_y \leq i_y \leq u_y\f$
305  and \f$l_z \leq i_z \leq u_z\f$.
306 
307  @{
308  */
309 #ifndef SWIG
310 #ifndef IMP_DOXYGEN
311 
312  typedef internal::GridIndexIterator<ExtendedGridIndexD<D>, ItHelper>
313  IndexIterator;
314 
315 #else
316  class IndexIterator;
317 #endif
318  IndexIterator indexes_begin(const ExtendedGridIndexD<D> &lb,
319  const ExtendedGridIndexD<D> &ub) const {
320  ExtendedGridIndexD<D> eub = ub.get_offset(1, 1, 1);
321  if (lb == ub) {
322  return IndexIterator();
323  } else {
324  IMP_INTERNAL_CHECK(internal::get_is_non_empty(lb, eub), "empty range");
325  return IndexIterator(lb, eub, ItHelper(this));
326  }
327  }
328  IndexIterator indexes_end(const ExtendedGridIndexD<D> &,
329  const ExtendedGridIndexD<D> &) const {
330  // IMP_INTERNAL_CHECK(lb <= ub, "empty range");
331  return IndexIterator();
332  }
333 #endif
334 
335  Vector<GridIndexD<D> > get_indexes(
336  const ExtendedGridIndexD<D> &lb, const ExtendedGridIndexD<D> &ub) const {
337  return Vector<GridIndexD<D> >(indexes_begin(lb, ub),
338  indexes_end(lb, ub));
339  }
340  /** @} */
341 
342  template <class Functor, class Grid>
343  Functor apply(const Grid &g, Functor f) const {
344  for (typename Data::const_iterator it = data_.begin(); it != data_.end();
345  ++it) {
346  f(g, it->first, g.get_center(it->first));
347  }
348  return f;
349  }
350 #ifndef IMP_DOXYGEN
351  //! Return the index which has no lower index in each coordinate
352  ExtendedGridIndexD<D> get_minimum_extended_index() const {
353  IMP_USAGE_CHECK(!data_.empty(), "No voxels in grid.");
354  GridIndexD<D> reti = data_.begin()->first;
355  ExtendedGridIndexD<D> ret(reti.begin(), reti.end());
356  for (typename Data::const_iterator it = data_.begin(); it != data_.end();
357  ++it) {
358  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
359  ret.access_data().get_data()[i] = std::min(ret[i], it->first[i]);
360  }
361  }
362  return ret;
363  }
364  //! Return the index that has no higher index in each coordinate
365  ExtendedGridIndexD<D> get_maximum_extended_index() const {
366  IMP_USAGE_CHECK(!data_.empty(), "No voxels in grid.");
367  GridIndexD<D> reti = data_.begin()->first;
368  ExtendedGridIndexD<D> ret(reti.begin(), reti.end());
369  for (typename Data::const_iterator it = data_.begin(); it != data_.end();
370  ++it) {
371  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
372  ret.access_data().get_data()[i] = std::min(ret[i], it->first[i]);
373  }
374  }
375  return ret;
376  }
377 #endif
378 };
379 IMPALGEBRA_END_NAMESPACE
380 
381 #endif /* IMPALGEBRA_GRID_STORAGES_H */
Basic types used by IMP.
#define IMP_SHOWABLE_INLINE(Name, how_to_show)
Declare the methods needed by an object that can be printed.
#define IMP_IF_CHECK(level)
Execute the code block if a certain level checks are on.
Definition: check_macros.h:104
#define IMP_FAILURE(message)
A runtime failure for IMP.
Definition: check_macros.h:72
A class to represent a voxel grid.
Represent a real cell in a grid (one within the bounding box)
Definition: grid_indexes.h:162
An index in an infinite grid on space.
Definition: grid_indexes.h:30
#define IMP_BRACKET(Value, Index, bounds_check_expr, expr)
Implement operator[] and at() for C++, and getitem for Python.
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Definition: check_macros.h:139
Macros to handle array indexing.
#define IMP_COPY_CONSTRUCTOR(Name, Base)
Use a copy_from method to create a copy constructor and operator=.
#define IMP_UNUSED(variable)
bool get_has_index(const ExtendedGridIndexD< D > &i) const
Return true if the voxel has been added.
A bounding box in D dimensions.
GridIndexD< D > get_index(const ExtendedGridIndexD< D > &i) const
requires get_has_index(i) is true.
A class for storing lists of IMP items.
Simple 3D vector class.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
A class to represent a voxel grid.
ParticleIndexes get_indexes(const ParticlesTemp &ps)
Get the indexes from a list of particles.
GridIndexD< D > add_voxel(const ExtendedGridIndexD< D > &i, const VT &gi)
Add a voxel to the storage, this voxel will now have a GridIndex3D.