IMP  2.1.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, class Storage,
76  // swig needs this for some reason
77  class Value, class EmbeddingT = DefaultEmbeddingD<D> >
78 class GridD : public Storage, public EmbeddingT, public GeometricPrimitiveD<D> {
79  private:
81 #ifndef IMP_DOXYGEN
82  protected:
83  struct GetVoxel {
84  mutable This *home_;
85  GetVoxel(This *home) : home_(home) {}
86  typedef Value &result_type;
87  typedef const GridIndexD<D> &argument_type;
88  result_type operator()(argument_type i) const {
89  // std::cout << i << std::endl;
90  return home_->operator[](i);
91  }
92  };
93 
94  struct ConstGetVoxel {
95  const This *home_;
96  ConstGetVoxel(const This *home) : home_(home) {}
97  typedef const Value &result_type;
98  typedef const GridIndexD<D> &argument_type;
99  result_type operator()(argument_type i) const {
100  // std::cout << i << std::endl;
101  return home_->operator[](i);
102  }
103  };
104 
105  VectorD<D> get_sides(const Ints &ns, const BoundingBoxD<D> &bb) const {
106  VectorD<D> ret = bb.get_corner(1);
107  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
108  ret[i] -= bb.get_corner(0)[i];
109  ret[i] /= ns[i];
110  }
111  return ret;
112  }
113  template <class NS>
114  Ints get_ns(const NS &ds, const BoundingBoxD<D> &bb) const {
115  Ints dd(ds.size());
116  for (unsigned int i = 0; i < ds.size(); ++i) {
117  IMP_USAGE_CHECK(ds[i] > 0,
118  "Number of voxels cannot be 0 on dimension: " << i);
119  double bside = bb.get_corner(1)[i] - bb.get_corner(0)[i];
120  double d = bside / ds[i];
121  double cd = std::ceil(d);
122  dd[i] = std::max<int>(1, static_cast<int>(cd));
123  }
124  return dd;
125  }
126 #endif
127  public:
128  typedef EmbeddingT Embedding;
129  typedef VectorD<D> Vector;
130  /** Create a grid from a bounding box and the counts in each direction.
131  */
132  GridD(const Ints counts, const BoundingBoxD<D> &bb,
133  Value default_value = Value())
134  : Storage(counts, default_value),
135  Embedding(bb.get_corner(0), get_sides(counts, bb)) {
136  IMP_USAGE_CHECK(D == 3, "Only in 3D");
137  }
138  /** Create a grid from a bounding box and the side of the cubical
139  voxel.
140  */
141  GridD(double side, const BoundingBoxD<D> &bb,
142  const Value &default_value = Value())
143  : Storage(get_ns(Floats(bb.get_dimension(), side), bb), default_value),
144  Embedding(bb.get_corner(0),
145  get_ones_vector_kd(bb.get_dimension(), side)) {
147  Storage::get_is_bounded(),
148  "This grid constructor can only be used with bounded grids.");
149  }
150 
151  GridD(const Storage &storage, const Embedding &embed)
152  : Storage(storage), Embedding(embed) {}
153  /** Construct a grid from the cubical voxel side and the origin.
154  */
155  GridD(double side, const VectorD<D> &origin,
156  const Value &default_value = Value())
157  : Storage(default_value),
158  Embedding(origin, get_ones_vector_kd(origin.get_dimension(), side)) {}
159  //! An empty, undefined grid.
160  GridD() : Storage(Value()) {}
161  /* \name Indexing
162  The vector must fall in a valid voxel to get and must be within
163  the volume of the grid to set.
164  @{
165  */
166  IMP_BRACKET(
167  Value, VectorD<D>,
168  Storage::get_has_index(Embedding::get_extended_index(i)),
169  return Storage::operator[](get_index(Embedding::get_extended_index(i))));
170 /** @} */
171 
172 #ifdef SWIG
173  const Value &__getitem__(const GridIndexD<D> &i) const;
174  void __setitem__(const GridIndexD<D> &i, const Value &vt);
175 #else
176  using Storage::__getitem__;
177  using Storage::__setitem__;
178  using Storage::operator[];
179 #endif
180  // ! Add a voxel to a sparse grid.
181  GridIndexD<D> add_voxel(const VectorD<D> &pt, const Value &vt) {
182  IMP_USAGE_CHECK(!Storage::get_is_dense(),
183  "add_voxel() only works on sparse grids.");
184  ExtendedGridIndexD<D> ei = Embedding::get_extended_index(pt);
185  return Storage::add_voxel(ei, vt);
186  }
187 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
188  Value &get_voxel_always(const VectorD<D> &pt) {
189  ExtendedGridIndexD<D> ei = Embedding::get_extended_index(pt);
190  return Storage::get_voxel_always(ei);
191  }
192  const Value &get_value_always(const VectorD<D> &pt) const {
193  ExtendedGridIndexD<D> ei = Embedding::get_extended_index(pt);
194  return Storage::get_value_always(ei);
195  }
196 #endif
197 #ifndef SWIG
198  using Storage::get_has_index;
199  using Storage::get_index;
200  using Storage::add_voxel;
201 #else
202  bool get_has_index(const ExtendedGridIndexD<D> &i) const;
203  GridIndexD<D> get_index(const ExtendedGridIndexD<D> &i) const;
204  GridIndexD<D> add_voxel(const ExtendedGridIndexD<D> &i, const Value &vt);
205 #endif
206  //! Convert an index back to an extended index
208  return ExtendedGridIndexD<D>(index.begin(), index.end());
209  }
210 #ifndef SWIG
211  using Embedding::get_extended_index;
212 #else
213  ExtendedGridIndexD<D> get_extended_index(const VectorD<D> &i) const;
214 #endif
215 
217  ExtendedGridIndexD<D> min = Storage::get_minimum_extended_index();
218  ExtendedGridIndexD<D> max = Storage::get_maximum_extended_index();
219  return get_bounding_box(min) + get_bounding_box(max);
220  }
221 #ifndef SWIG
223 #else
224  BoundingBoxD<D> get_bounding_box(const ExtendedGridIndexD<D> &i) const;
225  BoundingBoxD<D> get_bounding_box(const GridIndexD<D> &i) const;
226 #endif
227 
228  //! Change the bounding box but not the grid or contents
229  /** The origin is set to corner 0 of the new bounding box and the grid
230  voxels are resized as needed.
231  */
233  Floats nuc(bb3.get_dimension());
234  for (unsigned int i = 0; i < bb3.get_dimension(); ++i) {
235  double side = bb3.get_corner(1)[i] - bb3.get_corner(0)[i];
236  IMP_USAGE_CHECK(side > 0, "Can't have flat grid");
237  nuc[i] = side / Storage::get_number_of_voxels(i);
238  }
239  Embedding::set_unit_cell(VectorD<D>(nuc.begin(), nuc.end()));
240  Embedding::set_origin(bb3.get_corner(0));
241  }
242 
243  /** \name Get nearest
244  If the point is in the bounding box of the grid, this is the index
245  of the voxel containing the point,
246  otherwise it is the closest one in the bounding box. This can only be
247  used with bounded grids, right now.
248  @{
249  */
250  GridIndexD<D> get_nearest_index(const VectorD<D> &pt) const {
252  Storage::get_is_dense(), "get_nearest_index "
253  << "only works on dense grids.");
254  ExtendedGridIndexD<D> ei = get_nearest_extended_index(pt);
255  return get_index(ei);
256  }
257  ExtendedGridIndexD<D> get_nearest_extended_index(const VectorD<D> &pt) const {
259  Storage::get_is_bounded(), "get_nearest_index "
260  << "only works on bounded grids.");
261  ExtendedGridIndexD<D> ei = Embedding::get_extended_index(pt);
262  for (unsigned int i = 0; i < pt.get_dimension(); ++i) {
263  ei.access_data().get_data()[i] = std::max(0, ei[i]);
264  ei.access_data().get_data()[i] =
265  std::min<int>(Storage::get_number_of_voxels(i) - 1, ei[i]);
266  }
267  return ei;
268  }
269 /** @} */
270 
271 /** \name Voxel iterators
272 
273  These iterators go through a range of voxels in the grid. These voxels
274  include any that touch or are contained in the shape passed to the
275  begin/end calls.
276  @{
277  */
278 #ifndef SWIG
279 #ifdef IMP_DOXYGEN
280  class VoxelIterator;
281  class VoxelConstIterator;
282 #else
283  typedef boost::transform_iterator<GetVoxel, typename Storage::IndexIterator>
284  VoxelIterator;
285  typedef boost::transform_iterator<
286  ConstGetVoxel, typename Storage::IndexIterator> VoxelConstIterator;
287 #endif
288  VoxelIterator voxels_begin(const BoundingBoxD<D> &bb) {
289  return VoxelIterator(indexes_begin(bb), GetVoxel(this));
290  }
291  VoxelIterator voxels_end(const BoundingBoxD<D> &bb) {
292  // ExtendedIndex lb= get_extended_index(bb.get_corner(0));
293  // ExtendedIndex ub= get_extended_index(bb.get_corner(1));
294  return VoxelIterator(indexes_end(bb), GetVoxel(this));
295  }
296 
297  VoxelConstIterator voxels_begin(const BoundingBoxD<D> &bb) const {
298  return VoxelConstIterator(indexes_begin(bb), ConstGetVoxel(this));
299  }
300  VoxelConstIterator voxels_end(const BoundingBoxD<D> &bb) const {
301  return VoxelConstIterator(indexes_end(bb), ConstGetVoxel(this));
302  }
303  using Storage::indexes_begin;
304  using Storage::indexes_end;
305  typename Storage::IndexIterator indexes_begin(
306  const BoundingBoxD<D> &bb) const {
307  ExtendedGridIndexD<3> lb = get_extended_index(bb.get_corner(0));
308  ExtendedGridIndexD<3> ub = get_extended_index(bb.get_corner(1));
309  return Storage::indexes_begin(lb, ub);
310  }
311  typename Storage::IndexIterator indexes_end(const BoundingBoxD<D> &) const {
312  // ExtendedIndex lb= get_extended_index(bb.get_corner(0));
313  // ExtendedIndex ub= get_extended_index(bb.get_corner(1));
314  return Storage::indexes_end(ExtendedGridIndexD<3>(),
315  ExtendedGridIndexD<3>());
316  }
317 #endif
318  /** @} */
319  /** \name Apply
320  In C++, iterating through all the voxels can be slow, and it
321  can be faster to use functional programming to apply a
322  function to each voxel. The passed apply function takes three
323  arguments, the grid, the Grid::Index of the voxel and the
324  Grid::Vector for the center of the voxel.
325  @{ */
326  template <class Functor>
327  Functor apply(const Functor &f) const {
328  return Storage::apply(*this, f);
329  }
330  /** @} */
331 };
332 
333 template <int D, class Storage,
334  // swig needs this for some reason
335  class Value, class Embedding>
336 inline BoundingBoxD<D> get_bounding_box(
337  const GridD<D, Storage, Value, Embedding> &g) {
338  return g.get_bounding_box();
339 }
340 
341 IMPALGEBRA_END_NAMESPACE
342 
343 #endif /* IMPALGEBRA_GRID_D_H */
Basic types used by IMP.
Basic types used by IMP.
GridD()
An empty, undefined grid.
Definition: GridD.h:160
const VectorD< D > & get_corner(unsigned int i) const
For 0 return lower corner and 1 upper corner.
Definition: BoundingBoxD.h:123
GridD(double side, const VectorD< D > &origin, const Value &default_value=Value())
Definition: GridD.h:155
Ints get_index(const kernel::ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
A class to represent a voxel grid.
A class to represent a voxel grid.
Represent a real cell in a grid (one within the bounding box)
Definition: grid_indexes.h:136
GridD(double side, const BoundingBoxD< D > &bb, const Value &default_value=Value())
Definition: GridD.h:141
An index in an infinite grid on space.
Definition: grid_indexes.h:28
A voxel grid in d-dimensional space space.
Definition: GridD.h:78
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
VectorD< D > get_ones_vector_kd(unsigned int Di, double v=1)
Return a vector of ones (or another constant)
Definition: VectorD.h:481
void set_bounding_box(const BoundingBoxD< D > &bb3)
Change the bounding box but not the grid or contents.
Definition: GridD.h:232
A Cartesian vector in D-dimensions.
Definition: VectorD.h:48
ExtendedGridIndexD< D > get_extended_index(const GridIndexD< D > &index) const
Convert an index back to an extended index.
Definition: GridD.h:207
A bounding box in D dimensions.
A class for storing lists of IMP items.
An axis-aligned bounding box.
Definition: BoundingBoxD.h:26
Simple 3D vector class.
GridD(const Ints counts, const BoundingBoxD< D > &bb, Value default_value=Value())
Definition: GridD.h:132
#define IMP_BRACKET(Value, Index, bounds_check_expr, expr)
Declare an efficient stl-compatible map.
BoundingBoxD< D > get_bounding_box(const BoundingBoxD< D > &g)
Definition: BoundingBoxD.h:160