Added function to fit a line to a Vec3List and to fit an axis to

an alpha helix (actually fits a cylinder). I also added the axis fitting
on the CoordGroup level.

modules/geom/pymod/export_vec3.cc

@@ -92,5 +92,7 @@ void export_Vec3()
     .add_property("center", &Vec3List::GetCenter)
     .add_property("inertia", &Vec3List::GetInertia)
     .add_property("principal_axes", &Vec3List::GetPrincipalAxes)
+    .def("GetODRLine", &Vec3List::GetODRLine)
+    .def("FitCylinder", &Vec3List::FitCylinder)
   ;
 }

modules/geom/src/composite3_op.cc

@@ -178,5 +178,50 @@ bool IsInSphere(const Sphere& s, const Vec3& v){
   return Length(s.GetOrigin()-v)<=s.GetRadius();
 }
 
+Line3 Vec3List::FitCylinder(){
+  Line3 axis=this->GetODRLine(),axis_old;
+  Real radius,res_sum_old,res_sum,delta=0.01,prec=0.001,err;
+  unsigned long n_step=1000, n_res=this->size();
+  Vec3 v,gradient;
+  radius=0.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) {
+    res_sum+=pow(Distance(axis,(*i))-radius,2.);
+  }
+  unsigned long k=0;
+  err=2.0*prec;
+  while (err>prec and k<n_step) {
+    res_sum_old=res_sum;
+    axis_old=axis;
+    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);
+      for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+        res_sum+=pow(Distance(axis,(*i))-radius,2.);
+      }
+      gradient[j]=res_sum-res_sum_old;
+    }
+    gradient=Normalize(gradient);
+    axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()-delta*gradient);
+    res_sum=0.0;
+    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+      res_sum+=pow(Distance(axis,(*i))-radius,2.);
+    }
+    err=fabs((res_sum-res_sum_old)/float(n_res));
+    k++;
+  }
+  return axis;
+}
+
 } // ns
 

modules/geom/src/vec3.cc

@@ -68,5 +68,12 @@ Vec3 Vec3List::GetCenter() const
   return center/=this->size();
 }
 
+Line3 Vec3List::GetODRLine()
+{
+  Vec3 center=this->GetCenter();
+  Vec3 direction=this->GetPrincipalAxes().GetRow(2);
+  return Line3(center,center+direction);
+}
+
 
 }

modules/geom/src/vec3.hh

@@ -34,7 +34,7 @@ namespace geom {
 // fw decl
 class Vec2;
 class Vec4;
-
+class Line3;
 
 /// \brief Three dimensional vector class, using Real precision.
 class DLLEXPORT_OST_GEOM Vec3:
@@ -193,6 +193,7 @@ inline std::ostream& operator<<(std::ostream& os, const Vec3& v)
 #include <ost/geom/vec2.hh>
 #include <ost/geom/vec4.hh>
 #include <ost/geom/mat3.hh>
+#include <ost/geom/composite3.hh>
 
 namespace geom {
 
@@ -216,6 +217,9 @@ public:
   Vec3 GetCenter() const;
   
   Mat3 GetPrincipalAxes() const;
+  Line3 GetODRLine();
+  Line3 FitCylinder();
+
 };
 
 

modules/mol/alg/pymod/export_trajectory_analysis.cc

@@ -46,4 +46,6 @@ void export_TrajectoryAnalysis()
   def("AnalyzeMinDistance", &AnalyzeMinDistance, (arg("traj"), arg("view1"), arg("view2"), arg("stride")=1));
   def("AnalyzeMinDistanceBetwCenterOfMassAndView", &AnalyzeMinDistanceBetwCenterOfMassAndView, (arg("traj"), arg("view_cm"), arg("view_atoms"), arg("stride")=1));
   def("AnalyzeAromaticRingInteraction", &AnalyzeAromaticRingInteraction, (arg("traj"), arg("view_ring1"), arg("view_ring2"), arg("stride")=1));
+  def("AnalyzeAlphaHelixAxis", &AnalyzeAlphaHelixAxis, (arg("traj"), arg("protein_segment"), arg("directions"), arg("centers"), arg("stride")=1));
+  def("AnalyzeBestFitLine", &AnalyzeBestFitLine, (arg("traj"), arg("protein_segment"), arg("directions"), arg("centers"), arg("stride")=1));
 }

modules/mol/alg/src/trajectory_analysis.cc

@@ -230,4 +230,51 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
   }
   return dist;  
   }
+
+  void AnalyzeAlphaHelixAxis(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions,
+                             geom::Vec3List& centers, unsigned int stride)
+  {
+    CheckHandleValidity(traj);
+    if (prot_seg.GetAtomCount()==0){
+      throw std::runtime_error("EntityView is empty");
+    }
+    std::vector<unsigned long> indices_ca;
+    geom::Line3 axis;
+    directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
+    centers.reserve(ceil(traj.GetFrameCount()/float(stride)));
+    GetCaIndices(prot_seg, indices_ca);
+    for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+      CoordFramePtr frame=traj.GetFrame(i);
+      axis=frame->FitCylinder(indices_ca);
+      directions.push_back(axis.GetDirection());
+      centers.push_back(axis.GetOrigin());
+    }
+    return;
+  }
+ 
+  //std::vector<geom::Line3> AnalyzeBestFitLine(const CoordGroupHandle& traj, const EntityView& prot_seg,
+  //                                               unsigned int stride)
+  void AnalyzeBestFitLine(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions,
+                          geom::Vec3List& centers, unsigned int stride)
+  {
+    CheckHandleValidity(traj);
+    if (prot_seg.GetAtomCount()==0){
+      throw std::runtime_error("EntityView is empty");
+    }
+    std::vector<unsigned long> indices_ca;
+    geom::Line3 axis;
+    directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
+    centers.reserve(ceil(traj.GetFrameCount()/float(stride)));
+    GetCaIndices(prot_seg, indices_ca);
+    for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+      CoordFramePtr frame=traj.GetFrame(i);
+      axis=frame->GetODRLine(indices_ca);
+      directions.push_back(axis.GetDirection());
+      centers.push_back(axis.GetOrigin());
+    }
+    return;
+  }
+ 
+
+
 }}} //ns

