IMP  2.0.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= i[2]
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 "
65  << i);
66  return ii;
67  }
68 #ifndef SWIG
69  struct NonDefault {
70  VT default_;
71  NonDefault(const VT &def): default_(def){}
72  template <class P>
73  bool operator()(const P &def) const {
74  return def.second != default_;
75  }
76  };
77 #endif
78  unsigned int get_extent() const {
79  return extent_;
80  }
81  void copy_from(const DenseGridStorageD &o) {
82  default_= o.default_;
83  extent_= o.extent_;
84  data_.reset(new VT[extent_]);
85  std::copy(o.data_.get(), o.data_.get()+o.extent_,
86  data_.get());
88  }
89  void set_number_of_voxels(Ints dims) {
90  extent_=1;
91  for (unsigned int i=0; i < dims.size(); ++i) {
92  extent_*=dims[i];
93  }
94  data_.reset(new VT[extent_]);
95  std::fill(data_.get(), data_.get()+get_extent(), default_);
96  //BoundedGridRangeD<D>::set_number_of_voxels(dims);
97  }
98  public:
100  typedef VT Value;
101  DenseGridStorageD(const Ints &counts, const VT &default_value=VT()):
102  BoundedGridRangeD<D>(counts),
103  default_(default_value) {
104  set_number_of_voxels(counts);
105  }
106  IMP_BRACKET(VT, GridIndexD<D>, true, return data_[index(i)]);
107  /** \name Direct indexing
108  One can directly access a particular voxel based on its index
109  in the array of all voxels. This can be faster than using a
110  GridIndexD.
111  @{
112  */
113  IMP_BRACKET(VT, unsigned int, i<extent_, return data_[i]);
114  /** @}
115  */
116 #ifndef IMP_DOXYGEN
117  DenseGridStorageD(const VT &default_value=VT()):
118  extent_(0), default_(default_value) {
119  }
120  static bool get_is_dense() {
121  return true;
122  }
123 #endif
125 #ifndef SWIG
126  const VT* get_raw_data() const {
127  return data_.get();
128  }
129  VT* get_raw_data() {
130  return data_.get();
131  }
132 #endif
133 
134 #ifndef IMP_DOXYGEN
135  GridIndexD<D> add_voxel(const ExtendedGridIndexD<D>&, const VT&) {
136  IMP_FAILURE("Cannot add voxels to dense grid");
137  }
138 #endif
139 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
140  VT &get_voxel_always(const ExtendedGridIndexD<D> &i) {
141  GridIndexD<D> gi(i.begin(), i.end());
142  return operator[](gi);
143  }
144  const VT &get_value_always(const ExtendedGridIndexD<D> &i) const {
145  GridIndexD<D> gi(i.begin(), i.end());
146  return operator[](gi);
147  }
148 #endif
149 
150  /** \name All voxel iterators
151  The value type is VT.
152  @{
153  */
154 #ifndef SWIG
155  typedef VT* AllVoxelIterator;
156  typedef const VT* AllVoxelConstIterator;
157  AllVoxelIterator all_voxels_begin() {
158  return data_.get();
159  }
160  AllVoxelIterator all_voxels_end() {
161  return data_.get()+get_extent();
162  }
163  AllVoxelConstIterator all_voxels_begin() const {
164  return data_.get();
165  }
166  AllVoxelConstIterator all_voxels_end() const {
167  return data_.get()+get_extent();
168  }
169 #endif
170  base::Vector<VT> get_all_voxels() const {
171  return base::Vector<VT>(all_voxels_begin(),
172  all_voxels_end());
173  }
174  /** @} */
175  template <class Functor, class Grid>
176  Functor apply(const Grid &g, const Functor &fi) const {
177  Functor f=fi;
178  typename Grid::ExtendedIndex lb(Ints(g.get_dimension(),
179  0));
180  typename Grid::ExtendedIndex ub(BoundedGridRangeD<D>::get_end_index());
181  typename Grid::Vector corner= g.get_bounding_box().get_corner(0);
182  typename Grid::Vector cell= g.get_unit_cell();
183  typename Grid::Index cur;
184  typename Grid::Vector center;
185  internal::GridApplier<Functor, Grid, D-1>::apply(g, lb, ub,
186  corner, cell, cur,
187  center,f);
188  return f;
189  }
190 };
191 
192 
193 
194 
195 
196 
197 
198 
199 
200 
201 
202 /** Store a grid as a sparse set of voxels (only the voxels which have
203  been added are actually stored). The
204  get_has_index() functions allow one to tell if a voxel has been added.
205  \unstable{SparseGridStorageD}
206 
207  Base should be one of BoundedGridRangeD or UnboundedGridRangeD.
208  \see Grid3D
209 */
210 template <int D, class VT, class Base,
211  class Map=typename IMP::base::map<GridIndexD<D>, VT> >
212 class SparseGridStorageD: public Base {
213  typedef Map Data;
214  struct GetIndex {
215  typedef GridIndexD<D> result_type;
216  typedef typename Data::const_iterator::value_type argument_type;
217  template <class T>
218  GridIndexD<D> operator()(const T&t) const {
219  return t.first;
220  }
221  };
222  struct ItHelper {
223  const SparseGridStorageD<D, VT, Base> *stor_;
224  ItHelper(const SparseGridStorageD<D, VT, Base> *stor): stor_(stor){}
225  bool get_is_good(const ExtendedGridIndexD<D> &ei) {
226  /*std::cout << "Checking " << ei << " getting "
227  << stor_->get_has_index(ei) << std::endl;*/
228  return stor_->get_has_index(ei);
229  }
230  typedef GridIndexD<D> ReturnType;
231  ReturnType get_return(const ExtendedGridIndexD<D> &ei) const {
232  return stor_->get_index(ei);
233  }
234  ItHelper(): stor_(nullptr){}
235  };
236 
237  Data data_;
238  VT default_;
239  public:
240  typedef VT Value;
241  SparseGridStorageD(const Ints &counts,
242  const VT &default_value): Base(counts),
243  default_(default_value) {
244  }
245  IMP_SHOWABLE_INLINE(SparseGridStorage3D, out << "Sparse grid with "
246  << data_.size() << " cells set");
247  //! Add a voxel to the storage, this voxel will now have a GridIndex3D
249  IMP_USAGE_CHECK(Base::get_has_index(i), "Out of grid domain "
250  << i);
251  GridIndexD<D> ret(i.begin(), i.end());
252  data_[ret]=gi;
253  return ret;
254  }
255  void remove_voxel(const GridIndexD<D>& gi) {
256  data_.erase(gi);
257  }
258 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
259  SparseGridStorageD(const VT &def): default_(def) {
260  }
261  static bool get_is_dense() {
262  return false;
263  }
264  using Base::get_number_of_voxels;
265 #endif
266  unsigned int get_number_of_voxels() const {
267  return data_.size();
268  }
269  //! Return true if the voxel has been added
270  bool get_has_index(const ExtendedGridIndexD<D>&i) const {
271  return data_.find(GridIndexD<D>(i.begin(), i.end())) != data_.end();
272  }
273  //! requires get_has_index(i) is true.
275  IMP_USAGE_CHECK(get_has_index(i), "Index is not a valid "
276  << "voxel " << i);
277  return GridIndexD<D>(i.begin(), i.end());
278  }
279 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
280  VT &get_voxel_always(const ExtendedGridIndexD<D> &i) {
281  GridIndexD<D> gi(i.begin(), i.end());
282  typename Map::iterator it= data_.find(gi);
283  if (it == data_.end()) {
284  return data_.insert(std::make_pair(gi, default_)).first->second;
285  } else {
286  return it->second;
287  }
288  }
289  const VT &get_value_always(const ExtendedGridIndexD<D> &i) const {
290  GridIndexD<D> gi(i.begin(), i.end());
291  typename Map::const_iterator it= data_.find(gi);
292  if (it == data_.end()) {
293  return default_;
294  } else {
295  return it->second;
296  }
297  }
298 #endif
299  /** \name Operator []
300  Operator[] isn't very useful at the moment as it can only
301  be used with a cell which has already been set. This
302  behavior/the existence of these functions is likely to change.
303  @{
304  */
305  IMP_BRACKET(VT, GridIndexD<D>, true,
306  IMP_USAGE_CHECK(data_.find(i) != data_.end(),
307  "Invalid index " << i);
308  return data_.find(i)->second);
309  /** @} */
310 
311  /** \name Iterators through set cells
312  Iterate through the voxels which have been set. The value
313  type is a pair of GridIndex3D and VT.
314  @{
315  */
316 #ifndef SWIG
317  typedef typename Data::const_iterator AllConstIterator;
318  AllConstIterator all_begin() const {
319  return data_.begin();
320  }
321  AllConstIterator all_end() const {
322  return data_.end();
323  }
324 #endif
325 
326  base::Vector<GridIndexD<D> > get_all_indexes() const {
327  return base::Vector<GridIndexD<D> >
328  (boost::make_transform_iterator(all_begin(), GetIndex()),
329  boost::make_transform_iterator(all_end(), GetIndex()));
330  }
331  /** @} */
332 
333 
334 
335  /** \name Index Iterators
336 
337  Iterate through a range of actual indexes. The value
338  type for the iterator is an GridIndex3D.
339 
340  The range is defined by a pair of indexes. It includes
341  all indexes in the axis aligned box defined by lb
342  as the lower corner and the second as the ub. That is, if
343  lb is \f$(l_x, l_y, l_z)\f$ and ub is \f$(u_x, u_y, u_z)\f$,
344  then the range includes all
345  indexes \f$(i_x, i_y, i_z)\f$ such that \f$l_x \leq i_x \leq u_x\f$,
346  \f$l_y \leq i_y \leq u_y\f$
347  and \f$l_z \leq i_z \leq u_z\f$.
348 
349  @{
350  */
351 #ifndef SWIG
352 #ifndef IMP_DOXYGEN
353 
354  typedef internal::GridIndexIterator<ExtendedGridIndexD<D>,
355  ItHelper > IndexIterator;
356 
357 #else
358 class IndexIterator;
359 #endif
360 IndexIterator indexes_begin(const ExtendedGridIndexD<D>& lb,
361  const ExtendedGridIndexD<D>& ub) const {
362  ExtendedGridIndexD<D> eub=ub.get_offset(1,1,1);
363  if (lb == ub) {
364  return IndexIterator();
365  } else {
366  IMP_INTERNAL_CHECK(internal::get_is_non_empty(lb, eub),
367  "empty range");
368  return IndexIterator(lb, eub, ItHelper(this));
369  }
370 }
371 IndexIterator indexes_end(const ExtendedGridIndexD<D>&,
372  const ExtendedGridIndexD<D>&) const {
373  //IMP_INTERNAL_CHECK(lb <= ub, "empty range");
374  return IndexIterator();
375 }
376 #endif
377 
378 base::Vector<GridIndexD<D> >
379 get_indexes(const ExtendedGridIndexD<D>& lb,
380  const ExtendedGridIndexD<D>& ub) const {
381  return base::Vector<GridIndexD<D> >(indexes_begin(lb, ub),
382  indexes_end(lb, ub));
383 }
384 /** @} */
385 
386 
387 template <class Functor, class Grid>
388 Functor apply(const Grid &g, Functor f) const {
389  for (typename Data::const_iterator it= data_.begin();
390  it != data_.end(); ++it) {
391  f(g, it->first, g.get_center(it->first));
392  }
393  return f;
394 }
395 
396 
397 //! Return the index which has no lower index in each coordinate
399  IMP_USAGE_CHECK(!data_.empty(), "No voxels in grid.");
400  GridIndexD<D> reti=data_.begin()->first;
401  ExtendedGridIndexD<D> ret(reti.begin(), reti.end());
402  for (typename Data::const_iterator it= data_.begin();
403  it != data_.end(); ++it) {
404  for (unsigned int i=0; i< ret.get_dimension(); ++i) {
405  ret.access_data().get_data()[i]
406  = std::min(ret[i], it->first[i]);
407  }
408  }
409  return ret;
410 }
411 //! Return the index that has no higher index in each coordinate
413  IMP_USAGE_CHECK(!data_.empty(), "No voxels in grid.");
414  GridIndexD<D> reti=data_.begin()->first;
415  ExtendedGridIndexD<D> ret(reti.begin(), reti.end());
416  for (typename Data::const_iterator it= data_.begin();
417  it != data_.end(); ++it) {
418  for (unsigned int i=0; i< ret.get_dimension(); ++i) {
419  ret.access_data().get_data()[i]
420  = std::min(ret[i], it->first[i]);
421  }
422  }
423  return ret;
424 }
425 };
426 IMPALGEBRA_END_NAMESPACE
427 
428 
429 #endif /* IMPALGEBRA_GRID_STORAGES_H */