IMP  2.0.1
The Integrative Modeling Platform
fft_based_rigid_fitting.h
Go to the documentation of this file.
1 /**
2  * \file IMP/multifit/fft_based_rigid_fitting.h
3  * \brief FFT based fitting
4  *
5  * Copyright 2007-2013 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPMULTIFIT_FFT_BASED_RIGID_FITTING_H
10 #define IMPMULTIFIT_FFT_BASED_RIGID_FITTING_H
11 
12 #include "fftw3.h"
13 #include <IMP/multifit/multifit_config.h>
14 #include <IMP/atom/Hierarchy.h>
15 #include <IMP/base_types.h>
16 #include <IMP/multifit/internal/FFTWGrid.h>
17 #include <IMP/multifit/FittingSolutionRecord.h>
18 #include <IMP/multifit/internal/FFTWPlan.h>
19 #include <IMP/em/DensityMap.h>
20 #include <IMP/em/SampledDensityMap.h>
21 #include <IMP/em/KernelParameters.h>
22 #include <IMP/multifit/internal/fft_fitting_utils.h>
23 #include <boost/scoped_array.hpp>
24 
25 IMPMULTIFIT_BEGIN_NAMESPACE
26 
27 class IMPMULTIFITEXPORT FFTFittingOutput : public base::Object {
28 public:
29  FFTFittingOutput() : base::Object("FFTFittingOutput%1%") {}
30 
31  IMP_OBJECT_INLINE(FFTFittingOutput,
32  { out << best_fits_.size() << " final fits; "
33  << best_trans_per_rot_.size()
34  << " best translations per rotation"; }, {});
35 public:
36  FittingSolutionRecords best_fits_; //final fits
37  FittingSolutionRecords best_trans_per_rot_;
38 };
39 
40 class IMPMULTIFITEXPORT FFTFitting : public base::Object {
41  IMP_OBJECT_INLINE(FFTFitting, {IMP_UNUSED(out);}, {});
42 
43  protected:
44  //logging
45  FittingSolutionRecords best_trans_per_rot_log_;
46  //map parameters
47  algebra::Transformation3D cen_trans_;
48  atom::Hierarchy high_mol_;
49  double resolution_;
50  unsigned int nx_,ny_,nz_; //extent
51  unsigned long nvox_; //total number of voxels
52  unsigned int nx_half_,ny_half_,nz_half_;// half of the map extent
53  double spacing_; //map voxel size
54  double origx_,origy_,origz_;// map origin
55  internal::FFTWGrid<double> low_map_data_; // low resolution map
56  Pointer<em::DensityMap> low_map_;
57  Pointer<em::SampledDensityMap> sampled_map_;//sampled from protein
58  internal::FFTWGrid<double> sampled_map_data_,fftw_r_grid_mol_;
59  // high resolution map
60  internal::FFTWGrid<double> reversed_fftw_data_;
61  boost::scoped_array<double> kernel_filter_;
62  unsigned int kernel_filter_ext_;
63  boost::scoped_array<double> gauss_kernel_; // low-pass (Gaussian) kernel
64  unsigned int gauss_kernel_ext_; //Gaussian kernel extent
65  unsigned long gauss_kernel_nvox_; //Gaussian kernel number of voxels
66  boost::scoped_array<double> filtered_kernel_; //filtered low-pass kernel
67  // unsigned long filtered_kernel_nvox_;
68  //number of voxels in the filtered kernel
69  unsigned filtered_kernel_ext_; //filtered low-pass kernel extent
70 
71  double sampled_norm_, asmb_norm_; //normalization factors for both maps
72  algebra::Vector3D map_cen_;
73  //FFT variables
74  unsigned long fftw_nvox_r2c_; /* FFTW real to complex voxel count */
75  unsigned long fftw_nvox_c2r_; /* FFTW complex to real voxel count */
76  internal::FFTWGrid<fftw_complex> fftw_grid_lo_,fftw_grid_hi_;
77  internal::FFTWPlan fftw_plan_forward_lo_,fftw_plan_forward_hi_;
78  internal::FFTWPlan fftw_plan_reverse_hi_;
79  double fftw_scale_; //eq to 1./nvox_
80 
81  //molecule to fit
82  atom::Hierarchy orig_mol_;
83  atom::Hierarchy orig_mol_copy_;//keep a copy in the original transformation
84  atom::Hierarchy copy_mol_;//rotated mol because
85  //we use an alternative rotating mechanism
86  core::RigidBody orig_rb_;
87  int num_angle_per_voxel_;
88  int num_fits_reported_;
89  double low_cutoff_;
90  int corr_mode_;
91  algebra::Vector3D orig_cen_;
92  //paddding
93  double fftw_pad_factor_; // grid size expansion factor for FFT padding
94  unsigned int fftw_zero_padding_extent_[3]; // padding extent
95  unsigned margin_ignored_in_conv_[3]; // margin that can be ignored
96  //in convolution
97  internal::RotScoresVec fits_hash_; //stores best fits
98  multifit::FittingSolutionRecords best_fits_;
99  internal::FFTScores fft_scores_;
100  unsigned int inside_num_;
101  unsigned int inside_num_flipped_;
102  internal::FFTScores fft_scores_flipped_;
103  // algebra::Rotation3Ds rots_;
104  multifit::internal::EulerAnglesList rots_;
105 
106  void prepare_probe(atom::Hierarchy mol2fit);
107  void prepare_lowres_map (em::DensityMap *dmap);
108  void prepare_kernels();
109  void copy_density_data(em::DensityMap *dmap,double *data_array);
110  void get_unwrapped_index (int wx, int wy, int wz,
111  int &ix, int &iy, int &iz);
112  void prepare_poslist_flipped (em::DensityMap *dmap);
113  void prepare_poslist(em::DensityMap *dmap);
114  void pad_resolution_map();
115  em::DensityMap* crop_margin(em::DensityMap *in_map);
116  // void fftw_translational_search(const algebra::Rotation3D &rot,int i);
117  void fftw_translational_search(const multifit::internal::EulerAngles &rot,
118  int i);
119  //! Detect the top fits
120  FittingSolutionRecords detect_top_fits(
121  const internal::RotScoresVec &rot_scores,
122  bool cluster_fits, double max_translation,
123  double max_clustering_trans,
124  double max_clustering_rotation);
125  public:
126  FFTFitting() : base::Object("FFTFitting%1%") {}
127  //! Fit a molecule inside its density
128  /**
129  \param[in] dmap the density map to fit into
130  \param[in] density_threshold voxels below this value will be treated as 0
131  \param[in] mol2fit the molecule to fit. The molecule has to be a rigid body
132  \param[in] angle_sampling_interval_rad Sample every internal angles
133  \param[in] num_fits_to_report number of top fits to report
134  \param[in] cluster_fits if true the fits are clustered.
135  Not recommended for refinement mode
136  \param[in] num_angle_per_voxel number of rotations to save per voxel
137  \param[in] max_clustering_translation cluster transformations whose
138  translational distance is lower than the parameter
139  \param[in] max_clustering_angle cluster transformations whose
140  rotational distance is lower than the parameter
141  \param[in] angles_filename a file containing angles to sample.
142  if not specified, all angles are sampled
143  */
144  FFTFittingOutput *do_global_fitting(em::DensityMap *dmap,
145  double density_threshold,
146  atom::Hierarchy mol2fit,
147  double angle_sampling_interval_rad,
148  int num_fits_to_report,
149  double max_clustering_translation,
150  double max_clustering_angle,
151  bool cluster_fits=true,
152  int num_angle_per_voxel=1,
153  const std::string &angles_filename="");
154  //! Locally fit a molecule inside its density
155  /**
156  \param[in] dmap the density map to fit into
157  \param[in] density_threshold voxels below this value will be treated as 0
158  \param[in] mol2fit the molecule to fit. The molecule has to be a rigid body
159  \param[in] angle_sampling_interval_rad sample the mol
160  within the range of +- this angle
161  \param[in] max_translation sample the mol within +-
162  this translation is all directions
163  \param[in] cluster_fits if true the fits are clustered.
164  Not recommended for refinement mode
165  \param[in] num_angle_per_voxel number of rotations to save per voxel
166  \param[in] max_clustering_translation cluster transformations whose
167  translational distance is lower than the parameter
168  \param[in] max_clustering_angle cluster transformations whose
169  rotational distance is lower than the parameter
170  \param[in] angles_filename a file containing angles to sample.
171  if not specified, all angles are sampled
172  */
173  FFTFittingOutput *do_local_fitting(em::DensityMap *dmap,
174  double density_threshold,
175  atom::Hierarchy mol2fit,
176  double angle_sampling_interval_rad,
177  double max_angle_sampling_rad,
178  double max_translation,
179  int num_fits_to_report,
180  bool cluster_fits,
181  int num_angle_per_voxel,
182  double max_clustering_translation,
183  double max_clustering_rotation,
184  const std::string &angles_filename="");
185 };
186 
187 //! FFT fit of a molecule in the density
188 /**
189 \param[in] mol2fit a rigid body molecule to fit
190 \param[in] dmap the map to fit into
191 \param[in] angle_sampling_interval_rad sampling internal in radians
192  */
193 IMPMULTIFITEXPORT
194 multifit::FittingSolutionRecords fft_based_rigid_fitting(
195  atom::Hierarchy mol2fit,
196  em::DensityMap *dmap,
197  double density_threshold,
198  double angle_sampling_interval_rad);
199 
200 IMPMULTIFIT_END_NAMESPACE
201 
202 #endif /* IMPMULTIFIT_FFT_BASED_RIGID_FITTING_H */