00001 #ifndef CCTBX_MAPTBX_EIGHT_POINT_INTERPOLATION_H
00002 #define CCTBX_MAPTBX_EIGHT_POINT_INTERPOLATION_H
00003
00004 #include <cctbx/coordinates.h>
00005 #include <scitbx/math/modulo.h>
00006 #include <scitbx/array_family/accessors/c_grid_padded.h>
00007 #include <scitbx/math/utils.h>
00008 #include <scitbx/math/floating_point_epsilon.h>
00009 #include <cctbx/crystal/direct_space_asu.h>
00010
00011 namespace cctbx { namespace maptbx {
00012
00013 template <typename IndexType,
00014 typename FloatType,
00015 typename SignedIntType=long>
00016 class get_corner
00017 {
00018 public:
00019 typedef typename IndexType::value_type iv_t;
00020
00021 get_corner() {}
00022
00023 get_corner(
00024 IndexType const& grid_n,
00025 fractional<FloatType> const& x_frac)
00026 {
00027 for(std::size_t i=0;i<3;i++) {
00028 FloatType xn = x_frac[i] * static_cast<FloatType>(grid_n[i]);
00029 SignedIntType ixn = scitbx::math::float_int_conversions<
00030 FloatType, SignedIntType>::ifloor(xn);
00031 i_grid[i] = scitbx::math::mod_positive(
00032 ixn, static_cast<SignedIntType>(grid_n[i]));
00033 weights_[i][1] = xn - static_cast<FloatType>(ixn);
00034 weights_[i][0] = 1 - weights_[i][1];
00035 }
00036 }
00037
00038 get_corner(
00039 scitbx::mat3<FloatType> const& gridding_matrix,
00040 scitbx::vec3<FloatType> const& site_cart)
00041 {
00042 scitbx::vec3<FloatType> grid_float = gridding_matrix * site_cart;
00043 for(std::size_t i=0;i<3;i++) {
00044 SignedIntType ixn = scitbx::math::float_int_conversions<
00045 FloatType, SignedIntType>::ifloor(grid_float[i]);
00046 i_grid[i] = ixn;
00047 weights_[i][1] = grid_float[i] - static_cast<FloatType>(ixn);
00048 weights_[i][0] = 1 - weights_[i][1];
00049 }
00050 }
00051
00052 get_corner(
00053 crystal::direct_space_asu::asu_mappings<FloatType> & am,
00054 IndexType const& grid_n,
00055 fractional<FloatType> const& x_frac)
00056 {
00057 cartesian<FloatType> const & x_cart = am.process(x_frac).mappings().back()[0].mapped_site();
00058 fractional<FloatType> new_x_frac = am.unit_cell().fractionalize(x_cart);
00059 FloatType epsilon = scitbx::math::floating_point_epsilon<FloatType>::get() * 10;
00060 for ( std::size_t i=0; i<3; ++i )
00061 if ( std::abs(new_x_frac[i]) < epsilon )
00062 new_x_frac[i] = 0;
00063 for(std::size_t i=0;i<3;i++) {
00064 FloatType xn = new_x_frac[i] * static_cast<FloatType>(grid_n[i]);
00065 SignedIntType ixn = scitbx::math::float_int_conversions<
00066 FloatType, SignedIntType>::ifloor(xn);
00067 i_grid[i] = ixn;
00068 weights_[i][1] = xn - static_cast<FloatType>(ixn);
00069 weights_[i][0] = 1 - weights_[i][1];
00070 }
00071 }
00072
00073 FloatType
00074 weight(iv_t s0, iv_t s1, iv_t s2) const
00075 {
00076 return weights_[0][s0] * weights_[1][s1] * weights_[2][s2];
00077 }
00078
00079 IndexType
00080 closest_grid_point(IndexType const& grid_n) const
00081 {
00082 IndexType result = i_grid;
00083 for(std::size_t i=0;i<3;i++) {
00084 if (weights_[i][1] > weights_[i][0]) {
00085 result[i]++;
00086 if (result[i] == grid_n[i]) {
00087 result[i] = 0;
00088 }
00089 }
00090 }
00091 return result;
00092 }
00093
00094 IndexType i_grid;
00095
00096 protected:
00097 FloatType weights_[3][2];
00098 };
00099
00100 template <typename FloatType>
00101 FloatType
00102 eight_point_interpolation(
00103 af::const_ref<FloatType, af::c_grid_padded<3> > const& map,
00104 fractional<FloatType> const& x_frac)
00105 {
00106 typedef af::c_grid_padded<3>::index_type index_t;
00107 typedef typename index_t::value_type iv_t;
00108 index_t const& grid_n = map.accessor().focus();
00109 get_corner<index_t, FloatType> corner(grid_n, x_frac);
00110 FloatType result = 0;
00111 for(iv_t s0=0;s0<2;s0++) { iv_t i0 = (corner.i_grid[0] + s0) % grid_n[0];
00112 for(iv_t s1=0;s1<2;s1++) { iv_t i1 = (corner.i_grid[1] + s1) % grid_n[1];
00113 for(iv_t s2=0;s2<2;s2++) { iv_t i2 = (corner.i_grid[2] + s2) % grid_n[2];
00114 result += map(i0,i1,i2) * corner.weight(s0,s1,s2);
00115 }}}
00116 return result;
00117 }
00118
00119 template <typename FloatType>
00120 typename af::c_grid_padded<3>::index_type
00121 closest_grid_point(
00122 af::flex_grid<> const& map_accessor,
00123 fractional<FloatType> const& x_frac)
00124 {
00125 af::c_grid_padded<3> c_grid(map_accessor);
00126 typedef af::c_grid_padded<3>::index_type index_t;
00127 index_t const& grid_n = c_grid.focus();
00128 return get_corner<index_t, FloatType>(grid_n, x_frac)
00129 .closest_grid_point(grid_n);
00130 }
00131
00132 template <typename FloatType>
00133 FloatType
00134 non_crystallographic_eight_point_interpolation(
00135 af::const_ref<FloatType, af::flex_grid<> > const& map,
00136 scitbx::mat3<FloatType> const& gridding_matrix,
00137 scitbx::vec3<FloatType> const& site_cart,
00138 bool allow_out_of_bounds=false,
00139 FloatType const& out_of_bounds_substitute_value=0)
00140 {
00141 CCTBX_ASSERT(map.accessor().nd() == 3);
00142 typedef typename af::flex_grid<>::index_type index_t;
00143 typedef typename index_t::value_type iv_t;
00144 index_t map_index(3, 0);
00145 get_corner<index_t, FloatType> corner(gridding_matrix, site_cart);
00146 for(unsigned i=0;i<3;i++) {
00147 if(corner.i_grid[i] < map.accessor().origin()[i] ||
00148 corner.i_grid[i] >= map.accessor().focus()[i]-1) {
00149 if(!allow_out_of_bounds) {
00150 throw error("non_crystallographic_eight_point_interpolation:"
00151 " point required for interpolation is out of bounds.");
00152 }
00153 else {
00154 return out_of_bounds_substitute_value;
00155 }
00156 }
00157 }
00158 FloatType result = 0;
00159 for(iv_t s0=0;s0<2;s0++) { map_index[0] = (corner.i_grid[0] + s0);
00160 for(iv_t s1=0;s1<2;s1++) { map_index[1] = (corner.i_grid[1] + s1);
00161 for(iv_t s2=0;s2<2;s2++) { map_index[2] = (corner.i_grid[2] + s2);
00162 result += map(map_index) * corner.weight(s0,s1,s2);
00163 }}}
00164 return result;
00165 }
00166
00167 template <typename FloatType>
00168 FloatType
00169 asu_eight_point_interpolation(
00170 af::const_ref<FloatType, af::flex_grid<> > const& map,
00171 crystal::direct_space_asu::asu_mappings<FloatType> & am,
00172 fractional<FloatType> const& site_cart)
00173 {
00174 CCTBX_ASSERT(map.accessor().nd() == 3);
00175 typedef typename af::flex_grid<>::index_type index_t;
00176 typedef typename index_t::value_type iv_t;
00177 index_t map_index(3, 0);
00178 index_t const& grid_n = map.accessor().focus();
00179 get_corner<index_t, FloatType> corner(am, grid_n, site_cart);
00180 FloatType epsilon = scitbx::math::floating_point_epsilon<FloatType>::get() * 10;
00181 FloatType result = 0;
00182 for(iv_t s0=0;s0<2;s0++) {
00183 map_index[0] = (corner.i_grid[0] + s0);
00184 for(iv_t s1=0;s1<2;s1++) {
00185 map_index[1] = (corner.i_grid[1] + s1);
00186 for(iv_t s2=0;s2<2;s2++) {
00187 map_index[2] = (corner.i_grid[2] + s2);
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197 if ( ! map.accessor().is_valid_index(map_index) ) {
00198 fractional<FloatType> lmap;
00199 for ( std::size_t i=0; i<3; ++i ) {
00200 lmap[i] = static_cast<FloatType>(map_index[i]) / grid_n[i];
00201 }
00202 cartesian<FloatType> const & xmap = am.process(lmap).mappings().back()[0].mapped_site();
00203 fractional<FloatType> nxmap = am.unit_cell().fractionalize(xmap);
00204 for ( std::size_t i=0; i<3; ++i ) {
00205 if ( std::abs(nxmap[i]) < epsilon )
00206 nxmap[i] = 0;
00207 map_index[i] = scitbx::math::float_int_conversions<FloatType, long>::
00208 ifloor(nxmap[i] * static_cast<FloatType>(grid_n[i]));
00209 }
00210 }
00211 result += map(map_index) * corner.weight(s0,s1,s2);
00212 }
00213 }
00214 }
00215 return result;
00216 }
00217
00218 }}
00219
00220 #endif // CCTBX_MAPTBX_EIGHT_POINT_INTERPOLATION_H