001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021
022package org.biojava.nbio.structure.symmetry.core;
023
024import org.biojava.nbio.structure.symmetry.geometry.SuperPosition;
025
026import javax.vecmath.AxisAngle4d;
027import javax.vecmath.Matrix4d;
028import javax.vecmath.Point3d;
029import javax.vecmath.Vector3d;
030import java.util.ArrayList;
031import java.util.Arrays;
032import java.util.List;
033
034/**
035 *
036 * @author Peter
037 */
038public class C2RotationSolver implements QuatSymmetrySolver {
039        private Subunits subunits = null;
040        private QuatSymmetryParameters parameters = null;
041        private Vector3d centroid = new Vector3d();
042        private Matrix4d centroidInverse = new Matrix4d();
043
044        private RotationGroup rotations = new RotationGroup();
045
046
047        public C2RotationSolver(Subunits subunits, QuatSymmetryParameters parameters) {
048                if (subunits.getSubunitCount() != 2) {
049                        throw new IllegalArgumentException("C2RotationSolver can only be applied to cases with 2 centers");
050                }
051                this.subunits = subunits;
052                this.parameters = parameters;
053        }
054
055        @Override
056        public RotationGroup getSymmetryOperations() {
057                if (rotations.getOrder() == 0) {
058                        solve();
059                }
060                return rotations;
061        }
062
063        private void solve() {
064                initialize();
065                Vector3d trans = new Vector3d(subunits.getCentroid());
066                trans.negate();
067                List<Point3d[]> traces = subunits.getTraces();
068
069//              Point3d[] x = SuperPosition.clonePoint3dArray(traces.get(0));
070//              SuperPosition.center(x);
071//              Point3d[] y = SuperPosition.clonePoint3dArray(traces.get(1));
072//              SuperPosition.center(y);
073
074                Point3d[] x = SuperPosition.clonePoint3dArray(traces.get(0));
075                SuperPosition.translate(new Point3d(trans), x);
076                Point3d[] y = SuperPosition.clonePoint3dArray(traces.get(1));
077                SuperPosition.translate(new Point3d(trans), y);
078
079                AxisAngle4d axisAngle = new AxisAngle4d();
080
081                Matrix4d transformation = SuperPosition.superposeAtOrigin(x, y, axisAngle);
082
083                // if rmsd or angle deviation is above threshold, stop
084                double angleThresholdRadians = Math.toRadians(parameters.getAngleThreshold());
085                double deltaAngle = Math.abs(Math.PI-axisAngle.angle);
086
087                if (deltaAngle > angleThresholdRadians) {
088                        rotations.setC1(subunits.getSubunitCount());
089                        return;
090                }
091
092                // add unit operation
093                addEOperation();
094
095                // add C2 operation
096                int fold = 2;
097                combineWithTranslation(transformation);
098                List<Integer> permutation = Arrays.asList(1,0);
099                QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation);
100                scores.setRmsdCenters(0.0); // rmsd for superposition of two subunits centers is zero by definition
101
102                if (scores.getRmsd() > parameters.getRmsdThreshold() || deltaAngle > angleThresholdRadians) {
103                        rotations.setC1(subunits.getSubunitCount());
104                        return;
105                }
106
107                Rotation symmetryOperation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
108                rotations.addRotation(symmetryOperation);
109        }
110
111        private void addEOperation() {
112                List<Integer> permutation = Arrays.asList(new Integer[]{0,1});
113                Matrix4d transformation = new Matrix4d();
114                transformation.setIdentity();
115                combineWithTranslation(transformation);
116                AxisAngle4d axisAngle = new AxisAngle4d();
117                QuatSymmetryScores scores = new QuatSymmetryScores();
118                int fold = 1; // ??
119                Rotation rotation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
120                rotations.addRotation(rotation);
121        }
122
123        /**
124         * Adds translational component to rotation matrix
125         * @param rotTrans
126         * @param rotation
127         * @return
128         */
129        private void combineWithTranslation(Matrix4d rotation) {
130                rotation.setTranslation(centroid);
131                rotation.mul(rotation, centroidInverse);
132        }
133
134        private Rotation createSymmetryOperation(List<Integer> permutation, Matrix4d transformation, AxisAngle4d axisAngle, int fold, QuatSymmetryScores scores) {
135                Rotation s = new Rotation();
136                s.setPermutation(new ArrayList<Integer>(permutation));
137                s.setTransformation(new Matrix4d(transformation));
138                s.setAxisAngle(new AxisAngle4d(axisAngle));
139                s.setFold(fold);
140                s.setScores(scores);
141                return s;
142        }
143
144        private void initialize() {
145                // translation to centered coordinate system
146                centroid = new Vector3d(subunits.getCentroid());
147           // translation back to original coordinate system
148                Vector3d reverse = new Vector3d(centroid);
149                reverse.negate();
150                centroidInverse.set(reverse);
151//        // On LINUX there seems to be a bug with vecmath, and element m33 is zero. Here we make sure it's 1.
152                centroidInverse.setElement(3, 3, 1.0);
153        }
154
155}