From 8e84db9719d7a9213f4414dfc64c3d1e815bbe27 Mon Sep 17 00:00:00 2001
From: Niklaus Johner <nij2003@med.cornell.edu>
Date: Tue, 27 Sep 2011 15:23:34 -0400
Subject: [PATCH] AnalyzeBestFitPlane function added to trajectory_analysis.
Added a function to extract the best fit plane from a trajectory,
and on the CoordFrame and Vec3List level.
---
modules/geom/src/vec3.cc | 8 +++++-
modules/geom/src/vec3.hh | 2 ++
.../alg/pymod/export_trajectory_analysis.cc | 1 +
modules/mol/alg/src/trajectory_analysis.cc | 25 +++++++++++++++++--
modules/mol/alg/src/trajectory_analysis.hh | 1 +
modules/mol/base/src/coord_frame.cc | 10 ++++++++
modules/mol/base/src/coord_frame.hh | 1 +
7 files changed, 45 insertions(+), 3 deletions(-)
diff --git a/modules/geom/src/vec3.cc b/modules/geom/src/vec3.cc
index bcc62d42a..ea07f3cca 100644
--- a/modules/geom/src/vec3.cc
+++ b/modules/geom/src/vec3.cc
@@ -75,5 +75,11 @@ Line3 Vec3List::GetODRLine()
return Line3(center,center+direction);
}
-
+Plane Vec3List::GetODRPlane()
+{
+ Vec3 origin=this->GetCenter();
+ Vec3 normal=this->GetPrincipalAxes().GetRow(0);
+ return Plane(origin,normal);
+}
+
}
diff --git a/modules/geom/src/vec3.hh b/modules/geom/src/vec3.hh
index a2d1d3cea..ac3458837 100644
--- a/modules/geom/src/vec3.hh
+++ b/modules/geom/src/vec3.hh
@@ -35,6 +35,7 @@ namespace geom {
class Vec2;
class Vec4;
class Line3;
+class Plane;
/// \brief Three dimensional vector class, using Real precision.
class DLLEXPORT_OST_GEOM Vec3:
@@ -218,6 +219,7 @@ public:
Mat3 GetPrincipalAxes() const;
Line3 GetODRLine();
+ Plane GetODRPlane();
Line3 FitCylinder(const Vec3 initial_direction, const Vec3 center);
};
diff --git a/modules/mol/alg/pymod/export_trajectory_analysis.cc b/modules/mol/alg/pymod/export_trajectory_analysis.cc
index 5671dbca0..7e42d3eee 100644
--- a/modules/mol/alg/pymod/export_trajectory_analysis.cc
+++ b/modules/mol/alg/pymod/export_trajectory_analysis.cc
@@ -48,6 +48,7 @@ void export_TrajectoryAnalysis()
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));
+ def("AnalyzeBestFitPlane", &AnalyzeBestFitPlane, (arg("traj"), arg("protein_segment"), arg("normals"), arg("origins"), arg("stride")=1));
def("AnalyzeHelicity", &AnalyzeHelicity, (arg("traj"), arg("protein_segment"), arg("stride")=1));
def("CreateMeanStructure", &CreateMeanStructure, (arg("traj"), arg("selection"), arg("from")=0, arg("to")=-1, arg("stride")=1 ));
}
diff --git a/modules/mol/alg/src/trajectory_analysis.cc b/modules/mol/alg/src/trajectory_analysis.cc
index c7222e599..c5cc63c68 100644
--- a/modules/mol/alg/src/trajectory_analysis.cc
+++ b/modules/mol/alg/src/trajectory_analysis.cc
@@ -272,7 +272,7 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
geom::Line3 axis;
directions.reserve(ceil(traj.GetFrameCount()/float(stride)));
centers.reserve(ceil(traj.GetFrameCount()/float(stride)));
- GetCaIndices(prot_seg, indices_ca);
+ GetIndices(prot_seg, indices_ca);
for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
CoordFramePtr frame=traj.GetFrame(i);
axis=frame->GetODRLine(indices_ca);
@@ -281,7 +281,28 @@ std::vector<Real> AnalyzeAromaticRingInteraction(const CoordGroupHandle& traj, c
}
return;
}
-
+
+ void AnalyzeBestFitPlane(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& normals,
+ geom::Vec3List& origins, unsigned int stride)
+ {
+ CheckHandleValidity(traj);
+ if (prot_seg.GetAtomCount()==0){
+ throw std::runtime_error("EntityView is empty");
+ }
+ std::vector<unsigned long> indices_ca;
+ geom::Plane best_plane;
+ normals.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ origins.reserve(ceil(traj.GetFrameCount()/float(stride)));
+ GetIndices(prot_seg, indices_ca);
+ for (size_t i=0; i<traj.GetFrameCount(); i+=stride) {
+ CoordFramePtr frame=traj.GetFrame(i);
+ best_plane=frame->GetODRPlane(indices_ca);
+ normals.push_back(best_plane.GetNormal());
+ origins.push_back(best_plane.GetOrigin());
+ }
+ return;
+ }
+
std::vector<Real> AnalyzeHelicity(const CoordGroupHandle& traj, const EntityView& prot_seg,
unsigned int stride)
{
diff --git a/modules/mol/alg/src/trajectory_analysis.hh b/modules/mol/alg/src/trajectory_analysis.hh
index d2ba49b6e..501638925 100644
--- a/modules/mol/alg/src/trajectory_analysis.hh
+++ b/modules/mol/alg/src/trajectory_analysis.hh
@@ -45,6 +45,7 @@ namespace ost { namespace mol { namespace alg {
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);
+ void DLLEXPORT_OST_MOL_ALG AnalyzeBestFitPlane(const CoordGroupHandle& traj, const EntityView& prot_seg, geom::Vec3List& normals, geom::Vec3List& origins, unsigned int stride=1);
EntityHandle DLLEXPORT_OST_MOL_ALG CreateMeanStructure(const CoordGroupHandle& traj, const EntityView& selection, int from=0, int to=-1, unsigned int stride=1);
std::vector<Real> DLLEXPORT_OST_MOL_ALG AnalyzeHelicity(const CoordGroupHandle& traj, const EntityView& prot_seg, unsigned int stride=1);
}}}//ns
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index bb70ac62a..c6d13aa78 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -274,6 +274,16 @@ namespace ost { namespace mol {
}
return atoms_pos_list.GetODRLine();
}
+
+ geom::Plane CoordFrame::GetODRPlane(std::vector<unsigned long>& indices_ca){
+ //Returns the normal to the best fit plane 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.GetODRPlane();
+ }
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
diff --git a/modules/mol/base/src/coord_frame.hh b/modules/mol/base/src/coord_frame.hh
index db9645dc6..dbd51c29b 100644
--- a/modules/mol/base/src/coord_frame.hh
+++ b/modules/mol/base/src/coord_frame.hh
@@ -61,6 +61,7 @@ public:
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::Plane GetODRPlane(std::vector<unsigned long>& indices_ca);
geom::Line3 FitCylinder(std::vector<unsigned long>& indices_ca);
Real GetAlphaHelixContent(std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_c,
std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n);
--
GitLab