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.geometry.CalcPoint;
025import org.biojava.nbio.structure.geometry.UnitQuaternions;
026
027import javax.vecmath.AxisAngle4d;
028import javax.vecmath.Matrix4d;
029import javax.vecmath.Point3d;
030import javax.vecmath.Quat4d;
031import javax.vecmath.Vector3d;
032
033import java.util.ArrayList;
034import java.util.Arrays;
035import java.util.List;
036
037/**
038 *
039 * @author Peter
040 */
041public class C2RotationSolver implements QuatSymmetrySolver {
042        private QuatSymmetrySubunits subunits;
043        private QuatSymmetryParameters parameters;
044        private Vector3d centroid = new Vector3d();
045        private Matrix4d centroidInverse = new Matrix4d();
046
047        private RotationGroup rotations = new RotationGroup();
048
049
050        public C2RotationSolver(QuatSymmetrySubunits subunits, QuatSymmetryParameters parameters) {
051                if (subunits.getSubunitCount() != 2) {
052                        throw new IllegalArgumentException("C2RotationSolver can only be applied to cases with 2 centers");
053                }
054                this.subunits = subunits;
055                this.parameters = parameters;
056        }
057
058        @Override
059        public RotationGroup getSymmetryOperations() {
060                if (rotations.getOrder() == 0) {
061                        solve();
062                }
063                return rotations;
064        }
065
066        private void solve() {
067                initialize();
068                Vector3d trans = new Vector3d(subunits.getCentroid());
069                trans.negate();
070                List<Point3d[]> traces = subunits.getTraces();
071
072//              Point3d[] x = SuperPosition.clonePoint3dArray(traces.get(0));
073//              SuperPosition.center(x);
074//              Point3d[] y = SuperPosition.clonePoint3dArray(traces.get(1));
075//              SuperPosition.center(y);
076
077                Point3d[] x = CalcPoint.clonePoint3dArray(traces.get(0));
078                CalcPoint.translate(trans, x);
079                Point3d[] y = CalcPoint.clonePoint3dArray(traces.get(1));
080                CalcPoint.translate(trans, y);
081
082                // TODO implement this piece of code using at origin superposition
083                Quat4d quat = UnitQuaternions.relativeOrientation(
084                                x, y);
085                AxisAngle4d axisAngle = new AxisAngle4d();
086                Matrix4d transformation = new Matrix4d();
087
088                transformation.set(quat);
089                axisAngle.set(quat);
090
091                Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
092                if (axis.lengthSquared() < 1.0E-6) {
093                        axisAngle.x = 0;
094                        axisAngle.y = 0;
095                        axisAngle.z = 1;
096                        axisAngle.angle = 0;
097                } else {
098                        axis.normalize();
099                        axisAngle.x = axis.x;
100                        axisAngle.y = axis.y;
101                        axisAngle.z = axis.z;
102                }
103
104                CalcPoint.transform(transformation, y);
105
106                // if rmsd or angle deviation is above threshold, stop
107                double angleThresholdRadians = Math.toRadians(parameters.getAngleThreshold());
108                double deltaAngle = Math.abs(Math.PI-axisAngle.angle);
109
110                if (deltaAngle > angleThresholdRadians) {
111                        rotations.setC1(subunits.getSubunitCount());
112                        return;
113                }
114
115                // add unit operation
116                addEOperation();
117
118                // add C2 operation
119                int fold = 2;
120                combineWithTranslation(transformation);
121                List<Integer> permutation = Arrays.asList(1,0);
122                QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation);
123                scores.setRmsdCenters(0.0); // rmsd for superposition of two subunits centers is zero by definition
124
125                if (scores.getRmsd() > parameters.getRmsdThreshold() || deltaAngle > angleThresholdRadians) {
126                        rotations.setC1(subunits.getSubunitCount());
127                        return;
128                }
129
130                Rotation symmetryOperation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
131                rotations.addRotation(symmetryOperation);
132        }
133
134        private void addEOperation() {
135                List<Integer> permutation = Arrays.asList(0,1);
136                Matrix4d transformation = new Matrix4d();
137                transformation.setIdentity();
138                combineWithTranslation(transformation);
139                AxisAngle4d axisAngle = new AxisAngle4d();
140                QuatSymmetryScores scores = new QuatSymmetryScores();
141                int fold = 1; // ??
142                Rotation rotation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
143                rotations.addRotation(rotation);
144        }
145
146        /**
147         * Adds translational component to rotation matrix
148         * @param rotation
149         * @return
150         */
151        private void combineWithTranslation(Matrix4d rotation) {
152                rotation.setTranslation(centroid);
153                rotation.mul(rotation, centroidInverse);
154        }
155
156        private Rotation createSymmetryOperation(List<Integer> permutation, Matrix4d transformation, AxisAngle4d axisAngle, int fold, QuatSymmetryScores scores) {
157                Rotation s = new Rotation();
158                s.setPermutation(new ArrayList<Integer>(permutation));
159                s.setTransformation(new Matrix4d(transformation));
160                s.setAxisAngle(new AxisAngle4d(axisAngle));
161                s.setFold(fold);
162                s.setScores(scores);
163                return s;
164        }
165
166        private void initialize() {
167                // translation to centered coordinate system
168                centroid = new Vector3d(subunits.getCentroid());
169           // translation back to original coordinate system
170                Vector3d reverse = new Vector3d(centroid);
171                reverse.negate();
172                centroidInverse.set(reverse);
173//        // On LINUX there seems to be a bug with vecmath, and element m33 is zero. Here we make sure it's 1.
174                centroidInverse.setElement(3, 3, 1.0);
175        }
176
177}