Modified the Vec3List::FitCylinder function to also return the radius

of the fitted cylinder. The input was also changed from taking both a direction
and center as initial guess to just the direction. The center is guessed from
the geometrical center of the atoms and fitted during the optimization
procedure.

modules/geom/pymod/export_vec3.cc

@@ -21,9 +21,11 @@
 #include <ost/geom/vec3.hh>
 #include <ost/geom/geom.hh>
 #include <ost/geom/export_helper/vector.hh>
+#include <ost/export_helper/pair_to_tuple_conv.hh>
 
 using namespace boost::python;
 
+
 const Real Vec3_getitem(const geom::Vec3& v, int i) {
   return v.At(i);
 }
@@ -116,4 +118,6 @@ void export_Vec3()
     .def("GetODRLine", &Vec3List::GetODRLine)
     .def("FitCylinder", &Vec3List::FitCylinder)
   ;
+  to_python_converter<std::pair<geom::Line3, Real>,
+  ost::PairToTupleConverter<geom::Line3, Real> >();
 }

modules/geom/src/vec3.cc

@@ -91,76 +91,96 @@ Plane Vec3List::GetODRPlane() const
   Vec3 normal=this->GetPrincipalAxes().GetRow(0);
   return Plane(origin,normal);
 }
-  
-Line3 Vec3List::FitCylinder(const Vec3& initial_direction, const Vec3& center) const
-{
-  Line3 axis=Line3(center,center+initial_direction), axis_old;
-  Real radius,res_sum_old,res_sum,delta_0=0.01,prec=0.0000001,err,norm,delta;
-  unsigned long n_step=1000, n_res=this->size();
-  Vec3 v,gradient;
-  
-  radius=0.0;
-  delta=delta_0;
-  for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-    radius+=geom::Distance(axis,(*i));
-  }
-  radius/=Real(n_res);
-  res_sum=0.0;
-  for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-    Real r=Distance(axis,(*i))-radius;
-    res_sum+=r*r;
-  }
-  unsigned long k=0;
-  err=2.0*prec;
-  while (err>prec && k<n_step) {
-    res_sum_old=res_sum;
-    axis_old=axis;
+
+std::pair<Line3, Real> Vec3List::FitCylinder(const Vec3& initial_direction) const
+  { 
+    Vec3 center=this->GetCenter();
+    Line3 axis=Line3(center,center+initial_direction), axis_old;
+    Real radius,res_sum_old,res_sum,delta_0=0.01,prec=0.0000001,err,norm,delta;
+    unsigned long n_step=1000, n_res=this->size();
+    Vec3 v,gradient_dir,gradient_center;
+    
     radius=0.0;
-    if (k>50) {
-      delta=delta_0*50.0*50.0/(k*k);
-    }
+    delta=delta_0;
     for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-      radius+=Distance(axis,(*i));
+      radius+=geom::Distance(axis,(*i));
     }
     radius/=Real(n_res);
-    for (int j=0; j!=3; ++j){
-      res_sum=0.0;
-      v=Vec3(0.0,0.0,0.0);
-      v[j]=delta;
-      axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()+v);
+    res_sum=0.0;
+    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+      Real r=Distance(axis,(*i))-radius;
+      res_sum+=r*r;
+    }
+    unsigned long k=0;
+    err=2.0*prec;
+    while (err>prec && k<n_step) {
+      res_sum_old=res_sum;
+      axis_old=axis;
+      //radius=0.0;
+      if (k>50) {
+        delta=delta_0*50.0*50.0/(k*k);
+      }
+      //for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+      //  radius+=Distance(axis,(*i));
+      //}
+      radius/=Real(n_res);
+      for (int j=0; j!=3; ++j){
+        res_sum=0.0;
+        v=Vec3(0.0,0.0,0.0);
+        v[j]=delta;
+        axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()+v);
+        radius=0.0;
+        for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+          radius+=Distance(axis,(*i));
+        }
+        radius/=Real(n_res);
+        for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+          Real r=Distance(axis,(*i))-radius;
+          res_sum+=r*r;
+        }
+        gradient_dir[j]=(res_sum-res_sum_old)/delta;
+      }
+      norm=Dot(gradient_dir,gradient_dir);
+      if (norm>1.) {
+        gradient_dir=Normalize(gradient_dir);
+      }
+      for (int j=0; j!=3; ++j){
+        res_sum=0.0;
+        v=Vec3(0.0,0.0,0.0);
+        v[j]=delta;
+        axis=Line3(axis_old.GetOrigin()+v,axis_old.GetOrigin()+axis_old.GetDirection()+v);
+        radius=0.0;
+        for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+          radius+=Distance(axis,(*i));
+        }
+        radius/=Real(n_res);
+        for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+          Real r=Distance(axis,(*i))-radius;
+          res_sum+=r*r;
+        }
+        gradient_center[j]=(res_sum-res_sum_old)/delta;
+      }
+      norm=Dot(gradient_center,gradient_center);
+      if (norm>1.) {
+        gradient_center=Normalize(gradient_center);
+      }      
+      axis=Line3(axis_old.GetOrigin()-50*delta*gradient_center,axis_old.GetOrigin()-50*delta*gradient_center+axis_old.GetDirection()-delta*gradient_dir);
       radius=0.0;
       for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
         radius+=Distance(axis,(*i));
       }
       radius/=Real(n_res);
