IMP  2.0.0
The Integrative Modeling Platform
GridD.h
Go to the documentation of this file.
1 /**
2  * \file IMP/algebra/GridD.h \brief A class to represent a voxel grid.
3  *
4  * Copyright 2007-2013 IMP Inventors. All rights reserved.
5  *
6  */
7 
8 #ifndef IMPALGEBRA_GRID_D_H
9 #define IMPALGEBRA_GRID_D_H
10 
11 #include <IMP/algebra/algebra_config.h>
12 
13 #include <IMP/base/types.h>
14 #include "grid_embeddings.h"
15 #include "grid_indexes.h"
16 #include "Vector3D.h"
17 #include "BoundingBoxD.h"
18 #include "GeometricPrimitiveD.h"
19 #include <boost/iterator/transform_iterator.hpp>
20 #include <IMP/base/map.h>
21 #include <IMP/base/Vector.h>
22 
23 #include <limits>
24 
25 IMPALGEBRA_BEGIN_NAMESPACE
26 
27 //! A voxel grid in d-dimensional space space.
28 /**
29  First some terminology:
30  - a voxel is the data stored at a given location is space
31  - an Index is a way of identifying a particular voxel. That is, given
32  an index, it is easy to get the voxel, but not vice-versa
33  - an ExtendedIndex identifies a particular region in space, but
34  it may not have a corresponding voxel (if it is outside of the
35  region the grid is built on or if that voxel has not yet been
36  added to the sparse grid).
37 
38  \imp provides support for a variety of spatial grids. The grid support in
39  C++ is implemented by combining several different layers to specify
40  what capabilities are desired. These layers are:
41  - Data: any type of data can be stored in a voxel of the grid
42  - Boundedness: By using UnboundedGridStorage3D or BoundedGridStorage3D,
43  one can choose whether you want a grid over a finite region of space
44  or over the whole space.
45  - Storage: by choosing SparseGridStorage3D or DenseGridStorage3D, you can
46  choose whether you want to store all voxels or only a subset of the
47  voxels. The former is faster and more compact when most of the voxels are
48  used, the latter when only a few are used (say <1/4).]
49  - Geometry: The Grid3D class itself provides a geometric layer, mapping
50  Vector3D objects into voxels in the grid.
51 
52  These are implemented as mix-ins, so each layer provides a set of accessible
53  functionality as methods/types in the final class.
54 
55  \par Basic operations
56  Creating a grid with a given cell size and upper and lower
57  bounds
58  \code
59  BoundingBox3D bb(Vector3D(10,10,10), Vector3D(100,100,100));
60  typedef Grid3D<Ints> Grid;
61  Grid grid(5, bb, 0.0);
62  \endcode
63 
64  Iterate over the set of voxels incident on a bounding box:
65  \code
66  BoundingBoxD<3> bb(Vector3D(20.2,20.3,20.5), Vector3D(31.3,32.5,38.9));
67  for (Grid::IndexIterator it= grid.voxels_begin(bb);
68  it != grid.voxels_end(bb); ++it) {
69  it->push_back(1);
70  }
71  \endcode
72  \see DenseGridStorage3D
73  \see SparseGridStorageD
74 */
75 template <int D,
76  class Storage,
77  // swig needs this for some reason
78  class Value,
79  class EmbeddingT=DefaultEmbeddingD<D> >
80 class GridD: public Storage, public EmbeddingT,
81  public GeometricPrimitiveD<D>
82 {
83  private:
85 #ifndef IMP_DOXYGEN
86  protected:
87  struct GetVoxel {
88  mutable This *home_;
89  GetVoxel(This *home): home_(home) {}
90  typedef Value& result_type;
91  typedef const GridIndexD<D>& argument_type;
92  result_type operator()(argument_type i) const {
93  //std::cout << i << std::endl;
94  return home_->operator[](i);
95  }
96  };
97 
98  struct ConstGetVoxel {
99  const This *home_;
100  ConstGetVoxel(const This *home): home_(home) {}
101  typedef const Value& result_type;
102  typedef const GridIndexD<D>& argument_type;
103  result_type operator()(argument_type i) const {
104  //std::cout << i << std::endl;
105  return home_->operator[](i);
106  }
107  };
108 
109  Floats get_sides(const Ints &ns,
110  const BoundingBoxD<D> &bb) const {
111  Floats ret(bb.get_dimension());
112  for (unsigned int i=0; i< ret.size(); ++i) {
113  ret[i]= (bb.get_corner(1)[i]-bb.get_corner(0)[i])/ns[i];
114  }
115  return ret;
116  }
117  template <class NS>
118  Ints get_ns(const NS &ds,
119  const BoundingBoxD<D> &bb) const {
120  Ints dd(ds.size());
121  for (unsigned int i=0; i< ds.size(); ++i ) {
122  IMP_USAGE_CHECK(ds[i]>0, "Number of voxels cannot be 0 on dimension: "
123  << i);
124  double bside= bb.get_corner(1)[i]- bb.get_corner(0)[i];
125  double d= bside/ds[i];
126  double cd= std::ceil(d);
127  dd[i]= std::max<int>(1, static_cast<int>(cd));
128  }
129  return dd;
130  }
131 #endif
132  public:
133  typedef EmbeddingT Embedding;
134  typedef VectorD<D> Vector;
135  /** Create a grid from a bounding box and the counts in each direction.
136  */
137  GridD(const Ints counts,
138  const BoundingBoxD<D> &bb,
139  Value default_value=Value()):
140  Storage(counts, default_value),
141  Embedding(bb.get_corner(0), get_sides(counts, bb)) {
142  IMP_USAGE_CHECK(D==3, "Only in 3D");
143  }
144  /** Create a grid from a bounding box and the side of the cubical
145  voxel.
146  */
147  GridD(double side,
148  const BoundingBoxD<D> &bb,
149  const Value& default_value=Value()):
150  Storage(get_ns(Floats(bb.get_dimension(), side), bb), default_value),
151  Embedding(bb.get_corner(0),
152  VectorD<D>(Floats(bb.get_dimension(), side))){
153  IMP_USAGE_CHECK(Storage::get_is_bounded(),
154  "This grid constructor can only be used with bounded grids.");
155  }
156 
157  GridD(const Storage&storage, const Embedding &embed): Storage(storage),
158  Embedding(embed){
159  }
160  /** Construct a grid from the cubical voxel side and the origin.
161  */
162  GridD(double side,
163  const VectorD<D> &origin,
164  const Value& default_value= Value()):
165  Storage(default_value), Embedding(origin,
166  VectorD<D>(Floats(origin.get_dimension(),
167  side))){
168  }
169  //! An empty, undefined grid.
170  GridD(): Storage(Value()){
171  }
172  /* \name Indexing
173  The vector must fall in a valid voxel to get and must be within
174  the volume of the grid to set.
175  @{
176  */
177  IMP_BRACKET(Value, VectorD<D>,
178  Storage::get_has_index(Embedding::get_extended_index(i)),
179  return Storage::operator[](get_index(Embedding
180  ::get_extended_index(i))));
181  /** @} */
182 
183 #ifdef SWIG
184  const Value& __getitem__(const GridIndexD<D> &i) const;
185  void __setitem__(const GridIndexD<D> &i, const Value &vt);
186 #else
187  using Storage::__getitem__;
188  using Storage::__setitem__;
189  using Storage::operator[];
190 #endif
191  // ! Add a voxel to a sparse grid.
192  GridIndexD<D> add_voxel(const VectorD<D>& pt, const Value &vt) {
193  IMP_USAGE_CHECK(!Storage::get_is_dense(),
194  "add_voxel() only works on sparse grids.");
195  ExtendedGridIndexD<D> ei= Embedding::get_extended_index(pt);
196  return Storage::add_voxel(ei, vt);
197  }
198 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
199  Value &get_voxel_always(const VectorD<D>& pt) {
200  ExtendedGridIndexD<D> ei= Embedding::get_extended_index(pt);
201  return Storage::get_voxel_always(ei);
202  }
203  const Value &
204  get_value_always(const VectorD<D>& pt) const {
205  ExtendedGridIndexD<D> ei= Embedding::get_extended_index(pt);
206  return Storage::get_value_always(ei);
207  }
208 #endif
209 #ifndef SWIG
210  using Storage::get_has_index;
211  using Storage::get_index;
212  using Storage::add_voxel;
213 #else
214  bool get_has_index(const ExtendedGridIndexD<D>&i) const;
215  GridIndexD<D> get_index(const ExtendedGridIndexD<D> &i) const;
216  GridIndexD<D> add_voxel(const ExtendedGridIndexD<D> &i,
217  const Value &vt);
218 #endif
219  //! Convert an index back to an extended index
221  return ExtendedGridIndexD<D>(index.begin(), index.end());
222  }
223 #ifndef SWIG
224  using Embedding::get_extended_index;
225 #else
226  ExtendedGridIndexD<D> get_extended_index(const VectorD<D> &i) const;
227 #endif
228 
229  BoundingBoxD<D> get_bounding_box() const {
230  ExtendedGridIndexD<D> min= Storage::get_minimum_extended_index();
231  ExtendedGridIndexD<D> max= Storage::get_maximum_extended_index();
232  return get_bounding_box(min)+get_bounding_box(max);
233  }
234 #ifndef SWIG
235  using Embedding::get_bounding_box;
236 #else
237  BoundingBoxD<D> get_bounding_box(const ExtendedGridIndexD<D> &i) const;
238  BoundingBoxD<D> get_bounding_box(const GridIndexD<D> &i) const;
239 #endif
240 
241  //! Change the bounding box but not the grid or contents
242  /** The origin is set to corner 0 of the new bounding box and the grid
243  voxels are resized as needed.
244  */
246  Floats nuc(bb3.get_dimension());
247  for (unsigned int i=0; i< bb3.get_dimension(); ++i) {
248  double side= bb3.get_corner(1)[i]- bb3.get_corner(0)[i];
249  IMP_USAGE_CHECK(side>0, "Can't have flat grid");
250  nuc[i]= side/Storage::get_number_of_voxels(i);
251  }
252  Embedding::set_unit_cell(VectorD<D>(nuc.begin(), nuc.end()));
253  Embedding::set_origin(bb3.get_corner(0));
254  }
255 
256  /** \name Get nearest
257  If the point is in the bounding box of the grid, this is the index
258  of the voxel containing the point,
259  otherwise it is the closest one in the bounding box. This can only be
260  used with bounded grids, right now.
261  @{
262  */
263  GridIndexD<D> get_nearest_index(const VectorD<D>& pt) const {
264  IMP_USAGE_CHECK(Storage::get_is_dense(), "get_nearest_index "
265  << "only works on dense grids.");
266  ExtendedGridIndexD<D> ei= get_nearest_extended_index(pt);
267  return get_index(ei);
268  }
269  ExtendedGridIndexD<D>
270  get_nearest_extended_index(const VectorD<D>& pt) const {
271  IMP_USAGE_CHECK(Storage::get_is_bounded(), "get_nearest_index "
272  << "only works on bounded grids.");
273  ExtendedGridIndexD<D> ei= Embedding::get_extended_index(pt);
274  boost::scoped_array<int> is(new int[pt.get_dimension()]);
275  for (unsigned int i=0; i< pt.get_dimension(); ++i) {
276  is[i]= std::max(0, ei[i]);
277  is[i]= std::min<int>(Storage::get_number_of_voxels(i)-1, is[i]);
278  }
279  return ExtendedGridIndexD<D>(is.get(), is.get()+pt.get_dimension());
280  }
281  /** @} */
282 
283  /** \name Voxel iterators
284 
285  These iterators go through a range of voxels in the grid. These voxels
286  include any that touch or are contained in the shape passed to the
287  begin/end calls.
288  @{
289  */
290 #ifndef SWIG
291 #ifdef IMP_DOXYGEN
292  class VoxelIterator;
293  class VoxelConstIterator;
294 #else
295  typedef boost::transform_iterator<GetVoxel, typename Storage::IndexIterator>
296  VoxelIterator;
297  typedef boost::transform_iterator<ConstGetVoxel,
298  typename Storage::IndexIterator>
299  VoxelConstIterator;
300 #endif
301  VoxelIterator voxels_begin(const BoundingBoxD<D> &bb) {
302  return VoxelIterator(indexes_begin(bb), GetVoxel(this));
303  }
304  VoxelIterator voxels_end(const BoundingBoxD<D> &bb) {
305  //ExtendedIndex lb= get_extended_index(bb.get_corner(0));
306  //ExtendedIndex ub= get_extended_index(bb.get_corner(1));
307  return VoxelIterator(indexes_end(bb),
308  GetVoxel(this));
309  }
310 
311  VoxelConstIterator voxels_begin(const BoundingBoxD<D> &bb) const {
312  return VoxelConstIterator(indexes_begin(bb),
313  ConstGetVoxel(this));
314  }
315  VoxelConstIterator voxels_end(const BoundingBoxD<D> &bb) const {
316  return VoxelConstIterator(indexes_end(bb),
317  ConstGetVoxel(this));
318  }
319  using Storage::indexes_begin;
320  using Storage::indexes_end;
321  typename Storage::IndexIterator
322  indexes_begin(const BoundingBoxD<D> &bb) const {
323  ExtendedGridIndexD<3> lb= get_extended_index(bb.get_corner(0));
324  ExtendedGridIndexD<3> ub= get_extended_index(bb.get_corner(1));
325  return Storage::indexes_begin(lb, ub);
326  }
327  typename Storage::IndexIterator indexes_end(const BoundingBoxD<D> &) const {
328  //ExtendedIndex lb= get_extended_index(bb.get_corner(0));
329  //ExtendedIndex ub= get_extended_index(bb.get_corner(1));
330  return Storage::indexes_end(ExtendedGridIndexD<3>(),
331  ExtendedGridIndexD<3>());
332  }
333 #endif
334  /** @} */
335  /** \name Apply
336  In C++, iterating through all the voxels can be slow, and it
337  can be faster to use functional programming to apply a
338  function to each voxel. The passed apply function takes three
339  arguments, the grid, the Grid::Index of the voxel and the
340  Grid::Vector for the center of the voxel.
341  @{ */
342  template <class Functor>
343  Functor apply(const Functor &f) const {
344  return Storage::apply(*this, f);
345  }
346  /** @} */
347 };
348 
349 
350 
351 template <int D,
352  class Storage,
353  // swig needs this for some reason
354  class Value,
355  class Embedding>
356 inline BoundingBoxD<D> get_bounding_box(const
357  GridD<D, Storage, Value,
358  Embedding> &g) {
359  return g.get_bounding_box();
360 }
361 
362 IMPALGEBRA_END_NAMESPACE
363 
364 
365 #endif /* IMPALGEBRA_GRID_D_H */