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