00001 #ifndef CCTBX_MAPTBX_GRIDDING_H
00002 #define CCTBX_MAPTBX_GRIDDING_H
00003
00004 #include <cctbx/sgtbx/search_symmetry.h>
00005 #include <scitbx/fftpack/gridding.h>
00006 #include <scitbx/array_family/tiny_algebra.h>
00007
00008 namespace cctbx { namespace maptbx {
00009
00010 namespace detail {
00011 static const af::tiny<int, 3> one_one_one(1,1,1);
00012 }
00013
00014 template <typename IndexValueType>
00015 af::tiny<IndexValueType, 3>
00016 determine_gridding(
00017 uctbx::unit_cell const& unit_cell,
00018 double d_min,
00019 double resolution_factor,
00020 af::tiny<IndexValueType, 3> const& mandatory_factors=detail::one_one_one,
00021 IndexValueType max_prime=5,
00022 bool assert_shannon_sampling=true)
00023 {
00024 CCTBX_ASSERT(d_min > 0);
00025 CCTBX_ASSERT(resolution_factor > 0);
00026 if (assert_shannon_sampling) {
00027 CCTBX_ASSERT(resolution_factor <= 0.5);
00028 }
00029 af::tiny<IndexValueType, 3>
00030 grid(unit_cell.max_miller_indices(d_min * 2 * resolution_factor));
00031 grid *= IndexValueType(2);
00032 grid += IndexValueType(1);
00033 return scitbx::fftpack::adjust_gridding_array(
00034 grid, max_prime, mandatory_factors);
00035 }
00036
00037 template <typename IndexValueType>
00038 af::tiny<IndexValueType, 3>
00039 determine_gridding(
00040 uctbx::unit_cell const& unit_cell,
00041 double d_min,
00042 double resolution_factor,
00043 sgtbx::search_symmetry_flags const& symmetry_flags,
00044 sgtbx::space_group_type const& space_group_type,
00045 af::tiny<IndexValueType, 3> const& mandatory_factors=detail::one_one_one,
00046 IndexValueType max_prime=5,
00047 bool assert_shannon_sampling=true)
00048 {
00049 typedef IndexValueType i_v_t;
00050 sgtbx::structure_seminvariants ss;
00051 sgtbx::space_group sub_space_group;
00052 af::tiny<i_v_t, 3> mandatory_factors_ = mandatory_factors;
00053 if (symmetry_flags.use_seminvariants()) {
00054 ss = sgtbx::structure_seminvariants(space_group_type.group());
00055 mandatory_factors_ = ss.refine_gridding(mandatory_factors_);
00056 sub_space_group = sgtbx::search_symmetry(
00057 symmetry_flags, space_group_type, ss).subgroup();
00058 }
00059 else {
00060 sub_space_group = sgtbx::search_symmetry(
00061 symmetry_flags, space_group_type).subgroup();
00062 }
00063 mandatory_factors_ = sub_space_group.refine_gridding(mandatory_factors_);
00064 af::tiny<i_v_t, 3>
00065 grid = determine_gridding(
00066 unit_cell, d_min, resolution_factor,
00067 mandatory_factors_, max_prime, assert_shannon_sampling);
00068 std::size_t best_size = 0;
00069 af::tiny<i_v_t, 3> best_grid(0,0,0);
00070 i_v_t g_limit = af::max(grid) + 1;
00071 af::tiny<i_v_t, 3> loop;
00072 for(loop[0]=grid[0];loop[0]<g_limit;loop[0]+=mandatory_factors_[0])
00073 for(loop[1]=grid[1];loop[1]<g_limit;loop[1]+=mandatory_factors_[1])
00074 for(loop[2]=grid[2];loop[2]<g_limit;loop[2]+=mandatory_factors_[2]) {
00075 af::tiny<i_v_t, 3> trial_grid = scitbx::fftpack::adjust_gridding_array(
00076 loop, max_prime, mandatory_factors_);
00077 if (symmetry_flags.use_seminvariants()) {
00078 trial_grid = ss.refine_gridding(trial_grid);
00079 }
00080 trial_grid = sub_space_group.refine_gridding(trial_grid);
00081 CCTBX_ASSERT(scitbx::fftpack::adjust_gridding_array(
00082 trial_grid, max_prime, mandatory_factors_).all_eq(trial_grid));
00083 if (best_size == 0 && trial_grid.all_eq(grid)) {
00084 return grid;
00085 }
00086 std::size_t trial_size = 1;
00087 for(std::size_t i=0;i<3;i++) trial_size *= trial_grid[i];
00088 CCTBX_ASSERT(trial_size > 0);
00089 if (best_size == 0 || trial_size < best_size) {
00090 best_grid = trial_grid;
00091 best_size = trial_size;
00092 }
00093 }
00094 return best_grid;
00095 }
00096
00097 }}
00098
00099 #endif // CCTBX_MAPTBX_GRIDDING_H
