added function to efficiently superpose coord group frames

modules/mol/alg/molalg.rst

@@ -22,4 +22,22 @@
   naming of the atoms is ambigous. For these residues, the overlap of both 
   possible solutions to the fixed atoms, that is, everything that is not 
   ambigous is calculated. The solution that gives higher overlap is then used to 
-  calculate the actual overlap score.
\ No newline at end of file
+  calculate the actual overlap score.
+  
+
+.. function:: SuperposeFrames(frames, sel, from=0, to=-1, ref=-1)
+
+  This function superposes the frames of the given coord group and returns them
+  as a new coord group.
+  
+  :param frames: The source coord group.
+  :param from: index of the first frame
+  :param to: index of the last frame plus one. If set to -1, the value is set to 
+     the number of frames in the coord group
+  :param sel: An entity view containing the selection of atoms to be used for     
+    superposition. If set to an invalid view, all atoms in the coord group are 
+    used.
+  :param ref: The index of the reference frame to use for superposition. If set 
+     to -1, the each frame is superposed to the previous frame.
+     
+  :returns: A newly created coord group containing the superposed frames.
\ No newline at end of file

modules/mol/alg/pymod/wrap_mol_alg.cc

@@ -23,7 +23,7 @@
 #include <boost/python.hpp>
 #include <ost/config.hh>
 #include <ost/mol/alg/local_dist_test.hh>
-
+#include <ost/mol/alg/superpose_frames.hh>
 using namespace boost::python;
 
 void export_svdSuperPose();
@@ -40,4 +40,7 @@ BOOST_PYTHON_MODULE(_mol_alg)
   #endif
   
   def("LocalDistTest", &ost::mol::alg::LocalDistTest);
+  def("SuperposeFrames", &ost::mol::alg::SuperposeFrames, 
+      (arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0, 
+       arg("end")=-1, arg("ref")=-1));  
 }

modules/mol/alg/src/CMakeLists.txt

@@ -3,12 +3,14 @@ set(OST_MOL_ALG_HEADERS
   module_config.hh
   sec_structure_segments.hh
   local_dist_test.hh
+  superpose_frames.hh
 )
 
 set(OST_MOL_ALG_SOURCES
   svd_superpose.cc
   sec_structure_segments.cc
   local_dist_test.cc
+  superpose_frames.cc
 )
 
 set(MOL_ALG_DEPS mol)

modules/mol/alg/src/superpose_frames.cc

