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.align.util; 022 023import java.io.StringWriter; 024import java.util.Locale; 025 026import javax.vecmath.AxisAngle4d; 027import javax.vecmath.Matrix3d; 028import javax.vecmath.Matrix4d; 029import javax.vecmath.Vector3d; 030 031import org.biojava.nbio.structure.Atom; 032import org.biojava.nbio.structure.AtomImpl; 033import org.biojava.nbio.structure.Calc; 034import org.biojava.nbio.structure.StructureException; 035import org.biojava.nbio.structure.align.model.AFPChain; 036import org.biojava.nbio.structure.contact.Pair; 037import org.biojava.nbio.structure.geometry.Matrices; 038import org.biojava.nbio.structure.jama.Matrix; 039 040/** 041 * Calculates the rotation axis for an alignment 042 * 043 * <p>A superposition of two structures is generally represented as a rotation 044 * matrix plus a translation vector. However, it can also be represented as an 045 * axis of rotation plus some translation. 046 * 047 * <p>This class calculates the rotation axis and stores it as four properties: 048 * <ul><li>A unit vector parallel to the rotation axis ({@link #getRotationAxis()}) 049 * <li>The angle of rotation ({@link #getAngle()}) 050 * <li>A point on the rotation axis ({@link #getRotationPos()}) 051 * <li>Some translation parallel to the axis ({@link #getScrewTranslation()}) 052 * </ul> 053 * 054 * <p>The axis of rotation is poorly defined and numerically unstable for small 055 * angles. Therefore it's direction is left as null for angles less than 056 * {@link #MIN_ANGLE}. 057 * 058 * @author Spencer Bliven 059 * 060 */ 061public final class RotationAxis { 062 063 /** 064 * Minimum angle to calculate rotation axes. 5 degrees. 065 */ 066 static final double MIN_ANGLE = Math.toRadians(5.); 067 068 private double theta; 069 private Atom rotationAxis; // axis of rotation 070 private Atom rotationPos; // a point on the axis of rotation 071 private Atom screwTranslation; //translation parallel to the axis of rotation 072 private Atom otherTranslation; // translation perpendicular to the axis of rotation 073 074 /** 075 * The rotation angle 076 * @return the angle, in radians 077 */ 078 public double getAngle() { 079 return theta; 080 } 081 082 /** 083 * Get a unit vector along the rotation axis 084 * @return rotationAxis 085 */ 086 public Atom getRotationAxis() { 087 return rotationAxis; 088 } 089 090 /** 091 * Returns the rotation axis and angle in a single javax.vecmath.AxisAngle4d object 092 * @return 093 */ 094 public AxisAngle4d getAxisAngle4d() { 095 return new AxisAngle4d(rotationAxis.getX(),rotationAxis.getY(),rotationAxis.getZ(),theta); 096 } 097 098 /** 099 * Get a position on the rotation axis. 100 * 101 * Specifically, project the origin onto the rotation axis 102 * @return rotationPos 103 */ 104 public Atom getRotationPos() { 105 return rotationPos; 106 } 107 108 /** 109 * Get the component of translation parallel to the axis of rotation 110 * @return 111 */ 112 public Atom getScrewTranslation() { 113 return screwTranslation; 114 } 115 116 public Vector3d getVector3dScrewTranslation() { 117 return new Vector3d(screwTranslation.getX(),screwTranslation.getY(),screwTranslation.getZ()); 118 } 119 120 public double getTranslation() { 121 return Calc.amount(screwTranslation); 122 } 123 124 /** 125 * Calculate the rotation axis for the first block of an AFPChain 126 * @param afpChain 127 * @throws StructureException 128 * @throws NullPointerException if afpChain does not contain a valid rotation matrix and shift vector 129 */ 130 public RotationAxis(AFPChain afpChain) throws StructureException { 131 if(afpChain.getAlnLength() < 1) { 132 throw new StructureException("No aligned residues"); 133 } 134 init(afpChain.getBlockRotationMatrix()[0],afpChain.getBlockShiftVector()[0]); 135 } 136 137 /** 138 * Create a rotation axis from a vector, a point, and an angle. 139 * 140 * The result will be a pure rotation, with no screw component. 141 * @param axis A vector parallel to the axis of rotation 142 * @param pos A point on the axis of rotation 143 * @param theta The angle to rotate (radians) 144 */ 145 public RotationAxis(Atom axis, Atom pos, double theta) { 146 this.rotationAxis = Calc.unitVector(axis); 147 this.rotationPos = (Atom) pos.clone(); 148 this.theta = theta; 149 this.screwTranslation = new AtomImpl(); //zero 150 this.otherTranslation = null; //deprecated 151 } 152 153 /** 154 * Determine the location of the rotation axis based on a rotation matrix and a translation vector 155 * @param rotation 156 * @param translation 157 */ 158 public RotationAxis(Matrix rotation, Atom translation) { 159 init(rotation, translation); 160 } 161 162 /** 163 * Create a rotation axis from a Matrix4d containing a rotational 164 * component and a translational component. 165 * 166 * @param transform 167 */ 168 public RotationAxis(Matrix4d transform) { 169 170 Matrix rot = Matrices.getRotationJAMA(transform); 171 Atom transl = Calc.getTranslationVector(transform); 172 init(rot,transl); 173 } 174 175 /** 176 * Get the rotation matrix corresponding to this axis 177 * @return A 3x3 rotation matrix 178 */ 179 public Matrix getRotationMatrix() { 180 return getRotationMatrix(theta); 181 } 182 183 /** 184 * Get the rotation matrix corresponding to a rotation about this axis 185 * @param theta The amount to rotate 186 * @return A 3x3 rotation matrix 187 */ 188 public Matrix getRotationMatrix(double theta) { 189 if( rotationAxis == null) { 190 // special case for pure translational axes 191 return Matrix.identity(3, 3); 192 } 193 double x = rotationAxis.getX(); 194 double y = rotationAxis.getY(); 195 double z = rotationAxis.getZ(); 196 double cos = Math.cos(theta); 197 double sin = Math.sin(theta); 198 double com = 1 - cos; 199 return new Matrix(new double[][] { 200 {com*x*x + cos, com*x*y+sin*z, com*x*z+-sin*y}, 201 {com*x*y-sin*z, com*y*y+cos, com*y*z+sin*x}, 202 {com*x*z+sin*y, com*y*z-sin*x, com*z*z+cos}, 203 }); 204 } 205 206 /** 207 * Returns the rotation order o that gives the lowest value of {@code |2PI / o - theta}, 208 * given that the value is strictly lower than {@code threshold}, for orders {@code o=1,...,maxOrder}. 209 */ 210 public int guessOrderFromAngle(double threshold, int maxOrder) { 211 double bestDelta = threshold; 212 int bestOrder = 1; 213 for (int order = 2; order < maxOrder; order++) { 214 double delta = Math.abs(2 * Math.PI / order - theta); 215 if (delta < bestDelta) { 216 bestOrder = order; 217 bestDelta = delta; 218 } 219 } 220 return bestOrder; 221 } 222 223 224 /** 225 * Returns a matrix that describes both rotation and translation. 226 */ 227 public Matrix getFullMatrix() { 228 return null; // TODO, easy 229 } 230 231 /** 232 * Initialize variables 233 * 234 * @param rotation 235 * @param translation 236 */ 237 private void init(Matrix rotation, Atom translation) { 238 if(rotation.getColumnDimension() != 3 || rotation.getRowDimension() != 3) { 239 throw new IllegalArgumentException("Expected 3x3 rotation matrix"); 240 } 241 242 243 // Calculate angle 244 double c = (rotation.trace()-1)/2.0; //=cos(theta) 245 // c is sometimes slightly out of the [-1,1] range due to numerical instabilities 246 if( -1-1e-8 < c && c < -1 ) c = -1; 247 if( 1+1e-8 > c && c > 1 ) c = 1; 248 if( -1 > c || c > 1 ) { 249 throw new IllegalArgumentException("Input matrix is not a valid rotation matrix."); 250 } 251 this.theta = Math.acos(c); 252 253 if(theta < MIN_ANGLE) { 254 calculateTranslationalAxis(rotation,translation); 255 } else { 256 calculateRotationalAxis(rotation, translation, c); 257 } 258 } 259 260 /** 261 * Calculate the rotation axis for the normal case, where there is a 262 * significant rotation angle 263 * @param rotation 264 * @param translation 265 * @param c 266 */ 267 private void calculateRotationalAxis(Matrix rotation, Atom translation, 268 double c) { 269 // Calculate magnitude of rotationAxis components, but not signs 270 double sum=0; 271 double[] rotAx = new double[3]; 272 for(int i=0;i<3;i++) { 273 rotAx[i] = Math.sqrt(rotation.get(i, i)-c); 274 sum+=rotAx[i]*rotAx[i]; 275 } 276 for(int i=0;i<3;i++) { 277 rotAx[i] /= Math.sqrt(sum); 278 } 279 280 // Now determine signs 281 double d0 = rotation.get(2,1)-rotation.get(1,2); //=2u[0]*sin(theta) 282 double d1 = rotation.get(0,2)-rotation.get(2,0); //=2u[1]*sin(theta) 283 double d2 = rotation.get(1,0)-rotation.get(0,1); //=2u[2]*sin(theta) 284 285 double s12 = rotation.get(2,1)+rotation.get(1,2); //=2*u[1]*u[2]*(1-cos(theta)) 286 double s02 = rotation.get(0,2)+rotation.get(2,0); //=2*u[0]*u[2]*(1-cos(theta)) 287 double s01 = rotation.get(1,0)+rotation.get(0,1); //=2*u[0]*u[1]*(1-cos(theta)) 288 289 //Only use biggest d for the sign directly, for numerical stability around 180deg 290 if( Math.abs(d0) < Math.abs(d1) ) { // not d0 291 if( Math.abs(d1) < Math.abs(d2) ) { //d2 292 if(d2>=0){ //u[2] positive 293 if( s02 < 0 ) rotAx[0] = -rotAx[0]; 294 if( s12 < 0 ) rotAx[1] = -rotAx[1]; 295 } else { //u[2] negative 296 rotAx[2] = -rotAx[2]; 297 if( s02 >= 0 ) rotAx[0] = -rotAx[0]; 298 if( s12 >= 0 ) rotAx[1] = -rotAx[1]; 299 } 300 } else { //d1 301 if(d1>=0) {//u[1] positive 302 if( s01 < 0) rotAx[0] = -rotAx[0]; 303 if( s12 < 0) rotAx[2] = -rotAx[2]; 304 } else { //u[1] negative 305 rotAx[1] = -rotAx[1]; 306 if( s01 >= 0) rotAx[0] = -rotAx[0]; 307 if( s12 >= 0) rotAx[2] = -rotAx[2]; 308 } 309 } 310 } else { // not d1 311 if( Math.abs(d0) < Math.abs(d2) ) { //d2 312 if(d2>=0){ //u[2] positive 313 if( s02 < 0 ) rotAx[0] = -rotAx[0]; 314 if( s12 < 0 ) rotAx[1] = -rotAx[1]; 315 } else { //u[2] negative 316 rotAx[2] = -rotAx[2]; 317 if( s02 >= 0 ) rotAx[0] = -rotAx[0]; 318 if( s12 >= 0 ) rotAx[1] = -rotAx[1]; 319 } 320 } else { //d0 321 if(d0>=0) { //u[0] positive 322 if( s01 < 0 ) rotAx[1] = -rotAx[1]; 323 if( s02 < 0 ) rotAx[2] = -rotAx[2]; 324 } else { //u[0] negative 325 rotAx[0] = -rotAx[0]; 326 if( s01 >= 0 ) rotAx[1] = -rotAx[1]; 327 if( s02 >= 0 ) rotAx[2] = -rotAx[2]; 328 } 329 } 330 } 331 332 rotationAxis = new AtomImpl(); 333 rotationAxis.setCoords(rotAx); 334 335 // Calculate screw = (rotationAxis dot translation)*u 336 double dotProduct = Calc.scalarProduct(rotationAxis, translation); 337 338 screwTranslation = Calc.scale(rotationAxis, dotProduct); 339 otherTranslation = Calc.subtract(translation, screwTranslation); 340 341 Atom hypot = Calc.vectorProduct(otherTranslation,rotationAxis); 342 Calc.scaleEquals(hypot,.5/Math.tan(theta/2.0)); 343 344 // Calculate rotation axis position 345 rotationPos = Calc.scaleAdd(.5,otherTranslation, hypot); 346 } 347 348 /** 349 * Handle cases with small angles of rotation 350 * @param rotation 351 * @param translation 352 */ 353 private void calculateTranslationalAxis(Matrix rotation, Atom translation) { 354 // set axis parallel to translation 355 rotationAxis = Calc.scale(translation, 1./Calc.amount(translation)); 356 357 // position is undefined 358 rotationPos = null; 359 360 screwTranslation = translation; 361 otherTranslation = new AtomImpl(); 362 otherTranslation.setCoords(new double[] {0,0,0}); 363 } 364 365 /** 366 * Returns a Jmol script which will display the axis of rotation. This 367 * consists of a cyan arrow along the axis, plus an arc showing the angle 368 * of rotation. 369 * <p> 370 * As the rotation angle gets smaller, the axis of rotation becomes poorly 371 * defined and would need to get farther and farther away from the protein. 372 * This is not particularly useful, so we arbitrarily draw it parallel to 373 * the translation and omit the arc. 374 * @param atoms Some atoms from the protein, used for determining the bounds 375 * of the axis. 376 * 377 * @return The Jmol script, suitable for calls to 378 * {@link org.biojava.nbio.structure.align.gui.jmol.StructureAlignmentJmol#evalString() jmol.evalString()} 379 */ 380 public String getJmolScript(Atom[] atoms){ 381 return getJmolScript(atoms, 0); 382 } 383 384 /** 385 * Find a segment of the axis that covers the specified set of atoms. 386 * <p> 387 * Projects the input atoms onto the rotation axis and returns the bounding 388 * points. 389 * <p> 390 * In the case of a pure translational axis, the axis location is undefined 391 * so the center of mass will be used instead. 392 * @param atoms 393 * @return two points defining the axis segment 394 */ 395 public Pair<Atom> getAxisEnds(Atom[] atoms) { 396 // Project each Atom onto the rotation axis to determine limits 397 double min, max; 398 min = max = Calc.scalarProduct(rotationAxis,atoms[0]); 399 for(int i=1;i<atoms.length;i++) { 400 double prod = Calc.scalarProduct(rotationAxis,atoms[i]); 401 if(prod<min) min = prod; 402 if(prod>max) max = prod; 403 } 404 double uLen = Calc.scalarProduct(rotationAxis,rotationAxis);// Should be 1, but double check 405 min/=uLen; 406 max/=uLen; 407 408 // Project the origin onto the axis. If the axis is undefined, use the center of mass 409 Atom axialPt; 410 if(rotationPos == null) { 411 Atom center = Calc.centerOfMass(atoms); 412 413 // Project center onto the axis 414 Atom centerOnAxis = Calc.scale(rotationAxis, Calc.scalarProduct(center, rotationAxis)); 415 416 // Remainder is projection of origin onto axis 417 axialPt = Calc.subtract(center, centerOnAxis); 418 419 } else { 420 axialPt = rotationPos; 421 } 422 423 // Find end points of the rotation axis to display 424 Atom axisMin = (Atom) axialPt.clone(); 425 Calc.scaleAdd(min, rotationAxis, axisMin); 426 Atom axisMax = (Atom) axialPt.clone(); 427 Calc.scaleAdd(max, rotationAxis, axisMax); 428 429 return new Pair<>(axisMin, axisMax); 430 } 431 /** 432 * Returns a Jmol script which will display the axis of rotation. This 433 * consists of a cyan arrow along the axis, plus an arc showing the angle 434 * of rotation. 435 * <p> 436 * As the rotation angle gets smaller, the axis of rotation becomes poorly 437 * defined and would need to get farther and farther away from the protein. 438 * This is not particularly useful, so we arbitrarily draw it parallel to 439 * the translation and omit the arc. 440 * @param atoms Some atoms from the protein, used for determining the bounds 441 * of the axis. 442 * @param axisID in case of representing more than one axis in the same jmol 443 * panel, indicate the ID number. 444 * 445 * @return The Jmol script, suitable for calls to 446 * {@link org.biojava.nbio.structure.align.gui.jmol.StructureAlignmentJmol#evalString() jmol.evalString()} 447 */ 448 public String getJmolScript(Atom[] atoms, int axisID){ 449 final double width=.5;// width of JMol object 450 final String axisColor = "yellow"; //axis color 451 final String screwColor = "orange"; //screw translation color 452 453 Pair<Atom> endPoints = getAxisEnds(atoms); 454 Atom axisMin = endPoints.getFirst(); 455 Atom axisMax = endPoints.getSecond(); 456 457 StringWriter result = new StringWriter(); 458 459 // set arrow heads to a reasonable length 460 result.append("set defaultDrawArrowScale 2.0;"); 461 462 // draw axis of rotation 463 result.append( 464 String.format(Locale.US, "draw ID rot"+axisID+" CYLINDER {%f,%f,%f} {%f,%f,%f} WIDTH %f COLOR %s ;", 465 axisMin.getX(),axisMin.getY(),axisMin.getZ(), 466 axisMax.getX(),axisMax.getY(),axisMax.getZ(), width, axisColor )); 467 468 // draw screw component 469 boolean positiveScrew = Math.signum(rotationAxis.getX()) == Math.signum(screwTranslation.getX()); 470 if( positiveScrew ) { 471 // screw is in the same direction as the axis 472 result.append( String.format(Locale.US, 473 "draw ID screw"+axisID+" VECTOR {%f,%f,%f} {%f,%f,%f} WIDTH %f COLOR %s ;", 474 axisMax.getX(),axisMax.getY(),axisMax.getZ(), 475 screwTranslation.getX(),screwTranslation.getY(),screwTranslation.getZ(), 476 width, screwColor )); 477 } else { 478 // screw is in the opposite direction as the axis 479 result.append( String.format(Locale.US, 480 "draw ID screw"+axisID+" VECTOR {%f,%f,%f} {%f,%f,%f} WIDTH %f COLOR %s ;", 481 axisMin.getX(),axisMin.getY(),axisMin.getZ(), 482 screwTranslation.getX(),screwTranslation.getY(),screwTranslation.getZ(), 483 width, screwColor )); 484 } 485 486 // draw angle of rotation 487 if(rotationPos != null) { 488 result.append(System.getProperty("line.separator")); 489 result.append(String.format(Locale.US, "draw ID rotArc"+axisID+" ARC {%f,%f,%f} {%f,%f,%f} {0,0,0} {0,%f,%d} SCALE 500 DIAMETER %f COLOR %s;", 490 axisMin.getX(),axisMin.getY(),axisMin.getZ(), 491 axisMax.getX(),axisMax.getY(),axisMax.getZ(), 492 Math.toDegrees(theta), 493 positiveScrew ? 0 : 1 , // draw at the opposite end from the screw arrow 494 width, axisColor )); 495 } 496 497 return result.toString(); 498 } 499 500 /** 501 * Projects a given point onto the axis of rotation 502 * @param point 503 * @return An atom which lies on the axis, or null if the RotationAxis is purely translational 504 */ 505 public Atom getProjectedPoint(Atom point) { 506 if(rotationPos == null) { 507 // translation only 508 return null; 509 } 510 511 Atom localPoint = Calc.subtract(point, rotationPos); 512 double dot = Calc.scalarProduct(localPoint, rotationAxis); 513 514 Atom localProjected = Calc.scale(rotationAxis, dot); 515 Atom projected = Calc.add(localProjected, rotationPos); 516 return projected; 517 } 518 519 /** 520 * Get the distance from a point to the axis of rotation 521 * @param point 522 * @return The distance to the axis, or NaN if the RotationAxis is purely translational 523 */ 524 public double getProjectedDistance(Atom point) { 525 Atom projected = getProjectedPoint(point); 526 if( projected == null) { 527 // translation only 528 return Double.NaN; 529 } 530 531 532 return Calc.getDistance(point, projected); 533 534 } 535 536 public void rotate(Atom[] atoms, double theta) { 537 Matrix rot = getRotationMatrix(theta); 538 if(rotationPos == null) { 539 // Undefined rotation axis; do nothing 540 return; 541 } 542 Atom negPos; 543 544 negPos = Calc.invert(rotationPos); 545 546 for(Atom a: atoms) { 547 Calc.shift(a, negPos); 548 } 549 Calc.rotate(atoms, rot); 550 for(Atom a: atoms) { 551 Calc.shift(a, rotationPos); 552 } 553 } 554 555 /** 556 * Calculate the rotation angle for a structure 557 * @param afpChain 558 * @return The rotation angle, in radians 559 * @throws StructureException If the alignment doesn't contain any blocks 560 * @throws NullPointerException If the alignment doesn't have a rotation matrix set 561 */ 562 public static double getAngle(AFPChain afpChain) throws StructureException { 563 if(afpChain.getBlockNum() < 1) { 564 throw new StructureException("No aligned residues"); 565 } 566 Matrix rotation = afpChain.getBlockRotationMatrix()[0]; 567 568 if(rotation == null) { 569 throw new NullPointerException("AFPChain does not contain a rotation matrix"); 570 } 571 return getAngle(rotation); 572 } 573 /** 574 * Calculate the rotation angle for a given matrix 575 * @param rotation Rotation matrix 576 * @return The angle, in radians 577 */ 578 public static double getAngle(Matrix rotation) { 579 double c = (rotation.trace()-1)/2.0; //=cos(theta) 580 // c is sometimes slightly out of the [-1,1] range due to numerical instabilities 581 if( -1-1e-8 < c && c < -1 ) c = -1; 582 if( 1+1e-8 > c && c > 1 ) c = 1; 583 if( -1 > c || c > 1 ) { 584 throw new IllegalArgumentException("Input matrix is not a valid rotation matrix."); 585 } 586 return Math.acos(c); 587 } 588 589 /** 590 * 591 * @return If the rotation axis is well defined, rather than purely translational 592 */ 593 public boolean isDefined() { 594 return rotationPos != null; 595 } 596 597 /** 598 * Quickly compute the rotation angle from a rotation matrix. 599 * @param transform 4D transformation matrix. Translation components are ignored. 600 * @return Angle, from 0 to PI 601 */ 602 public static double getAngle(Matrix4d transform) { 603 // Calculate angle 604 double c = (transform.m00 + transform.m11 + transform.m22 - 1)/2.0; //=cos(theta) 605 // c is sometimes slightly out of the [-1,1] range due to numerical instabilities 606 if( -1-1e-8 < c && c < -1 ) c = -1; 607 if( 1+1e-8 > c && c > 1 ) c = 1; 608 if( -1 > c || c > 1 ) { 609 throw new IllegalArgumentException("Input matrix is not a valid rotation matrix."); 610 } 611 return Math.acos(c); 612 } 613 614 /** 615 * Quickly compute the rotation angle from a rotation matrix. 616 * @param transform 3D rotation matrix 617 * @return Angle, from 0 to PI 618 */ 619 public static double getAngle(Matrix3d transform) { 620 // Calculate angle 621 double c = (transform.m00 + transform.m11 + transform.m22 - 1)/2.0; //=cos(theta) 622 // c is sometimes slightly out of the [-1,1] range due to numerical instabilities 623 if( -1-1e-8 < c && c < -1 ) c = -1; 624 if( 1+1e-8 > c && c > 1 ) c = 1; 625 if( -1 > c || c > 1 ) { 626 throw new IllegalArgumentException("Input matrix is not a valid rotation matrix."); 627 } 628 return Math.acos(c); 629 } 630}