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