diff --git a/modules/mol/alg/molalg.rst b/modules/mol/alg/molalg.rst
index 00735f6def328ed899de4005fb2f803cd23f0206..2daf3ca7597bca27d86b9d86d53e7c67a33b5bb6 100644
--- a/modules/mol/alg/molalg.rst
+++ b/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
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 240a029a42dd523e8a4dfab72be40fe77bd74fa7..67f4fbf9e13454d105b0f69844297b06a1508309 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/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));
}
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index 17dc9acbf2a5d32ac17e62b0e50af18f35122b2f..e613a9c735cff3588598ae612eec5ff0ce805202 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/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)
diff --git a/modules/mol/alg/src/superpose_frames.cc b/modules/mol/alg/src/superpose_frames.cc
new file mode 100644
index 0000000000000000000000000000000000000000..560293c166dcec0f31986493dcea0518d1e1929f
--- /dev/null
+++ b/modules/mol/alg/src/superpose_frames.cc
@@ -0,0 +1,170 @@
+//------------------------------------------------------------------------------
+// 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;
+}
+
+}}}
diff --git a/modules/mol/alg/src/superpose_frames.hh b/modules/mol/alg/src/superpose_frames.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9416c5a9236d9bcd090eccedc24b8655ee15a2a5
--- /dev/null
+++ b/modules/mol/alg/src/superpose_frames.hh
@@ -0,0 +1,34 @@
+//------------------------------------------------------------------------------
+// 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