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