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 */
021
022package org.biojava.nbio.structure.geometry;
023
024import javax.vecmath.Matrix3d;
025import javax.vecmath.Matrix4d;
026import javax.vecmath.Point3d;
027import javax.vecmath.Vector3d;
028
029import org.slf4j.Logger;
030import org.slf4j.LoggerFactory;
031
032/**
033 * Implementation of the Quaternion-Based Characteristic Polynomial algorithm
034 * for RMSD and Superposition calculations.
035 * <p>
036 * Usage:
037 * <p>
038 * The input consists of 2 Point3d arrays of equal length. The input coordinates
039 * are not changed.
040 * 
041 * <pre>
042 *    Point3d[] x = ...
043 *    Point3d[] y = ...
044 *    SuperPositionQCP qcp = new SuperPositionQCP();
045 *    qcp.set(x, y);
046 * </pre>
047 * <p>
048 * or with weighting factors [0 - 1]]
049 * 
050 * <pre>
051 *    double[] weights = ...
052 *    qcp.set(x, y, weights);
053 * </pre>
054 * <p>
055 * For maximum efficiency, create a SuperPositionQCP object once and reuse it.
056 * <p>
057 * A. Calculate rmsd only
058 * 
059 * <pre>
060 * double rmsd = qcp.getRmsd();
061 * </pre>
062 * <p>
063 * B. Calculate a 4x4 transformation (rotation and translation) matrix
064 * 
065 * <pre>
066 * Matrix4d rottrans = qcp.getTransformationMatrix();
067 * </pre>
068 * <p>
069 * C. Get transformated points (y superposed onto the reference x)
070 * 
071 * <pre>
072 * Point3d[] ySuperposed = qcp.getTransformedCoordinates();
073 * </pre>
074 * <p>
075 * Citations:
076 * <p>
077 * Liu P, Agrafiotis DK, & Theobald DL (2011) Reply to comment on: "Fast
078 * determination of the optimal rotation matrix for macromolecular
079 * superpositions." Journal of Computational Chemistry 32(1):185-186.
080 * [http://dx.doi.org/10.1002/jcc.21606]
081 * <p>
082 * Liu P, Agrafiotis DK, & Theobald DL (2010) "Fast determination of the optimal
083 * rotation matrix for macromolecular superpositions." Journal of Computational
084 * Chemistry 31(7):1561-1563. [http://dx.doi.org/10.1002/jcc.21439]
085 * <p>
086 * Douglas L Theobald (2005) "Rapid calculation of RMSDs using a
087 * quaternion-based characteristic polynomial." Acta Crystallogr A
088 * 61(4):478-480. [http://dx.doi.org/10.1107/S0108767305015266 ]
089 * <p>
090 * This is an adoption of the original C code QCProt 1.4 (2012, October 10) to
091 * Java. The original C source code is available from
092 * http://theobald.brandeis.edu/qcp/ and was developed by
093 * <p>
094 * Douglas L. Theobald Department of Biochemistry MS 009 Brandeis University 415
095 * South St Waltham, MA 02453 USA
096 * <p>
097 * dtheobald@brandeis.edu
098 * <p>
099 * Pu Liu Johnson & Johnson Pharmaceutical Research and Development, L.L.C. 665
100 * Stockton Drive Exton, PA 19341 USA
101 * <p>
102 * pliu24@its.jnj.com
103 * <p>
104 * 
105 * @author Douglas L. Theobald (original C code)
106 * @author Pu Liu (original C code)
107 * @author Peter Rose (adopted to Java)
108 * @author Aleix Lafita (adopted to Java)
109 */
110public final class SuperPositionQCP extends SuperPositionAbstract {
111
112        private static final Logger logger = LoggerFactory.getLogger(SuperPositionQCP.class);
113
114        private double evec_prec = 1E-6;
115        private double eval_prec = 1E-11;
116
117        private Point3d[] x;
118        private Point3d[] y;
119
120        private double[] weight;
121        private double wsum;
122
123        private Point3d[] xref;
124        private Point3d[] yref;
125        private Point3d xtrans;
126        private Point3d ytrans;
127
128        private double e0;
129        private Matrix3d rotmat = new Matrix3d();
130        private Matrix4d transformation = new Matrix4d();
131        private double rmsd = 0;
132        private double Sxy, Sxz, Syx, Syz, Szx, Szy;
133        private double SxxpSyy, Szz, mxEigenV, SyzmSzy, SxzmSzx, SxymSyx;
134        private double SxxmSyy, SxypSyx, SxzpSzx;
135        private double Syy, Sxx, SyzpSzy;
136        private boolean rmsdCalculated = false;
137        private boolean transformationCalculated = false;
138        private boolean centered = false;
139
140        /**
141         * Default constructor for the quaternion based superposition algorithm.
142         * 
143         * @param centered
144         *            true if the point arrays are centered at the origin (faster),
145         *            false otherwise
146         */
147        public SuperPositionQCP(boolean centered) {
148                super(centered);
149        }
150
151        /**
152         * Constructor with option to set the precision values.
153         * 
154         * @param centered
155         *            true if the point arrays are centered at the origin (faster),
156         *            false otherwise
157         * @param evec_prec
158         *            required eigenvector precision
159         * @param eval_prec
160         *            required eigenvalue precision
161         */
162        public SuperPositionQCP(boolean centered, double evec_prec, double eval_prec) {
163                super(centered);
164                this.evec_prec = evec_prec;
165                this.eval_prec = eval_prec;
166        }
167
168        /**
169         * Sets the two input coordinate arrays. These input arrays must be of equal
170         * length. Input coordinates are not modified.
171         * 
172         * @param x
173         *            3d points of reference coordinate set
174         * @param y
175         *            3d points of coordinate set for superposition
176         */
177        private void set(Point3d[] x, Point3d[] y) {
178                this.x = x;
179                this.y = y;
180                rmsdCalculated = false;
181                transformationCalculated = false;
182        }
183
184        /**
185         * Sets the two input coordinate arrays and weight array. All input arrays
186         * must be of equal length. Input coordinates are not modified.
187         * 
188         * @param x
189         *            3d points of reference coordinate set
190         * @param y
191         *            3d points of coordinate set for superposition
192         * @param weight
193         *            a weight in the inclusive range [0,1] for each point
194         */
195        private void set(Point3d[] x, Point3d[] y, double[] weight) {
196                this.x = x;
197                this.y = y;
198                this.weight = weight;
199                rmsdCalculated = false;
200                transformationCalculated = false;
201        }
202
203        /**
204         * Return the RMSD of the superposition of input coordinate set y onto x.
205         * Note, this is the fasted way to calculate an RMSD without actually
206         * superposing the two sets. The calculation is performed "lazy", meaning
207         * calculations are only performed if necessary.
208         * 
209         * @return root mean square deviation for superposition of y onto x
210         */
211        private double getRmsd() {
212                if (!rmsdCalculated) {
213                        calcRmsd(x, y);
214                        rmsdCalculated = true;
215                }
216                return rmsd;
217        }
218
219        /**
220         * Weighted superposition.
221         * 
222         * @param fixed
223         * @param moved
224         * @param weight
225         *            array of weigths for each equivalent point position
226         * @return
227         */
228        public Matrix4d weightedSuperpose(Point3d[] fixed, Point3d[] moved, double[] weight) {
229                set(moved, fixed, weight);
230                getRotationMatrix();
231                if (!centered) {
232                        calcTransformation();
233                } else {
234                        transformation.set(rotmat);
235                }
236                return transformation;
237        }
238
239        private Matrix3d getRotationMatrix() {
240                getRmsd();
241                if (!transformationCalculated) {
242                        calcRotationMatrix();
243                        transformationCalculated = true;
244                }
245                return rotmat;
246        }
247
248        /**
249         * Calculates the RMSD value for superposition of y onto x. This requires
250         * the coordinates to be precentered.
251         * 
252         * @param x
253         *            3d points of reference coordinate set
254         * @param y
255         *            3d points of coordinate set for superposition
256         */
257        private void calcRmsd(Point3d[] x, Point3d[] y) {
258                if (centered) {
259                        innerProduct(y, x);
260                } else {
261                        // translate to origin
262                        xref = CalcPoint.clonePoint3dArray(x);
263                        xtrans = CalcPoint.centroid(xref);
264                        logger.debug("x centroid: " + xtrans);
265                        xtrans.negate();
266                        CalcPoint.translate(new Vector3d(xtrans), xref);
267
268                        yref = CalcPoint.clonePoint3dArray(y);
269                        ytrans = CalcPoint.centroid(yref);
270                        logger.debug("y centroid: " + ytrans);
271                        ytrans.negate();
272                        CalcPoint.translate(new Vector3d(ytrans), yref);
273                        innerProduct(yref, xref);
274                }
275                calcRmsd(wsum);
276        }
277
278        /**
279         * Superposition coords2 onto coords1 -- in other words, coords2 is rotated,
280         * coords1 is held fixed
281         */
282        private void calcTransformation() {
283
284                // transformation.set(rotmat,new Vector3d(0,0,0), 1);
285                transformation.set(rotmat);
286                // long t2 = System.nanoTime();
287                // System.out.println("create transformation: " + (t2-t1));
288                // System.out.println("m3d -> m4d");
289                // System.out.println(transformation);
290
291                // combine with x -> origin translation
292                Matrix4d trans = new Matrix4d();
293                trans.setIdentity();
294                trans.setTranslation(new Vector3d(xtrans));
295                transformation.mul(transformation, trans);
296                // System.out.println("setting xtrans");
297                // System.out.println(transformation);
298
299                // combine with origin -> y translation
300                ytrans.negate();
301                Matrix4d transInverse = new Matrix4d();
302                transInverse.setIdentity();
303                transInverse.setTranslation(new Vector3d(ytrans));
304                transformation.mul(transInverse, transformation);
305                // System.out.println("setting ytrans");
306                // System.out.println(transformation);
307        }
308
309        /**
310         * Calculates the inner product between two coordinate sets x and y
311         * (optionally weighted, if weights set through
312         * {@link #set(Point3d[], Point3d[], double[])}). It also calculates an
313         * upper bound of the most positive root of the key matrix.
314         * http://theobald.brandeis.edu/qcp/qcprot.c
315         * 
316         * @param coords1
317         * @param coords2
318         * @return
319         */
320        private void innerProduct(Point3d[] coords1, Point3d[] coords2) {
321                double x1, x2, y1, y2, z1, z2;
322                double g1 = 0.0, g2 = 0.0;
323
324                Sxx = 0;
325                Sxy = 0;
326                Sxz = 0;
327                Syx = 0;
328                Syy = 0;
329                Syz = 0;
330                Szx = 0;
331                Szy = 0;
332                Szz = 0;
333
334                if (weight != null) {
335                        wsum = 0;
336                        for (int i = 0; i < coords1.length; i++) {
337
338                                wsum += weight[i];
339
340                                x1 = weight[i] * coords1[i].x;
341                                y1 = weight[i] * coords1[i].y;
342                                z1 = weight[i] * coords1[i].z;
343
344                                g1 += x1 * coords1[i].x + y1 * coords1[i].y + z1 * coords1[i].z;
345
346                                x2 = coords2[i].x;
347                                y2 = coords2[i].y;
348                                z2 = coords2[i].z;
349
350                                g2 += weight[i] * (x2 * x2 + y2 * y2 + z2 * z2);
351
352                                Sxx += (x1 * x2);
353                                Sxy += (x1 * y2);
354                                Sxz += (x1 * z2);
355
356                                Syx += (y1 * x2);
357                                Syy += (y1 * y2);
358                                Syz += (y1 * z2);
359
360                                Szx += (z1 * x2);
361                                Szy += (z1 * y2);
362                                Szz += (z1 * z2);
363                        }
364                } else {
365                        for (int i = 0; i < coords1.length; i++) {
366                                g1 += coords1[i].x * coords1[i].x + coords1[i].y * coords1[i].y + coords1[i].z * coords1[i].z;
367                                g2 += coords2[i].x * coords2[i].x + coords2[i].y * coords2[i].y + coords2[i].z * coords2[i].z;
368
369                                Sxx += coords1[i].x * coords2[i].x;
370                                Sxy += coords1[i].x * coords2[i].y;
371                                Sxz += coords1[i].x * coords2[i].z;
372
373                                Syx += coords1[i].y * coords2[i].x;
374                                Syy += coords1[i].y * coords2[i].y;
375                                Syz += coords1[i].y * coords2[i].z;
376
377                                Szx += coords1[i].z * coords2[i].x;
378                                Szy += coords1[i].z * coords2[i].y;
379                                Szz += coords1[i].z * coords2[i].z;
380                        }
381                        wsum = coords1.length;
382                }
383
384                e0 = (g1 + g2) * 0.5;
385        }
386
387        private int calcRmsd(double len) {
388                double Sxx2 = Sxx * Sxx;
389                double Syy2 = Syy * Syy;
390                double Szz2 = Szz * Szz;
391
392                double Sxy2 = Sxy * Sxy;
393                double Syz2 = Syz * Syz;
394                double Sxz2 = Sxz * Sxz;
395
396                double Syx2 = Syx * Syx;
397                double Szy2 = Szy * Szy;
398                double Szx2 = Szx * Szx;
399
400                double SyzSzymSyySzz2 = 2.0 * (Syz * Szy - Syy * Szz);
401                double Sxx2Syy2Szz2Syz2Szy2 = Syy2 + Szz2 - Sxx2 + Syz2 + Szy2;
402
403                double c2 = -2.0 * (Sxx2 + Syy2 + Szz2 + Sxy2 + Syx2 + Sxz2 + Szx2 + Syz2 + Szy2);
404                double c1 = 8.0 * (Sxx * Syz * Szy + Syy * Szx * Sxz + Szz * Sxy * Syx - Sxx * Syy * Szz - Syz * Szx * Sxy
405                                - Szy * Syx * Sxz);
406
407                SxzpSzx = Sxz + Szx;
408                SyzpSzy = Syz + Szy;
409                SxypSyx = Sxy + Syx;
410                SyzmSzy = Syz - Szy;
411                SxzmSzx = Sxz - Szx;
412                SxymSyx = Sxy - Syx;
413                SxxpSyy = Sxx + Syy;
414                SxxmSyy = Sxx - Syy;
415
416                double Sxy2Sxz2Syx2Szx2 = Sxy2 + Sxz2 - Syx2 - Szx2;
417
418                double c0 = Sxy2Sxz2Syx2Szx2 * Sxy2Sxz2Syx2Szx2
419                                + (Sxx2Syy2Szz2Syz2Szy2 + SyzSzymSyySzz2) * (Sxx2Syy2Szz2Syz2Szy2 - SyzSzymSyySzz2)
420                                + (-(SxzpSzx) * (SyzmSzy) + (SxymSyx) * (SxxmSyy - Szz))
421                                                * (-(SxzmSzx) * (SyzpSzy) + (SxymSyx) * (SxxmSyy + Szz))
422                                + (-(SxzpSzx) * (SyzpSzy) - (SxypSyx) * (SxxpSyy - Szz))
423                                                * (-(SxzmSzx) * (SyzmSzy) - (SxypSyx) * (SxxpSyy + Szz))
424                                + (+(SxypSyx) * (SyzpSzy) + (SxzpSzx) * (SxxmSyy + Szz))
425                                                * (-(SxymSyx) * (SyzmSzy) + (SxzpSzx) * (SxxpSyy + Szz))
426                                + (+(SxypSyx) * (SyzmSzy) + (SxzmSzx) * (SxxmSyy - Szz))
427                                                * (-(SxymSyx) * (SyzpSzy) + (SxzmSzx) * (SxxpSyy - Szz));
428
429                mxEigenV = e0;
430
431                int i;
432                for (i = 1; i < 51; ++i) {
433                        double oldg = mxEigenV;
434                        double x2 = mxEigenV * mxEigenV;
435                        double b = (x2 + c2) * mxEigenV;
436                        double a = b + c1;
437                        double delta = ((a * mxEigenV + c0) / (2.0 * x2 * mxEigenV + b + a));
438                        mxEigenV -= delta;
439
440                        if (Math.abs(mxEigenV - oldg) < Math.abs(eval_prec * mxEigenV))
441                                break;
442                }
443
444                if (i == 50) {
445                        logger.warn(String.format("More than %d iterations needed!", i));
446                } else {
447                        logger.info(String.format("%d iterations needed!", i));
448                }
449
450                /*
451                 * the fabs() is to guard against extremely small, but *negative*
452                 * numbers due to floating point error
453                 */
454                rmsd = Math.sqrt(Math.abs(2.0 * (e0 - mxEigenV) / len));
455
456                return 1;
457        }
458
459        private int calcRotationMatrix() {
460                double a11 = SxxpSyy + Szz - mxEigenV;
461                double a12 = SyzmSzy;
462                double a13 = -SxzmSzx;
463                double a14 = SxymSyx;
464                double a21 = SyzmSzy;
465                double a22 = SxxmSyy - Szz - mxEigenV;
466                double a23 = SxypSyx;
467                double a24 = SxzpSzx;
468                double a31 = a13;
469                double a32 = a23;
470                double a33 = Syy - Sxx - Szz - mxEigenV;
471                double a34 = SyzpSzy;
472                double a41 = a14;
473                double a42 = a24;
474                double a43 = a34;
475                double a44 = Szz - SxxpSyy - mxEigenV;
476                double a3344_4334 = a33 * a44 - a43 * a34;
477                double a3244_4234 = a32 * a44 - a42 * a34;
478                double a3243_4233 = a32 * a43 - a42 * a33;
479                double a3143_4133 = a31 * a43 - a41 * a33;
480                double a3144_4134 = a31 * a44 - a41 * a34;
481                double a3142_4132 = a31 * a42 - a41 * a32;
482                double q1 = a22 * a3344_4334 - a23 * a3244_4234 + a24 * a3243_4233;
483                double q2 = -a21 * a3344_4334 + a23 * a3144_4134 - a24 * a3143_4133;
484                double q3 = a21 * a3244_4234 - a22 * a3144_4134 + a24 * a3142_4132;
485                double q4 = -a21 * a3243_4233 + a22 * a3143_4133 - a23 * a3142_4132;
486
487                double qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4;
488
489                /*
490                 * The following code tries to calculate another column in the adjoint
491                 * matrix when the norm of the current column is too small. Usually this
492                 * commented block will never be activated. To be absolutely safe this
493                 * should be uncommented, but it is most likely unnecessary.
494                 */
495                if (qsqr < evec_prec) {
496                        q1 = a12 * a3344_4334 - a13 * a3244_4234 + a14 * a3243_4233;
497                        q2 = -a11 * a3344_4334 + a13 * a3144_4134 - a14 * a3143_4133;
498                        q3 = a11 * a3244_4234 - a12 * a3144_4134 + a14 * a3142_4132;
499                        q4 = -a11 * a3243_4233 + a12 * a3143_4133 - a13 * a3142_4132;
500                        qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4;
501
502                        if (qsqr < evec_prec) {
503                                double a1324_1423 = a13 * a24 - a14 * a23, a1224_1422 = a12 * a24 - a14 * a22;
504                                double a1223_1322 = a12 * a23 - a13 * a22, a1124_1421 = a11 * a24 - a14 * a21;
505                                double a1123_1321 = a11 * a23 - a13 * a21, a1122_1221 = a11 * a22 - a12 * a21;
506
507                                q1 = a42 * a1324_1423 - a43 * a1224_1422 + a44 * a1223_1322;
508                                q2 = -a41 * a1324_1423 + a43 * a1124_1421 - a44 * a1123_1321;
509                                q3 = a41 * a1224_1422 - a42 * a1124_1421 + a44 * a1122_1221;
510                                q4 = -a41 * a1223_1322 + a42 * a1123_1321 - a43 * a1122_1221;
511                                qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4;
512
513                                if (qsqr < evec_prec) {
514                                        q1 = a32 * a1324_1423 - a33 * a1224_1422 + a34 * a1223_1322;
515                                        q2 = -a31 * a1324_1423 + a33 * a1124_1421 - a34 * a1123_1321;
516                                        q3 = a31 * a1224_1422 - a32 * a1124_1421 + a34 * a1122_1221;
517                                        q4 = -a31 * a1223_1322 + a32 * a1123_1321 - a33 * a1122_1221;
518                                        qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4;
519
520                                        if (qsqr < evec_prec) {
521                                                /*
522                                                 * if qsqr is still too small, return the identity
523                                                 * matrix.
524                                                 */
525                                                rotmat.setIdentity();
526
527                                                return 0;
528                                        }
529                                }
530                        }
531                }
532
533                double normq = Math.sqrt(qsqr);
534                q1 /= normq;
535                q2 /= normq;
536                q3 /= normq;
537                q4 /= normq;
538
539                logger.debug("q: " + q1 + " " + q2 + " " + q3 + " " + q4);
540
541                double a2 = q1 * q1;
542                double x2 = q2 * q2;
543                double y2 = q3 * q3;
544                double z2 = q4 * q4;
545
546                double xy = q2 * q3;
547                double az = q1 * q4;
548                double zx = q4 * q2;
549                double ay = q1 * q3;
550                double yz = q3 * q4;
551                double ax = q1 * q2;
552
553                rotmat.m00 = a2 + x2 - y2 - z2;
554                rotmat.m01 = 2 * (xy + az);
555                rotmat.m02 = 2 * (zx - ay);
556
557                rotmat.m10 = 2 * (xy - az);
558                rotmat.m11 = a2 - x2 + y2 - z2;
559                rotmat.m12 = 2 * (yz + ax);
560
561                rotmat.m20 = 2 * (zx + ay);
562                rotmat.m21 = 2 * (yz - ax);
563                rotmat.m22 = a2 - x2 - y2 + z2;
564
565                return 1;
566        }
567
568        @Override
569        public double getRmsd(Point3d[] fixed, Point3d[] moved) {
570                set(moved, fixed);
571                return getRmsd();
572        }
573
574        @Override
575        public Matrix4d superpose(Point3d[] fixed, Point3d[] moved) {
576                set(moved, fixed);
577                getRotationMatrix();
578                if (!centered) {
579                        calcTransformation();
580                } else {
581                        transformation.set(rotmat);
582                }
583                return transformation;
584        }
585
586        /**
587         * @param fixed
588         * @param moved
589         * @param weight
590         *            array of weigths for each equivalent point position
591         * @return weighted RMSD.
592         */
593        public double getWeightedRmsd(Point3d[] fixed, Point3d[] moved, double[] weight) {
594                set(moved, fixed, weight);
595                return getRmsd();
596        }
597
598        /**
599         * The QCP method can be used as a two-step calculation: first compute the
600         * RMSD (fast) and then compute the superposition.
601         * 
602         * This method assumes that the RMSD of two arrays of points has been
603         * already calculated using {@link #getRmsd(Point3d[], Point3d[])} method
604         * and calculates the transformation of the same two point arrays.
605         * 
606         * @param fixed
607         * @param moved
608         * @return transformation matrix as a Matrix4d to superpose moved onto fixed
609         *         point arrays
610         */
611        public Matrix4d superposeAfterRmsd() {
612
613                if (!rmsdCalculated) {
614                        throw new IllegalStateException("The RMSD was not yet calculated. Use the superpose() method instead.");
615                }
616
617                getRotationMatrix();
618                if (!centered) {
619                        calcTransformation();
620                } else {
621                        transformation.set(rotmat);
622                }
623                return transformation;
624        }
625
626}