IMP  2.0.1
The Integrative Modeling Platform
grid_utility.h
Go to the documentation of this file.
1 /**
2  * \file IMP/algebra/grid_utility.h \brief A class to represent a voxel grid.
3  *
4  * Copyright 2007-2013 IMP Inventors. All rights reserved.
5  *
6  */
7 
8 #ifndef IMPALGEBRA_GRID_UTILITY_H
9 #define IMPALGEBRA_GRID_UTILITY_H
10 
11 #include <IMP/algebra/algebra_config.h>
12 #include "GridD.h"
13 #include "internal/grid_interpolation.h"
14 #include "grid_indexes.h"
15 #include "internal/grid_3d_impl.h"
16 IMPALGEBRA_BEGIN_NAMESPACE
17 
18 
19 
20 /** Get the value from the grid with linear interpolation. Values outside the
21  bounding box are snapped to the bounding box (effectively extending the
22  boundary values out to infinity).
23 */
24 template <int D, class Storage, class Value, class Embedding>
25 inline Value get_linearly_interpolated(const GridD<D, Storage, Value,
26  Embedding> &g,
27  const VectorD<D> &pt) {
28  base::Vector<VectorD<D> > corners=internal::get_interpolation_corners(g, pt);
29  Floats values= internal::get_interpolation_values(g, corners);
30  VectorD<D> fraction= internal::get_interpolation_fraction(g, pt);
31  return internal::get_interpolation_value(values, fraction);
32 }
33 
34 
35 
36 //! Use trilinear interpolation to compute a smoothed value at v
37 /** The voxel values are assumed to be at the center of the voxel
38  and the passed outside value is used for voxels outside the
39  grid. The type Voxel must support get_linearly_interpolated().
40  \see get_linearly_interpolated()
41  \relatesalso GridD
42 */
43 template <class Storage, class Embedding>
44 inline const typename Storage::Value
45 get_trilinearly_interpolated(const GridD<3, Storage,
46  typename Storage::Value, Embedding> &g,
47  const Vector3D &v,
48  const typename Storage::Value& outside=0) {
49  // trilirp in z, y, x
50  const Vector3D halfside= g.get_unit_cell()*.5;
51  const Vector3D bottom_sample= g.get_bounding_box().get_corner(0)+halfside;
52  const Vector3D top_sample= g.get_bounding_box().get_corner(1)-halfside;
53  for (unsigned int i=0; i< 3; ++i){
54  if (v[i] < bottom_sample[i]
55  || v[i] >= top_sample[i]) {
56  //std::cout << v << " was rejected." << std::endl;
57  return outside;
58  }
59  }
60  using namespace internal::trilep_helpers;
61  int ivox[3];
62  algebra::Vector3D r;
63  internal::trilep_helpers::compute_voxel(g, v, ivox, r);
64  typename Storage::Value is[4];
65  for (unsigned int i=0; i< 4; ++i) {
66  // operator >> has high precidence compared. Go fig.
67  unsigned int bx= ((i&2) >> 1);
68  unsigned int by= (i&1);
69  IMP_INTERNAL_CHECK((bx==0 || bx==1) && (by==0 || by==1),
70  "Logic error in trilerp");
71  is[i]=get_linearly_interpolated(1-r[2],
72  get_value(g, ivox[0]+bx, ivox[1]+by,
73  ivox[2], outside),
74  get_value(g, ivox[0]+bx, ivox[1]+by,
75  ivox[2]+1U, outside));
76  }
77  typename Storage::Value js[2];
78  for (unsigned int i=0; i< 2; ++i) {
79  js[i]= get_linearly_interpolated(1-r[1], is[i*2], is[i*2+1]);
80  }
81  return get_linearly_interpolated(1-r[0], js[0], js[1]);
82 }
83 
84 IMPALGEBRA_END_NAMESPACE
85 
86 
87 #include "internal/grid_3d_impl.h"
88 
89 
90 /** Iterate over each voxel in grid. The voxel index is
91  GridIndexD<3> voxel_index and the coordinates of the center is
92  Vector3D voxel_center and the index of the voxel is
93  loop_voxel_index.
94  \relatesalso Grid3D
95  */
96 #define IMP_GRID3D_FOREACH_VOXEL(g, action) \
97  { \
98  unsigned int next_loop_voxel_index=0; \
99  const IMP::algebra::Vector3D macro_map_unit_cell=g.get_unit_cell(); \
100  const int macro_map_nx=g.get_number_of_voxels(0); \
101  const int macro_map_ny=g.get_number_of_voxels(1); \
102  const int macro_map_nz=g.get_number_of_voxels(2); \
103  const IMP::algebra::Vector3D macro_map_origin \
104  =g.get_origin(); \
105  IMP::algebra::GridIndexD<3> voxel_index; \
106  int *voxel_index_data=voxel_index.access_data().get_data(); \
107  IMP::algebra::Vector3D voxel_center; \
108  for (voxel_index_data[0]=0; \
109  voxel_index_data[0]< macro_map_nx; \
110  ++voxel_index_data[0]) { \
111  voxel_center[0]= macro_map_origin[0] \
112  +(voxel_index_data[0]+.5) \
113  *macro_map_unit_cell[0]; \
114  for (voxel_index_data[1]=0; voxel_index_data[1]< macro_map_ny; \
115  ++voxel_index_data[1]) { \
116  voxel_center[1]= macro_map_origin[1] \
117  +(voxel_index_data[1]+.5)*macro_map_unit_cell[1]; \
118  for (voxel_index_data[2]=0; voxel_index_data[2]< macro_map_nz; \
119  ++voxel_index_data[2]) { \
120  voxel_center[2]= macro_map_origin[2] \
121  +(voxel_index_data[2]+.5)*macro_map_unit_cell[2]; \
122  unsigned int loop_voxel_index=next_loop_voxel_index; \
123  IMP_UNUSED(loop_voxel_index); \
124  ++next_loop_voxel_index; \
125  {action}; \
126  } \
127  } \
128  } \
129  } \
130 
131 
132 
133 /** Iterate over each voxel in a subset of the grid that are less than
134  center. The voxel index is unsigned int voxel_index[3]. Use this if,
135  for example you want to find nearby pairs of voxels once each.
136  \relatesalso Grid3D
137 */
138 #define IMP_GRID3D_FOREACH_SMALLER_EXTENDED_INDEX_RANGE(grid, center, \
139  lower_corner, \
140  upper_corner, \
141  action) \
142  { int voxel_index[3]; \
143  IMP_USAGE_CHECK(lower_corner <= upper_corner, "Inverted range " \
144  << lower_corner << " " << upper_corner); \
145  IMP_USAGE_CHECK(lower_corner <= center, "Center not in range " \
146  << lower_corner << " " << center); \
147  IMP_USAGE_CHECK(center <= upper_corner, "Center not in range " \
148  << center << upper_corner); \
149  for (voxel_index[0]=lower_corner[0]; \
150  voxel_index[0] <= upper_corner[0]; ++voxel_index[0]) { \
151  if (voxel_index[0] > center[0]) break; \
152  for (voxel_index[1]=lower_corner[1]; \
153  voxel_index[1] <= upper_corner[1]; ++voxel_index[1]) { \
154  if (voxel_index[0] == center[0] \
155  && voxel_index[1] > center[1]) break; \
156  for (voxel_index[2]=lower_corner[2]; \
157  voxel_index[2] <= upper_corner[2]; ++voxel_index[2]) { \
158  if (voxel_index[0] == center[0] && voxel_index[1] == center[1]\
159  && voxel_index[2] >= center[2]) break; \
160  { action} \
161  } \
162  } \
163  } \
164  }
165 #endif /* IMPALGEBRA_GRID_UTILITY_H */