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