Source code

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.symmetry.axis;
022
023import org.biojava.nbio.structure.geometry.CalcPoint;
024import org.biojava.nbio.structure.geometry.MomentsOfInertia;
025import org.biojava.nbio.structure.symmetry.core.Helix;
026import org.biojava.nbio.structure.symmetry.core.HelixLayers;
027import org.biojava.nbio.structure.symmetry.core.QuatSymmetryResults;
028import org.biojava.nbio.structure.symmetry.core.QuatSymmetrySubunits;
029import org.slf4j.Logger;
030import org.slf4j.LoggerFactory;
031
032import javax.vecmath.*;
033
034import java.util.*;
035
036public class HelixAxisAligner extends AxisAligner {
037
038        private static final Logger logger = LoggerFactory
039                        .getLogger(HelixAxisAligner.class);
040
041
042        private static final Vector3d Y_AXIS = new Vector3d(0,1,0);
043        private static final Vector3d Z_AXIS = new Vector3d(0,0,1);
044
045        private QuatSymmetrySubunits subunits = null;
046        private HelixLayers helixLayers = null;
047
048        private Matrix4d transformationMatrix = new Matrix4d();
049        private Matrix4d reverseTransformationMatrix = new Matrix4d();
050        private Vector3d referenceVector = new Vector3d();
051        private Vector3d principalRotationVector = new Vector3d();
052        private Vector3d[] principalAxesOfInertia = null;
053        private List<List<Integer>> alignedOrbits = null;
054
055        private Vector3d minBoundary = new Vector3d();
056        private Vector3d maxBoundary = new Vector3d();
057        private double xzRadiusMax = Double.MIN_VALUE;
058
059        boolean modified = true;
060
061        public HelixAxisAligner(QuatSymmetryResults results) {
062                this.subunits = new QuatSymmetrySubunits(results.getSubunitClusters());
063                this.helixLayers = results.getHelixLayers();
064                if (subunits == null) {
065                        throw new IllegalArgumentException("HelixAxisTransformation: Subunits are null");
066                } else if (helixLayers == null) {
067                        throw new IllegalArgumentException("HelixAxisTransformation: HelixLayers is null");
068                } else if (subunits.getSubunitCount() == 0) {
069                        throw new IllegalArgumentException("HelixAxisTransformation: Subunits is empty");
070                } else if (helixLayers.size() == 0) {
071                        throw new IllegalArgumentException("HelixAxisTransformation: HelixLayers is empty");
072                }
073        }
074
075
076        /* (non-Javadoc)
077         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getTransformation()
078         */
079        @Override
080        public String getSymmetry() {
081                run();
082                return "H";
083        }
084
085        /* (non-Javadoc)
086         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getTransformation()
087         */
088        @Override
089        public Matrix4d getTransformation() {
090                run();
091                return transformationMatrix;
092        }
093
094        /* (non-Javadoc)
095         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getRotationMatrix()
096         */
097        @Override
098        public Matrix3d getRotationMatrix() {
099                run();
100                Matrix3d m = new Matrix3d();
101                transformationMatrix.getRotationScale(m);
102                return m;
103        }
104
105        /* (non-Javadoc)
106         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getReverseTransformation()
107         */
108        @Override
109        public Matrix4d getReverseTransformation() {
110                run();
111                return reverseTransformationMatrix;
112        }
113
114        /* (non-Javadoc)
115         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getPrincipalRotationAxis()
116         */
117        @Override
118        public Vector3d getPrincipalRotationAxis() {
119                run();
120                return principalRotationVector;
121        }
122
123        /* (non-Javadoc)
124         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getRotationReferenceAxis()
125         */
126        @Override
127        public Vector3d getRotationReferenceAxis() {
128                run();
129                return referenceVector;
130        }
131
132        /* (non-Javadoc)
133         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getPrincipalAxesOfInertia()
134         */
135        @Override
136        public Vector3d[] getPrincipalAxesOfInertia() {
137                run();
138                return principalAxesOfInertia;
139        }
140
141        /* (non-Javadoc)
142         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getDimension()
143         */
144        @Override
145        public Vector3d getDimension() {
146                run();
147                Vector3d dimension = new Vector3d();
148                dimension.sub(maxBoundary, minBoundary);
149                dimension.scale(0.5);
150                return dimension;
151        }
152
153        /* (non-Javadoc)
154         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getXYRadius()
155         */
156        @Override
157        public double getRadius() {
158                run();
159                return xzRadiusMax;
160        }
161
162        /* (non-Javadoc)
163         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getGeometicCenterTransformation()
164         */
165        @Override
166        public Matrix4d getGeometicCenterTransformation() {
167                run();
168
169                Matrix4d geometricCentered = new Matrix4d(reverseTransformationMatrix);
170                geometricCentered.setTranslation(new Vector3d(getGeometricCenter()));
171
172                return geometricCentered;
173        }
174
175        /* (non-Javadoc)
176         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getGeometricCenter()
177         */
178        @Override
179        public Point3d getGeometricCenter() {
180                run();
181
182                Point3d geometricCenter = new Point3d();
183                Vector3d translation = new Vector3d();
184//              reverseTransformationMatrix.get(translation);
185
186                // TODO does this apply to the helic case?
187                // calculate adjustment around z-axis and transform adjustment to
188                //  original coordinate frame with the reverse transformation
189
190//              Vector3d corr = new Vector3d(0,minBoundary.y+getDimension().y, 0);
191//              reverseTransformationMatrix.transform(corr);
192//              geometricCenter.set(corr);
193
194                reverseTransformationMatrix.transform(translation);
195                geometricCenter.add(translation);
196                return geometricCenter;
197        }
198
199        /* (non-Javadoc)
200         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getCentroid()
201         */
202        @Override
203        public Point3d getCentroid() {
204                return new Point3d(subunits.getCentroid());
205        }
206
207        /* (non-Javadoc)
208         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getSubunits()
209         */
210        @Override
211        public QuatSymmetrySubunits getSubunits() {
212                return subunits;
213        }
214
215        public HelixLayers getHelixLayers() {
216                run();
217                return helixLayers;
218        }
219
220        /* (non-Javadoc)
221         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getOrbits()
222         */
223        @Override
224        public List<List<Integer>> getOrbits() {
225                run();
226                return alignedOrbits;
227        }
228
229        /**
230         * @return
231         */
232
233        private void run () {
234                if (modified) {
235                        calcPrincipalRotationVector();
236                        calcPrincipalAxes();
237                        calcReferenceVector();
238                        calcTransformation();
239                        // orient helix along Y axis by rotating 90 degrees around X-axis
240                        transformationMatrix = reorientHelix(0);
241
242                        calcReverseTransformation();
243                        calcBoundaries();
244                        calcAlignedOrbits();
245                        calcCenterOfRotation();
246
247                        // orient helix along Y axis by rotating 90 degrees around X-axis
248//                      transformationMatrix = reorientHelix(0);
249//                      calcReverseTransformation();
250
251                        modified = false;
252                }
253        }
254
255        public Point3d calcCenterOfRotation() {
256                List<Integer> line = getLongestLayerLine();
257
258                // can't determine center of rotation if there are only 2 points
259                // TODO does this ever happen??
260                if (line.size() < 3) {
261                        return subunits.getCentroid();
262                }
263
264                Point3d centerOfRotation = new Point3d();
265                List<Point3d> centers = subunits.getOriginalCenters();
266
267                // calculate helix mid points for each set of 3 adjacent subunits
268                for (int i = 0; i < line.size()-2; i++) {
269                        Point3d p1 = new Point3d(centers.get(line.get(i)));
270                        Point3d p2 = new Point3d(centers.get(line.get(i+1)));
271                        Point3d p3 = new Point3d(centers.get(line.get(i+2)));
272                        transformationMatrix.transform(p1);
273                        transformationMatrix.transform(p2);
274                        transformationMatrix.transform(p3);
275                        centerOfRotation.add(getMidPoint(p1, p2, p3));
276                }
277
278                // average over all midpoints to find best center of rotation
279                centerOfRotation.scale(1/(line.size()-2));
280                // since helix is aligned along the y-axis, with an origin at y = 0, place the center of rotation there
281                centerOfRotation.y = 0;
282                // transform center of rotation to the original coordinate frame
283                reverseTransformationMatrix.transform(centerOfRotation);
284//              System.out.println("center of rotation: " + centerOfRotation);
285                return centerOfRotation;
286        }
287
288        private List<Integer> getLongestLayerLine() {
289                int len = 0;
290                int index = 0;
291
292                Helix helix = helixLayers.getByLargestContacts();
293                List<List<Integer>> layerLines = helix.getLayerLines();
294                for (int i = 0; i < layerLines.size(); i++) {
295                        if (layerLines.get(i).size() > len) {
296                                len = layerLines.get(i).size();
297                                index = i;
298                        }
299                }
300                return layerLines.get(index);
301        }
302
303        /**
304         * Return a midpoint of a helix, calculated from three positions
305         * of three adjacent subunit centers.
306         * @param p1 center of first subunit
307         * @param p2 center of second subunit
308         * @param p3 center of third subunit
309         * @return midpoint of helix
310         */
311        private Point3d getMidPoint(Point3d p1, Point3d p2, Point3d p3) {
312                Vector3d v1 = new Vector3d();
313                v1.sub(p1, p2);
314                Vector3d v2 = new Vector3d();
315                v2.sub(p3, p2);
316                Vector3d v3 = new Vector3d();
317                v3.add(v1);
318                v3.add(v2);
319                v3.normalize();
320
321                // calculat the total distance between to subunits
322                double dTotal = v1.length();
323                // calculate the rise along the y-axis. The helix axis is aligned with y-axis,
324                // therfore, the rise between subunits is the y-distance
325                double rise = p2.y - p1.y;
326                // use phythagorean theoremm to calculate chord length between two subunit centers
327                double chord = Math.sqrt(dTotal*dTotal - rise*rise);
328//              System.out.println("Chord d: " + dTotal + " rise: " + rise + "chord: " + chord);
329                double angle = helixLayers.getByLargestContacts().getAxisAngle().angle;
330
331                // using the axis angle and the chord length, we can calculate the radius of the helix
332                // http://en.wikipedia.org/wiki/Chord_%28geometry%29
333                double radius = chord/Math.sin(angle/2)/2; // can this go to zero?
334//              System.out.println("Radius: " + radius);
335
336                // project the radius onto the vector that points toward the helix axis
337                v3.scale(radius);
338                v3.add(p2);
339//              System.out.println("Angle: " + Math.toDegrees(helixLayers.getByLowestAngle().getAxisAngle().angle));
340                Point3d cor = new Point3d(v3);
341                return cor;
342        }
343
344        private Matrix4d reorientHelix(int index) {
345                Matrix4d matrix = new Matrix4d();
346                matrix.setIdentity();
347                matrix.setRotation(new AxisAngle4d(1,0,0,Math.PI/2*(index+1)));
348                matrix.mul(transformationMatrix);
349        return matrix;
350        }
351
352        /**
353         * Returns a list of orbits (an orbit is a cyclic permutation of subunit indices that are related
354         * by a rotation around the principal rotation axis) ordered from the +z direction to the -z direction (z-depth).
355         * Within an orbit, subunit indices are ordered with a clockwise rotation around the z-axis.
356         * @return list of orbits ordered by z-depth
357         */
358        private void calcAlignedOrbits() {
359                Map<Double, List<Integer>> depthMap = new TreeMap<Double, List<Integer>>();
360                double[] depth = getSubunitZDepth();
361                alignedOrbits = calcOrbits();
362
363                // calculate the mean depth of orbit along z-axis
364                for (List<Integer> orbit: alignedOrbits) {
365                        // calculate the mean depth along z-axis for each orbit
366                        double meanDepth = 0;
367                        for (int subunit: orbit) {
368                                meanDepth += depth[subunit];
369                        }
370                        meanDepth /= orbit.size();
371
372                        if (depthMap.get(meanDepth) != null) {
373                                meanDepth += 0.01;
374                        }
375//                      System.out.println("calcAlignedOrbits: " + meanDepth + " orbit: " + orbit);
376                        depthMap.put(meanDepth, orbit);
377                }
378
379                // now fill orbits back into list ordered by depth
380                alignedOrbits.clear();
381                for (List<Integer> orbit: depthMap.values()) {
382                        // order subunit in a clockwise rotation around the z-axis
383                        /// starting at the 12 O-clock position (+y position)
384                        // TODO how should this be aligned??
385        //              alignWithReferenceAxis(orbit);
386                        alignedOrbits.add(orbit);
387                }
388        }
389
390
391        private void calcTransformation() {
392                calcTransformationBySymmetryAxes();
393                // make sure this value is zero. On Linux this value is 0.
394                transformationMatrix.setElement(3, 3, 1.0);
395        }
396
397        private void calcReverseTransformation() {
398                reverseTransformationMatrix.invert(transformationMatrix);
399        }
400
401        private void calcTransformationBySymmetryAxes() {
402                Vector3d[] axisVectors = new Vector3d[2];
403                axisVectors[0] = new Vector3d(principalRotationVector);
404                axisVectors[1] = new Vector3d(referenceVector);
405
406                //  y,z axis centered at the centroid of the subunits
407                Vector3d[] referenceVectors = new Vector3d[2];
408                referenceVectors[0] = new Vector3d(Z_AXIS);
409                referenceVectors[1] = new Vector3d(Y_AXIS);
410
411                transformationMatrix = alignAxes(axisVectors, referenceVectors);
412
413                // combine with translation
414                Matrix4d combined = new Matrix4d();
415                combined.setIdentity();
416                Vector3d trans = new Vector3d(subunits.getCentroid());
417                trans.negate();
418                combined.setTranslation(trans);
419                transformationMatrix.mul(combined);
420
421                // for helical geometry, set a canonical view for the Z direction
422                calcZDirection();
423        }
424
425        /**
426         * Returns a transformation matrix that rotates refPoints to match
427         * coordPoints
428         * @param refPoints the points to be aligned
429         * @param referenceVectors
430         * @return
431         */
432        private Matrix4d alignAxes(Vector3d[] axisVectors, Vector3d[] referenceVectors) {
433                Matrix4d m1 = new Matrix4d();
434                AxisAngle4d a = new AxisAngle4d();
435                Vector3d axis = new Vector3d();
436
437                // calculate rotation matrix to rotate refPoints[0] into coordPoints[0]
438                Vector3d v1 = new Vector3d(axisVectors[0]);
439                Vector3d v2 = new Vector3d(referenceVectors[0]);
440                double dot = v1.dot(v2);
441                if (Math.abs(dot) < 0.999) {
442                        axis.cross(v1,v2);
443                        axis.normalize();
444                        a.set(axis, v1.angle(v2));
445                        m1.set(a);
446                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
447                        m1.setElement(3,  3, 1.0);
448                } else if (dot > 0) {
449                        // parallel axis, nothing to do -> identity matrix
450                        m1.setIdentity();
451                } else if (dot < 0) {
452                        // anti-parallel axis, flip around x-axis
453                        m1.set(flipX());
454                }
455
456                // apply transformation matrix to all refPoints
457                m1.transform(axisVectors[0]);
458                m1.transform(axisVectors[1]);
459
460                // calculate rotation matrix to rotate refPoints[1] into coordPoints[1]
461                v1 = new Vector3d(axisVectors[1]);
462                v2 = new Vector3d(referenceVectors[1]);
463                Matrix4d m2 = new Matrix4d();
464                dot = v1.dot(v2);
465                if (Math.abs(dot) < 0.999) {
466                        axis.cross(v1,v2);
467                        axis.normalize();
468                        a.set(axis, v1.angle(v2));
469                        m2.set(a);
470                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
471                        m2.setElement(3,  3, 1.0);
472                } else if (dot > 0) {
473                        // parallel axis, nothing to do -> identity matrix
474                        m2.setIdentity();
475                } else if (dot < 0) {
476                        // anti-parallel axis, flip around z-axis
477                        m2.set(flipZ());
478                }
479
480                // apply transformation matrix to all refPoints
481                m2.transform(axisVectors[0]);
482                m2.transform(axisVectors[1]);
483
484                // combine the two rotation matrices
485                m2.mul(m1);
486
487                // the RMSD should be close to zero
488                Point3d[] axes = new Point3d[2];
489                axes[0] = new Point3d(axisVectors[0]);
490                axes[1] = new Point3d(axisVectors[1]);
491                Point3d[] ref = new Point3d[2];
492                ref[0] = new Point3d(referenceVectors[0]);
493                ref[1] = new Point3d(referenceVectors[1]);
494                if (CalcPoint.rmsd(axes, ref) > 0.1) {
495                        logger.warn("AxisTransformation: axes alignment is off. RMSD: "
496                                        + CalcPoint.rmsd(axes, ref));
497                }
498
499                return m2;
500        }
501
502        private void calcPrincipalAxes() {
503                MomentsOfInertia moi = new MomentsOfInertia();
504
505                for (Point3d[] list: subunits.getTraces()) {
506                        for (Point3d p: list) {
507                                moi.addPoint(p, 1.0);
508                        }
509                }
510                principalAxesOfInertia = moi.getPrincipalAxes();
511        }
512
513        /**
514         * Calculates the min and max boundaries of the structure after it has been
515         * transformed into its canonical orientation.
516         */
517        private void calcBoundaries() {
518                minBoundary.x = Double.MAX_VALUE;
519                maxBoundary.x = Double.MIN_VALUE;
520                minBoundary.y = Double.MAX_VALUE;
521                maxBoundary.x = Double.MIN_VALUE;
522                minBoundary.z = Double.MAX_VALUE;
523                maxBoundary.z = Double.MIN_VALUE;
524                xzRadiusMax = Double.MIN_VALUE;
525
526                Point3d probe = new Point3d();
527
528                for (Point3d[] list: subunits.getTraces()) {
529                        for (Point3d p: list) {
530                                probe.set(p);
531                                transformationMatrix.transform(probe);
532
533                                minBoundary.x = Math.min(minBoundary.x, probe.x);
534                                maxBoundary.x = Math.max(maxBoundary.x, probe.x);
535                                minBoundary.y = Math.min(minBoundary.y, probe.y);
536                                maxBoundary.y = Math.max(maxBoundary.y, probe.y);
537                                minBoundary.z = Math.min(minBoundary.z, probe.z);
538                                maxBoundary.z = Math.max(maxBoundary.z, probe.z);
539                                xzRadiusMax = Math.max(xzRadiusMax, Math.sqrt(probe.x*probe.x + probe.z * probe.z));
540                        }
541                }
542//              System.out.println("MinBoundary: " + minBoundary);
543//              System.out.println("MaxBoundary: " + maxBoundary);
544//              System.out.println("zxRadius: " + xzRadiusMax);
545        }
546
547        /*
548         * Modifies the rotation part of the transformation axis for
549         * a Cn symmetric complex, so that the narrower end faces the
550         * viewer, and the wider end faces away from the viewer. Example: 3LSV
551         */
552        private void calcZDirection() {
553                calcBoundaries();
554
555                // if the longer part of the structure faces towards the back (-z direction),
556                // rotate around y-axis so the longer part faces the viewer (+z direction)
557                if (Math.abs(minBoundary.z) > Math.abs(maxBoundary.z)) {
558                        Matrix4d rot = flipY();
559                        rot.mul(transformationMatrix);
560                        transformationMatrix.set(rot);
561                }
562        }
563
564        private double[] getSubunitZDepth() {
565                int n = subunits.getSubunitCount();
566                double[] depth = new double[n];
567                Point3d probe = new Point3d();
568
569                // transform subunit centers into z-aligned position and calculate
570                // z-coordinates (depth) along the z-axis.
571                for (int i = 0; i < n; i++) {
572                        Point3d p= subunits.getCenters().get(i);
573                        probe.set(p);
574                        transformationMatrix.transform(probe);
575                        depth[i] = probe.z;
576                }
577                return depth;
578        }
579
580        /**
581         * Returns a list of list of subunit ids that form an "orbit", i.e. they
582         * are transformed into each other during a rotation around the principal symmetry axis (z-axis)
583         * @return
584         */
585        private List<List<Integer>> calcOrbits() {
586                int n = subunits.getSubunitCount();
587
588                List<List<Integer>> orbits = new ArrayList<List<Integer>>();
589                for (int i = 0; i < n; i++) {
590                        orbits.add(Collections.singletonList(i));
591                }
592
593                return orbits;
594        }
595
596        /**
597         * Returns a vector along the principal rotation axis for the
598         * alignment of structures along the z-axis
599         * @return principal rotation vector
600         */
601        private void calcPrincipalRotationVector() {
602//              AxisAngle4d axisAngle = helixLayers.getByLowestAngle().getAxisAngle();
603                AxisAngle4d axisAngle = helixLayers.getByLargestContacts().getAxisAngle();
604                principalRotationVector = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
605        }
606
607        /**
608         * Returns a vector perpendicular to the principal rotation vector
609         * for the alignment of structures in the xy-plane
610         * @return reference vector
611         */
612        private void calcReferenceVector() {
613                referenceVector = getReferenceAxisCylic();
614
615                if (referenceVector == null) {
616                        logger.warn("no reference vector found. Using y-axis.");
617                        referenceVector = new Vector3d(Y_AXIS);
618                }
619                // make sure reference vector is perpendicular principal roation vector
620                referenceVector = orthogonalize(principalRotationVector, referenceVector);
621        }
622
623        private Vector3d orthogonalize(Vector3d vector1, Vector3d vector2) {
624                double dot = vector1.dot(vector2);
625                Vector3d ref = new Vector3d(vector2);
626//              System.out.println("p.r: " + dot);
627//              System.out.println("Orig refVector: " + referenceVector);
628                if (dot < 0) {
629                        vector2.negate();
630                }
631                if (Math.abs(dot) < 0.00001) {
632                        logger.info("HelixAxisAligner: reference axis parallel");
633                }
634                vector2.cross(vector1, vector2);
635//              System.out.println("Intermed. refVector: " + vector2);
636                vector2.normalize();
637//              referenceVector.cross(referenceVector, principalRotationVector);
638                vector2.cross(vector1, vector2);
639                vector2.normalize();
640                if (ref.dot(vector2) < 0) {
641                        vector2.negate();
642                }
643//              System.out.println("Mod. refVector: " + vector2);
644                return vector2;
645        }
646        /**
647         * Returns the default reference vector for the alignment of Cn structures
648         * @return
649         */
650        private Vector3d getReferenceAxisCylic() {
651                // get principal axis vector that is perpendicular to the principal
652                // rotation vector
653                Vector3d vmin = null;
654                double dotMin = 1.0;
655                for (Vector3d v: principalAxesOfInertia) {
656                        if (Math.abs(principalRotationVector.dot(v)) < dotMin) {
657                                dotMin = Math.abs(principalRotationVector.dot(v));
658                                vmin = new Vector3d(v);
659                        }
660                }
661                if (principalRotationVector.dot(vmin) < 0) {
662                        vmin.negate();
663                }
664
665                return vmin;
666        }
667
668        private static Matrix4d flipX() {
669                Matrix4d rot = new Matrix4d();
670                rot.m00 = 1;
671                rot.m11 = -1;
672                rot.m22 = -1;
673                rot.m33 = 1;
674                return rot;
675        }
676
677        private static Matrix4d flipY() {
678                Matrix4d rot = new Matrix4d();
679                rot.m00 = -1;
680                rot.m11 = 1;
681                rot.m22 = -1;
682                rot.m33 = 1;
683                return rot;
684        }
685
686        private static Matrix4d flipZ() {
687                Matrix4d rot = new Matrix4d();
688                rot.m00 = -1;
689                rot.m11 = -1;
690                rot.m22 = 1;
691                rot.m33 = 1;
692                return rot;
693        }
694
695}