8 #ifndef IMPALGEBRA_GRID_UTILITY_H
9 #define IMPALGEBRA_GRID_UTILITY_H
11 #include <IMP/algebra/algebra_config.h>
13 #include "internal/grid_interpolation.h"
15 #include "internal/grid_3d_impl.h"
16 IMPALGEBRA_BEGIN_NAMESPACE
22 template <
int D,
class Storage,
class Value,
class Embedding>
25 return internal::get_interpolation_value(g, pt);
35 template <
class Storage,
class Embedding>
38 const Vector3D &v,
const typename Storage::Value &outside = 0) {
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]) {
49 using namespace internal::trilep_helpers;
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) {
56 unsigned int bx = ((i & 2) >> 1);
57 unsigned int by = (i & 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));
64 typename Storage::Value js[2];
65 for (
unsigned int i = 0; i < 2; ++i) {
71 IMPALGEBRA_END_NAMESPACE
73 #include "internal/grid_3d_impl.h"
81 #define IMP_GRID3D_FOREACH_VOXEL(g, action) \
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; \
120 #define IMP_GRID3D_FOREACH_SMALLER_EXTENDED_INDEX_RANGE( \
121 grid, center, lower_corner, upper_corner, action) \
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]) \
A class to represent a voxel grid.
A class to represent a voxel grid.
A voxel grid in d-dimensional space.
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Base class for a simple primitive-like type.
A Cartesian vector in D-dimensions.
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.
Value get_linearly_interpolated(const GridD< D, Storage, Value, Embedding > &g, const VectorD< D > &pt)
Get the value from the grid with linear interpolation.