/*
 *                    BioJava development code
 *
 * This code may be freely distributed and modified under the
 * terms of the GNU Lesser General Public Licence.  This should
 * be distributed with the code.  If you do not have a copy,
 * see:
 *
 *      http://www.gnu.org/copyleft/lesser.html
 *
 * Copyright for this code is held jointly by the individual
 * authors.  These should be listed in @author doc comments.
 *
 * For more information on the BioJava project and its aims,
 * or to join the biojava-l mailing list, visit the home page
 * at:
 *
 *      http://www.biojava.org/
 *
 */
package org.biojava.nbio.structure.symmetry.core;

import org.biojava.nbio.structure.symmetry.geometry.MomentsOfInertia;
import org.biojava.nbio.structure.symmetry.geometry.SuperPosition;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import javax.vecmath.*;
import java.util.*;

public class RotationAxisAligner extends AxisAligner{
        
        private static final Logger LOGGER = LoggerFactory.getLogger(RotationAxisAligner.class);
        
        private static final Vector3d X_AXIS = new Vector3d(1,0,0);
        private static final Vector3d Y_AXIS = new Vector3d(0,1,0);
        private static final Vector3d Z_AXIS = new Vector3d(0,0,1);

        private Subunits subunits = null;
        private RotationGroup rotationGroup = null;

        private Matrix4d transformationMatrix = new Matrix4d();
        private Matrix4d reverseTransformationMatrix = new Matrix4d();
        private Vector3d referenceVector = new Vector3d();
        private Vector3d principalRotationVector = new Vector3d();
        private Vector3d[] principalAxesOfInertia = null;
         List<List<Integer>> alignedOrbits = null;

        private Vector3d minBoundary = new Vector3d();
        private Vector3d maxBoundary = new Vector3d();
        private double xyRadiusMax = Double.MIN_VALUE;

        boolean modified = true;

        public RotationAxisAligner(QuatSymmetryResults results) {
                this.subunits = results.getSubunits();
                this.rotationGroup = results.getRotationGroup();

                if (subunits == null) {
                        throw new IllegalArgumentException("AxisTransformation: Subunits are null");
                } else if (rotationGroup == null) {
                        throw new IllegalArgumentException("AxisTransformation: RotationGroup is null");
                } else if (subunits.getSubunitCount() == 0) {
                        throw new IllegalArgumentException("AxisTransformation: Subunits is empty");
                } else if (rotationGroup.getOrder() == 0) {
                        throw new IllegalArgumentException("AxisTransformation: RotationGroup is empty");
                }
        }

        /* (non-Javadoc)
         * @see org.biojava.nbio.structure.quaternary.core.AxisAligner#getTransformation()
         */
        @Override
        public String getSymmetry() {
                run();
                return rotationGroup.getPointGroup();
        }

        @Override
        public Matrix4d getTransformation() {
                run();
                return transformationMatrix;
        }

        @Override
        public Matrix3d getRotationMatrix() {
                run();
                Matrix3d m = new Matrix3d();
                transformationMatrix.getRotationScale(m);
                return m;
        }

        @Override
        public Matrix4d getReverseTransformation() {
                run();
                return reverseTransformationMatrix;
        }

        @Override
        public Vector3d getPrincipalRotationAxis() {
                run();
                return principalRotationVector;
        }

        @Override
        public Vector3d getRotationReferenceAxis() {
                run();
                return referenceVector;
        }

        @Override
        public Vector3d[] getPrincipalAxesOfInertia() {
                run();
                return principalAxesOfInertia;
        }

        @Override
        public Vector3d getDimension() {
                run();
                Vector3d dimension = new Vector3d();
                dimension.sub(maxBoundary, minBoundary);
                dimension.scale(0.5);
                return dimension;
        }

        /**
         * Returns the radius for drawing the minor rotation axis in the xy-plane
         * @return double radius in xy-plane
         */
        @Override
        public double getRadius() {
                run();
                return xyRadiusMax;
        }

