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