001/*
002 *                  BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 * Created on Dec 4, 2005
021 *
022 */
023package org.biojava.nbio.structure.geometry;
024
025import javax.vecmath.Matrix4d;
026import javax.vecmath.Point3d;
027import javax.vecmath.Vector3d;
028
029import org.biojava.nbio.structure.jama.Matrix;
030import org.biojava.nbio.structure.jama.SingularValueDecomposition;
031
032/**
033 * A class that calculates the superposition between two sets of points using an
034 * SVD Matrix Decomposition. It was introduced by Wolfgang Kabsch, hence the
035 * alternative name Kabsh algorithm. Inspired by the biopython SVDSuperimposer
036 * class.
037 *
038 * @author Andreas Prlic
039 * @author Aleix Lafita
040 * @since 1.5
041 * @version %I% %G%
042 *
043 */
044public class SuperPositionSVD extends SuperPositionAbstract {
045
046        /**
047         * Constructor for the SVD superposition algorithm.
048         *
049         * @param centered
050         *            true if the point arrays are centered at the origin (faster),
051         *            false otherwise
052         */
053        public SuperPositionSVD(boolean centered) {
054                super(centered);
055        }
056
057        @Override
058        public Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
059
060                checkInput(fixed, moved);
061
062                Point3d cena = CalcPoint.centroid(fixed);
063                Point3d cenb = CalcPoint.centroid(moved);
064
065                double[][] centAcoords = new double[][] { { cena.x, cena.y, cena.z } };
066                Matrix centroidA = new Matrix(centAcoords);
067
068                double[][] centBcoords = new double[][] { { cenb.x, cenb.y, cenb.z } };
069                Matrix centroidB = new Matrix(centBcoords);
070
071                // center at centroid
072
073                cena.negate();
074                cenb.negate();
075
076                Point3d[] ats1 = CalcPoint.clonePoint3dArray(fixed);
077                CalcPoint.translate(new Vector3d(cena), ats1);
078
079                Point3d[] ats2 = CalcPoint.clonePoint3dArray(moved);
080                CalcPoint.translate(new Vector3d(cenb), ats2);
081
082                double[][] coordSet1 = new double[ats1.length][3];
083                double[][] coordSet2 = new double[ats2.length][3];
084
085                // copy the atoms into the internal coords;
086                for (int i = 0; i < ats1.length; i++) {
087                        coordSet1[i] = new double[3];
088                        ats1[i].get(coordSet1[i]);
089                        coordSet2[i] = new double[3];
090                        ats2[i].get(coordSet2[i]);
091                }
092
093                // now this is the bridge to the Jama package:
094                Matrix a = new Matrix(coordSet1);
095                Matrix b = new Matrix(coordSet2);
096
097                // # correlation matrix
098
099                Matrix b_trans = b.transpose();
100                Matrix corr = b_trans.times(a);
101
102                SingularValueDecomposition svd = corr.svd();
103
104                Matrix u = svd.getU();
105                // v is alreaady transposed ! difference to numermic python ...
106                Matrix vt = svd.getV();
107
108                Matrix vt_orig = (Matrix) vt.clone();
109                Matrix u_transp = u.transpose();
110
111                Matrix rot_nottrans = vt.times(u_transp);
112                Matrix rot = rot_nottrans.transpose();
113
114                // check if we have found a reflection
115
116                double det = rot.det();
117
118                if (det < 0) {
119                        vt = vt_orig.transpose();
120                        vt.set(2, 0, (0 - vt.get(2, 0)));
121                        vt.set(2, 1, (0 - vt.get(2, 1)));
122                        vt.set(2, 2, (0 - vt.get(2, 2)));
123
124                        Matrix nv_transp = vt.transpose();
125                        rot_nottrans = nv_transp.times(u_transp);
126                        rot = rot_nottrans.transpose();
127
128                }
129
130                Matrix cb_tmp = centroidB.times(rot);
131                Matrix tran = centroidA.minus(cb_tmp);
132
133                return Matrices.getTransformation(rot, tran);
134
135        }
136
137        @Override
138        public double getRmsd(Point3d[] x, Point3d[] y) {
139                Point3d[] yclone = CalcPoint.clonePoint3dArray(y);
140                superposeAndTransform(x, yclone);
141                return CalcPoint.rmsd(x, yclone);
142        }
143
144}