IMP  2.1.0
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-2013 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>
14 #include <IMP/base/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 <IMP/base/map.h>
22 #include <IMP/base/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  base::Vector<VT> get_all_voxels() const {
154  return base::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 IMP::base::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(
268  VT, GridIndexD<D>, true,
269  IMP_USAGE_CHECK(data_.find(i) != data_.end(), "Invalid index " << 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 #endif
283 
284  base::Vector<GridIndexD<D> > get_all_indexes() const {
285  return base::Vector<GridIndexD<D> >(
286  boost::make_transform_iterator(all_begin(), GetIndex()),
287  boost::make_transform_iterator(all_end(), GetIndex()));
288  }
289 /** @} */
290 
291 /** \name Index Iterators
292 
293  Iterate through a range of actual indexes. The value
294  type for the iterator is an GridIndex3D.
295 
296  The range is defined by a pair of indexes. It includes
297  all indexes in the axis aligned box defined by lb
298  as the lower corner and the second as the ub. That is, if
299  lb is \f$(l_x, l_y, l_z)\f$ and ub is \f$(u_x, u_y, u_z)\f$,
300  then the range includes all
301  indexes \f$(i_x, i_y, i_z)\f$ such that \f$l_x \leq i_x \leq u_x\f$,
302  \f$l_y \leq i_y \leq u_y\f$
303  and \f$l_z \leq i_z \leq u_z\f$.
304 
305  @{
306  */
307 #ifndef SWIG
308 #ifndef IMP_DOXYGEN
309 
310  typedef internal::GridIndexIterator<ExtendedGridIndexD<D>, ItHelper>
311  IndexIterator;
312 
313 #else
314  class IndexIterator;
315 #endif
316  IndexIterator indexes_begin(const ExtendedGridIndexD<D> &lb,
317  const ExtendedGridIndexD<D> &ub) const {
318  ExtendedGridIndexD<D> eub = ub.get_offset(1, 1, 1);
319  if (lb == ub) {
320  return IndexIterator();
321  } else {
322  IMP_INTERNAL_CHECK(internal::get_is_non_empty(lb, eub), "empty range");
323  return IndexIterator(lb, eub, ItHelper(this));
324  }
325  }
326  IndexIterator indexes_end(const ExtendedGridIndexD<D> &,
327  const ExtendedGridIndexD<D> &) const {
328  // IMP_INTERNAL_CHECK(lb <= ub, "empty range");
329  return IndexIterator();
330  }
331 #endif
332 
333  base::Vector<GridIndexD<D> > get_indexes(
334  const ExtendedGridIndexD<D> &lb, const ExtendedGridIndexD<D> &ub) const {
335  return base::Vector<GridIndexD<D> >(indexes_begin(lb, ub),
336  indexes_end(lb, ub));
337  }
338  /** @} */
339 
340  template <class Functor, class Grid>
341  Functor apply(const Grid &g, Functor f) const {
342  for (typename Data::const_iterator it = data_.begin(); it != data_.end();
343  ++it) {
344  f(g, it->first, g.get_center(it->first));
345  }
346  return f;
347  }
348 
349  //! Return the index which has no lower index in each coordinate
351  IMP_USAGE_CHECK(!data_.empty(), "No voxels in grid.");
352  GridIndexD<D> reti = data_.begin()->first;
353  ExtendedGridIndexD<D> ret(reti.begin(), reti.end());
354  for (typename Data::const_iterator it = data_.begin(); it != data_.end();
355  ++it) {
356  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
357  ret.access_data().get_data()[i] = std::min(ret[i], it->first[i]);
358  }
359  }
360  return ret;
361  }
362  //! Return the index that has no higher index in each coordinate
364  IMP_USAGE_CHECK(!data_.empty(), "No voxels in grid.");
365  GridIndexD<D> reti = data_.begin()->first;
366  ExtendedGridIndexD<D> ret(reti.begin(), reti.end());
367  for (typename Data::const_iterator it = data_.begin(); it != data_.end();
368  ++it) {
369  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
370  ret.access_data().get_data()[i] = std::min(ret[i], it->first[i]);
371  }
372  }
373  return ret;
374  }
375 };
376 IMPALGEBRA_END_NAMESPACE
377 
378 #endif /* IMPALGEBRA_GRID_STORAGES_H */
Basic types used by IMP.
ParticleIndexes get_indexes(const ParticlesTemp &ps)
ExtendedGridIndexD< D > get_maximum_extended_index() const
Return the index that has no higher index in each coordinate.
A class to represent a voxel grid.
Represent a real cell in a grid (one within the bounding box)
Definition: grid_indexes.h:136
An index in an infinite grid on space.
Definition: grid_indexes.h:28
#define IMP_UNUSED(variable)
#define IMP_SHOWABLE_INLINE(Name, how_to_show)
Declare the methods needed by an object that can be printed.
Various general useful macros for IMP.
#define IMP_COPY_CONSTRUCTOR(Name, Base)
Use a copy_from method to create a copy constructor and operator=.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
bool get_has_index(const ExtendedGridIndexD< D > &i) const
Return true if the voxel has been added.
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
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.
#define IMP_IF_CHECK(level)
Execute the code block if a certain level checks are on.
Simple 3D vector class.
A class to represent a voxel grid.
#define IMP_FAILURE(message)
A runtime failure for IMP.
ExtendedGridIndexD< D > get_minimum_extended_index() const
Return the index which has no lower index in each coordinate.
#define IMP_BRACKET(Value, Index, bounds_check_expr, expr)
GridIndexD< D > add_voxel(const ExtendedGridIndexD< D > &i, const VT &gi)
Add a voxel to the storage, this voxel will now have a GridIndex3D.
Declare an efficient stl-compatible map.