diff --git a/CMakeLists.txt b/CMakeLists.txt
index 413959a6b913767e97ea83e7a8a06720d4235952..fea7a6b8f88ff46894ae742808e5e2ec75d15aa0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -154,8 +154,6 @@ else()
set(_UBUNTU_LAYOUT OFF)
endif()
-add_definitions(-DEIGEN2_SUPPORT)
-
if (COMPOUND_LIB)
set(_COMP_LIB "${COMPOUND_LIB}")
if (NOT IS_ABSOLUTE "${COMPOUND_LIB}")
diff --git a/modules/geom/src/vec3.cc b/modules/geom/src/vec3.cc
index a36db7122a72a7f043a6a3cc160788987fd8088c..d55ab5d2ca8ec784ee4e637a40b4de2e6d200b46 100644
--- a/modules/geom/src/vec3.cc
+++ b/modules/geom/src/vec3.cc
@@ -60,7 +60,7 @@ Mat3 Vec3List::GetPrincipalAxes() const
Mat3 inertia=this->GetInertia();
EMat3 inertia_mat(inertia.Data());
- Eigen::SVD<EMat3> svd(inertia_mat);
+ Eigen::JacobiSVD<EMat3> svd(inertia_mat,Eigen::ComputeThinU);
EMat3 rot=svd.matrixU();
Mat3 axes(rot.data());
return axes;
diff --git a/modules/gfx/src/impl/cartoon_renderer.cc b/modules/gfx/src/impl/cartoon_renderer.cc
index f9eee510679686554c1f7bfac9fd6550f7e3f837..00894483df3e7c0e3433463810fb403a80a84a46 100644
--- a/modules/gfx/src/impl/cartoon_renderer.cc
+++ b/modules/gfx/src/impl/cartoon_renderer.cc
@@ -23,10 +23,7 @@
#include "cartoon_renderer.hh"
-#include <Eigen/Core>
-#include <Eigen/Array>
#include <Eigen/SVD>
-#include <Eigen/LU>
#include <ost/gfx/entity.hh>
#include <ost/gfx/impl/tabulated_trig.hh>
@@ -138,7 +135,7 @@ namespace {
A.row(i)=to_eigen(points[i]-cen);
}
- Eigen::SVD<EMatX> svd(A);
+ Eigen::JacobiSVD<EMatX> svd(A,Eigen::ComputeThinV);
EMatX V=svd.matrixV();
geom::Vec3 ax(V(0,0),V(1,0),V(2,0));
diff --git a/modules/mol/alg/src/superpose_frames.cc b/modules/mol/alg/src/superpose_frames.cc
index f5cbba61073416b7c47dd05393c9104f4ee0b0b2..8c7c789a7dcbf664262e05703469bc3cb5bc84c8 100644
--- a/modules/mol/alg/src/superpose_frames.cc
+++ b/modules/mol/alg/src/superpose_frames.cc
@@ -17,10 +17,7 @@
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------
-#include <Eigen/Core>
-#include <Eigen/Array>
-#include <Eigen/SVD>
-#include <Eigen/LU>
+#include <Eigen/Eigenvalues>
#include <ost/message.hh>
#include <ost/mol/mol.hh>
#include <ost/mol/alg/superpose_frames.hh>
@@ -100,7 +97,7 @@ void AddSuperposedFrame(CoordGroupHandle& superposed, EMatX3& ref_mat,EMatX3& re
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);
+ Eigen::JacobiSVD<EMat3> svd(frame_centered*ref_centered,Eigen::ComputeThinU | Eigen::ComputeThinV);
EMat3 matrixVT=svd.matrixV().transpose();
//determine rotation
Real detv=matrixVT.determinant();
diff --git a/modules/mol/alg/src/svd_superpose.cc b/modules/mol/alg/src/svd_superpose.cc
index fd789b867c8f5f1383319b04628dbe2452c4868f..c8c068c0be576bab51592b701f83c07d598670f0 100644
--- a/modules/mol/alg/src/svd_superpose.cc
+++ b/modules/mol/alg/src/svd_superpose.cc
@@ -21,10 +21,8 @@
#include <iostream>
#include <boost/bind.hpp>
-#include <Eigen/Core>
-#include <Eigen/Array>
#include <Eigen/SVD>
-#include <Eigen/LU>
+#include <Eigen/Geometry>
#include <ost/base.hh>
#include <ost/geom/vec3.hh>
@@ -94,7 +92,7 @@ Real calc_rmsd_for_ematx(const EMatX& atoms1,
const EMatX& atoms2,
const EMat4& transformation)
{
- EMatX transformed_atoms1 = EMatX::Zero(atoms1.rows(), 3);
+ EMatX transformed = EMatX::Zero(atoms1.rows(), 3);
EMatX vector = EMatX::Zero(4,1);
EMatX transformed_vector = EMatX::Zero(4,1);
@@ -103,17 +101,12 @@ Real calc_rmsd_for_ematx(const EMatX& atoms1,
for(int i=0;i<atoms1.rows();++i){
vector.block<3,1>(0,0)=atoms1.block<1,3>(i,0).transpose();
transformed_vector = transformation*vector;
- transformed_atoms1.block<1,3>(i,0)=transformed_vector.block<3,1>(0,0).transpose();
+ transformed.block<1,3>(i,0)=transformed_vector.block<3,1>(0,0).transpose();
}
- EMatX diff = EMatX::Zero(atoms1.rows(),atoms1.cols());
- EMatX squared_dist = EMatX::Zero(atoms1.rows(),1);
-
- diff = transformed_atoms1-atoms2;
- squared_dist = (diff.cwise()*diff).rowwise().sum();
-
- return sqrt(squared_dist.sum()/squared_dist.rows());
-
+ transformed = transformed - atoms2;
+ transformed = transformed.array()*transformed.array();
+ return sqrt(transformed.rowwise().sum().sum()/atoms1.rows());
}
Real CalculateRMSD(const mol::EntityView& ev1,
@@ -204,7 +197,8 @@ public:
EMatX TransformEMatX(const EMatX& mat, const EMat4& transformation) const;
- std::vector<int> CreateMatchingSubsets(const EMatX& atoms, const EMatX& atoms_ref, Real distance_threshold) const;
+ void CreateMatchingSubsets(std::vector<int>& matching_subset, const EMatX& atoms,
+ const EMatX& atoms_ref, Real distance_threshold) const;
SuperpositionResult IterativeMinimize(int ncycles, Real distance_threshold) const;
@@ -287,7 +281,7 @@ SuperpositionResult MeanSquareMinimizerImpl::IterativeMinimize(int max_cycles, R
for(;cycles<max_cycles;++cycles){
transformed_atoms = this->TransformEMatX(atoms1_, transformation_matrix);
- matching_indices = this->CreateMatchingSubsets(transformed_atoms, atoms2_, distance_threshold);
+ this->CreateMatchingSubsets(matching_indices, transformed_atoms, atoms2_, distance_threshold);
if(matching_indices.size()<3){
std::stringstream ss;
@@ -312,7 +306,7 @@ SuperpositionResult MeanSquareMinimizerImpl::IterativeMinimize(int max_cycles, R
diff = transformation_matrix_old-transformation_matrix;
- if(diff.cwise().abs().sum()<0.0001){
+ if(diff.array().abs().sum()<0.0001){
++cycles;
break;
}
@@ -348,24 +342,25 @@ EMatX MeanSquareMinimizerImpl::TransformEMatX(const EMatX& mat, const EMat4& tra
return transformed_mat;
}
-std::vector<int> MeanSquareMinimizerImpl::CreateMatchingSubsets(const EMatX& atoms, const EMatX& atoms_ref, Real distance_threshold) const{
+void MeanSquareMinimizerImpl::CreateMatchingSubsets(std::vector<int>& matching_subset,
+ const EMatX& atoms, const EMatX& atoms_ref,
+ Real distance_threshold) const{
- std::vector<int> return_vec;
+ matching_subset.clear();
EMatX diff = EMatX::Zero(atoms.rows(),atoms.cols());
EMatX dist = EMatX::Zero(atoms.rows(),1);
diff = atoms-atoms_ref;
- dist = (diff.cwise()*diff).rowwise().sum();
- dist = dist.cwise().sqrt();
+ dist = (diff.array()*diff.array()).rowwise().sum();
+
+ Real squared_dist_threshold = distance_threshold * distance_threshold;
for(int i=0; i<dist.rows(); ++i){
- if(dist(i,0) <= distance_threshold){
- return_vec.push_back(i);
+ if(dist(i,0) <= squared_dist_threshold){
+ matching_subset.push_back(i);
}
}
-
- return return_vec;
}
@@ -381,7 +376,7 @@ SuperpositionResult MeanSquareMinimizerImpl::Minimize(const EMatX& atoms, const
EMatX atoms_ref_shifted = MatrixShiftedBy(atoms_ref, avg_ref).transpose();
// determine rotational component
- Eigen::SVD<EMat3> svd(atoms_ref_shifted*atoms_shifted);
+ Eigen::JacobiSVD<EMat3> svd(atoms_ref_shifted*atoms_shifted,Eigen::ComputeThinU | Eigen::ComputeThinV);
EMatX matrixVT=svd.matrixV().transpose();
//determine rotation
diff --git a/modules/mol/base/src/bounding_box.cc b/modules/mol/base/src/bounding_box.cc
index 4b70f78172ecdfde25af3d612944823423b46138..85c40ed1acbb0832893da76f57af96815a3829cc 100644
--- a/modules/mol/base/src/bounding_box.cc
+++ b/modules/mol/base/src/bounding_box.cc
@@ -22,7 +22,7 @@
*/
#include <limits>
-#include <Eigen/QR>
+#include <Eigen/Eigenvalues>
#include <ost/mol/mol.hh>
#include <ost/mol/bounding_box.hh>
namespace ost { namespace mol {