+      res_sum=0.0;
       for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
         Real r=Distance(axis,(*i))-radius;
         res_sum+=r*r;
       }
-      gradient[j]=(res_sum-res_sum_old)/delta;
-    }
-    norm=Dot(gradient,gradient);
-    if (norm>1.) {
-      gradient=Normalize(gradient);
+      err=fabs((res_sum-res_sum_old)/float(n_res));
+      k++;
     }
-    axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()-delta*gradient);
-    radius=0.0;
-    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-      radius+=Distance(axis,(*i));
+    if (err>prec) {
+      std::cout<<"axis fitting did not converge"<<std::endl;
     }
-    radius/=Real(n_res);
-    res_sum=0.0;
-    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-      Real r=Distance(axis,(*i))-radius;
-      res_sum+=r*r;
-    }
-    err=fabs((res_sum-res_sum_old)/float(n_res));
-    k++;
+    return std::make_pair(axis,radius);
   }
-  if (err>prec) {
-    std::cout<<"axis fitting did not converge"<<std::endl;
-  }
-  return axis;
-}
-
 }

modules/geom/src/vec3.hh

@@ -316,10 +316,13 @@ public:
   Plane GetODRPlane() const;
 
   //This function fits a cylinder to the positions in Vec3List
-  //It takes as argument an initial guess for the direction and the geometric
-  //center of the atoms. The center is not changed during optimisation as the
-  //best fitting cylinder can be shown to have its axis pass through the geometric center
-  Line3 FitCylinder(const Vec3& initial_direction, const Vec3& center) const;
+  //It takes as argument an initial guess for the direction.
+  //The center is set to the geometric centero of the atoms
+  //and is not changed during optimisation as the best fitting cylinder
+  //can be shown to have its axis pass through the geometric center
+  //It returns a pair containing a line3, giving the direction of the Cylinder
+  //and a Real containing the radius.
+  std::pair<Line3, Real> FitCylinder(const Vec3& initial_direction) const;
 };
 } // ns geom
 

modules/mol/alg/pymod/structure_analysis.py

@@ -127,7 +127,7 @@ def CalculateHelixAxis(sele1):
     return
   eh=sele1.GetHandle()
   f=GetFrameFromEntity(eh)
-  return f.FitCylinder(sele1)
+  return f.FitCylinder(sele1)[0]
 
 
 def CalculateDistanceDifferenceMatrix(sele1,sele2):

modules/mol/alg/src/trajectory_analysis.cc

@@ -242,6 +242,7 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
       throw Error("EntityView is empty");
     }
     std::vector<unsigned long> indices_ca;
+    std::pair<geom::Line3, Real> cylinder;
     geom::Line3 axis;
     Real sign;
     directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
@@ -250,7 +251,8 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
     unsigned int n_atoms=indices_ca.size();
     for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
       CoordFramePtr frame=traj.GetFrame(i);
