use of de.bioforscher.jstructure.mathematics.Transformation in project jstructure by JonStargaryen.
the class SVDSuperimposer method align.
@Override
public StructureAlignmentResult align(AtomContainer reference, AtomContainer query) {
AtomContainer originalReference = reference;
AtomContainer originalCandidate = query;
Pair<GroupContainer, GroupContainer> atomContainerPair = AbstractAlignmentAlgorithm.comparableGroupContainerPair(reference, query, minimalSetOfAtomNames, maximalSetOfAtomNames);
reference = atomContainerPair.getLeft();
query = atomContainerPair.getRight();
// calculate centroids
double[] centroid1 = reference.calculate().centroid().getValue();
double[] centroid2 = query.calculate().centroid().getValue();
// center atoms
reference.calculate().center();
query.calculate().center();
// compose covariance matrix and calculate SVD
RealMatrix matrix1 = convertToMatrix(reference);
RealMatrix matrix2 = convertToMatrix(query);
RealMatrix covariance = matrix2.transpose().multiply(matrix1);
SingularValueDecomposition svd = new SingularValueDecomposition(covariance);
// R = (V * U')'
RealMatrix ut = svd.getU().transpose();
RealMatrix rotationMatrix = svd.getV().multiply(ut).transpose();
// check if reflection
if (new LUDecomposition(rotationMatrix).getDeterminant() < 0) {
RealMatrix v = svd.getV().transpose();
v.setEntry(2, 0, (0 - v.getEntry(2, 0)));
v.setEntry(2, 1, (0 - v.getEntry(2, 1)));
v.setEntry(2, 2, (0 - v.getEntry(2, 2)));
rotationMatrix = v.transpose().multiply(ut).transpose();
}
double[][] rotation = rotationMatrix.getData();
// calculate translation
double[] translation = LinearAlgebra.on(centroid1).subtract(LinearAlgebra.on(centroid2).multiply(rotation)).getValue();
logger.trace("rotation matrix\n{}\ntranslation vector\n{}", Arrays.deepToString(rotationMatrix.getData()), Arrays.toString(translation));
/* transform 2nd atom select - employ neutral translation (3D vector of zeros), because the atoms are already
* centered and calculate RMSD */
query.calculate().transform(new Transformation(rotation));
double rmsd = calculateRmsd(reference, query);
// return alignment
return new StructureAlignmentResult(originalReference, originalCandidate, query, rmsd, translation, rotation);
}
Aggregations