IMP logo
IMP Reference Guide  2.19.0
The Integrative Modeling Platform
grid_embeddings.h
Go to the documentation of this file.
1 /**
2  * \file IMP/algebra/grid_embeddings.h
3  * \brief A class to represent a voxel grid.
4  *
5  * Copyright 2007-2022 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPALGEBRA_GRID_EMBEDDINGS_H
10 #define IMPALGEBRA_GRID_EMBEDDINGS_H
11 
12 #include <IMP/algebra/algebra_config.h>
13 
14 #include <IMP/types.h>
15 #include "grid_indexes.h"
16 #include "Vector3D.h"
17 #include "BoundingBoxD.h"
18 #include <boost/iterator/transform_iterator.hpp>
19 #include <cereal/access.hpp>
20 #include <IMP/Vector.h>
21 #include <IMP/check_macros.h>
22 #include <IMP/exception.h>
23 
24 #include <limits>
25 
26 IMPALGEBRA_BEGIN_NAMESPACE
27 
28 //! Embed a grid as an evenly spaced axis aligned grid.
29 template <int D>
31  VectorD<D> origin_;
32  VectorD<D> unit_cell_;
33  // inverse
34  VectorD<D> inverse_unit_cell_;
35 
36  friend class cereal::access;
37 
38  template<class Archive> void save(Archive &ar) const {
39  ar(origin_, unit_cell_);
40  }
41 
42  template<class Archive> void load(Archive &ar) {
43  ar(origin_, unit_cell_);
44  set_unit_cell(unit_cell_); // sets inverse
45  }
46 
47  template <class O>
48  VectorD<D> get_elementwise_product(VectorD<D> v0, const O &v1) const {
49  for (unsigned int i = 0; i < get_dimension(); ++i) {
50  v0[i] *= v1[i];
51  }
52  return v0;
53  }
54  template <class O>
55  VectorD<D> get_uniform_offset(const O &v0, double o) const {
56  Floats ret(get_dimension());
57  for (unsigned int i = 0; i < get_dimension(); ++i) {
58  ret[i] = v0[i] + o;
59  }
60  return VectorD<D>(ret.begin(), ret.end());
61  }
62  void initialize_from_box(Ints ns, const BoundingBoxD<D> &bb) {
63  Floats nuc(bb.get_dimension());
64  for (unsigned int i = 0; i < bb.get_dimension(); ++i) {
65  double side = bb.get_corner(1)[i] - bb.get_corner(0)[i];
66  IMP_USAGE_CHECK(side > 0, "Can't have flat grid");
67  nuc[i] = side / ns[i];
68  }
69  set_unit_cell(VectorD<D>(nuc.begin(), nuc.end()));
70  set_origin(bb.get_corner(0));
71  }
72 
73  public:
74  DefaultEmbeddingD(const VectorD<D> &origin, const VectorD<D> &cell) {
75  set_origin(origin);
76  set_unit_cell(cell);
77  }
79  void set_origin(const VectorD<D> &o) { origin_ = o; }
80  const VectorD<D> get_origin() const { return origin_; }
81  unsigned int get_dimension() const { return get_origin().get_dimension(); }
82  void set_unit_cell(const VectorD<D> &o) {
83  unit_cell_ = o;
84  Floats iuc(o.get_dimension());
85  for (unsigned int i = 0; i < get_dimension(); ++i) {
86  iuc[i] = 1.0 / unit_cell_[i];
87  }
88  inverse_unit_cell_ = VectorD<D>(iuc.begin(), iuc.end());
89  }
90 #ifndef IMP_DOXYGEN
91  //! Return the vector (1/u[0], 1/u[1], 1/u[2])
92  const VectorD<D> &get_inverse_unit_cell() const { return inverse_unit_cell_; }
93 #endif
94  //! Return the unit cell, relative to the origin.
95  /** That is, the unit cell is
96  \code
97  BoundingBoxD<D>(get_zeros_vector_d<D>(),get_unit_cell());
98  \endcode
99  */
100  const VectorD<D> &get_unit_cell() const { return unit_cell_; }
101  //! Return the index that would contain the voxel if the grid extended there
102  /** For example vectors below the "lower left" corner of the grid have
103  indexes with all negative components. This operation will always
104  succeed.
105  */
108  origin_.get_dimension());
109  for (unsigned int i = 0; i < get_dimension(); ++i) {
110  double d = o[i] - origin_[i];
111  double fi = d * inverse_unit_cell_[i];
112  ret.access_data().get_data()[i] = static_cast<int>(std::floor(fi));
113  }
114  return ret;
115  }
116  GridIndexD<D> get_index(const VectorD<D> &o) const {
118  origin_.get_dimension());
119  for (unsigned int i = 0; i < get_dimension(); ++i) {
120  double d = o[i] - origin_[i];
121  double fi = d * inverse_unit_cell_[i];
122  ret.access_data().get_data()[i] = static_cast<int>(std::floor(fi));
123  }
124  return ret;
125  }
126  /** \name Center
127  Return the coordinates of the center of the voxel.
128  @{
129  */
130  VectorD<D> get_center(const ExtendedGridIndexD<D> &ei) const {
131  return origin_ +
132  get_elementwise_product(get_unit_cell(), get_uniform_offset(ei, .5));
133  }
134  VectorD<D> get_center(const GridIndexD<D> &ei) const {
135  return origin_ +
136  get_elementwise_product(get_unit_cell(), get_uniform_offset(ei, .5));
137  }
138  /** @} */
139 
140  /** \name Bounding box
141  Return the bounding box of the voxel.
142  @{
143  */
144  BoundingBoxD<D> get_bounding_box(const ExtendedGridIndexD<D> &ei) const {
145  return BoundingBoxD<D>(
146  origin_ + get_elementwise_product(unit_cell_, ei),
147  origin_ +
148  get_elementwise_product(unit_cell_, get_uniform_offset(ei, 1)));
149  }
150  BoundingBoxD<D> get_bounding_box(const GridIndexD<D> &ei) const {
151  return BoundingBoxD<D>(
152  origin_ + get_elementwise_product(unit_cell_, ei),
153  origin_ +
154  get_elementwise_product(unit_cell_, get_uniform_offset(ei, 1)));
155  }
156  /** @} */
157  IMP_SHOWABLE_INLINE(DefaultEmbeddingD, out << "origin: " << origin_
158  << " unit cell: " << unit_cell_);
159 };
160 
161 #if !defined(IMP_DOXYGEN)
162 typedef DefaultEmbeddingD<1> DefaultEmbedding1D;
163 typedef Vector<DefaultEmbedding1D> DefaultEmbedding1Ds;
164 
165 typedef DefaultEmbeddingD<2> DefaultEmbedding2D;
166 typedef Vector<DefaultEmbedding2D> DefaultEmbedding2Ds;
167 
168 typedef DefaultEmbeddingD<3> DefaultEmbedding3D;
169 typedef Vector<DefaultEmbedding3D> DefaultEmbedding3Ds;
170 
171 typedef DefaultEmbeddingD<4> DefaultEmbedding4D;
172 typedef Vector<DefaultEmbedding4D> DefaultEmbedding4Ds;
173 
174 typedef DefaultEmbeddingD<5> DefaultEmbedding5D;
175 typedef Vector<DefaultEmbedding5D> DefaultEmbedding5Ds;
176 
177 typedef DefaultEmbeddingD<6> DefaultEmbedding6D;
178 typedef Vector<DefaultEmbedding6D> DefaultEmbedding6Ds;
179 
180 typedef DefaultEmbeddingD<-1> DefaultEmbeddingKD;
181 typedef Vector<DefaultEmbeddingKD> DefaultEmbeddingKDs;
182 #endif
183 
184 /** Embedding of a grid as log-evenly spaced axis aligned grid.*/
185 template <int D>
187  VectorD<D> origin_;
188  VectorD<D> unit_cell_;
189  VectorD<D> base_;
190 
191  friend class cereal::access;
192 
193  template<class Archive> void serialize(Archive &ar) {
194  ar(origin_, unit_cell_, base_);
195  }
196 
197  template <class O>
198  VectorD<D> get_coordinates(const O &index) const {
199  VectorD<D> ret = origin_;
200  for (unsigned int i = 0; i < unit_cell_.get_dimension(); ++i) {
201  IMP_USAGE_CHECK(index[i] >= 0, "Out of range index in log graph.'");
202  if (base_[i] != 1) {
203  IMP_USAGE_CHECK(index[i] >= 0,
204  "Log grid axis must have positive index.");
205  ret[i] += unit_cell_[i] * (1.0 - std::pow(base_[i], index[i])) /
206  (1.0 - base_[i]);
207  } else {
208  ret[i] += unit_cell_[i] * index[i];
209  }
210  }
211  return ret;
212  }
213  template <class O>
214  VectorD<D> get_uniform_offset(const O &v0, double o) const {
215  Floats ret(get_dimension());
216  for (unsigned int i = 0; i < get_dimension(); ++i) {
217  ret[i] = v0[i] + o;
218  }
219  return VectorD<D>(ret.begin(), ret.end());
220  }
221  void initialize_from_box(Ints ns, const BoundingBoxD<D> &bb) {
222  Floats nuc(bb.get_dimension());
223  for (unsigned int i = 0; i < bb.get_dimension(); ++i) {
224  double side = bb.get_corner(1)[i] - bb.get_corner(0)[i];
225  IMP_USAGE_CHECK(side > 0, "Can't have flat grid");
226  nuc[i] = side / ns[i];
227  }
228  set_unit_cell(VectorD<D>(nuc.begin(), nuc.end()));
229  set_origin(bb.get_corner(0));
230  }
231 
232  public:
233  LogEmbeddingD(const VectorD<D> &origin, const VectorD<D> &cell,
234  const VectorD<D> &base) {
235  set_origin(origin);
236  set_unit_cell(cell, base);
237  }
238  /** Embedding of a grid as a log-evenly distributed axis-aligned grid
239  over the bounding box bb.
240 
241  @param bb the bounding box in which the grid is embedded
242  @param bases bases[i] is a positive log base used for the grid
243  spacing in dimension i. Set base[i] to 1 in order
244  to create a standard evenly spaced grid along
245  dimension i.
246  @param counts counts[i] is the number of discrete points in dimension i
247  @param bound_centers if true, then the bounding box tightly bounds
248  the centers of the voxels, not their extents.
249  */
250  LogEmbeddingD(const BoundingBoxD<D> &bb, const VectorD<D> &bases,
251  const Ints &counts, bool bound_centers = false) {
252  set_origin(bb.get_corner(0));
253  VectorD<D> cell = bb.get_corner(0);
254  for (unsigned int i = 0; i < bases.get_dimension(); ++i) {
255  IMP_ALWAYS_CHECK(bases[i] > 0,
256  "LogEmbedding base #" << i << " cannot be negative",
258  // cell[i](1-base[i]^counts[i])/(1-base[i])=width[i]
259  if (bases[i] != 1) {
260  cell[i] = (bb.get_corner(1)[i] - bb.get_corner(0)[i]) * (bases[i] - 1) /
261  (std::pow(bases[i], counts[i]) - 1.0);
262  } else {
263  cell[i] = (bb.get_corner(1)[i] - bb.get_corner(0)[i]) / counts[i];
264  }
266  .9 * cell[i] < bb.get_corner(1)[i] - bb.get_corner(0)[i],
267  "Too large a cell side");
268  IMP_INTERNAL_CHECK(cell[i] > 0, "Non-positive cell side");
269  }
270  set_unit_cell(cell, bases);
271  if (bound_centers) {
272  VectorD<D> lower_corner = get_center(GridIndexD<D>(
273  typename GridIndexD<D>::Filled(), bases.get_dimension(), 0));
274  VectorD<D> upper_corner = get_coordinates(
275  get_uniform_offset(GridIndexD<D>(counts.begin(), counts.end()), -.5));
276  VectorD<D> extents = upper_corner - lower_corner;
277  VectorD<D> uc = cell;
278  VectorD<D> orig = uc;
279  for (unsigned int i = 0; i < uc.get_dimension(); ++i) {
280  uc[i] =
281  (bb.get_corner(1)[i] - bb.get_corner(0)[i]) / extents[i] * uc[i];
282  if (base_[i] == 1) {
283  orig[i] = bb.get_corner(0)[i] - .5 * uc[i];
284  } else {
285  /*orig[i]+ uc[i]*(1.0-std::pow(base_[i], index[i]))
286  /(1.0-base_[i])==bb[i]*/
287  orig[i] = bb.get_corner(0)[i] -
288  uc[i] * (1.0 - std::pow(bases[i], .5)) / (1.0 - bases[i]);
289  }
290  }
291  set_origin(orig);
292  set_unit_cell(uc, bases);
293  }
294  }
295  LogEmbeddingD(const VectorD<D> &, const VectorD<D> &) {
296  IMP_FAILURE("not supported");
297  }
298  LogEmbeddingD() {}
299  void set_origin(const VectorD<D> &o) { origin_ = o; }
300  const VectorD<D> get_origin() const { return origin_; }
301  unsigned int get_dimension() const { return get_origin().get_dimension(); }
302  void set_unit_cell(const VectorD<D> &o, const VectorD<D> &base) {
303  unit_cell_ = o;
304  base_ = base;
305  }
306  void set_unit_cell(const VectorD<D> &o) { unit_cell_ = o; }
307  //! Return the unit cell, relative to the origin.
308  /** That is, the unit cell is
309  \code
310  BoundingBoxD<D>(get_zeros_vector_d<D>(),get_unit_cell());
311  \endcode
312  */
313  const VectorD<D> &get_unit_cell() const { return unit_cell_; }
314  //! Return the index that would contain the voxel if the grid extended there
315  /** For example vectors below the "lower left" corner of the grid have
316  indexes with all negative components. This operation will always
317  succeed.
318  */
321  origin_.get_dimension());
322  for (unsigned int i = 0; i < get_dimension(); ++i) {
323  double d = o[i] - origin_[i];
324  // cache everything
325  double fi = d / unit_cell_[i];
326  double li = std::log(fi) / std::log(base_[i]);
327  ret.access_data().get_data()[i] = static_cast<int>(std::floor(li));
328  }
329  return ret;
330  }
331  GridIndexD<D> get_index(const VectorD<D> &o) const {
332  ExtendedGridIndexD<D> ei = get_extended_index(o);
333  return GridIndexD<D>(ei.begin(), ei.end());
334  }
335  /** \name Center
336  Return the coordinates of the center of the voxel.
337  @{
338  */
339  VectorD<D> get_center(const ExtendedGridIndexD<D> &ei) const {
340  return get_coordinates(get_uniform_offset(ei, .5));
341  }
342  VectorD<D> get_center(const GridIndexD<D> &ei) const {
343  return get_coordinates(get_uniform_offset(ei, .5));
344  }
345  /** @} */
346 
347  /** \name Bounding box
348  Return the bounding box of the voxel.
349  @{
350  */
351  BoundingBoxD<D> get_bounding_box(const ExtendedGridIndexD<D> &ei) const {
352  return BoundingBoxD<D>(get_coordinates(ei),
353  get_coordinates(get_uniform_offset(ei, 1)));
354  }
355  BoundingBoxD<D> get_bounding_box(const GridIndexD<D> &ei) const {
356  return get_bounding_box(ExtendedGridIndexD<D>(ei.begin(), ei.end()));
357  }
358  /** @} */
359  IMP_SHOWABLE_INLINE(LogEmbeddingD, out << "origin: " << origin_
360  << " base: " << base_);
361 };
362 
363 #if !defined(IMP_DOXYGEN)
364 typedef LogEmbeddingD<1> LogEmbedding1D;
365 typedef Vector<LogEmbedding1D> LogEmbedding1Ds;
366 
367 typedef LogEmbeddingD<2> LogEmbedding2D;
368 typedef Vector<LogEmbedding2D> LogEmbedding2Ds;
369 
370 typedef LogEmbeddingD<3> LogEmbedding3D;
371 typedef Vector<LogEmbedding3D> LogEmbedding3Ds;
372 
373 typedef LogEmbeddingD<4> LogEmbedding4D;
374 typedef Vector<LogEmbedding4D> LogEmbedding4Ds;
375 
376 typedef LogEmbeddingD<5> LogEmbedding5D;
377 typedef Vector<LogEmbedding5D> LogEmbedding5Ds;
378 
379 typedef LogEmbeddingD<6> LogEmbedding6D;
380 typedef Vector<LogEmbedding6D> LogEmbedding6Ds;
381 
382 typedef LogEmbeddingD<-1> LogEmbeddingKD;
383 typedef Vector<LogEmbeddingKD> LogEmbeddingKDs;
384 #endif
385 
386 IMPALGEBRA_END_NAMESPACE
387 
388 #endif /* IMPALGEBRA_GRID_EMBEDDINGS_H */
Basic types used by IMP.
#define IMP_SHOWABLE_INLINE(Name, how_to_show)
Declare the methods needed by an object that can be printed.
VectorD< D > get_elementwise_product(const algebra::VectorD< D > &a, const algebra::VectorD< D > &b)
Return the vector that is the elementwise product of the two.
Definition: VectorD.h:461
const VectorD< D > & get_corner(unsigned int i) const
For 0 return lower corner and for 1, the upper corner.
Definition: BoundingBoxD.h:140
#define IMP_FAILURE(message)
A runtime failure for IMP.
Definition: check_macros.h:72
A class to represent a voxel grid.
Represent a real cell in a grid (one within the bounding box)
Definition: grid_indexes.h:170
An index in an infinite grid on space.
Definition: grid_indexes.h:31
ExtendedGridIndexD< D > get_extended_index(const VectorD< D > &o) const
Return the index that would contain the voxel if the grid extended there.
Exception definitions and assertions.
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Definition: check_macros.h:139
A Cartesian vector in D-dimensions.
Definition: VectorD.h:56
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
ExtendedGridIndexD< D > get_extended_index(const VectorD< D > &o) const
Return the index that would contain the voxel if the grid extended there.
const VectorD< D > & get_unit_cell() const
Return the unit cell, relative to the origin.
A bounding box in D dimensions.
A class for storing lists of IMP items.
An axis-aligned bounding box.
Definition: BoundingBoxD.h:28
Helper macros for throwing and handling exceptions.
Embed a grid as an evenly spaced axis aligned grid.
Simple 3D vector class.
const VectorD< D > & get_unit_cell() const
Return the unit cell, relative to the origin.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
const double * get_data() const
Return a pointer to the data stored.
Definition: VectorBaseD.h:236
#define IMP_ALWAYS_CHECK(condition, message, exception_name)
Throw an exception if a check fails.
Definition: check_macros.h:61
An exception for an invalid value being passed to IMP.
Definition: exception.h:136
LogEmbeddingD(const BoundingBoxD< D > &bb, const VectorD< D > &bases, const Ints &counts, bool bound_centers=false)