        /**
         * Returns a transformation matrix transform polyhedra for Cn structures.
         * The center in this matrix is the geometric center, rather then the centroid.
         * In Cn structures those are usually not the same.
         * @return
         */
        @Override
        public Matrix4d getGeometicCenterTransformation() {
                run();

                Matrix4d geometricCentered = new Matrix4d(reverseTransformationMatrix);
                geometricCentered.setTranslation(new Vector3d(getGeometricCenter()));

                return geometricCentered;
        }

        /**
         * Returns the geometric center of polyhedron. In the case of the Cn
         * point group, the centroid and geometric center are usually not
         * identical.
         * @return
         */
        @Override
        public Point3d getGeometricCenter() {
                run();

                Point3d geometricCenter = new Point3d();
                Vector3d translation = new Vector3d();
                reverseTransformationMatrix.get(translation);

                // calculate adjustment around z-axis and transform adjustment to
                //  original coordinate frame with the reverse transformation
                if (rotationGroup.getPointGroup().startsWith("C")) {
                        Vector3d corr = new Vector3d(0,0, minBoundary.z+getDimension().z);
                        reverseTransformationMatrix.transform(corr);
                        geometricCenter.set(corr);
                }

                geometricCenter.add(translation);
                return geometricCenter;
        }

        @Override
        public Point3d getCentroid() {
                return new Point3d(subunits.getCentroid());
        }

        @Override
        public Subunits getSubunits() {
                return subunits;
        }

        public RotationGroup getRotationGroup() {
                return rotationGroup;
        }

        @Override
        public List<List<Integer>> getOrbits() {
                return alignedOrbits;
        }

        /**
         * @return
         */

        private void run () {
                if (modified) {
                        calcPrincipalRotationVector();
                        calcPrincipalAxes();
                        // initial alignment with draft reference axis
                        calcReferenceVector();
                        calcTransformation();

                        // refine ref. axis for cyclic and dihedral systems
                        if ((rotationGroup.getPointGroup().startsWith("C") &&
                                        !rotationGroup.getPointGroup().startsWith("C2")) ||
                                (rotationGroup.getPointGroup().startsWith("D") &&
                                                        !rotationGroup.getPointGroup().startsWith("D2"))
                                        ) {
                                refineReferenceVector();
                                calcTransformation();
                        }
                        calcReverseTransformation();
                        calcBoundaries();
                        if (! rotationGroup.getPointGroup().equals("Helical")) {
                                calcAlignedOrbits();
                        }
                        modified = false;
                }
        }
        /**
         * Returns a list of orbits (an orbit is a cyclic permutation of subunit indices that are related
         * by a rotation around the principal rotation axis) ordered from the +z direction to the -z direction (z-depth).
         * Within an orbit, subunit indices are ordered with a clockwise rotation around the z-axis.
         * @return list of orbits ordered by z-depth
         */
        private void calcAlignedOrbits() {
                Map<Double, List<Integer>> depthMap = new TreeMap<Double, List<Integer>>();
                double[] depth = getSubunitZDepth();
                alignedOrbits = calcOrbits();

                // calculate the mean depth of orbit along z-axis
                for (List<Integer> orbit: alignedOrbits) {
                        // calculate the mean depth along z-axis for each orbit
                        double meanDepth = 0;
                        for (int subunit: orbit) {
                                meanDepth += depth[subunit];
                        }
                        meanDepth /= orbit.size();

                        if (depthMap.get(meanDepth) != null) {
                                // System.out.println("Conflict in depthMap");
                                meanDepth += 0.01;
                        }
                        depthMap.put(meanDepth, orbit);
                }

                // now fill orbits back into list ordered by depth
                alignedOrbits.clear();
                for (List<Integer> orbit: depthMap.values()) {
                        // order subunit in a clockwise rotation around the z-axis
                        /// starting at the 12 O-clock position (+y position)
                        alignWithReferenceAxis(orbit);
                        alignedOrbits.add(orbit);
                }
        }

