IMP  2.0.1
The Integrative Modeling Platform
rigid_fitting.h
Go to the documentation of this file.
1 /**
2  * \file IMP/em/rigid_fitting.h
3  * \brief preforms rigid fitting between a set of particles and a density map
4  *
5  * Copyright 2007-2013 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPEM_RIGID_FITTING_H
9 #define IMPEM_RIGID_FITTING_H
10 
11 #include <IMP/core/MonteCarlo.h>
13 #include <IMP/RestraintSet.h>
15 #include <IMP/VersionInfo.h>
16 #include <IMP/Particle.h>
17 #include <IMP/Model.h>
18 #include "DensityMap.h"
19 #include "FitRestraint.h"
20 #include <IMP/em/em_config.h>
21 #include <IMP/OptimizerState.h>
22 #include <IMP/core/rigid_bodies.h>
23 #include <IMP/ScoreState.h>
24 #include <algorithm>
25 IMPEM_BEGIN_NAMESPACE
26 
27 
28 //! A simple list of fitting solutions.
29 /**
30  \see local_rigid_fitting_around_point
31  \see local_rigid_fitting_around_points
32  \see local_rigid_fitting_grid_search
33  \see compute_fitting_scores
34  */
35 class IMPEMEXPORT FittingSolutions {
36 typedef std::pair<algebra::Transformation3D,Float> FittingSolution;
37 struct sort_by_cc
38 {
39  bool operator()(const FittingSolution &s1, const FittingSolution & s2) const
40  {
41  return s1.second < s2.second;
42  }
43 };
44 public:
46  //! Get the number of solutions in the set
47  inline int get_number_of_solutions() const {return fs_.size();}
48  //! Get the score of the i'th solution
49  /**
50  \return the i'th transformation, or throw an exception
51  if the index is out of range
52  */
53  inline algebra::Transformation3D get_transformation(unsigned int i) const {
54  IMP_USAGE_CHECK(i<fs_.size(),"The index requested ("<<
55  i<<") in get_transformation is our of range ("<<
56  fs_.size()<<")"<<std::endl);
57  return fs_[i].first;
58  }
59  //! Get the score of the i'th solution
60  /**
61  \return the i'th score, or throw an exception
62  if the index is out of range
63  */
64  inline Float get_score(unsigned int i) const {
65  IMP_USAGE_CHECK(i<fs_.size(),"The index requested ("<<
66  i<<") in get_transformation is out of range ("<<
67  fs_.size()<<")"<<std::endl);
68  return fs_[i].second;
69  }
70  void set_score(unsigned int i,Float score) {
71  IMP_USAGE_CHECK(i<fs_.size(),"The index requested ("<<
72  i<<") in get_transformation is out of range ("<<
73  fs_.size()<<")"<<std::endl);
74  fs_[i].second=score;
75  }
76  //! Add a solution to the fitting solution set
77  void add_solution(const algebra::Transformation3D &t,Float score);
78  //! Sort solutions by cross-correlation scores
79  void sort(bool reverse=false);
80  //! Multiply each transformation (T) by t,
81  //!such that the new transformation are T*t
83  for(unsigned int i=0;i<fs_.size();i++) fs_[i].first=fs_[i].first*t;
84  }
85  inline algebra::Transformation3Ds get_transformations() const {
87  for(unsigned int i=0;i<fs_.size();i++) all_ts.push_back(fs_[i].first);
88  return all_ts;
89  }
90  void show(std::ostream& out=std::cout) const {
91  for(std::vector<FittingSolution>::const_iterator it = fs_.begin();
92  it != fs_.end(); it++) {
93  out<<"("<<it->first<<" , "<<it->second<<")"<<std::endl;
94  }
95  }
96 protected:
97  std::vector<FittingSolution> fs_;
98 };
100 
101 //! Local rigid fitting of a rigid body around a center point
102 /**
103 \brief Fit a set of particles to a density map around an anchor point.
104  The fitting is assessed using the cross-correaltion score.
105  The optimization is a standard MC/CG procedure.
106  The function returns a list of solutions sortedo the cross-correlation
107  score.
108 \note The returned cross-correlation score is 1-cc, as we usually want to
109  minimize a scroing function. Thus a score of 1 means no-correlation
110  and a score of 0. is perfect correlation.
111 \note The input rigid body should be also IMP::atom::Hierarchy
112 \param[in] p The rigid body to fit
113 \param[in] refiner Refiner to yield rigid body members
114 \param[in] weight_key The weight key of the particles in the rigid body
115 \param[in] dmap The density map to fit to
116 \param[in] anchor_centroid The point to fit the particles around
117 \param[in] display_log If provided, then intermediate states
118  in during the optimization procedure are printed
119 \param[in] number_of_optimization_runs number of Monte Carlo optimizations
120 \param[in] number_of_mc_steps number of steps in a Monte Carlo optimization
121 \param[in] number_of_cg_steps number of Conjugate Gradients steps in
122  a Monte Carlo step
123 \param[in] max_translation maximum translation step in a MC optimization step
124 \param[in] max_rotation maximum rotation step in a single MC optimization step
125 \param[in] fast if true the density map of the rigid body is not resampled
126  but transformed at each iteration of the optimization
127 \return the refined fitting solutions
128 */
130  Particle *p, Refiner *refiner,
131  const FloatKey &weight_key,
132  DensityMap *dmap, const algebra::Vector3D &anchor_centroid,
133  OptimizerStates display_log,
134  Int number_of_optimization_runs = 5, Int number_of_mc_steps = 10,
135  Int number_of_cg_steps=100,
136  Float max_translation=2., Float max_rotation=.3,bool fast=false);
137 
138 //! Local rigid fitting of a rigid body
139 /**
140 \brief Fit a set of particles to a density map around their centroid.
141  The fitting is assessed using the cross-correaltion score.
142  The optimization is a standard MC/CG procedure.
143  The function returns a list of solutions sortedo the cross-correlation
144  score.
145 \note The returned cross-correlation score is 1-cc, as we usually want to
146  minimize a scroing function. Thus a score of 1 means no-correlation
147  and a score of 0. is perfect correlation.
148 \note The input rigid body should be also IMP::atom::Hierarchy
149 \param[in] p The root of the hierarchy to fit
150 \param[in] refiner The refiner to get the leaves of the particle
151 \param[in] weight_key The weight key of the particles in the rigid body
152 \param[in] dmap The density map to fit to
153 \param[in] display_log If provided, then intermediate states
154  in during the optimization procedure are printed
155 \param[in] number_of_optimization_runs number of Monte Carlo optimizations
156 \param[in] number_of_mc_steps number of steps in a Monte Carlo optimization
157 \param[in] number_of_cg_steps number of Conjugate Gradients steps in
158  a Monte Carlo step
159 \param[in] max_translation maximum translation step in a MC optimization step
160 \param[in] max_rotation maximum rotation step in radians in a single
161  MC optimization step
162 \param[in] fast if true the density map of the rigid body is not resampled
163  but transformed at each iteration of the optimization
164 \return the refined fitting solutions
165 */
166 
168  Particle *p, Refiner *refiner,
169  const FloatKey &weight_key,
170  DensityMap *dmap,
171  OptimizerStates display_log,
172  Int number_of_optimization_runs = 5, Int number_of_mc_steps = 10,
173  Int number_of_cg_steps=100,
174  Float max_translation=2., Float max_rotation=.3,
175  bool fast=true) {
176  IMP_LOG_VERBOSE("Start: local_rigid_fitting\n");
177  algebra::Vector3D rb_cen=
178  IMP::core::get_centroid(core::XYZs(refiner->get_refined(p)));
179  IMP_LOG_VERBOSE("centroid is:"<<rb_cen<<"\n");
181  p, refiner,weight_key, dmap,
182  rb_cen,display_log,
183  number_of_optimization_runs, number_of_mc_steps,
184  number_of_cg_steps, max_translation, max_rotation,fast);
185 }
186 
187 
188 //! Local rigid fitting of a rigid body around a set of center points
189 /**
190 \brief Fit a set of particles to a density map around each of the input points.
191  The function apply local_rigid_fitting_around_point around each center.
192 \note The input rigid body should be also IMP::atom::Hierarchy
193 \param[in] p The rigid body to fit
194 \param[in] refiner Refiner to yield rigid body members
195 \param[in] wei_key The weight key of the particles in the rigid body
196 \param[in] dmap The density map to fit to
197 \param[in] anchor_centroids The points to fit the particles around
198 \param[in] display_log If provided, then intermediate states
199  in during the optimization procedure are printed
200 \param[in] number_of_optimization_runs number of Monte Carlo optimizations
201 \param[in] number_of_mc_steps number of steps in a Monte Carlo optimization
202 \param[in] number_of_cg_steps number of Conjugate Gradients steps in
203  a Monte Carlo step
204 \param[in] max_translation maximum translation step in a MC optimization step
205 \param[in] max_rotation maximum rotation step in a single MC optimization step
206 \return the refined fitting solutions
207 */
208 IMPEMEXPORT FittingSolutions local_rigid_fitting_around_points(
209  Particle *p,Refiner *refiner,
210  const FloatKey &wei_key,
211  DensityMap *dmap, const algebra::Vector3Ds &anchor_centroids,
212  OptimizerStates display_log,
213  Int number_of_optimization_runs = 5, Int number_of_mc_steps = 10,
214  Int number_of_cg_steps=100,
215  Float max_translation=2., Float max_rotation=.3);
216 
217 
218 //! Local grid search rigid fitting
219 /**
220 \brief Fit a set of particles to a density map around their centroid.
221  The fitting is assessed using the cross-correaltion score.
222  The optimization is a grid search
223 \note The transformations are not clustered.
224 \note The returned cross-correlation score is 1-cc, as we usually want to
225  minimize a scroing function. Thus a score of 1 means no-correlation
226  and a score of 0. is perfect correlation.
227 \param[in] ps The particles to be fitted (treated rigid)
228 \param[in] wei_key The weight key of the particles in the rigid body
229 \param[in] dmap The density map to fit to
230 \param[in] max_voxels_translation Sample translations within
231  -max_voxel_translation to max_voxel_translation
232 \param[in] translation_step The translation sampling step
233 \param[in] max_angle_in_radians Sample rotations with +- max_angle_in_radians
234  around the current rotation
235 \param[in] number_of_rotations The number of rotations to sample
236 \return the refined fitting solutions
237 */
238 IMPEMEXPORT FittingSolutions local_rigid_fitting_grid_search(
239  const ParticlesTemp &ps,
240  const FloatKey &wei_key,
241  DensityMap *dmap,
242  Int max_voxels_translation=2,
243  Int translation_step=1,
244  Float max_angle_in_radians = 0.174,
245  Int number_of_rotations = 100);
246 
247 
248 //! Compute fitting scores for a given set of rigid transformations
249 /**
250 \brief Score how well a set of particles fit to the map in various
251  rigid transformations.
252 \param[in] ps The particles to be fitted (treated rigid)
253 \param[in] em_map The density map to fit to
254 \param[in] wei_key The weight key of the particles in the rigid body
255 \param[in] fast_version If fast_version is used the sampled density map of the
256  input particles (ps) is not resampled for each
257  transformation but instead the corresponding grid
258  is rotated. This option significantly improves the
259  running times but the returned scores are less accurate
260 \param[in] transformations A set of rigid transformations
261 \param[in] fast_version if true the density map of each transformation
262  is interpolated
263 \param[in] local_score if true a local cross correlation score is used
264 \return The scored fitting solutions
265 \note the function assumes the density map holds its density
266  */
267 IMPEMEXPORT FittingSolutions compute_fitting_scores(const ParticlesTemp &ps,
268  DensityMap *em_map,
269  const algebra::Transformation3Ds &transformations,
270  bool fast_version=false, bool local_score=false,
271  const FloatKey &wei_key=atom::Mass::get_mass_key());
272 
273 
274 //! Compute fitting scores for a given set of rigid transformations
275 /**
276 \brief Score how well a rigid body fits to the map
277 \param[in] em_map The density map to fit to
278 \param[in] rb The rigid body
279 \param[in] refiner The rigid body refiner
280 \param[in] transformations Transformations of the rigid body
281 \return The scored fitting solutions
282 \note the function assumes the density map holds its density
283  */
285  DensityMap *em_map,
286  core::RigidBody rb,Refiner *refiner,
287  const algebra::Transformation3Ds& transformations) {
288  return compute_fitting_scores(refiner->get_refined(rb),em_map,
289  transformations,true);
290 }
291 
292 
293 
294 //! Compute fitting scores for a given set of rigid transformations
295 /**
296 \brief Scores how well a set of particles fit a map
297 \param[in] ps The particles to be fitted
298 \param[in] em_map The density map to fit to
299 \param[in] wei_key The weight key of the particles in the rigid body
300 \note the function assumes the density map holds its density
301  */
302 IMPEMEXPORT Float compute_fitting_score(const ParticlesTemp &ps,
303  DensityMap *em_map,
304  FloatKey wei_key=atom::Mass::get_mass_key());
305 
306 
307 IMPEM_END_NAMESPACE
308 #endif /* IMPEM_RIGID_FITTING_H */