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
24 template <
int D,
class Storage,
class Value,
class Embedding>
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);
43 template <
class Storage,
class Embedding>
44 inline const typename Storage::Value
46 typename Storage::Value, Embedding> &g,
48 const typename Storage::Value& outside=0) {
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]) {
60 using namespace internal::trilep_helpers;
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) {
67 unsigned int bx= ((i&2) >> 1);
68 unsigned int by= (i&1);
70 "Logic error in trilerp");
72 get_value(g, ivox[0]+bx, ivox[1]+by,
74 get_value(g, ivox[0]+bx, ivox[1]+by,
75 ivox[2]+1U, outside));
77 typename Storage::Value js[2];
78 for (
unsigned int i=0; i< 2; ++i) {
84 IMPALGEBRA_END_NAMESPACE
87 #include "internal/grid_3d_impl.h"
96 #define IMP_GRID3D_FOREACH_VOXEL(g, action) \
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 \
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; \
138 #define IMP_GRID3D_FOREACH_SMALLER_EXTENDED_INDEX_RANGE(grid, center, \
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; \