        /**
         * Returns an ordered list of subunit ids (orbit) in such a way that the subunit
         * indices start at the 12 o-clock (+y axis) and proceed in a clockwise direction
         * to the 11 o-clock position to close the "orbit".
         *
         * @param orbit list of subunit indices that are transformed into each other by a rotation
         * @return list of subunit indices ordered in a clockwise direction
         */

        private List<Integer> alignWithReferenceAxis(List<Integer> orbit) {
                int n = subunits.getSubunitCount();
                int fold = rotationGroup.getRotation(0).getFold();
                if (fold < 2) {
                        return orbit;
                }
                Vector3d probe = new Vector3d();

                double dotMin1 = Double.MIN_VALUE;
                double dotMin2 = Double.MIN_VALUE;
                int index1 = 0;
                int index2 = 0;
                Vector3d Y1 = new Vector3d(0,1,0);
                Vector3d Y2 = new Vector3d(0,1,0);
                Matrix3d m = new Matrix3d();
                double angle = -2*Math.PI/fold;
                m.rotZ(0.1*angle); // add small offset, since two subunits may be equidistant to the y-axis
                m.transform(Y1);
                m.rotZ(1.1*angle);
                m.transform(Y2);
                // transform subunit centers into z-aligned position and calculate
                // width in xy direction.
                for (int i: orbit) {
                        Point3d p = subunits.getCenters().get(i);
                        probe.set(p);
                        transformationMatrix.transform(probe);
                        // find subunit that lines up with y-axis
                        double dot1 = Y1.dot(probe);
                        if (dot1 > dotMin1) {
                                dotMin1 = dot1;
                                index1 = i;
                        }
                        // find next subunit (rotated by one fold around z-axis - clockwise)
                        double dot2 = Y2.dot(probe);
                        if (dot2 > dotMin2) {
                                dotMin2 = dot2;
                                index2 = i;
                        }
                }
//              System.out.println("Index1/2: " + index1 + " - " + index2);
//              System.out.println("Orbit0: " + orbit);
                // order subunit indices in a clockwise orientation around the z-axis
                // bring subunit into position 0
                for (int i = 0; i < n; i++) {
                        if (orbit.get(0) == index1) {
                                break;
                        }
                        Collections.rotate(orbit,1);
                }
//              System.out.println("Orbit1: " + orbit);
                // bring second subunit  onto position 1
                if (orbit.get(1) == index2) {
                        return orbit;
                }
                Collections.reverse(orbit.subList(1,  orbit.size()));
                if (orbit.get(1) != index2) {
                        LOGGER.warn("alignWithReferenceAxis failed");
                }
//              System.out.println("Orbit2: " + orbit);
                return orbit;
        }



        private void calcTransformation() {
                if (rotationGroup.getPointGroup().equals("C1")) {
                        calcTransformationByInertiaAxes();
                } else {
                        calcTransformationBySymmetryAxes();
                }
                // make sure this value is zero. On Linux this value is 0.
                transformationMatrix.setElement(3, 3, 1.0);
        }

        private void calcReverseTransformation() {
                reverseTransformationMatrix.invert(transformationMatrix);
        }

        private void calcTransformationBySymmetryAxes() {
                Vector3d[] axisVectors = new Vector3d[2];
                axisVectors[0] = new Vector3d(principalRotationVector);
                axisVectors[1] = new Vector3d(referenceVector);

                //  y,z axis centered at the centroid of the subunits
                Vector3d[] referenceVectors = new Vector3d[2];
                referenceVectors[0] = new Vector3d(Z_AXIS);
                referenceVectors[1] = new Vector3d(Y_AXIS);

                transformationMatrix = alignAxes(axisVectors, referenceVectors);

                // combine with translation
                Matrix4d combined = new Matrix4d();
                combined.setIdentity();
                Vector3d trans = new Vector3d(subunits.getCentroid());
                trans.negate();
                combined.setTranslation(trans);
                transformationMatrix.mul(combined);

                // for cyclic geometry, set a canonical view for the Z direction
                if (rotationGroup.getPointGroup().startsWith("C")) {
                        calcZDirection();
                }
        }

