IMP logo
IMP Reference Guide  develop.29b24c6ad6,2021/05/15
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-2021 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/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/Vector.h>
21 #include <IMP/showable_macros.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.
29 /**
30  First some terminology:
31  - a voxel is the data stored at a given location in 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(bb.get_dimension());
119  for (unsigned int i = 0; i < bb.get_dimension(); ++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 
145  /** Create a grid from a bounding box and the side of the cubical
146  voxel.
147 
148  This constructor works for all bounded grids.
149  */
150  GridD(double side, const BoundingBoxD<D> &bb,
151  const Value &default_value = Value())
152  : Storage(get_ns(Floats(bb.get_dimension(), side), bb), default_value),
153  Embedding(bb.get_corner(0),
154  get_ones_vector_kd(bb.get_dimension(), side)) {
156  Storage::get_is_bounded(),
157  "This grid constructor can only be used with bounded grids.");
158  }
159 
160  /** Create a grid from a bounding box and the sides of the variable
161  size voxel.
162 
163  This constructor works for all bounded grids.
164  */
165  GridD(const VectorD<D> &sides, const BoundingBoxD<D> &bb,
166  const Value &default_value = Value())
167  : Storage(get_ns(sides.get_coordinates(), bb), default_value),
168  Embedding(bb.get_corner(0), sides) {
170  Storage::get_is_bounded(),
171  "This grid constructor can only be used with bounded grids.");
172  }
173 
174  /** Advanced constructor if you want to explicitly init storage
175  and embedding. */
176  GridD(const Storage &storage, const Embedding &embed)
177  : Storage(storage), Embedding(embed) {}
178 
179  /** Construct a grid from the cubical voxel side and the origin.
180 
181  This constructor is only valid with unbounded (hence sparse)
182  grids.
183  */
184  GridD(double side, const VectorD<D> &origin,
185  const Value &default_value = Value())
186  : Storage(default_value),
187  Embedding(origin, get_ones_vector_kd(origin.get_dimension(), side)) {}
188 
189  /** Construct a grid from the variable voxel side and the origin.
190 
191  This constructor is only valid with unbounded (hence sparse)
192  grids.
193  */
194  GridD(const VectorD<D> &sides, const VectorD<D> &origin,
195  const Value &default_value = Value())
196  : Storage(default_value),
197  Embedding(origin, sides) {}
198 
199  //! An empty, undefined grid.
200  /** Make sure you initialize it before you try to use it. */
201  GridD() : Storage(Value()) {}
202 
203  IMP_SHOWABLE_INLINE(GridD, out << "Grid");
204 
205  /* \name Indexing
206  The vector must fall in a valid voxel to get and must be within
207  the volume of the grid to set.
208  @{
209  */
210  IMP_BRACKET(
211  Value, VectorD<D>,
212  Storage::get_has_index(Embedding::get_extended_index(i)),
213  return Storage::operator[](get_index(Embedding::get_extended_index(i))));
214 /** @} */
215 
216 #ifdef SWIG
217  const Value &__getitem__(const GridIndexD<D> &i) const;
218  void __setitem__(const GridIndexD<D> &i, const Value &vt);
219 #else
220  using Storage::__getitem__;
221  using Storage::__setitem__;
222  using Storage::operator[];
223 #endif
224  // ! Add a voxel to a sparse grid.
225  GridIndexD<D> add_voxel(const VectorD<D> &pt, const Value &vt) {
226  IMP_USAGE_CHECK(!Storage::get_is_dense(),
227  "add_voxel() only works on sparse grids.");
228  ExtendedGridIndexD<D> ei = Embedding::get_extended_index(pt);
229  return Storage::add_voxel(ei, vt);
230  }
231 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
232  Value &get_voxel_always(const VectorD<D> &pt) {
233  ExtendedGridIndexD<D> ei = Embedding::get_extended_index(pt);
234  return Storage::get_voxel_always(ei);
235  }
236  const Value &get_value_always(const VectorD<D> &pt) const {
237  ExtendedGridIndexD<D> ei = Embedding::get_extended_index(pt);
238  return Storage::get_value_always(ei);
239  }
240 #endif
241 #ifndef SWIG
242  using Storage::get_has_index;
243  using Storage::get_index;
244  using Storage::add_voxel;
245 #endif
246  //! Convert an index back to an extended index
248  return ExtendedGridIndexD<D>(index.begin(), index.end());
249  }
250 #ifndef SWIG
251  using Embedding::get_extended_index;
252 #endif
253 
254  BoundingBoxD<D> get_bounding_box() const {
255  ExtendedGridIndexD<D> min = Storage::get_minimum_extended_index();
256  ExtendedGridIndexD<D> max = Storage::get_maximum_extended_index();
257  return get_bounding_box(min) + get_bounding_box(max);
258  }
259 #ifndef SWIG
260  using Embedding::get_bounding_box;
261 #else
262  BoundingBoxD<D> get_bounding_box(const ExtendedGridIndexD<D> &i) const;
263  BoundingBoxD<D> get_bounding_box(const GridIndexD<D> &i) const;
264 #endif
265 
266  //! Change the bounding box but not the grid or contents
267  /** The origin is set to corner 0 of the new bounding box and the grid
268  voxels are resized as needed.
269  */
271  Floats nuc(bb3.get_dimension());
272  for (unsigned int i = 0; i < bb3.get_dimension(); ++i) {
273  double side = bb3.get_corner(1)[i] - bb3.get_corner(0)[i];
274  IMP_USAGE_CHECK(side > 0, "Can't have flat grid");
275  nuc[i] = side / Storage::get_number_of_voxels(i);
276  }
277  Embedding::set_unit_cell(VectorD<D>(nuc.begin(), nuc.end()));
278  Embedding::set_origin(bb3.get_corner(0));
279  }
280 
281  /** \name Get nearest
282  If the point is in the bounding box of the grid, this is the index
283  of the voxel containing the point,
284  otherwise it is the closest one in the bounding box. This can only be
285  used with bounded grids, right now.
286  @{
287  */
288  GridIndexD<D> get_nearest_index(const VectorD<D> &pt) const {
290  Storage::get_is_dense(), "get_nearest_index "
291  << "only works on dense grids.");
292  ExtendedGridIndexD<D> ei = get_nearest_extended_index(pt);
293  return get_index(ei);
294  }
295  ExtendedGridIndexD<D> get_nearest_extended_index(const VectorD<D> &pt) const {
297  Storage::get_is_bounded(), "get_nearest_index "
298  << "only works on bounded grids.");
299  ExtendedGridIndexD<D> ei = Embedding::get_extended_index(pt);
300  for (unsigned int i = 0; i < pt.get_dimension(); ++i) {
301  ei.access_data().get_data()[i] = std::max(0, ei[i]);
302  ei.access_data().get_data()[i] =
303  std::min<int>(Storage::get_number_of_voxels(i) - 1, ei[i]);
304  }
305  return ei;
306  }
307 /** @} */
308 
309 /** \name Voxel iterators
310 
311  These iterators go through a range of voxels in the grid. These voxels
312  include any that touch or are contained in the shape passed to the
313  begin/end calls.
314  @{
315  */
316 #ifndef SWIG
317 #ifdef IMP_DOXYGEN
318  class VoxelIterator;
319  class VoxelConstIterator;
320 #else
321  typedef boost::transform_iterator<GetVoxel, typename Storage::IndexIterator>
322  VoxelIterator;
323  typedef boost::transform_iterator<
324  ConstGetVoxel, typename Storage::IndexIterator> VoxelConstIterator;
325 #endif
326  VoxelIterator voxels_begin(const BoundingBoxD<D> &bb) {
327  return VoxelIterator(indexes_begin(bb), GetVoxel(this));
328  }
329  VoxelIterator voxels_end(const BoundingBoxD<D> &bb) {
330  // ExtendedIndex lb= get_extended_index(bb.get_corner(0));
331  // ExtendedIndex ub= get_extended_index(bb.get_corner(1));
332  return VoxelIterator(indexes_end(bb), GetVoxel(this));
333  }
334 
335  VoxelConstIterator voxels_begin(const BoundingBoxD<D> &bb) const {
336  return VoxelConstIterator(indexes_begin(bb), ConstGetVoxel(this));
337  }
338  VoxelConstIterator voxels_end(const BoundingBoxD<D> &bb) const {
339  return VoxelConstIterator(indexes_end(bb), ConstGetVoxel(this));
340  }
341  using Storage::indexes_begin;
342  using Storage::indexes_end;
343  typename Storage::IndexIterator indexes_begin(const BoundingBoxD<D> &bb)
344  const {
345  ExtendedGridIndexD<D> lb = get_extended_index(bb.get_corner(0));
346  ExtendedGridIndexD<D> ub = get_extended_index(bb.get_corner(1));
347  return Storage::indexes_begin(lb, ub);
348  }
349  typename Storage::IndexIterator indexes_end(const BoundingBoxD<D> &) const {
350  // ExtendedIndex lb= get_extended_index(bb.get_corner(0));
351  // ExtendedIndex ub= get_extended_index(bb.get_corner(1));
352  return Storage::indexes_end(ExtendedGridIndexD<D>(),
353  ExtendedGridIndexD<D>());
354  }
355 
356  typedef internal::GridIndexIterator<
357  ExtendedGridIndexD<D>,
358  internal::AllItHelp<ExtendedGridIndexD<D>, ExtendedGridIndexD<D> > >
359  ExtendedIndexIterator;
360  ExtendedIndexIterator extended_indexes_begin(const BoundingBoxD<D> &bb)
361  const {
362  ExtendedGridIndexD<D> lb = get_extended_index(bb.get_corner(0));
363  ExtendedGridIndexD<D> ub = get_extended_index(bb.get_corner(1));
364  ExtendedGridIndexD<D> eub = ub.get_offset(1, 1, 1);
365  return ExtendedIndexIterator(lb, eub);
366  }
367  ExtendedIndexIterator extended_indexes_end(const BoundingBoxD<D> &) const {
368  // ExtendedIndex lb= get_extended_index(bb.get_corner(0));
369  // ExtendedIndex ub= get_extended_index(bb.get_corner(1));
370  return ExtendedIndexIterator();
371  }
372  using Storage::get_indexes;
373  typedef boost::iterator_range<typename Storage::IndexIterator> Indexes;
374  Indexes get_indexes(const BoundingBoxD<D> &bb) const {
375  return Indexes(indexes_begin(bb), indexes_end(bb));
376  }
377  typedef boost::iterator_range<typename Storage::AllIndexIterator> AllIndexes;
378  AllIndexes get_all_indexes() const {
379  return AllIndexes(Storage::all_indexes_begin(), Storage::all_indexes_end());
380  }
381  using Storage::get_extended_indexes;
382  typedef boost::iterator_range<ExtendedIndexIterator> ExtendedIndexes;
383  ExtendedIndexes get_extended_indexes(const BoundingBoxD<D> &bb) const {
384  return ExtendedIndexes(extended_indexes_begin(bb),
385  extended_indexes_end(bb));
386  }
387 #endif
388  /** @} */
389  /** \name Apply
390  In C++, iterating through all the voxels can be slow, and it
391  can be faster to use functional programming to apply a
392  function to each voxel. The passed apply function takes three
393  arguments, the grid, the Grid::Index of the voxel and the
394  Grid::Vector for the center of the voxel.
395  @{ */
396  template <class Functor>
397  Functor apply(const Functor &f) const {
398  return Storage::apply(*this, f);
399  }
400  /** @} */
401 };
402 
403 template <int D, class Storage,
404  // swig needs this for some reason
405  class Value, class Embedding>
406 inline BoundingBoxD<D> get_bounding_box(
407  const GridD<D, Storage, Value, Embedding> &g) {
408  return g.get_bounding_box();
409 }
410 
411 IMPALGEBRA_END_NAMESPACE
412 
413 #endif /* IMPALGEBRA_GRID_D_H */
Base class for geometric types.
Basic types used by IMP.
#define IMP_SHOWABLE_INLINE(Name, how_to_show)
Declare the methods needed by an object that can be printed.
const VectorD< D > & get_corner(unsigned int i) const
For 0 return lower corner and for 1, the upper corner.
Definition: BoundingBoxD.h:139
A class to represent a voxel grid.
GridD(double side, const VectorD< D > &origin, const Value &default_value=Value())
Definition: GridD.h:184
GridD(const VectorD< D > &sides, const VectorD< D > &origin, const Value &default_value=Value())
Definition: GridD.h:194
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.
Definition: GridD.h:79
GridD(const VectorD< D > &sides, const BoundingBoxD< D > &bb, const Value &default_value=Value())
Definition: GridD.h:165
#define IMP_BRACKET(Value, Index, bounds_check_expr, expr)
Implement operator[] and at() for C++, and getitem for Python.
VectorD< D > get_ones_vector_kd(unsigned int Di, double v=1)
Return a vector of ones (or another constant)
Definition: VectorD.h:309
GridD(const Storage &storage, const Embedding &embed)
Definition: GridD.h:176
Base class for a simple primitive-like type.
Definition: Value.h:23
Base class for geometric types.
A Cartesian vector in D-dimensions.
Definition: VectorD.h:52
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
GridD(double side, const BoundingBoxD< D > &bb, const Value &default_value=Value())
Definition: GridD.h:150
BoundingBoxD< 3 > get_bounding_box(const Cone3D &g)
Definition: Cone3D.h:64
A bounding box in D dimensions.
ExtendedGridIndexD< D > get_extended_index(const GridIndexD< D > &index) const
Convert an index back to an extended index.
Definition: GridD.h:247
GridD(const Ints counts, const BoundingBoxD< D > &bb, Value default_value=Value())
Definition: GridD.h:138
A class for storing lists of IMP items.
void set_bounding_box(const BoundingBoxD< D > &bb3)
Change the bounding box but not the grid or contents.
Definition: GridD.h:270
An axis-aligned bounding box.
Definition: BoundingBoxD.h:27
Simple 3D vector class.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
ParticleIndexes get_indexes(const ParticlesTemp &ps)
Get the indexes from a list of particles.
GridD()
An empty, undefined grid.
Definition: GridD.h:201
Macros to help with objects that can be printed to a stream.