IMP  2.2.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-2014 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, int &ix, int &iy, int &iz);
114  void prepare_poslist_flipped(em::DensityMap *dmap);
115  void prepare_poslist(em::DensityMap *dmap);
116  void pad_resolution_map();
117  em::DensityMap *crop_margin(em::DensityMap *in_map);
118  // void fftw_translational_search(const algebra::Rotation3D &rot,int i);
119  void fftw_translational_search(const multifit::internal::EulerAngles &rot,
120  int i);
121  //! Detect the top fits
122  FittingSolutionRecords detect_top_fits(
123  const internal::RotScoresVec &rot_scores, bool cluster_fits,
124  double max_translation, double max_clustering_trans,
125  double max_clustering_rotation);
126 
127  public:
128  FFTFitting() : base::Object("FFTFitting%1%") {}
129  //! Fit a molecule inside its density
130  /**
131  \param[in] dmap the density map to fit into
132  \param[in] density_threshold voxels below this value will be treated as 0
133  \param[in] mol2fit the molecule to fit. The molecule has to be a rigid body
134  \param[in] angle_sampling_interval_rad Sample every internal angles
135  \param[in] num_fits_to_report number of top fits to report
136  \param[in] cluster_fits if true the fits are clustered.
137  Not recommended for refinement mode
138  \param[in] num_angle_per_voxel number of rotations to save per voxel
139  \param[in] max_clustering_translation cluster transformations whose
140  translational distance is lower than the parameter
141  \param[in] max_clustering_angle cluster transformations whose
142  rotational distance is lower than the parameter
143  \param[in] angles_filename a file containing angles to sample.
144  if not specified, all angles are sampled
145  */
146  FFTFittingOutput *do_global_fitting(
147  em::DensityMap *dmap, double density_threshold, atom::Hierarchy mol2fit,
148  double angle_sampling_interval_rad, int num_fits_to_report,
149  double max_clustering_translation, double max_clustering_angle,
150  bool cluster_fits = true, int num_angle_per_voxel = 1,
151  const std::string &angles_filename = "");
152  //! Locally fit a molecule inside its density
153  /**
154  \param[in] dmap the density map to fit into
155  \param[in] density_threshold voxels below this value will be treated as 0
156  \param[in] mol2fit the molecule to fit. The molecule has to be a rigid body
157  \param[in] angle_sampling_interval_rad sample the mol
158  within the range of +- this angle
159  \param[in] max_angle_sampling_rad
160  \param[in] max_translation sample the mol within +-
161  this translation is all directions
162  \param[in] num_fits_to_report
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_rotation 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(
174  em::DensityMap *dmap, double density_threshold, atom::Hierarchy mol2fit,
175  double angle_sampling_interval_rad, double max_angle_sampling_rad,
176  double max_translation, int num_fits_to_report, bool cluster_fits,
177  int num_angle_per_voxel, double max_clustering_translation,
178  double max_clustering_rotation, const std::string &angles_filename = "");
179 };
180 
181 //! FFT fit of a molecule in the density
182 /**
183 \param[in] mol2fit a rigid body molecule to fit
184 \param[in] dmap the map to fit into
185 \param[in] density_threshold ignore density below this value
186 \param[in] angle_sampling_interval_rad sampling internal in radians
187  */
188 IMPMULTIFITEXPORT
189 multifit::FittingSolutionRecords fft_based_rigid_fitting(
190  atom::Hierarchy mol2fit, em::DensityMap *dmap, double density_threshold,
191  double angle_sampling_interval_rad);
192 
193 IMPMULTIFIT_END_NAMESPACE
194 
195 #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:94
IMP::atom::Hierarchy
The standard decorator for manipulating molecular structures.
Definition: atom/Hierarchy.h:198
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/Object.h:106
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::algebra::Vector3D
VectorD< 3 > Vector3D
Definition: VectorD.h:395
IMP::core::RigidBody
A decorator for a rigid body.
Definition: rigid_bodies.h:75
SampledDensityMap.h
Sampled density map.