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 */
021package org.biojava.nbio.structure.align.util;
022
023import java.io.StringWriter;
024
025import javax.vecmath.AxisAngle4d;
026import javax.vecmath.Matrix3d;
027import javax.vecmath.Matrix4d;
028import javax.vecmath.Vector3d;
029
030import org.biojava.nbio.structure.Atom;
031import org.biojava.nbio.structure.AtomImpl;
032import org.biojava.nbio.structure.Calc;
033import org.biojava.nbio.structure.StructureException;
034import org.biojava.nbio.structure.align.model.AFPChain;
035import org.biojava.nbio.structure.contact.Pair;
036import org.biojava.nbio.structure.geometry.Matrices;
037import org.biojava.nbio.structure.jama.Matrix;
038
039/**
040 * Calculates the rotation axis for an alignment
041 *
042 * <p>A superposition of two structures is generally represented as a rotation
043 * matrix plus a translation vector. However, it can also be represented as an
044 * axis of rotation plus some translation.
045 *
046 * <p>This class calculates the rotation axis and stores it as four properties:
047 * <ul><li>A unit vector parallel to the rotation axis ({@link #getRotationAxis()})
048 * <li>The angle of rotation ({@link #getAngle()})
049 * <li>A point on the rotation axis ({@link #getRotationPos()})
050 * <li>Some translation parallel to the axis ({@link #getScrewTranslation()})
051 * </ul>
052 *
053 * <p>The axis of rotation is poorly defined and numerically unstable for small
054 * angles. Therefore it's direction is left as null for angles less than
055 * {@link #MIN_ANGLE}.
056 *
057 * @author Spencer Bliven
058 *
059 */
060public final class RotationAxis {
061
062        /**
063         * Minimum angle to calculate rotation axes. 5 degrees.
064         */
065        static final double MIN_ANGLE = Math.toRadians(5.);
066
067        private double theta;
068        private Atom rotationAxis; // axis of rotation
069        private Atom rotationPos; // a point on the axis of rotation
070        private Atom screwTranslation; //translation parallel to the axis of rotation
071        private Atom otherTranslation; // translation perpendicular to the axis of rotation
072
073        /**
074         * The rotation angle
075         * @return the angle, in radians
076         */
077        public double getAngle() {
078                return theta;
079        }
080
081        /**
082         * Get a unit vector along the rotation axis
083         * @return rotationAxis
084         */
085        public Atom getRotationAxis() {
086                return rotationAxis;
087        }
088
089        /**
090         * Returns the rotation axis and angle in a single javax.vecmath.AxisAngle4d object
091         * @return
092         */
093        public AxisAngle4d getAxisAngle4d() {
094                return new AxisAngle4d(rotationAxis.getX(),rotationAxis.getY(),rotationAxis.getZ(),theta);
095        }
096
097        /**
098         * Get a position on the rotation axis.
099         *
100         * Specifically, project the origin onto the rotation axis
101         * @return rotationPos
102         */
103        public Atom getRotationPos() {
104                return rotationPos;
105        }
106
107        /**
108         * Get the component of translation parallel to the axis of rotation
109         * @return
110         */
111        public Atom getScrewTranslation() {
112                return screwTranslation;
113        }
114
115        public Vector3d getVector3dScrewTranslation() {
116                return new Vector3d(screwTranslation.getX(),screwTranslation.getY(),screwTranslation.getZ());
117        }
118
119        public double getTranslation() {
120                return Calc.amount(screwTranslation);
121        }
122
123        /**
124         * Calculate the rotation axis for the first block of an AFPChain
125         * @param afpChain
126         * @throws StructureException
127         * @throws NullPointerException if afpChain does not contain a valid rotation matrix and shift vector
128         */
129        public RotationAxis(AFPChain afpChain) throws StructureException {
130                if(afpChain.getAlnLength() < 1) {
131                        throw new StructureException("No aligned residues");
132                }
133                init(afpChain.getBlockRotationMatrix()[0],afpChain.getBlockShiftVector()[0]);
134        }
135
136        /**
137         * Create a rotation axis from a vector, a point, and an angle.
138         *
139         * The result will be a pure rotation, with no screw component.
140         * @param axis A vector parallel to the axis of rotation
141         * @param pos A point on the axis of rotation
142         * @param theta The angle to rotate (radians)
143         */
144        public RotationAxis(Atom axis, Atom pos, double theta) {
145                this.rotationAxis = Calc.unitVector(axis);
146                this.rotationPos = (Atom) pos.clone();
147                this.theta = theta;
148                this.screwTranslation = new AtomImpl(); //zero
149                this.otherTranslation = null; //deprecated
150        }
151
152        /**
153         * Determine the location of the rotation axis based on a rotation matrix and a translation vector
154         * @param rotation
155         * @param translation
156         */
157        public RotationAxis(Matrix rotation, Atom translation) {
158                init(rotation, translation);
159        }
160
161        /**
162         * Create a rotation axis from a Matrix4d containing a rotational
163         * component and a translational component.
164         *
165         * @param transform
166         */
167        public RotationAxis(Matrix4d transform) {
168
169                Matrix rot = Matrices.getRotationJAMA(transform);
170                Atom transl = Calc.getTranslationVector(transform);
171                init(rot,transl);
172        }
173
174        /**
175         * Get the rotation matrix corresponding to this axis
176         * @return A 3x3 rotation matrix
177         */
178        public Matrix getRotationMatrix() {
179                return getRotationMatrix(theta);
180        }
181
182        /**
183         * Get the rotation matrix corresponding to a rotation about this axis
184         * @param theta The amount to rotate
185         * @return A 3x3 rotation matrix
186         */
187        public Matrix getRotationMatrix(double theta) {
188                if( rotationAxis == null) {
189                        // special case for pure translational axes
190                        return Matrix.identity(3, 3);
191                }
192                double x = rotationAxis.getX();
193                double y = rotationAxis.getY();
194                double z = rotationAxis.getZ();
195                double cos = Math.cos(theta);
196                double sin = Math.sin(theta);
197                double com = 1 - cos;
198                return new Matrix(new double[][] {
199                                {com*x*x + cos, com*x*y+sin*z, com*x*z+-sin*y},
200                                {com*x*y-sin*z, com*y*y+cos, com*y*z+sin*x},
201                                {com*x*z+sin*y, com*y*z-sin*x, com*z*z+cos},
202                });
203        }
204
205        /**
206         * Returns the rotation order o that gives the lowest value of {@code |2PI / o - theta},
207         * given that the value is strictly lower than {@code threshold}, for orders {@code o=1,...,maxOrder}.
208         */
209        public int guessOrderFromAngle(double threshold, int maxOrder) {
210                double bestDelta = threshold;
211                int bestOrder = 1;
212                for (int order = 2; order < maxOrder; order++) {
213                        double delta = Math.abs(2 * Math.PI / order - theta);
214                        if (delta < bestDelta) {
215                                bestOrder = order;
216                                bestDelta = delta;
217                        }
218                }
219                return bestOrder;
220        }
221
222
223        /**
224         * Returns a matrix that describes both rotation and translation.
225         */
226        public Matrix getFullMatrix() {
227                return null; // TODO, easy
228        }
229
230        /**
231         * Initialize variables
232         *
233         * @param rotation
234         * @param translation
235         */
236        private void init(Matrix rotation, Atom translation) {
237                if(rotation.getColumnDimension() != 3 || rotation.getRowDimension() != 3) {
238                        throw new IllegalArgumentException("Expected 3x3 rotation matrix");
239                }
240
241
242                // Calculate angle
243                double c = (rotation.trace()-1)/2.0; //=cos(theta)
244                // c is sometimes slightly out of the [-1,1] range due to numerical instabilities
245                if( -1-1e-8 < c && c < -1 ) c = -1;
246                if( 1+1e-8 > c && c > 1 ) c = 1;
247                if( -1 > c || c > 1 ) {
248                        throw new IllegalArgumentException("Input matrix is not a valid rotation matrix.");
249                }
250                this.theta = Math.acos(c);
251
252                if(theta < MIN_ANGLE) {
253                        calculateTranslationalAxis(rotation,translation);
254                } else {
255                        calculateRotationalAxis(rotation, translation, c);
256                }
257        }
258
259        /**
260         * Calculate the rotation axis for the normal case, where there is a
261         * significant rotation angle
262         * @param rotation
263         * @param translation
264         * @param c
265         */
266        private void calculateRotationalAxis(Matrix rotation, Atom translation,
267                        double c) {
268                // Calculate magnitude of rotationAxis components, but not signs
269                double sum=0;
270                double[] rotAx = new double[3];
271                for(int i=0;i<3;i++) {
272                        rotAx[i] = Math.sqrt(rotation.get(i, i)-c);
273                        sum+=rotAx[i]*rotAx[i];
274                }
275                for(int i=0;i<3;i++) {
276                        rotAx[i] /= Math.sqrt(sum);
277                }
278
279                // Now determine signs
280                double d0 = rotation.get(2,1)-rotation.get(1,2); //=2u[0]*sin(theta)
281                double d1 = rotation.get(0,2)-rotation.get(2,0); //=2u[1]*sin(theta)
282                double d2 = rotation.get(1,0)-rotation.get(0,1); //=2u[2]*sin(theta)
283
284                double s12 = rotation.get(2,1)+rotation.get(1,2); //=2*u[1]*u[2]*(1-cos(theta))
285                double s02 = rotation.get(0,2)+rotation.get(2,0); //=2*u[0]*u[2]*(1-cos(theta))
286                double s01 = rotation.get(1,0)+rotation.get(0,1); //=2*u[0]*u[1]*(1-cos(theta))
287
288                //Only use biggest d for the sign directly, for numerical stability around 180deg
289                if( Math.abs(d0) < Math.abs(d1) ) { // not d0
290                        if( Math.abs(d1) < Math.abs(d2) ) { //d2
291                                if(d2>=0){ //u[2] positive
292                                        if( s02 < 0 ) rotAx[0] = -rotAx[0];
293                                        if( s12 < 0 ) rotAx[1] = -rotAx[1];
294                                } else { //u[2] negative
295                                        rotAx[2] = -rotAx[2];
296                                        if( s02 >= 0 ) rotAx[0] = -rotAx[0];
297                                        if( s12 >= 0 ) rotAx[1] = -rotAx[1];
298                                }
299                        } else { //d1
300                                if(d1>=0) {//u[1] positive
301                                        if( s01 < 0) rotAx[0] = -rotAx[0];
302                                        if( s12 < 0) rotAx[2] = -rotAx[2];
303                                } else { //u[1] negative
304                                        rotAx[1] = -rotAx[1];
305                                        if( s01 >= 0) rotAx[0] = -rotAx[0];
306                                        if( s12 >= 0) rotAx[2] = -rotAx[2];
307                                }
308                        }
309                } else { // not d1
310                        if( Math.abs(d0) < Math.abs(d2) ) { //d2
311                                if(d2>=0){ //u[2] positive
312                                        if( s02 < 0 ) rotAx[0] = -rotAx[0];
313                                        if( s12 < 0 ) rotAx[1] = -rotAx[1];
314                                } else { //u[2] negative
315                                        rotAx[2] = -rotAx[2];
316                                        if( s02 >= 0 ) rotAx[0] = -rotAx[0];
317                                        if( s12 >= 0 ) rotAx[1] = -rotAx[1];
318                                }
319                        } else { //d0
320                                if(d0>=0) { //u[0] positive
321                                        if( s01 < 0 ) rotAx[1] = -rotAx[1];
322                                        if( s02 < 0 ) rotAx[2] = -rotAx[2];
323                                } else { //u[0] negative
324                                        rotAx[0] = -rotAx[0];
325                                        if( s01 >= 0 ) rotAx[1] = -rotAx[1];
326                                        if( s02 >= 0 ) rotAx[2] = -rotAx[2];
327                                }
328                        }
329                }
330
331                rotationAxis = new AtomImpl();
332                rotationAxis.setCoords(rotAx);
333
334                // Calculate screw = (rotationAxis dot translation)*u
335                double dotProduct = Calc.scalarProduct(rotationAxis, translation);
336
337                screwTranslation = Calc.scale(rotationAxis, dotProduct);
338                otherTranslation = Calc.subtract(translation, screwTranslation);
339
340                Atom hypot = Calc.vectorProduct(otherTranslation,rotationAxis);
341                Calc.scaleEquals(hypot,.5/Math.tan(theta/2.0));
342
343                // Calculate rotation axis position
344                rotationPos = Calc.scaleAdd(.5,otherTranslation, hypot);
345        }
346
347        /**
348         * Handle cases with small angles of rotation
349         * @param rotation
350         * @param translation
351         */
352        private void calculateTranslationalAxis(Matrix rotation, Atom translation) {
353                // set axis parallel to translation
354                rotationAxis = Calc.scale(translation, 1./Calc.amount(translation));
355
356                // position is undefined
357                rotationPos = null;
358
359                screwTranslation = translation;
360                otherTranslation = new AtomImpl();
361                otherTranslation.setCoords(new double[] {0,0,0});
362        }
363
364        /**
365         * Returns a Jmol script which will display the axis of rotation. This
366         * consists of a cyan arrow along the axis, plus an arc showing the angle
367         * of rotation.
368         * <p>
369         * As the rotation angle gets smaller, the axis of rotation becomes poorly
370         * defined and would need to get farther and farther away from the protein.
371         * This is not particularly useful, so we arbitrarily draw it parallel to
372         * the translation and omit the arc.
373         * @param atoms Some atoms from the protein, used for determining the bounds
374         *        of the axis.
375         *
376         * @return The Jmol script, suitable for calls to
377         * {@link org.biojava.nbio.structure.align.gui.jmol.StructureAlignmentJmol#evalString() jmol.evalString()}
378         */
379        public String getJmolScript(Atom[] atoms){
380                return getJmolScript(atoms, 0);
381        }
382
383        /**
384         * Find a segment of the axis that covers the specified set of atoms.
385         * <p>
386         * Projects the input atoms onto the rotation axis and returns the bounding
387         * points.
388         * <p>
389         * In the case of a pure translational axis, the axis location is undefined
390         * so the center of mass will be used instead.
391         * @param atoms
392         * @return two points defining the axis segment
393         */
394        public Pair<Atom> getAxisEnds(Atom[] atoms) {
395                // Project each Atom onto the rotation axis to determine limits
396                double min, max;
397                min = max = Calc.scalarProduct(rotationAxis,atoms[0]);
398                for(int i=1;i<atoms.length;i++) {
399                        double prod = Calc.scalarProduct(rotationAxis,atoms[i]);
400                        if(prod<min) min = prod;
401                        if(prod>max) max = prod;
402                }
403                double uLen = Calc.scalarProduct(rotationAxis,rotationAxis);// Should be 1, but double check
404                min/=uLen;
405                max/=uLen;
406
407                // Project the origin onto the axis. If the axis is undefined, use the center of mass
408                Atom axialPt;
409                if(rotationPos == null) {
410                        Atom center = Calc.centerOfMass(atoms);
411
412                        // Project center onto the axis
413                        Atom centerOnAxis = Calc.scale(rotationAxis, Calc.scalarProduct(center, rotationAxis));
414
415                        // Remainder is projection of origin onto axis
416                        axialPt = Calc.subtract(center, centerOnAxis);
417
418                } else {
419                        axialPt = rotationPos;
420                }
421
422                // Find end points of the rotation axis to display
423                Atom axisMin = (Atom) axialPt.clone();
424                Calc.scaleAdd(min, rotationAxis, axisMin);
425                Atom axisMax = (Atom) axialPt.clone();
426                Calc.scaleAdd(max, rotationAxis, axisMax);
427
428                return new Pair<>(axisMin, axisMax);
429        }
430        /**
431         * Returns a Jmol script which will display the axis of rotation. This
432         * consists of a cyan arrow along the axis, plus an arc showing the angle
433         * of rotation.
434         * <p>
435         * As the rotation angle gets smaller, the axis of rotation becomes poorly
436         * defined and would need to get farther and farther away from the protein.
437         * This is not particularly useful, so we arbitrarily draw it parallel to
438         * the translation and omit the arc.
439         * @param atoms Some atoms from the protein, used for determining the bounds
440         *        of the axis.
441         * @param axisID in case of representing more than one axis in the same jmol
442         *                panel, indicate the ID number.
443         *
444         * @return The Jmol script, suitable for calls to
445         * {@link org.biojava.nbio.structure.align.gui.jmol.StructureAlignmentJmol#evalString() jmol.evalString()}
446         */
447        public String getJmolScript(Atom[] atoms, int axisID){
448                final double width=.5;// width of JMol object
449                final String axisColor = "yellow"; //axis color
450                final String screwColor = "orange"; //screw translation color
451
452                Pair<Atom> endPoints = getAxisEnds(atoms);
453                Atom axisMin = endPoints.getFirst();
454                Atom axisMax = endPoints.getSecond();
455
456                StringWriter result = new StringWriter();
457
458                // set arrow heads to a reasonable length
459                result.append("set defaultDrawArrowScale 2.0;");
460
461                // draw axis of rotation
462                result.append(
463                                String.format("draw ID rot"+axisID+" CYLINDER {%f,%f,%f} {%f,%f,%f} WIDTH %f COLOR %s ;",
464                                                axisMin.getX(),axisMin.getY(),axisMin.getZ(),
465                                                axisMax.getX(),axisMax.getY(),axisMax.getZ(), width, axisColor ));
466
467                // draw screw component
468                boolean positiveScrew = Math.signum(rotationAxis.getX()) == Math.signum(screwTranslation.getX());
469                if( positiveScrew ) {
470                        // screw is in the same direction as the axis
471                        result.append( String.format(
472                                        "draw ID screw"+axisID+" VECTOR {%f,%f,%f} {%f,%f,%f} WIDTH %f COLOR %s ;",
473                                        axisMax.getX(),axisMax.getY(),axisMax.getZ(),
474                                        screwTranslation.getX(),screwTranslation.getY(),screwTranslation.getZ(),
475                                        width, screwColor ));
476                } else {
477                        // screw is in the opposite direction as the axis
478                        result.append( String.format(
479                                        "draw ID screw"+axisID+" VECTOR {%f,%f,%f} {%f,%f,%f} WIDTH %f COLOR %s ;",
480                                        axisMin.getX(),axisMin.getY(),axisMin.getZ(),
481                                        screwTranslation.getX(),screwTranslation.getY(),screwTranslation.getZ(),
482                                        width, screwColor ));
483                }
484
485                // draw angle of rotation
486                if(rotationPos != null) {
487                        result.append(System.getProperty("line.separator"));
488                        result.append(String.format("draw ID rotArc"+axisID+" ARC {%f,%f,%f} {%f,%f,%f} {0,0,0} {0,%f,%d} SCALE 500 DIAMETER %f COLOR %s;",
489                                        axisMin.getX(),axisMin.getY(),axisMin.getZ(),
490                                        axisMax.getX(),axisMax.getY(),axisMax.getZ(),
491                                        Math.toDegrees(theta),
492                                        positiveScrew ? 0 : 1 , // draw at the opposite end from the screw arrow
493                                                        width, axisColor ));
494                }
495
496                return result.toString();
497        }
498
499        /**
500         * Projects a given point onto the axis of rotation
501         * @param point
502         * @return An atom which lies on the axis, or null if the RotationAxis is purely translational
503         */
504        public Atom getProjectedPoint(Atom point) {
505                if(rotationPos == null) {
506                        // translation only
507                        return null;
508                }
509
510                Atom localPoint = Calc.subtract(point, rotationPos);
511                double dot = Calc.scalarProduct(localPoint, rotationAxis);
512
513                Atom localProjected = Calc.scale(rotationAxis, dot);
514                Atom projected = Calc.add(localProjected, rotationPos);
515                return projected;
516        }
517
518        /**
519         * Get the distance from a point to the axis of rotation
520         * @param point
521         * @return The distance to the axis, or NaN if the RotationAxis is purely translational
522         */
523        public double getProjectedDistance(Atom point) {
524                Atom projected = getProjectedPoint(point);
525                if( projected == null) {
526                        // translation only
527                        return Double.NaN;
528                }
529
530
531                return Calc.getDistance(point, projected);
532
533        }
534
535        public void rotate(Atom[] atoms, double theta) {
536                Matrix rot = getRotationMatrix(theta);
537                if(rotationPos == null) {
538                        // Undefined rotation axis; do nothing
539                        return;
540                }
541                Atom negPos;
542
543                        negPos = Calc.invert(rotationPos);
544
545                for(Atom a: atoms) {
546                        Calc.shift(a, negPos);
547                }
548                Calc.rotate(atoms, rot);
549                for(Atom a: atoms) {
550                        Calc.shift(a, rotationPos);
551                }
552        }
553
554        /**
555         * Calculate the rotation angle for a structure
556         * @param afpChain
557         * @return The rotation angle, in radians
558         * @throws StructureException If the alignment doesn't contain any blocks
559         * @throws NullPointerException If the alignment doesn't have a rotation matrix set
560         */
561        public static double getAngle(AFPChain afpChain) throws StructureException {
562                if(afpChain.getBlockNum() < 1) {
563                        throw new StructureException("No aligned residues");
564                }
565                Matrix rotation = afpChain.getBlockRotationMatrix()[0];
566
567                if(rotation == null) {
568                        throw new NullPointerException("AFPChain does not contain a rotation matrix");
569                }
570                return getAngle(rotation);
571        }
572        /**
573         * Calculate the rotation angle for a given matrix
574         * @param rotation Rotation matrix
575         * @return The angle, in radians
576         */
577        public static double getAngle(Matrix rotation) {
578                double c = (rotation.trace()-1)/2.0; //=cos(theta)
579                // c is sometimes slightly out of the [-1,1] range due to numerical instabilities
580                if( -1-1e-8 < c && c < -1 ) c = -1;
581                if( 1+1e-8 > c && c > 1 ) c = 1;
582                if( -1 > c || c > 1 ) {
583                        throw new IllegalArgumentException("Input matrix is not a valid rotation matrix.");
584                }
585                return Math.acos(c);
586        }
587
588        /**
589         *
590         * @return If the rotation axis is well defined, rather than purely translational
591         */
592        public boolean isDefined() {
593                return rotationPos != null;
594        }
595
596        /**
597         * Quickly compute the rotation angle from a rotation matrix.
598         * @param transform 4D transformation matrix. Translation components are ignored.
599         * @return Angle, from 0 to PI
600         */
601        public static double getAngle(Matrix4d transform) {
602                // Calculate angle
603                double c = (transform.m00 + transform.m11 + transform.m22 - 1)/2.0; //=cos(theta)
604                // c is sometimes slightly out of the [-1,1] range due to numerical instabilities
605                if( -1-1e-8 < c && c < -1 ) c = -1;
606                if( 1+1e-8 > c && c > 1 ) c = 1;
607                if( -1 > c || c > 1 ) {
608                        throw new IllegalArgumentException("Input matrix is not a valid rotation matrix.");
609                }
610                return Math.acos(c);
611        }
612
613        /**
614         * Quickly compute the rotation angle from a rotation matrix.
615         * @param transform 3D rotation matrix
616         * @return Angle, from 0 to PI
617         */
618        public static double getAngle(Matrix3d transform) {
619                // Calculate angle
620                double c = (transform.m00 + transform.m11 + transform.m22 - 1)/2.0; //=cos(theta)
621                // c is sometimes slightly out of the [-1,1] range due to numerical instabilities
622                if( -1-1e-8 < c && c < -1 ) c = -1;
623                if( 1+1e-8 > c && c > 1 ) c = 1;
624                if( -1 > c || c > 1 ) {
625                        throw new IllegalArgumentException("Input matrix is not a valid rotation matrix.");
626                }
627                return Math.acos(c);
628        }
629}