#ifndef _TRANSFORM_H_ #define _TRANSFORM_H_ #include "Quaternion.h" #include "Pointlist.h" #include "Superpose.h" #include ////////////////////////////////////////////////////////////////////////// // CLASS PURPOSE: Holds the rotation Quaternion and translation offset // // for the optimal Transformation of one Pointlist to // // another (minimizing the RMS error). The class creates // // a Transform object and provides basic operations on // // that object // ////////////////////////////////////////////////////////////////////////// class Transform { private: Quaternion Rotation; // rotation Quaternion XYZpoint Offset; // translational offset, which is added public: static Transform identity; Transform(){(*this)=identity;}; // void constructor Transform(const Quaternion &Q, const XYZpoint &A) { Rotation = Q; Offset = A; } // constructor for setting Rotation and Offset directly Transform(const double matrix[3][3]); // convert a rotation matrix (such as is used by rotate_by_matrix) // into a transform (pure rotation, no translation) // matrix transforms point (1,0,0) into matrix[0][0], matrix[1][0], matrix[2][0] // Construct optimal transform from sublist of A to // sublist of B (minimizing weighted rmsd). // fix_start_length() is // used to interpret negative values of start and length. // (the size of Pointlist A is used for interpreting length). // The weights have the same subscripting as Pointlist A. Transform(const Pointlist &A, const Pointlist &B, const float *weight=NULL, int start_point_A=0, int length=-1, int start_point_B=0); // Same as Transform with weights, but using 0/1 weight Transform(const Pointlist &A, const Pointlist &B, const bool *present, int start_point_A=0, int length=-1, int start_point_B=0); // Construct a transformation that is a rotation of theta radians // about the axis through points B and C. // That is, if you have non-colinear points A,B,C,D, // applying the transformation to B or C leaves them unchanged, // applying the transformation to D increases the torsion angle(A,B,C,D) // by theta, but leaves distance(C,D) and bond angle(B,C,D) unchanged Transform(const XYZpoint&B, const XYZpoint&C, double theta); inline Quaternion& rotation(void) { return Rotation; } inline const Quaternion& rotation(void) const { return Rotation; } // convert the Rotation Quaternion to a rotation matrix (such as is used by rotate_by_matrix) // matrix transforms point (1,0,0) into matrix[0][0], matrix[1][0], matrix[2][0] void get_rotation_matrix(double matrix[3][3]) const; inline bool is_finite(void) const { return Rotation.is_finite() && Offset.is_finite(); } inline bool OK(void) const { assert(is_finite()); double rot_mag2=Rotation.mag2(); const double ERR_LIMIT = 1.e-4; assert( 1.-ERR_LIMIT < rot_mag2 && rot_mag2 < 1.+ERR_LIMIT); const double MAX_OFFSET2 = 1.e+40; double offset_mag2 = Offset.mag2(); assert (offset_mag2 < MAX_OFFSET2); return 1; } inline XYZpoint& offset(void) { return Offset; } inline const XYZpoint& offset(void) const { return Offset; } inline bool operator == (const Transform& t) const { return Rotation==t.Rotation && Offset==t.Offset; } // Modify transform to be an nmer-th root of unity with a positive // real part of the rotation vector void normalize_nth_root(int nmer); //////////////////////////////////////////////////////////////////////// // The following sets of 3 routines provide 3 versions of the inverse, // times and divide operators. (Also *= and /= operators) //////////////////////////////////////////////////////////////////////// Transform inv(void) const; // function to find inverse of a Transform // Return by value - copy overhead friend Transform& inverse(Transform& destn, const Transform& source); // Most efficient inverse routine, no copy overhead // in return. Requires dest, a and b to be created by // caller routine // Return by reference - no copy overhead Transform * inverse(void) const; // Computes inverse and returns a pointer to a new // tranform created to store it // Return by reference - no copy overhead //////////////////////////////////////////////////////////////////////// Transform operator*(const Transform &Q) const; // operator to multiply 2 Transforms // Return by value - copy overhead // a*b(P) = a*(b(P)) friend Transform& times(Transform &dest, const Transform& a, const Transform& b); // Most efficient product routine, no copy overhead // in return. Requires dest, a and b to be created by // caller routine // Return by reference - no copy overhead // dest may be same as a or b. Transform* times(const Transform* Q) const; // Computes product and returns a pointer to a new // tranform created to store it // Return by reference - no copy overhead Transform& operator *=(const Transform& Q); // computes product in place //////////////////////////////////////////////////////////////////////// Transform operator/(const Transform &Q) const; // operator to divide 2 Transforms i.e. perform a*inv(b) // Return by value - copy overhead // a/b(P) = a*(b.inv(P)) friend Transform& divide(Transform &dest, const Transform& a, const Transform& b); // Most efficient divide routine, no copy overhead // in return. Requires dest, a and b to be created by // caller routine // Return by reference - no copy overhead // dest may be same as a or b Transform* divide(const Transform* Q) const; // Computes quotient and returns a pointer to a new // tranform created to store it // Return by reference - no copy overhead Transform& operator /=(const Transform& Q); // computes (*this)*Q.inv(void) in place //////////////////////////////////////////////////////////////////////// void rotate(const XYZpoint &A, XYZpoint &B) const; // rotates XYZpoint A into location B , using // precomputed Rotation Quaternion void shift(const XYZpoint &A, XYZpoint &B) const; // translates XYZpoint A to location B, using // precomputed translational Offset void trans(const XYZpoint & source, XYZpoint &destn) const; // Transforms XYZpoint source into location destn by first // rotating and then translating, using precomputed // Rotation Quaternion and translational Offset //////////////////////////////////////////////////////////////////////// // For the Pointlist versions of the operators, // if number<0, then last point to Transform is A[A.NumPoints+number], // and we need to have startA<=A.NumPoints+number+1 // void rotate(const Pointlist &A, Pointlist &B, long int startA = 0, long int startB = 0, long int number = -1) const; // rotates ( A[startA]...A[startA+number-1] ) to // location ( B[startB]...B[startB+number-1] ) using // precomputed Rotation Quaternion. void shift(const Pointlist &A, Pointlist &B, long int startA = 0, long int startB = 0, long int number = -1) const; // translates( A[startA]...A[startA+number-1] ) to // location ( B[startB]...B[startB+number-1] ) using // precomputed translational Offset void trans(const Pointlist &A, Pointlist &B, long int startA = 0, long int startB = 0, long int number = -1) const; // Transforms ( A[startA]...A[startA+number-1] ) to // location ( B[startB]...B[startB+number-1] ) by first // rotating and then translating, using precomputed // Rotation Quaternion and translational Offset //////////////////////////////////////////////////////////////////////// void rotateIP(Pointlist &A, long int start = 0, long int number = -1) const; // rotates ( A[start]...A[start+number-1] ) INPLACE, using // precomputed Rotation Quaternion // If number = -1, rotate from Points[startA] to end of // array INPLACE void shiftIP(Pointlist &A, long int start = 0, long int number = -1) const; // translates ( A[start]...A[start+number-1] ) INPLACE, // using precomputed translational Offset // If number = -1, shift from Points[startA] to end of // array INPLACE void transIP(Pointlist &A, long int start = 0, long int number = -1) const ; // Transform ( A[start]...A[start+number-1] ) INPLACE, using // precomputed Rotation Quaternion and translational Offset // If number = -1, Transform from Points[startA] to end of // array INPLACE //////////////////////////////////////////////////////////////////////// // find the optimal rigid Transformation of Pointlist A to B // (with optional weighting), and store it in this void optimal_transform(const Pointlist &A, const Pointlist &B, const float *weight=0, int start_a=0, int length=-1, int start_b=0); // efficient computation of Rotation and Offset // for exactly 3 points void coplanar_trans(const Pointlist &A, const Pointlist &B, int start_a=0, int start_b=0); void xyplane(const XYZpoint& orig, const XYZpoint& x, const XYZpoint& y); // set this to a Transformation that takes orig->(0,0,0) // x->(>0, 0, 0) // y->(*,>0,0) friend ostream &operator<<(ostream &stream, const Transform& T); friend istream &operator>>(istream &stream, Transform &T); // I/O operators }; // INLINE PROCEDURES for greater efficiency inline Transform& inverse(Transform& destn, const Transform& source) { XYZpoint zero(0,0,0); source.Rotation.conjugate(destn.Rotation); destn.Offset = zero - destn.Rotation*source.Offset*source.Rotation; return destn; } inline Transform Transform::inv(void) const { Transform ret; return ::inverse(ret, *this); } inline Transform * Transform::inverse(void) const { Transform * ret = new Transform; ::inverse(*ret, *this); return ret; } inline Transform& times(Transform& destn, const Transform& a, const Transform& b) { // update offsets first in case destn is a or b. XYZpoint tmp; a.rotate(b.Offset, tmp); add(destn.Offset,a.Offset,tmp); destn.Rotation = a.Rotation * b.Rotation; return destn; } inline Transform Transform::operator*(const Transform& Q) const { Transform ret; return ::times(ret, *this, Q); } inline Transform * Transform::times(const Transform * Q) const { Transform * ret = new Transform; ::times(*ret, *this, *Q); return ret; } inline Transform & Transform::operator*=(const Transform& Q) { return ::times(*this, *this, Q); } inline Transform& divide(Transform& destn, const Transform& a, const Transform& b) { Quaternion conj; b.Rotation.conjugate(conj); destn.Rotation = a.Rotation * conj; XYZpoint tmp; destn.rotate(b.Offset, tmp); destn.Offset = a.Offset - tmp; return destn; } inline Transform Transform::operator/(const Transform& Q) const { Transform ret; return ::divide(ret, *this, Q); } inline Transform * Transform::divide(const Transform * Q) const { Transform * ret = new Transform; ::divide(*ret, *this, *Q); return ret; } inline Transform& Transform::operator/=(const Transform& Q) { return ::divide(*this, *this, Q); } inline void Transform::rotate(const XYZpoint &A, XYZpoint &B) const { #ifdef EXPLICIT_ALLOCATE_IN_ROTATE Quaternion temp2,conj,temp3; cross(temp2, Rotation, A); Rotation.conjugate(conj); cross(temp3, temp2, conj); temp3.axis_vector(B); #else ((Rotation*A)*Rotation.conjugate()).axis_vector(B); #endif } inline void Transform::shift(const XYZpoint &A, XYZpoint &B) const { add(B, Offset, A); } inline void Transform::trans(const XYZpoint &A, XYZpoint &B) const { rotate(A,B); B+= Offset; } // matrix transforms point (1,0,0) into matrix[0][0], matrix[1][0], matrix[2][0] inline void rotate_by_matrix(const double matrix[3][3], const XYZpoint & A, XYZpoint & B) { XYZpoint tmp; // use temporary to avoid problems if A and B same. tmp.X = matrix[0][0]*A.X + matrix[0][1]*A.Y + matrix[0][2]*A.Z; tmp.Y = matrix[1][0]*A.X + matrix[1][1]*A.Y + matrix[1][2]*A.Z; tmp.Z = matrix[2][0]*A.X + matrix[2][1]*A.Y + matrix[2][2]*A.Z; B = tmp; } // CHANGE LOG: // 22 Jan 2004 Kevin Karplus // Added constructor for rotation about an axis. // 24 Oct 2004 Kevin Karplus // Added operator == // 24 Jan 2005 Kevin Karplus // Moved normalize_nth_root from undertaker's Multimer.cc class. // Mon Feb 7 09:16:53 PST 2005 Kevin Karplus // Added constructor with bool array marking points that are present. // Tue Jun 13 09:11:52 PDT 2006 Kevin Karplus // Added conversion from rotation matrix to Transform // Tue Jun 13 09:25:21 PDT 2006 Kevin Karplus // Added get_rotation_matrix // Wed Jun 18 07:54:05 PDT 2008 Kevin Karplus // Added is_finite() test #endif