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