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