        private void calcTransformationByInertiaAxes() {
                Vector3d[] axisVectors = new Vector3d[2];
                axisVectors[0] = new Vector3d(principalAxesOfInertia[0]);
                axisVectors[1] = new Vector3d(principalAxesOfInertia[1]);

                Vector3d[] referenceVectors = new Vector3d[2];
                referenceVectors[0] = new Vector3d(Y_AXIS);
                referenceVectors[1] = new Vector3d(X_AXIS);

                // align inertia axes with y-x plane
                transformationMatrix = alignAxes(axisVectors, referenceVectors);

                // combine with translation
                Matrix4d translation = new Matrix4d();
                translation.setIdentity();
                Vector3d trans = new Vector3d(subunits.getCentroid());
                trans.negate();
                translation.setTranslation(trans);
                transformationMatrix.mul(translation);
        }

        /**
         * Returns a transformation matrix that rotates refPoints to match
         * coordPoints
         * @param refPoints the points to be aligned
         * @param referenceVectors
         * @return
         */
        private Matrix4d alignAxes(Vector3d[] axisVectors, Vector3d[] referenceVectors) {
                Matrix4d m1 = new Matrix4d();
                AxisAngle4d a = new AxisAngle4d();
                Vector3d axis = new Vector3d();

                // calculate rotation matrix to rotate refPoints[0] into coordPoints[0]
                Vector3d v1 = new Vector3d(axisVectors[0]);
                Vector3d v2 = new Vector3d(referenceVectors[0]);
                double dot = v1.dot(v2);
                if (Math.abs(dot) < 0.999) {
                        axis.cross(v1,v2);
                        axis.normalize();
                        a.set(axis, v1.angle(v2));
                        m1.set(a);
                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
                        m1.setElement(3,  3, 1.0);
                } else if (dot > 0) {
                        // parallel axis, nothing to do -> identity matrix
                        m1.setIdentity();
                } else if (dot < 0) {
                        // anti-parallel axis, flip around x-axis
                        m1.set(flipX());
                }

                // apply transformation matrix to all refPoints
                m1.transform(axisVectors[0]);
                m1.transform(axisVectors[1]);

                // calculate rotation matrix to rotate refPoints[1] into coordPoints[1]
                v1 = new Vector3d(axisVectors[1]);
                v2 = new Vector3d(referenceVectors[1]);
                Matrix4d m2 = new Matrix4d();
                dot = v1.dot(v2);
                if (Math.abs(dot) < 0.999) {
                        axis.cross(v1,v2);
                        axis.normalize();
                        a.set(axis, v1.angle(v2));
                        m2.set(a);
                        // make sure matrix element m33 is 1.0. It's 0 on Linux.
                        m2.setElement(3,  3, 1.0);
                } else if (dot > 0) {
                        // parallel axis, nothing to do -> identity matrix
                        m2.setIdentity();
                } else if (dot < 0) {
                        // anti-parallel axis, flip around z-axis
                        m2.set(flipZ());
                }

                // apply transformation matrix to all refPoints
                m2.transform(axisVectors[0]);
                m2.transform(axisVectors[1]);

                // combine the two rotation matrices
                m2.mul(m1);

                // the RMSD should be close to zero
                Point3d[] axes = new Point3d[2];
                axes[0] = new Point3d(axisVectors[0]);
                axes[1] = new Point3d(axisVectors[1]);
                Point3d[] ref = new Point3d[2];
                ref[0] = new Point3d(referenceVectors[0]);
                ref[1] = new Point3d(referenceVectors[1]);
                if (SuperPosition.rmsd(axes, ref) > 0.1) {
                        LOGGER.info("AxisTransformation: axes alignment is off. RMSD: " + SuperPosition.rmsd(axes, ref));
                }

                return m2;
        }

