cctbx: cctbx/maptbx/fft.h Source File

00001 #ifndef CCTBX_MAPTBX_FFT_H
00002 #define CCTBX_MAPTBX_FFT_H
00003 
00004 #include <cctbx/maptbx/structure_factors.h>
00005 #include <cctbx/maptbx/copy.h>
00006 #include <scitbx/fftpack/real_to_complex_3d.h>
00007 
00008 namespace cctbx { namespace maptbx {
00009 
00011 
00015   template <typename FloatType>
00016   af::versa<FloatType, af::c_grid_padded<3> >
00017   fft_to_real_map(
00018     sgtbx::space_group const& space_group,
00019     af::tiny<int, 3> const& n_real,
00020     af::const_ref<miller::index<> > const& miller_indices,
00021     af::const_ref<std::complex<FloatType> > const& data)
00022   {
00023     // determine padding
00024     scitbx::fftpack::real_to_complex_3d<FloatType> rfft(
00025       (af::c_grid<3>(n_real)));
00026     // allocate map and copy structure factors
00027     structure_factors::to_map<FloatType> to_map(
00028       space_group,
00029       /* anomalous_flag */ false,
00030       miller_indices,
00031       data,
00032       rfft.n_real(),
00033       af::c_grid_padded<3>(rfft.n_complex(), rfft.n_complex()),
00034       /* conjugate_flag */ true,
00035       /* treat_restricted */ true);
00036     // tedious but trivial conversion to suitable array reference
00037     af::ref<std::complex<FloatType>, af::c_grid<3> > complex_map(
00038       to_map.complex_map().begin(), af::c_grid<3>(rfft.n_complex()));
00039     // actual FFT
00040     rfft.backward(complex_map);
00041     // cast from complex to real
00042     return af::versa<FloatType, af::c_grid_padded<3> >(
00043       to_map.complex_map().handle(),
00044       af::c_grid_padded<3>(rfft.m_real(), rfft.n_real()));
00045   }
00046 
00048   template <typename FloatType>
00049   af::versa<FloatType, af::c_grid<3> >
00050   fft_to_real_map_unpadded(
00051     sgtbx::space_group const& space_group,
00052     af::tiny<int, 3> const& n_real,
00053     af::const_ref<miller::index<> > const& miller_indices,
00054     af::const_ref<std::complex<FloatType> > const& data)
00055   {
00056     af::versa<FloatType, af::c_grid_padded<3> > map = fft_to_real_map(
00057       space_group, n_real, miller_indices, data);
00058     unpad_in_place(map.begin(), map.accessor().all(), map.accessor().focus());
00059     return af::versa<FloatType, af::c_grid<3> >(
00060       map, af::c_grid<3>(map.accessor().focus()));
00061   }
00062 
00063 }} // namespace cctbx::maptbx
00064 
00065 #endif // CCTBX_MAPTBX_FFT_H