+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+//------------------------------------------------------------------------------
+
+#include <Eigen/Core>
+#include <Eigen/Array>
+#include <Eigen/SVD>
+#include <Eigen/LU>
+#include <ost/mol/mol.hh>
+#include <ost/mol/alg/superpose_frames.hh>
+
+namespace ost { namespace mol { namespace alg {
+
+
+typedef Eigen::Matrix<Real, 3, 3> EMat3;
+typedef Eigen::Matrix<Real, 1, 3> ECVec3;
+typedef Eigen::Matrix<Real, 3, 1> ERVec3;
+typedef Eigen::Matrix<Real, 3, Eigen::Dynamic> EMat3X;
+typedef Eigen::Matrix<Real, Eigen::Dynamic, 3> EMatX3;
+
+
+inline geom::Vec3 rvec_to_gvec(const ERVec3 &vec) {
+  return *reinterpret_cast<const geom::Vec3*>(&vec);
+}
+
+inline geom::Vec3 cvec_to_gvec(const ECVec3 &vec) {
+  return *reinterpret_cast<const geom::Vec3*>(&vec);
+}
+
+inline geom::Mat3 emat_to_gmat(const EMat3 &mat)
+{
+  return *reinterpret_cast<const geom::Mat3*>(&mat);
+}
+
+inline ERVec3 gvec_to_rvec(const geom::Vec3 &vec)
+{
+  return *reinterpret_cast<const ERVec3*>(&vec);
+}
+
+inline ECVec3 gvec_to_cvec(const geom::Vec3 &vec)
+{
+  return *reinterpret_cast<const ECVec3*>(&vec);
+}
+
+inline EMat3X row_sub(const EMat3X& m, const ERVec3& s)
+{
+  EMat3X r(m);
+  for (int col=0;  col<m.cols();++col) {
+    r.col(col)-=s;
+  }
+  return r;
+}
+
+inline EMatX3 col_sub(const EMatX3& m, const ECVec3& s)
+{
+  EMatX3 r(m);
+  for (int row=0;  row<m.rows();++row) {
+    r.row(row)-=s;
+  }
+  return r;
+}
+
+CoordGroupHandle SuperposeFrames(CoordGroupHandle cg, EntityView sel,
+                                 int begin, int end, int ref)
+{
+  int real_end=end==-1 ? cg.GetFrameCount() : end;
+  CoordFramePtr ref_frame;
+  std::vector<unsigned long> indices;
+  if (!sel.IsValid()) {
+    indices.reserve(cg.GetAtomCount());
+    for (size_t i=0;i<cg.GetAtomCount(); ++i) {
+      indices.push_back(i);
+    }
+  } else {
+    AtomViewList atoms=sel.GetAtomList();
+    indices.reserve(atoms.size());
+    for (AtomViewList::const_iterator i=atoms.begin(), 
+         e=atoms.end(); i!=e; ++i) {
+      indices.push_back(i->GetIndex());
+    }
+  }
+  CoordGroupHandle superposed=CreateCoordGroup(cg.GetEntity().GetAtomList());
+  EMatX3 ref_mat=EMatX3::Zero(indices.size(), 3);
+  EMatX3 ref_centered=EMatX3::Zero(indices.size(), 3);
+  EMat3X frame_centered=EMat3X::Zero(3, indices.size());
+  EMat3X frame_mat=EMat3X::Zero(3, indices.size());
+  ECVec3 ref_center;
+  ERVec3 frame_center;
+  if (ref!=-1) {
+    ref_frame=cg.GetFrame(ref);
+    for (size_t i=0; i<indices.size(); ++i) {
+      ref_mat.row(i)=gvec_to_cvec((*ref_frame)[indices[i]]);
+    }
+    ref_center=ref_mat.colwise().sum()/ref_mat.rows();
+    ref_centered=col_sub(ref_mat, ref_center);
+  } else {
+    ref_frame=cg.GetFrame(begin);
+  }
+  for (int i=begin; i<real_end; ++i) {
+    if (ref==i || (ref==-1 && i==begin)) {
+      superposed.AddFrame(*ref_frame);
+      ref_frame=superposed.GetFrame(superposed.GetFrameCount()-1);
+      continue;
+    }
+    if (ref==-1) {
+      for (size_t j=0; j<indices.size(); ++j) {
+        ref_mat.row(j)=gvec_to_cvec((*ref_frame)[indices[j]]);
+      }
+      ref_center=ref_mat.colwise().sum()/ref_mat.rows();
+      ref_centered=col_sub(ref_mat, ref_center);
+    }
+    CoordFramePtr frame=cg.GetFrame(i);
+    for (size_t j=0; j<indices.size(); ++j) {
+      frame_mat.col(j)=gvec_to_rvec((*frame)[indices[j]]);
+    }
+    std::vector<geom::Vec3> frame_data=*frame;
+    frame_center=frame_mat.rowwise().sum()/frame_mat.cols();
+    frame_centered=row_sub(frame_mat, frame_center);
+    //single value decomposition
+    Eigen::SVD<EMat3> svd(frame_centered*ref_centered);
+    EMat3 matrixVT=svd.matrixV().transpose();
+    //determine rotation
+    Real detv=matrixVT.determinant();
+    Real dett=svd.matrixU().determinant();
+    Real det=detv*dett;
+    EMat3 e_rot;
+    if (det<0) {
+      EMat3 tmat=EMat3::Identity();
+      tmat(2,2)=-1;
+      e_rot=(svd.matrixU()*tmat)*matrixVT;
+    }else{
+      e_rot=svd.matrixU()*matrixVT;
+    }
+    //  prepare rmsd calculation
+    geom::Vec3 shift=rvec_to_gvec(ref_center);
+    geom::Vec3 com_vec=-cvec_to_gvec(frame_center);
+    geom::Mat3 rot=emat_to_gmat(e_rot);
+    geom::Mat4 mat4_com, mat4_rot, mat4_shift;
+    mat4_rot.PasteRotation(rot);
+    mat4_shift.PasteTranslation(shift);
+    mat4_com.PasteTranslation(com_vec);
+    geom::Mat4 tf=geom::Mat4(mat4_shift*mat4_rot*mat4_com);
+    for (std::vector<geom::Vec3>::iterator c=frame_data.begin(), 
+         e2=frame_data.end(); c!=e2; ++c) {
+      *c=geom::Vec3(tf*geom::Vec4(*c));
+    }
+    superposed.AddFrame(frame_data);
+    if (ref==-1) {
+      ref_frame=superposed.GetFrame(superposed.GetFrameCount()-1);
+    }
+  }
+  return superposed;
+}
+
+}}}

modules/mol/alg/src/superpose_frames.hh

+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+//------------------------------------------------------------------------------
+#ifndef OST_MOL_ALG_SUPERPOSE_FRAMES_HH
+#define OST_MOL_ALG_SUPERPOSE_FRAMES_HH
+
+#include <ost/mol/coord_group.hh>
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/alg/module_config.hh>
+namespace ost { namespace mol { namespace alg {
+
+/// \brief returns a superposed version of coord group
+CoordGroupHandle DLLEXPORT_OST_MOL_ALG SuperposeFrames(CoordGroupHandle cg, 
+                                                       EntityView sel,
+                                                       int begin=0, int end=-1, 
+                                                       int ref=-1);
+}}}
+
+#endif