        private void calcPrincipalAxes() {
                MomentsOfInertia moi = new MomentsOfInertia();

                for (Point3d[] list: subunits.getTraces()) {
                        for (Point3d p: list) {
                                moi.addPoint(p, 1.0);
                        }
                }
                principalAxesOfInertia = moi.getPrincipalAxes();
        }

        /**
         * Calculates the min and max boundaries of the structure after it has been
         * transformed into its canonical orientation.
         */
        private void calcBoundaries() {
                minBoundary.x = Double.MAX_VALUE;
                maxBoundary.x = Double.MIN_VALUE;
                minBoundary.y = Double.MAX_VALUE;
                maxBoundary.x = Double.MIN_VALUE;
                minBoundary.z = Double.MAX_VALUE;
                maxBoundary.z = Double.MIN_VALUE;

                Point3d probe = new Point3d();

                for (Point3d[] list: subunits.getTraces()) {
                        for (Point3d p: list) {
                                probe.set(p);
                                transformationMatrix.transform(probe);

                                minBoundary.x = Math.min(minBoundary.x, probe.x);
                                maxBoundary.x = Math.max(maxBoundary.x, probe.x);
                                minBoundary.y = Math.min(minBoundary.y, probe.y);
                                maxBoundary.y = Math.max(maxBoundary.y, probe.y);
                                minBoundary.z = Math.min(minBoundary.z, probe.z);
                                maxBoundary.z = Math.max(maxBoundary.z, probe.z);
                                xyRadiusMax = Math.max(xyRadiusMax, Math.sqrt(probe.x*probe.x + probe.y * probe.y));
                        }
                }
        }

        /*
         * Modifies the rotation part of the transformation axis for
         * a Cn symmetric complex, so that the narrower end faces the
         * viewer, and the wider end faces away from the viewer. Example: 3LSV
         */
        private void calcZDirection() {
                calcBoundaries();

                // if the longer part of the structure faces towards the back (-z direction),
                // rotate around y-axis so the longer part faces the viewer (+z direction)
                if (Math.abs(minBoundary.z) > Math.abs(maxBoundary.z)) {
                        Matrix4d rot = flipY();
                        rot.mul(transformationMatrix);
                        transformationMatrix.set(rot);
                }
        }

        /**
         *
         */
        private List<List<Integer>> getOrbitsByXYWidth() {
                Map<Double, List<Integer>> widthMap = new TreeMap<Double, List<Integer>>();
                double[] width = getSubunitXYWidth();
                List<List<Integer>> orbits = calcOrbits();

                // calculate the mean width of orbits in XY-plane
                for (List<Integer> orbit: orbits) {
                        double meanWidth = 0;
                        for (int subunit: orbit) {
                                meanWidth += width[subunit];
                        }
                        meanWidth /= orbit.size();

                        if (widthMap.get(meanWidth) != null) {
                                meanWidth += 0.01;
                        }
                        widthMap.put(meanWidth, orbit);
                }

                // now fill orbits back into list ordered by width
                orbits.clear();
                for (List<Integer> orbit: widthMap.values()) {
                        orbits.add(orbit);
                }
                return orbits;
        }

        private double[] getSubunitXYWidth() {
                int n = subunits.getSubunitCount();
                double[] width = new double[n];
                Point3d probe = new Point3d();

                // transform subunit centers into z-aligned position and calculate
                // width in xy direction.
                for (int i = 0; i < n; i++) {
                        width[i] = Double.MIN_VALUE;
                        for (Point3d p: subunits.getTraces().get(i)) {
                                probe.set(p);
                                transformationMatrix.transform(probe);
                                width[i] = Math.max(width[i], Math.sqrt(probe.x*probe.x + probe.y*probe.y));
                        }
                }
                return width;
        }

        private double[] getSubunitZDepth() {
                int n = subunits.getSubunitCount();
                double[] depth = new double[n];
                Point3d probe = new Point3d();

                // transform subunit centers into z-aligned position and calculate
                // z-coordinates (depth) along the z-axis.
                for (int i = 0; i < n; i++) {
                        Point3d p= subunits.getCenters().get(i);
                        probe.set(p);
                        transformationMatrix.transform(probe);
                        depth[i] = probe.z;
                }
                return depth;
        }

