diff --git a/modules/geom/pymod/export_vec3.cc b/modules/geom/pymod/export_vec3.cc
index 0df8101549767a5c3a7af43bdde4bd33da3b9dd2..0bd91ba68991fef784ee76508e569a9139c2a6ec 100644
--- a/modules/geom/pymod/export_vec3.cc
+++ b/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> >();
}
diff --git a/modules/geom/src/vec3.cc b/modules/geom/src/vec3.cc
index 8909ce8b8648879741c01965ef13961614223e6b..a36db7122a72a7f043a6a3cc160788987fd8088c 100644
--- a/modules/geom/src/vec3.cc
+++ b/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;
-}
-
}
diff --git a/modules/geom/src/vec3.hh b/modules/geom/src/vec3.hh
index 58e63b89258294462a34e9ab28311ddde59a9670..1626bb2d5064df6b1873e799a16b1564143f5dee 100644
--- a/modules/geom/src/vec3.hh
+++ b/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
diff --git a/modules/mol/alg/pymod/structure_analysis.py b/modules/mol/alg/pymod/structure_analysis.py
index a9f600bca0a16ea181e2f9780c085208d7b8a7d6..3e013e4791a975f7bbb38f7f15021da9992292ac 100644
--- a/modules/mol/alg/pymod/structure_analysis.py
+++ b/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):
diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc
index 96dd5fabb0599b77fb68fa4d9ca8dea94a8de313..359530f7cf63f18778acb40507c251df01a81567 100644
--- a/modules/mol/alg/src/trajectory_analysis.cc
+++ b/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());
diff --git a/modules/mol/base/pymod/export_coord_frame.cc b/modules/mol/base/pymod/export_coord_frame.cc
index 7ea3c2f3ad83467e148489c8db779a2e07d27aee..70ccadd92054bab3e54931c3a362c63e57290619 100644
--- a/modules/mol/base/pymod/export_coord_frame.cc
+++ b/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;
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index ff04e00d8ba4660298d52f8e04e23d863b7a0309..49c294fcee0fc4ac784e79cd02935d4d4f66e2c6 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/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);
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index f319e1a3beb325a697ca72223dce688098fdb240..33d24d1d402999d3c0e24bc4afa0a587d047e377 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/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