From bae6f2fb87a54036f5da970107a81d5808a35624 Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Wed, 6 Sep 2017 22:03:31 +0200
Subject: [PATCH] Make sure, that the iterative superposition algorithms
converge
---
core/src/superpose.cc | 4 +++-
core/src/superpose.hh | 35 +++++++++++++++++++----------------
loop/src/backbone.cc | 14 ++++++++++----
3 files changed, 32 insertions(+), 21 deletions(-)
diff --git a/core/src/superpose.cc b/core/src/superpose.cc
index 1d4fad62..7eb50402 100644
--- a/core/src/superpose.cc
+++ b/core/src/superpose.cc
@@ -547,7 +547,9 @@ void RigidBlocks(EMatX3& pos_one, EMatX3& pos_two,
distance_thresh, current_indices,
false);
- unique_blocks[current_indices] = t;;
+ // only add if iterative superposition converged => number of indices
+ // must be larger than 3
+ if(current_indices.size() >= 3) unique_blocks[current_indices] = t;
}
for(std::map<std::vector<uint>, geom::Mat4>::iterator i = unique_blocks.begin();
diff --git a/core/src/superpose.hh b/core/src/superpose.hh
index 36b5a4b5..d9a558da 100644
--- a/core/src/superpose.hh
+++ b/core/src/superpose.hh
@@ -46,35 +46,38 @@ std::pair<geom::Mat4,Real> Superposition(EMatX3& pos_one, EMatX3& pos_two,
void FillRMSDMatrix(EMatX3List& position_list, Real** data);
// ITERATIVE MIN RMSD SUPERPOSITION ALGORITHMS
+// You'll always know what indices are actually used for the final
+// superposition. If there are less than 3, you can expect, that the
+// algorithm didn't converge.
// Get minimal possible RMSD between points defined in pos_one and pos_two
Real SuperposedRMSD(EMatX3& pos_one, EMatX3& pos_two,
- uint max_iterations,
- Real distance_thresh,
- std::vector<uint>& indices,
- bool apply_superposition = false);
+ uint max_iterations,
+ Real distance_thresh,
+ std::vector<uint>& indices,
+ bool apply_superposition = false);
// Get transformation to superpose points in pos_one onto the ones in pos_two
geom::Mat4 MinRMSDSuperposition(EMatX3& pos_one, EMatX3& pos_two,
- uint max_iterations,
- Real distance_thresh,
- std::vector<uint>& indices,
- bool apply_superposition = false);
+ uint max_iterations,
+ Real distance_thresh,
+ std::vector<uint>& indices,
+ bool apply_superposition = false);
// Get both things...
std::pair<geom::Mat4,Real> Superposition(EMatX3& pos_one, EMatX3& pos_two,
- uint max_iterations,
- Real distance_thresh,
- std::vector<uint>& indices,
- bool apply_superposition = false);
+ uint max_iterations,
+ Real distance_thresh,
+ std::vector<uint>& indices,
+ bool apply_superposition = false);
void RigidBlocks(EMatX3& pos_one, EMatX3& pos_two,
uint window_length,
- uint max_iterations,
- Real distance_thresh,
- std::vector<std::vector<uint> >& indices,
- std::vector<geom::Mat4>& transformations);
+ uint max_iterations,
+ Real distance_thresh,
+ std::vector<std::vector<uint> >& indices,
+ std::vector<geom::Mat4>& transformations);
// Fill row of given Eigen Matrix with 3 entries of Vec3
diff --git a/loop/src/backbone.cc b/loop/src/backbone.cc
index 8376a87c..0a04aa12 100644
--- a/loop/src/backbone.cc
+++ b/loop/src/backbone.cc
@@ -1131,11 +1131,17 @@ geom::Mat4 BackboneList::GetTransformIterative(const BackboneList& other,
BackboneList::ExtractEigenNCACPos(other, pos_two);
}
- std::vector<uint> empty_index_vec;
- return core::MinRMSDSuperposition(pos_one, pos_two, max_iterations,
- distance_thresh, empty_index_vec, false);
+ std::vector<uint> index_vec;
+ geom::Mat4 t = core::MinRMSDSuperposition(pos_one, pos_two, max_iterations,
+ distance_thresh, index_vec, false);
+ if(index_vec.size() >= 3) {
+ return t;
+ } else {
+ // iterative superposition did not converge, let's use normal superposition
+ // as fallback
+ return core::MinRMSDSuperposition(pos_one, pos_two);
+ }
}
-
///////////////////////////////////////////////////////////////////////////////
// BackboneList private
void BackboneList::ConstructBackboneList_(const String& sequence,
--
GitLab