modules/mol/alg/src/trajectory_analysis.hh

@@ -43,5 +43,7 @@ namespace ost { namespace mol { namespace alg {
   std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeMinDistance(const CoordGroupHandle& traj, const EntityView& view1, const EntityView& view2,unsigned int stride=1);
   std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeMinDistanceBetwCenterOfMassAndView(const CoordGroupHandle& traj, const EntityView& view_cm, const EntityView& view_atoms,unsigned int stride=1);
   std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, const EntityView& view_ring1, const EntityView& view_ring2,unsigned int stride=1);
+  void  DLLEXPORT_OST_MOL_ALG AnalyzeAlphaHelixAxis(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions, geom::Vec3List& centers, unsigned int stride=1);
+  void  DLLEXPORT_OST_MOL_ALG AnalyzeBestFitLine(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& directions, geom::Vec3List& centers, unsigned int stride=1);
 }}}//ns
 #endif

modules/mol/base/pymod/export_coord_frame.cc

@@ -46,6 +46,7 @@ void export_CoordFrame()
     .def("GetDistanceBetwCenterOfMass", get_dist_cm)
     .def("GetRMSD",get_rmsd)
     .def("GetMinDistance",get_min_dist)
+    .def("GetODRLine",&geom::Vec3List::GetODRLine)
   ;
   def("CreateCoordFrame",CreateCoordFrame);
 }

modules/mol/base/src/coord_frame.cc

@@ -84,6 +84,22 @@ namespace ost { namespace mol {
       indices.push_back(i->GetIndex());
     }
   }
+
+  void GetCaIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca){
+    mol::AtomView ca;
+    ResidueViewList residues=segment.GetResidueList();
+    indices_ca.reserve(residues.size());
+    //for (ResidueViewList::const_iterator res=residues.begin()+1,
+    //     e=residues.end(); res!=e-1; ++res){
+    //  if (!InSequence((*res).GetHandle(),(*(res+1)).GetHandle())) throw std::runtime_error("Residues are not in a continuous sequence");
+    //}
+    for (ResidueViewList::const_iterator res=residues.begin(),
+         e=residues.end(); res!=e; ++res) {
+      ca=res->FindAtom("CA");
+      CheckHandleValidity(ca);
+      indices_ca.push_back(ca.GetIndex());
+      }
+  }
   
   void GetMasses(const EntityView& sele, std::vector<Real>& masses){
   // This function returns a vector containing the atom masses of the atoms in an EntityView
@@ -218,5 +234,29 @@ namespace ost { namespace mol {
     GetIndicesAndMasses(view_cm,indices_cm,masses_cm);
     return this->GetMinDistBetwCenterOfMassAndView(indices_cm,masses_cm,indices_atoms);
   }
+
+  geom::Line3 CoordFrame::GetODRLine(std::vector<unsigned long>& indices_ca){
+  //Returns the best fit line to atoms with indices in indices_ca
+    geom::Vec3List atoms_pos_list;
+    atoms_pos_list.reserve(indices_ca.size());
+    for (std::vector<unsigned long>::const_iterator i1=indices_ca.begin(),e=indices_ca.end(); i1!=e; ++i1) {
+      atoms_pos_list.push_back((*this)[*i1]);
+    }
+    return atoms_pos_list.GetODRLine();
+  }
+ 
+  geom::Line3 CoordFrame::FitCylinder(std::vector<unsigned long>& indices_ca){
+  //Returns a line which is the axis of a fitted cylinder to the atoms with indices given in indices_ca
+    geom::Vec3List atoms_pos_list;
+    atoms_pos_list.reserve(indices_ca.size());
+    for (std::vector<unsigned long>::const_iterator i1=indices_ca.begin(),e=indices_ca.end(); i1!=e; ++i1) {
+      atoms_pos_list.push_back((*this)[*i1]);
+    }
+    return atoms_pos_list.FitCylinder();
+  }
+ 
+ 
+
+
 }}//ns
 

modules/mol/base/src/coord_frame.hh

@@ -60,13 +60,16 @@ public:
   Real GetMinDistBetwCenterOfMassAndView(std::vector<unsigned long>& indices_cm, std::vector<Real>& masses_cm,
                                          std::vector<unsigned long>& indices_atoms);
   Real GetMinDistBetwCenterOfMassAndView(const mol::EntityView& view_cm, const mol::EntityView& view_atoms);
+  geom::Line3 GetODRLine(std::vector<unsigned long>& indices_ca);
+  geom::Line3 FitCylinder(std::vector<unsigned long>& indices_ca);
 };
   
   void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices);
   void GetMasses(const EntityView& sele, std::vector<Real>& masses);
   void GetIndicesAndMasses(const EntityView& sele, std::vector<unsigned long>& indices,std::vector<Real>& masses);
   void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos);
-  
+  void GetCaIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca);
+
 typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
 typedef std::vector<CoordFramePtr> CoordFrameList;