        /**
         * Returns a list of list of subunit ids that form an "orbit", i.e. they
         * are transformed into each other during a rotation around the principal symmetry axis (z-axis)
         * @return
         */
        private List<List<Integer>> calcOrbits() {
                int n = subunits.getSubunitCount();
                int fold = rotationGroup.getRotation(0).getFold();

                List<List<Integer>> orbits = new ArrayList<List<Integer>>();
                boolean[] used = new boolean[n];
                Arrays.fill(used, false);

                for (int i = 0; i < n; i++) {
                        if (! used[i]) {
                                // determine the equivalent subunits
                                List<Integer> orbit = new ArrayList<Integer>(fold);
                                for (int j = 0; j < fold; j++) {
                                        List<Integer> permutation = rotationGroup.getRotation(j).getPermutation();
                                        orbit.add(permutation.get(i));
                                        used[permutation.get(i)] = true;
                                }
                                orbits.add(deconvolute(orbit));
                        }
                }
                return orbits;
        }

        private List<Integer> deconvolute(List<Integer> orbit) {
                if (rotationGroup.getOrder() < 2) {
                        return orbit;
                }
                List<Integer> p0 = rotationGroup.getRotation(0).getPermutation();
                List<Integer> p1 = rotationGroup.getRotation(1).getPermutation();
//              System.out.println("deconvolute");
//              System.out.println("Permutation0: " + p0);
//              System.out.println("Permutation1: " + p1);

                List<Integer> inRotationOrder = new ArrayList<Integer>(orbit.size());
                inRotationOrder.add(orbit.get(0));
                for (int i = 1; i < orbit.size(); i++) {
                        int current = inRotationOrder.get(i-1);
                        int index = p0.indexOf(current);
                        int next = p1.get(index);
                        if (!orbit.contains(next)) {
                                LOGGER.warn("deconvolute: inconsistency in permuation. Returning original order");
                                return orbit;
                        }
                        inRotationOrder.add(next);
                }
//              System.out.println("In order: " + inRotationOrder);
                return inRotationOrder;
        }

        /**
         * Returns a vector along the principal rotation axis for the
         * alignment of structures along the z-axis
         * @return principal rotation vector
         */
        private void calcPrincipalRotationVector() {
                Rotation rotation = rotationGroup.getRotation(0); // the rotation around the principal axis is the first rotation
                AxisAngle4d axisAngle = rotation.getAxisAngle();
                principalRotationVector = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
        }

        /**
         * Returns a vector perpendicular to the principal rotation vector
         * for the alignment of structures in the xy-plane
         * @return reference vector
         */
        private void calcReferenceVector() {
                referenceVector = null;
                if (rotationGroup.getPointGroup().startsWith("C")) {
                        referenceVector = getReferenceAxisCylic();
                } else if (rotationGroup.getPointGroup().startsWith("D")) {
                        referenceVector = getReferenceAxisDihedral();
                } else if (rotationGroup.getPointGroup().equals("T")) {
                        referenceVector = getReferenceAxisTetrahedral();
                } else if (rotationGroup.getPointGroup().equals("O")) {
                        referenceVector = getReferenceAxisOctahedral();
                } else if (rotationGroup.getPointGroup().equals("I")) {
                        referenceVector = getReferenceAxisIcosahedral();
                } else if (rotationGroup.getPointGroup().equals("Helical")) {
                        // TODO what should the reference vector be??
                        referenceVector = getReferenceAxisCylic();
                }

                if (referenceVector == null) {
                        LOGGER.warn("no reference vector found. Using y-axis.");
                        referenceVector = new Vector3d(Y_AXIS);
                }
                // make sure reference vector is perpendicular principal roation vector
                referenceVector = orthogonalize(principalRotationVector, referenceVector);
        }