-      axis=frame->FitCylinder(indices_ca);
+      cylinder=frame->FitCylinder(indices_ca);
+      axis=cylinder.first;
       sign=geom::Dot(axis.GetDirection(),(*frame)[indices_ca[n_atoms-1]]-axis.GetOrigin());
       sign=sign/fabs(sign);
       directions.push_back(sign*axis.GetDirection());

modules/mol/base/pymod/export_coord_frame.cc

@@ -43,7 +43,7 @@ Real (CoordFrame::*get_min_dist)(const mol::EntityView&, const mol::EntityView&)
 Real (CoordFrame::*get_alpha)(const mol::EntityView&) const = &CoordFrame::GetAlphaHelixContent;
 geom::Line3 (CoordFrame::*get_odr_line)(const mol::EntityView&) const = &CoordFrame::GetODRLine;
 geom::Plane (CoordFrame::*get_odr_plane)(const mol::EntityView&) const = &CoordFrame::GetODRPlane;
-geom::Line3 (CoordFrame::*fit_cylinder)(const mol::EntityView&) const = &CoordFrame::FitCylinder;
+std::pair<geom::Line3,Real> (CoordFrame::*fit_cylinder)(const mol::EntityView&) const = &CoordFrame::FitCylinder;
 // TODO: move to geom
 geom::Line3 (CoordFrame::*get_odr_line2)() const = &geom::Vec3List::GetODRLine;
 

modules/mol/base/src/coord_frame.cc

@@ -235,7 +235,7 @@ namespace ost { namespace mol {
     return this->GetODRPlane(indices);
   }  
   
-  geom::Line3 CoordFrame::FitCylinder(std::vector<unsigned long>& indices_ca) const {
+  std::pair<geom::Line3, Real> CoordFrame::FitCylinder(std::vector<unsigned long>& indices_ca) const {
     geom::Vec3List atoms_pos_list;
     int n_atoms=indices_ca.size();
     atoms_pos_list.reserve(n_atoms);
@@ -243,7 +243,7 @@ namespace ost { namespace mol {
       atoms_pos_list.push_back((*this)[*i1]);
     }
     //Initial guess
-    geom::Vec3 initial_axis=geom::Vec3(0.0,0.0,0.0),center=geom::Vec3(0.0,0.0,0.0);
+    geom::Vec3 initial_axis=geom::Vec3(0.0,0.0,0.0);
     if (n_atoms<5) {
       initial_axis=atoms_pos_list[n_atoms-1]-atoms_pos_list[0];
     }
@@ -252,14 +252,10 @@ namespace ost { namespace mol {
         initial_axis+=(*(i+4))-(*i);
       }
     }
-    for (geom::Vec3List::const_iterator i=atoms_pos_list.begin(),e=atoms_pos_list.end(); i!=e; ++i) {
-      center+=(*i);
-    }
-    center/=atoms_pos_list.size();
-    return atoms_pos_list.FitCylinder(initial_axis,center);
+    return atoms_pos_list.FitCylinder(initial_axis);
   }
 
-  geom::Line3 CoordFrame::FitCylinder(const mol::EntityView& view1) const {
+  std::pair<geom::Line3, Real> CoordFrame::FitCylinder(const mol::EntityView& view1) const {
     CheckHandleValidity(view1);
     std::vector<unsigned long> indices_ca;
     GetCaIndices(view1, indices_ca);

modules/mol/base/src/coord_frame.hh

@@ -182,10 +182,10 @@ public:
     It is assumed that we fit an alpha-helix and that the CA atoms are oredered sequentially
     This is used for the initial guess of the helix axis
   */
-  geom::Line3 FitCylinder(std::vector<unsigned long>& indices_ca) const;
+  std::pair<geom::Line3, Real> FitCylinder(std::vector<unsigned long>& indices_ca) const;
 
   //! see FitCylinder(std::vector<unsigned long>&)
-  geom::Line3 FitCylinder(const mol::EntityView& view1) const;
+  std::pair<geom::Line3, Real> FitCylinder(const mol::EntityView& view1) const;
 
   /*!
     Returns the percentage of residues in the EntityView (segment) that are in an alpha-helix