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}