        /**
         * Returns a normalized vector that represents a minor rotation axis, except
         * for Cn, this represents an axis orthogonal to the principal axis.
         * @return minor rotation axis
         */
        private void refineReferenceVector() {
                referenceVector = new Vector3d(Y_AXIS);
                if (rotationGroup.getPointGroup().startsWith("C")) {
                        referenceVector = getReferenceAxisCylicWithSubunitAlignment();
                } else if (rotationGroup.getPointGroup().startsWith("D")) {
                        referenceVector = getReferenceAxisDihedralWithSubunitAlignment();
                }

                referenceVector = orthogonalize(principalRotationVector, referenceVector);
        }

        private Vector3d orthogonalize(Vector3d vector1, Vector3d vector2) {
                double dot = vector1.dot(vector2);
                Vector3d ref = new Vector3d(vector2);
//              System.out.println("p.r: " + dot);
//              System.out.println("Orig refVector: " + referenceVector);
                if (dot < 0) {
                        vector2.negate();
                }
                vector2.cross(vector1, vector2);
//              System.out.println("Intermed. refVector: " + vector2);
                vector2.normalize();
//              referenceVector.cross(referenceVector, principalRotationVector);
                vector2.cross(vector1, vector2);
                vector2.normalize();
                if (ref.dot(vector2) < 0) {
                        vector2.negate();
                }
//              System.out.println("Mod. refVector: " + vector2);
                return vector2;
        }
        /**
         * Returns the default reference vector for the alignment of Cn structures
         * @return
         */
        private Vector3d getReferenceAxisCylic() {
                if (rotationGroup.getPointGroup().equals("C2")) {
                        Vector3d vr = new Vector3d(subunits.getOriginalCenters().get(0));
                        vr.sub(subunits.getCentroid());
                        vr.normalize();
                        return vr;
                }

                // get principal axis vector that is perpendicular to the principal
                // rotation vector
                Vector3d vmin = null;
                double dotMin = 1.0;
                for (Vector3d v: principalAxesOfInertia) {
                        if (Math.abs(principalRotationVector.dot(v)) < dotMin) {
                                dotMin = Math.abs(principalRotationVector.dot(v));
                                vmin = new Vector3d(v);
                        }
                }
                if (principalRotationVector.dot(vmin) < 0) {
                        vmin.negate();
                }

                return vmin;
        }


        /**
         * Returns a reference vector for the alignment of Cn structures.
         * @return reference vector
         */
        private Vector3d getReferenceAxisCylicWithSubunitAlignment() {
                if (rotationGroup.getPointGroup().equals("C2")) {
                        return referenceVector;
                }

                // find subunit that extends the most in the xy-plane
                List<List<Integer>> orbits = getOrbitsByXYWidth();
                // get the last orbit which is the widest
                List<Integer> widestOrbit = orbits.get(orbits.size()-1);
                List<Point3d> centers = subunits.getCenters();
                int subunit = widestOrbit.get(0);

                // calculate reference vector
                Vector3d refAxis = new Vector3d();
                refAxis.sub(centers.get(subunit), subunits.getCentroid());
                refAxis.normalize();
                return refAxis;
        }

        /**
         *
         */
        private Vector3d getReferenceAxisDihedralWithSubunitAlignment() {
                int maxFold = rotationGroup.getRotation(0).getFold();

                double minAngle = Double.MAX_VALUE;
                Vector3d refVec = null;

                Vector3d ref = getSubunitReferenceVector();

                for (int i = 0; i < rotationGroup.getOrder(); i++) {
                        if (rotationGroup.getRotation(i).getDirection() == 1 &&
                                        (rotationGroup.getRotation(i).getFold() < maxFold) ||
                                        rotationGroup.getPointGroup().equals("D2")) {

                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
                                v.normalize();

//                              System.out.println("Ref axis angle(+): " + Math.toDegrees(v.angle(ref)));
                                double angle =  v.angle(ref);
                                if (angle < minAngle) {
                                        minAngle = angle;
                                        refVec = v;
                                }
                                Vector3d vn = new Vector3d(v);
                                vn.negate();
//                              System.out.println("Ref axis angle(-): " + Math.toDegrees(vn.angle(ref)));
                                angle =  vn.angle(ref);
                                if (angle < minAngle) {
                                        minAngle = angle;
                                        refVec = vn;
                                }
                        }
                }
                refVec.normalize();
                return refVec;
        }

