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