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.QuatSymmetryResults;
026import org.biojava.nbio.structure.symmetry.core.Rotation;
027import org.biojava.nbio.structure.symmetry.core.RotationGroup;
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 RotationAxisAligner extends AxisAligner{
037
038        private static final Logger logger = LoggerFactory
039                        .getLogger(RotationAxisAligner.class);
040
041        private static final Vector3d X_AXIS = new Vector3d(1,0,0);
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 RotationGroup rotationGroup = 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         List<List<Integer>> alignedOrbits = null;
054
055        private Vector3d minBoundary = new Vector3d();
056        private Vector3d maxBoundary = new Vector3d();
057        private double xyRadiusMax = Double.MIN_VALUE;
058
059        boolean modified = true;
060
061        public RotationAxisAligner(QuatSymmetryResults results) {
062                this.subunits = new QuatSymmetrySubunits(results.getSubunitClusters());
063                this.rotationGroup = results.getRotationGroup();
064
065                if (subunits == null) {
066                        throw new IllegalArgumentException("AxisTransformation: Subunits are null");
067                } else if (rotationGroup == null) {
068                        throw new IllegalArgumentException("AxisTransformation: RotationGroup is null");
069                } else if (subunits.getSubunitCount() == 0) {
070                        throw new IllegalArgumentException("AxisTransformation: Subunits is empty");
071                } else if (rotationGroup.getOrder() == 0) {
072                        throw new IllegalArgumentException("AxisTransformation: RotationGroup is empty");
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 rotationGroup.getPointGroup();
083        }
084
085        @Override
086        public Matrix4d getTransformation() {
087                run();
088                return transformationMatrix;
089        }
090
091        @Override
092        public Matrix3d getRotationMatrix() {
093                run();
094                Matrix3d m = new Matrix3d();
095                transformationMatrix.getRotationScale(m);
096                return m;
097        }
098
099        @Override
100        public Matrix4d getReverseTransformation() {
101                run();
102                return reverseTransformationMatrix;
103        }
104
105        @Override
106        public Vector3d getPrincipalRotationAxis() {
107                run();
108                return principalRotationVector;
109        }
110
111        @Override
112        public Vector3d getRotationReferenceAxis() {
113                run();
114                return referenceVector;
115        }
116
117        @Override
118        public Vector3d[] getPrincipalAxesOfInertia() {
119                run();
120                return principalAxesOfInertia;
121        }
122
123        @Override
124        public Vector3d getDimension() {
125                run();
126                Vector3d dimension = new Vector3d();
127                dimension.sub(maxBoundary, minBoundary);
128                dimension.scale(0.5);
129                return dimension;
130        }
131
132        /**
133         * Returns the radius for drawing the minor rotation axis in the xy-plane
134         * @return double radius in xy-plane
135         */
136        @Override
137        public double getRadius() {
138                run();
139                return xyRadiusMax;
140        }
141
142        /**
143         * Returns a transformation matrix transform polyhedra for Cn structures.
144         * The center in this matrix is the geometric center, rather then the centroid.
145         * In Cn structures those are usually not the same.
146         * @return
147         */
148        @Override
149        public Matrix4d getGeometicCenterTransformation() {
150                run();
151
152                Matrix4d geometricCentered = new Matrix4d(reverseTransformationMatrix);
153                geometricCentered.setTranslation(new Vector3d(getGeometricCenter()));
154
155                return geometricCentered;
156        }
157
158        /**
159         * Returns the geometric center of polyhedron. In the case of the Cn
160         * point group, the centroid and geometric center are usually not
161         * identical.
162         * @return
163         */
164        @Override
165        public Point3d getGeometricCenter() {
166                run();
167
168                Point3d geometricCenter = new Point3d();
169                Vector3d translation = new Vector3d();
170                reverseTransformationMatrix.get(translation);
171
172                // calculate adjustment around z-axis and transform adjustment to
173                //  original coordinate frame with the reverse transformation
174                if (rotationGroup.getPointGroup().startsWith("C")) {
175                        Vector3d corr = new Vector3d(0,0, minBoundary.z+getDimension().z);
176                        reverseTransformationMatrix.transform(corr);
177                        geometricCenter.set(corr);
178                }
179
180                geometricCenter.add(translation);
181                return geometricCenter;
182        }
183
184        @Override
185        public Point3d getCentroid() {
186                return new Point3d(subunits.getCentroid());
187        }
188
189        @Override
190        public QuatSymmetrySubunits getSubunits() {
191                return subunits;
192        }
193
194        public RotationGroup getRotationGroup() {
195                return rotationGroup;
196        }
197
198        @Override
199        public List<List<Integer>> getOrbits() {
200                return alignedOrbits;
201        }
202
203        /**
204         * @return
205         */
206
207        private void run () {
208                if (modified) {
209                        calcPrincipalRotationVector();
210                        calcPrincipalAxes();
211                        // initial alignment with draft reference axis
212                        calcReferenceVector();
213                        calcTransformation();
214
215                        // refine ref. axis for cyclic and dihedral systems
216                        if ((rotationGroup.getPointGroup().startsWith("C") &&
217                                        !rotationGroup.getPointGroup().startsWith("C2")) ||
218                                (rotationGroup.getPointGroup().startsWith("D") &&
219                                                        !rotationGroup.getPointGroup().startsWith("D2"))
220                                        ) {
221                                refineReferenceVector();
222                                calcTransformation();
223                        }
224                        calcReverseTransformation();
225                        calcBoundaries();
226                        if (! rotationGroup.getPointGroup().equals("Helical")) {
227                                calcAlignedOrbits();
228                        }
229                        modified = false;
230                }
231        }
232        /**
233         * Returns a list of orbits (an orbit is a cyclic permutation of subunit indices that are related
234         * by a rotation around the principal rotation axis) ordered from the +z direction to the -z direction (z-depth).
235         * Within an orbit, subunit indices are ordered with a clockwise rotation around the z-axis.
236         * @return list of orbits ordered by z-depth
237         */
238        private void calcAlignedOrbits() {
239                Map<Double, List<Integer>> depthMap = new TreeMap<Double, List<Integer>>();
240                double[] depth = getSubunitZDepth();
241                alignedOrbits = calcOrbits();
242
243                // calculate the mean depth of orbit along z-axis
244                for (List<Integer> orbit: alignedOrbits) {
245                        // calculate the mean depth along z-axis for each orbit
246                        double meanDepth = 0;
247                        for (int subunit: orbit) {
248                                meanDepth += depth[subunit];
249                        }
250                        meanDepth /= orbit.size();
251
252                        if (depthMap.get(meanDepth) != null) {
253                                // System.out.println("Conflict in depthMap");
254                                meanDepth += 0.01;
255                        }
256                        depthMap.put(meanDepth, orbit);
257                }
258
259                // now fill orbits back into list ordered by depth
260                alignedOrbits.clear();
261                for (List<Integer> orbit: depthMap.values()) {
262                        // order subunit in a clockwise rotation around the z-axis
263                        /// starting at the 12 O-clock position (+y position)
264                        alignWithReferenceAxis(orbit);
265                        alignedOrbits.add(orbit);
266                }
267        }
268
269        /**
270         * Returns an ordered list of subunit ids (orbit) in such a way that the subunit
271         * indices start at the 12 o-clock (+y axis) and proceed in a clockwise direction
272         * to the 11 o-clock position to close the "orbit".
273         *
274         * @param orbit list of subunit indices that are transformed into each other by a rotation
275         * @return list of subunit indices ordered in a clockwise direction
276         */
277
278        private List<Integer> alignWithReferenceAxis(List<Integer> orbit) {
279                int n = subunits.getSubunitCount();
280                int fold = rotationGroup.getRotation(0).getFold();
281                if (fold < 2) {
282                        return orbit;
283                }
284                Vector3d probe = new Vector3d();
285
286                double dotMin1 = Double.MIN_VALUE;
287                double dotMin2 = Double.MIN_VALUE;
288                int index1 = 0;
289                int index2 = 0;
290                Vector3d Y1 = new Vector3d(0,1,0);
291                Vector3d Y2 = new Vector3d(0,1,0);
292                Matrix3d m = new Matrix3d();
293                double angle = -2*Math.PI/fold;
294                m.rotZ(0.1*angle); // add small offset, since two subunits may be equidistant to the y-axis
295                m.transform(Y1);
296                m.rotZ(1.1*angle);
297                m.transform(Y2);
298                // transform subunit centers into z-aligned position and calculate
299                // width in xy direction.
300                for (int i: orbit) {
301                        Point3d p = subunits.getCenters().get(i);
302                        probe.set(p);
303                        transformationMatrix.transform(probe);
304                        // find subunit that lines up with y-axis
305                        double dot1 = Y1.dot(probe);
306                        if (dot1 > dotMin1) {
307                                dotMin1 = dot1;
308                                index1 = i;
309                        }
310                        // find next subunit (rotated by one fold around z-axis - clockwise)
311                        double dot2 = Y2.dot(probe);
312                        if (dot2 > dotMin2) {
313                                dotMin2 = dot2;
314                                index2 = i;
315                        }
316                }
317//              System.out.println("Index1/2: " + index1 + " - " + index2);
318//              System.out.println("Orbit0: " + orbit);
319                // order subunit indices in a clockwise orientation around the z-axis
320                // bring subunit into position 0
321                for (int i = 0; i < n; i++) {
322                        if (orbit.get(0) == index1) {
323                                break;
324                        }
325                        Collections.rotate(orbit,1);
326                }
327//              System.out.println("Orbit1: " + orbit);
328                // bring second subunit  onto position 1
329                if (orbit.get(1) == index2) {
330                        return orbit;
331                }
332                Collections.reverse(orbit.subList(1,  orbit.size()));
333                if (orbit.get(1) != index2) {
334                        logger.warn("alignWithReferenceAxis failed");
335                }
336//              System.out.println("Orbit2: " + orbit);
337                return orbit;
338        }
339
340
341
342        private void calcTransformation() {
343                if (rotationGroup.getPointGroup().equals("C1")) {
344                        calcTransformationByInertiaAxes();
345                } else {
346                        calcTransformationBySymmetryAxes();
347                }
348                // make sure this value is zero. On Linux this value is 0.
349                transformationMatrix.setElement(3, 3, 1.0);
350        }
351
352        private void calcReverseTransformation() {
353                reverseTransformationMatrix.invert(transformationMatrix);
354        }
355
356        private void calcTransformationBySymmetryAxes() {
357                Vector3d[] axisVectors = new Vector3d[2];
358                axisVectors[0] = new Vector3d(principalRotationVector);
359                axisVectors[1] = new Vector3d(referenceVector);
360
361                //  y,z axis centered at the centroid of the subunits
362                Vector3d[] referenceVectors = new Vector3d[2];
363                referenceVectors[0] = new Vector3d(Z_AXIS);
364                referenceVectors[1] = new Vector3d(Y_AXIS);
365
366                transformationMatrix = alignAxes(axisVectors, referenceVectors);
367
368                // combine with translation
369                Matrix4d combined = new Matrix4d();
370                combined.setIdentity();
371                Vector3d trans = new Vector3d(subunits.getCentroid());
372                trans.negate();
373                combined.setTranslation(trans);
374                transformationMatrix.mul(combined);
375
376                // for cyclic geometry, set a canonical view for the Z direction
377                if (rotationGroup.getPointGroup().startsWith("C")) {
378                        calcZDirection();
379                }
380        }
381
382        private void calcTransformationByInertiaAxes() {
383                Vector3d[] axisVectors = new Vector3d[2];
384                axisVectors[0] = new Vector3d(principalAxesOfInertia[0]);
385                axisVectors[1] = new Vector3d(principalAxesOfInertia[1]);
386
387                Vector3d[] referenceVectors = new Vector3d[2];
388                referenceVectors[0] = new Vector3d(Y_AXIS);
389                referenceVectors[1] = new Vector3d(X_AXIS);
390
391                // align inertia axes with y-x plane
392                transformationMatrix = alignAxes(axisVectors, referenceVectors);
393
394                // combine with translation
395                Matrix4d translation = new Matrix4d();
396                translation.setIdentity();
397                Vector3d trans = new Vector3d(subunits.getCentroid());
398                trans.negate();
399                translation.setTranslation(trans);
400                transformationMatrix.mul(translation);
401        }
402
403        /**
404         * Returns a transformation matrix that rotates refPoints to match
405         * coordPoints
406         * @param refPoints the points to be aligned
407         * @param referenceVectors
408         * @return
409         */
410        private Matrix4d alignAxes(Vector3d[] axisVectors, Vector3d[] referenceVectors) {
411                Matrix4d m1 = new Matrix4d();
412                AxisAngle4d a = new AxisAngle4d();
413                Vector3d axis = new Vector3d();
414
415                // calculate rotation matrix to rotate refPoints[0] into coordPoints[0]
416                Vector3d v1 = new Vector3d(axisVectors[0]);
417                Vector3d v2 = new Vector3d(referenceVectors[0]);
418                double dot = v1.dot(v2);
419                if (Math.abs(dot) < 0.999) {
420                        axis.cross(v1,v2);
421                        axis.normalize();
422                        a.set(axis, v1.angle(v2));
423                        m1.set(a);
424                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
425                        m1.setElement(3,  3, 1.0);
426                } else if (dot > 0) {
427                        // parallel axis, nothing to do -> identity matrix
428                        m1.setIdentity();
429                } else if (dot < 0) {
430                        // anti-parallel axis, flip around x-axis
431                        m1.set(flipX());
432                }
433
434                // apply transformation matrix to all refPoints
435                m1.transform(axisVectors[0]);
436                m1.transform(axisVectors[1]);
437
438                // calculate rotation matrix to rotate refPoints[1] into coordPoints[1]
439                v1 = new Vector3d(axisVectors[1]);
440                v2 = new Vector3d(referenceVectors[1]);
441                Matrix4d m2 = new Matrix4d();
442                dot = v1.dot(v2);
443                if (Math.abs(dot) < 0.999) {
444                        axis.cross(v1,v2);
445                        axis.normalize();
446                        a.set(axis, v1.angle(v2));
447                        m2.set(a);
448                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
449                        m2.setElement(3,  3, 1.0);
450                } else if (dot > 0) {
451                        // parallel axis, nothing to do -> identity matrix
452                        m2.setIdentity();
453                } else if (dot < 0) {
454                        // anti-parallel axis, flip around z-axis
455                        m2.set(flipZ());
456                }
457
458                // apply transformation matrix to all refPoints
459                m2.transform(axisVectors[0]);
460                m2.transform(axisVectors[1]);
461
462                // combine the two rotation matrices
463                m2.mul(m1);
464
465                // the RMSD should be close to zero
466                Point3d[] axes = new Point3d[2];
467                axes[0] = new Point3d(axisVectors[0]);
468                axes[1] = new Point3d(axisVectors[1]);
469                Point3d[] ref = new Point3d[2];
470                ref[0] = new Point3d(referenceVectors[0]);
471                ref[1] = new Point3d(referenceVectors[1]);
472                if (CalcPoint.rmsd(axes, ref) > 0.1) {
473                        logger.warn("AxisTransformation: axes alignment is off. RMSD: "
474                                        + CalcPoint.rmsd(axes, ref));
475                }
476
477                return m2;
478        }
479
480        private void calcPrincipalAxes() {
481                MomentsOfInertia moi = new MomentsOfInertia();
482
483                for (Point3d[] list: subunits.getTraces()) {
484                        for (Point3d p: list) {
485                                moi.addPoint(p, 1.0);
486                        }
487                }
488                principalAxesOfInertia = moi.getPrincipalAxes();
489        }
490
491        /**
492         * Calculates the min and max boundaries of the structure after it has been
493         * transformed into its canonical orientation.
494         */
495        private void calcBoundaries() {
496                minBoundary.x = Double.MAX_VALUE;
497                maxBoundary.x = Double.MIN_VALUE;
498                minBoundary.y = Double.MAX_VALUE;
499                maxBoundary.x = Double.MIN_VALUE;
500                minBoundary.z = Double.MAX_VALUE;
501                maxBoundary.z = Double.MIN_VALUE;
502
503                Point3d probe = new Point3d();
504
505                for (Point3d[] list: subunits.getTraces()) {
506                        for (Point3d p: list) {
507                                probe.set(p);
508                                transformationMatrix.transform(probe);
509
510                                minBoundary.x = Math.min(minBoundary.x, probe.x);
511                                maxBoundary.x = Math.max(maxBoundary.x, probe.x);
512                                minBoundary.y = Math.min(minBoundary.y, probe.y);
513                                maxBoundary.y = Math.max(maxBoundary.y, probe.y);
514                                minBoundary.z = Math.min(minBoundary.z, probe.z);
515                                maxBoundary.z = Math.max(maxBoundary.z, probe.z);
516                                xyRadiusMax = Math.max(xyRadiusMax, Math.sqrt(probe.x*probe.x + probe.y * probe.y));
517                        }
518                }
519        }
520
521        /*
522         * Modifies the rotation part of the transformation axis for
523         * a Cn symmetric complex, so that the narrower end faces the
524         * viewer, and the wider end faces away from the viewer. Example: 3LSV
525         */
526        private void calcZDirection() {
527                calcBoundaries();
528
529                // if the longer part of the structure faces towards the back (-z direction),
530                // rotate around y-axis so the longer part faces the viewer (+z direction)
531                if (Math.abs(minBoundary.z) > Math.abs(maxBoundary.z)) {
532                        Matrix4d rot = flipY();
533                        rot.mul(transformationMatrix);
534                        transformationMatrix.set(rot);
535                }
536        }
537
538        /**
539         *
540         */
541        private List<List<Integer>> getOrbitsByXYWidth() {
542                Map<Double, List<Integer>> widthMap = new TreeMap<Double, List<Integer>>();
543                double[] width = getSubunitXYWidth();
544                List<List<Integer>> orbits = calcOrbits();
545
546                // calculate the mean width of orbits in XY-plane
547                for (List<Integer> orbit: orbits) {
548                        double meanWidth = 0;
549                        for (int subunit: orbit) {
550                                meanWidth += width[subunit];
551                        }
552                        meanWidth /= orbit.size();
553
554                        if (widthMap.get(meanWidth) != null) {
555                                meanWidth += 0.01;
556                        }
557                        widthMap.put(meanWidth, orbit);
558                }
559
560                // now fill orbits back into list ordered by width
561                orbits.clear();
562                for (List<Integer> orbit: widthMap.values()) {
563                        orbits.add(orbit);
564                }
565                return orbits;
566        }
567
568        private double[] getSubunitXYWidth() {
569                int n = subunits.getSubunitCount();
570                double[] width = new double[n];
571                Point3d probe = new Point3d();
572
573                // transform subunit centers into z-aligned position and calculate
574                // width in xy direction.
575                for (int i = 0; i < n; i++) {
576                        width[i] = Double.MIN_VALUE;
577                        for (Point3d p: subunits.getTraces().get(i)) {
578                                probe.set(p);
579                                transformationMatrix.transform(probe);
580                                width[i] = Math.max(width[i], Math.sqrt(probe.x*probe.x + probe.y*probe.y));
581                        }
582                }
583                return width;
584        }
585
586        private double[] getSubunitZDepth() {
587                int n = subunits.getSubunitCount();
588                double[] depth = new double[n];
589                Point3d probe = new Point3d();
590
591                // transform subunit centers into z-aligned position and calculate
592                // z-coordinates (depth) along the z-axis.
593                for (int i = 0; i < n; i++) {
594                        Point3d p= subunits.getCenters().get(i);
595                        probe.set(p);
596                        transformationMatrix.transform(probe);
597                        depth[i] = probe.z;
598                }
599                return depth;
600        }
601
602        /**
603         * Returns a list of list of subunit ids that form an "orbit", i.e. they
604         * are transformed into each other during a rotation around the principal symmetry axis (z-axis)
605         * @return
606         */
607        private List<List<Integer>> calcOrbits() {
608                int n = subunits.getSubunitCount();
609                int fold = rotationGroup.getRotation(0).getFold();
610
611                List<List<Integer>> orbits = new ArrayList<List<Integer>>();
612                boolean[] used = new boolean[n];
613                Arrays.fill(used, false);
614
615                for (int i = 0; i < n; i++) {
616                        if (! used[i]) {
617                                // determine the equivalent subunits
618                                List<Integer> orbit = new ArrayList<Integer>(fold);
619                                for (int j = 0; j < fold; j++) {
620                                        List<Integer> permutation = rotationGroup.getRotation(j).getPermutation();
621                                        orbit.add(permutation.get(i));
622                                        used[permutation.get(i)] = true;
623                                }
624                                orbits.add(deconvolute(orbit));
625                        }
626                }
627                return orbits;
628        }
629
630        private List<Integer> deconvolute(List<Integer> orbit) {
631                if (rotationGroup.getOrder() < 2) {
632                        return orbit;
633                }
634                List<Integer> p0 = rotationGroup.getRotation(0).getPermutation();
635                List<Integer> p1 = rotationGroup.getRotation(1).getPermutation();
636//              System.out.println("deconvolute");
637//              System.out.println("Permutation0: " + p0);
638//              System.out.println("Permutation1: " + p1);
639
640                List<Integer> inRotationOrder = new ArrayList<Integer>(orbit.size());
641                inRotationOrder.add(orbit.get(0));
642                for (int i = 1; i < orbit.size(); i++) {
643                        int current = inRotationOrder.get(i-1);
644                        int index = p0.indexOf(current);
645                        int next = p1.get(index);
646                        if (!orbit.contains(next)) {
647                                logger.warn("deconvolute: inconsistency in permuation. Returning original order");
648                                return orbit;
649                        }
650                        inRotationOrder.add(next);
651                }
652//              System.out.println("In order: " + inRotationOrder);
653                return inRotationOrder;
654        }
655
656        /**
657         * Returns a vector along the principal rotation axis for the
658         * alignment of structures along the z-axis
659         * @return principal rotation vector
660         */
661        private void calcPrincipalRotationVector() {
662                Rotation rotation = rotationGroup.getRotation(0); // the rotation around the principal axis is the first rotation
663                AxisAngle4d axisAngle = rotation.getAxisAngle();
664                principalRotationVector = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
665        }
666
667        /**
668         * Returns a vector perpendicular to the principal rotation vector
669         * for the alignment of structures in the xy-plane
670         * @return reference vector
671         */
672        private void calcReferenceVector() {
673                referenceVector = null;
674                if (rotationGroup.getPointGroup().startsWith("C")) {
675                        referenceVector = getReferenceAxisCylic();
676                } else if (rotationGroup.getPointGroup().startsWith("D")) {
677                        referenceVector = getReferenceAxisDihedral();
678                } else if (rotationGroup.getPointGroup().equals("T")) {
679                        referenceVector = getReferenceAxisTetrahedral();
680                } else if (rotationGroup.getPointGroup().equals("O")) {
681                        referenceVector = getReferenceAxisOctahedral();
682                } else if (rotationGroup.getPointGroup().equals("I")) {
683                        referenceVector = getReferenceAxisIcosahedral();
684                } else if (rotationGroup.getPointGroup().equals("Helical")) {
685                        // TODO what should the reference vector be??
686                        referenceVector = getReferenceAxisCylic();
687                }
688
689                if (referenceVector == null) {
690                        logger.warn("no reference vector found. Using y-axis.");
691                        referenceVector = new Vector3d(Y_AXIS);
692                }
693                // make sure reference vector is perpendicular principal roation vector
694                referenceVector = orthogonalize(principalRotationVector, referenceVector);
695        }
696
697        /**
698         * Returns a normalized vector that represents a minor rotation axis, except
699         * for Cn, this represents an axis orthogonal to the principal axis.
700         * @return minor rotation axis
701         */
702        private void refineReferenceVector() {
703                referenceVector = new Vector3d(Y_AXIS);
704                if (rotationGroup.getPointGroup().startsWith("C")) {
705                        referenceVector = getReferenceAxisCylicWithSubunitAlignment();
706                } else if (rotationGroup.getPointGroup().startsWith("D")) {
707                        referenceVector = getReferenceAxisDihedralWithSubunitAlignment();
708                }
709
710                referenceVector = orthogonalize(principalRotationVector, referenceVector);
711        }
712
713        private Vector3d orthogonalize(Vector3d vector1, Vector3d vector2) {
714                double dot = vector1.dot(vector2);
715                Vector3d ref = new Vector3d(vector2);
716//              System.out.println("p.r: " + dot);
717//              System.out.println("Orig refVector: " + referenceVector);
718                if (dot < 0) {
719                        vector2.negate();
720                }
721                vector2.cross(vector1, vector2);
722//              System.out.println("Intermed. refVector: " + vector2);
723                vector2.normalize();
724//              referenceVector.cross(referenceVector, principalRotationVector);
725                vector2.cross(vector1, vector2);
726                vector2.normalize();
727                if (ref.dot(vector2) < 0) {
728                        vector2.negate();
729                }
730//              System.out.println("Mod. refVector: " + vector2);
731                return vector2;
732        }
733        /**
734         * Returns the default reference vector for the alignment of Cn structures
735         * @return
736         */
737        private Vector3d getReferenceAxisCylic() {
738                if (rotationGroup.getPointGroup().equals("C2")) {
739                        Vector3d vr = new Vector3d(subunits.getOriginalCenters().get(0));
740                        vr.sub(subunits.getCentroid());
741                        vr.normalize();
742                        return vr;
743                }
744
745                // get principal axis vector that is perpendicular to the principal
746                // rotation vector
747                Vector3d vmin = null;
748                double dotMin = 1.0;
749                for (Vector3d v: principalAxesOfInertia) {
750                        if (Math.abs(principalRotationVector.dot(v)) < dotMin) {
751                                dotMin = Math.abs(principalRotationVector.dot(v));
752                                vmin = new Vector3d(v);
753                        }
754                }
755                if (principalRotationVector.dot(vmin) < 0) {
756                        vmin.negate();
757                }
758
759                return vmin;
760        }
761
762
763        /**
764         * Returns a reference vector for the alignment of Cn structures.
765         * @return reference vector
766         */
767        private Vector3d getReferenceAxisCylicWithSubunitAlignment() {
768                if (rotationGroup.getPointGroup().equals("C2")) {
769                        return referenceVector;
770                }
771
772                // find subunit that extends the most in the xy-plane
773                List<List<Integer>> orbits = getOrbitsByXYWidth();
774                // get the last orbit which is the widest
775                List<Integer> widestOrbit = orbits.get(orbits.size()-1);
776                List<Point3d> centers = subunits.getCenters();
777                int subunit = widestOrbit.get(0);
778
779                // calculate reference vector
780                Vector3d refAxis = new Vector3d();
781                refAxis.sub(centers.get(subunit), subunits.getCentroid());
782                refAxis.normalize();
783                return refAxis;
784        }
785
786        /**
787         *
788         */
789        private Vector3d getReferenceAxisDihedralWithSubunitAlignment() {
790                int maxFold = rotationGroup.getRotation(0).getFold();
791
792                double minAngle = Double.MAX_VALUE;
793                Vector3d refVec = null;
794
795                Vector3d ref = getSubunitReferenceVector();
796
797                for (int i = 0; i < rotationGroup.getOrder(); i++) {
798                        if (rotationGroup.getRotation(i).getDirection() == 1 &&
799                                        (rotationGroup.getRotation(i).getFold() < maxFold) ||
800                                        rotationGroup.getPointGroup().equals("D2")) {
801
802                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
803                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
804                                v.normalize();
805
806//                              System.out.println("Ref axis angle(+): " + Math.toDegrees(v.angle(ref)));
807                                double angle =  v.angle(ref);
808                                if (angle < minAngle) {
809                                        minAngle = angle;
810                                        refVec = v;
811                                }
812                                Vector3d vn = new Vector3d(v);
813                                vn.negate();
814//                              System.out.println("Ref axis angle(-): " + Math.toDegrees(vn.angle(ref)));
815                                angle =  vn.angle(ref);
816                                if (angle < minAngle) {
817                                        minAngle = angle;
818                                        refVec = vn;
819                                }
820                        }
821                }
822                refVec.normalize();
823                return refVec;
824        }
825
826        /**
827         *
828         */
829        private Vector3d getReferenceAxisDihedral() {
830                int maxFold = rotationGroup.getRotation(0).getFold();
831                // one exception: D2
832                if (maxFold == 2) {
833                        maxFold = 3;
834                }
835                // TODO how about D2, where minor axis = 2 = principal axis??
836                for (int i = 0; i < rotationGroup.getOrder(); i++) {
837                        if (rotationGroup.getRotation(i).getDirection() == 1 && rotationGroup.getRotation(i).getFold() < maxFold) {
838                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
839                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
840                                v.normalize();
841                                return v;
842                        }
843                }
844                return null;
845        }
846
847        private Vector3d getReferenceAxisTetrahedral() {
848                for (int i = 0; i < rotationGroup.getOrder(); i++) {
849                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
850                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
851                                double d = v.dot(principalRotationVector);
852                                if (rotationGroup.getRotation(i).getFold() == 3) {
853                                        // the dot product 0 is between to adjacent 3-fold axes
854                                        if (d > 0.3 && d < 0.9) {
855                                                return v;
856                                        }
857                                }
858                }
859                return null;
860        }
861
862        private Vector3d getReferenceAxisOctahedral() {
863                for (int i = 0; i < rotationGroup.getOrder(); i++) {
864                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
865                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
866                                double d = v.dot(principalRotationVector);
867                                if (rotationGroup.getRotation(i).getFold() == 4) {
868                                        // the dot product 0 is between to adjacent 4-fold axes
869                                        if (d > -0.1 && d < 0.1 ) {
870                                                return v;
871                                        }
872                                }
873                }
874                return null;
875        }
876
877        private Vector3d getReferenceAxisIcosahedral() {
878                for (int i = 0; i < rotationGroup.getOrder(); i++) {
879                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
880                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
881                                double d = v.dot(principalRotationVector);
882                                if (rotationGroup.getRotation(i).getFold() == 5) {
883                                        // the dot product of 0.447.. is between to adjacent 5-fold axes
884//                                      if (d > 0.447 && d < 0.448) {
885                                        if (d > 0.4 && d < 0.5) {
886                                                return v;
887                                        }
888                                }
889                }
890                return null;
891        }
892
893        private Vector3d getSubunitReferenceVector() {
894                        int n = subunits.getSubunitCount();
895                        Point3d probe = new Point3d();
896
897                        // transform subunit centers into z-aligned position and calculate
898                        // width in xy direction.
899                        double maxWidthSq = 0;
900                        Point3d ref = null;
901                        for (int i = 0; i < n; i++) {
902                                for (Point3d p: subunits.getTraces().get(i)) {
903                                        probe.set(p);
904                                        transformationMatrix.transform(probe);
905                                        double widthSq = probe.x*probe.x + probe.y*probe.y;
906                                        if (widthSq > maxWidthSq) {
907                                                maxWidthSq = widthSq;
908                                                ref = p;
909                                        }
910                                }
911                        }
912        //              System.out.println("width: " + maxWidthSq);
913                        Vector3d refVector = new Vector3d();
914                        refVector.sub(ref, subunits.getCentroid());
915                        refVector.normalize();
916                        return refVector;
917                }
918
919        private static Matrix4d flipX() {
920                Matrix4d rot = new Matrix4d();
921                rot.m00 = 1;
922                rot.m11 = -1;
923                rot.m22 = -1;
924                rot.m33 = 1;
925                return rot;
926        }
927
928        private static Matrix4d flipY() {
929                Matrix4d rot = new Matrix4d();
930                rot.m00 = -1;
931                rot.m11 = 1;
932                rot.m22 = -1;
933                rot.m33 = 1;
934                return rot;
935        }
936
937        private static Matrix4d flipZ() {
938                Matrix4d rot = new Matrix4d();
939                rot.m00 = -1;
940                rot.m11 = -1;
941                rot.m22 = 1;
942                rot.m33 = 1;
943                return rot;
944        }
945
946}