        /**
         *
         */
        private Vector3d getReferenceAxisDihedral() {
                int maxFold = rotationGroup.getRotation(0).getFold();
                // one exception: D2
                if (maxFold == 2) {
                        maxFold = 3;
                }
                // TODO how about D2, where minor axis = 2 = principal axis??
                for (int i = 0; i < rotationGroup.getOrder(); i++) {
                        if (rotationGroup.getRotation(i).getDirection() == 1 && rotationGroup.getRotation(i).getFold() < maxFold) {
                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
                                v.normalize();
                                return v;
                        }
                }
                return null;
        }

        private Vector3d getReferenceAxisTetrahedral() {
                for (int i = 0; i < rotationGroup.getOrder(); i++) {
                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
                                double d = v.dot(principalRotationVector);
                                if (rotationGroup.getRotation(i).getFold() == 3) {
                                        // the dot product 0 is between to adjacent 3-fold axes
                                        if (d > 0.3 && d < 0.9) {
                                                return v;
                                        }
                                }
                }
                return null;
        }

        private Vector3d getReferenceAxisOctahedral() {
                for (int i = 0; i < rotationGroup.getOrder(); i++) {
                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
                                double d = v.dot(principalRotationVector);
                                if (rotationGroup.getRotation(i).getFold() == 4) {
                                        // the dot product 0 is between to adjacent 4-fold axes
                                        if (d > -0.1 && d < 0.1 ) {
                                                return v;
                                        }
                                }
                }
                return null;
        }

        private Vector3d getReferenceAxisIcosahedral() {
                for (int i = 0; i < rotationGroup.getOrder(); i++) {
                                AxisAngle4d axisAngle = rotationGroup.getRotation(i).getAxisAngle();
                                Vector3d v = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
                                double d = v.dot(principalRotationVector);
                                if (rotationGroup.getRotation(i).getFold() == 5) {
                                        // the dot product of 0.447.. is between to adjacent 5-fold axes
//                                      if (d > 0.447 && d < 0.448) {
                                        if (d > 0.4 && d < 0.5) {
                                                return v;
                                        }
                                }
                }
                return null;
        }

        private Vector3d getSubunitReferenceVector() {
                        int n = subunits.getSubunitCount();
                        Point3d probe = new Point3d();

                        // transform subunit centers into z-aligned position and calculate
                        // width in xy direction.
                        double maxWidthSq = 0;
                        Point3d ref = null;
                        for (int i = 0; i < n; i++) {
                                for (Point3d p: subunits.getTraces().get(i)) {
                                        probe.set(p);
                                        transformationMatrix.transform(probe);
                                        double widthSq = probe.x*probe.x + probe.y*probe.y;
                                        if (widthSq > maxWidthSq) {
                                                maxWidthSq = widthSq;
                                                ref = p;
                                        }
                                }
                        }
        //              System.out.println("width: " + maxWidthSq);
                        Vector3d refVector = new Vector3d();
                        refVector.sub(ref, subunits.getCentroid());
                        refVector.normalize();
                        return refVector;
                }

        private static Matrix4d flipX() {
                Matrix4d rot = new Matrix4d();
                rot.m00 = 1;
                rot.m11 = -1;
                rot.m22 = -1;
                rot.m33 = 1;
                return rot;
        }

        private static Matrix4d flipY() {
                Matrix4d rot = new Matrix4d();
                rot.m00 = -1;
                rot.m11 = 1;
                rot.m22 = -1;
                rot.m33 = 1;
                return rot;
        }

        private static Matrix4d flipZ() {
                Matrix4d rot = new Matrix4d();
                rot.m00 = -1;
                rot.m11 = -1;
                rot.m22 = 1;
                rot.m33 = 1;
                return rot;
        }

}