matvec_utils.hh
00001 // matvec.hh 00002 00003 // Various utility routines for handling 3x3 matrices & vectors 00004 // Some routines using CCtbx routines are in cctbx_utils 00005 // 00006 // By default all are <double> unless there are special reasons 00007 // for eg float 00008 // (mainly for interface to external libraries such as CCP4) 00009 00010 // For various reasons, I am using up to 4 different representations 00011 // of 3x3 matrices (+ sgtbx::rt_mx 00012 // 00013 // 3x3 matrices are 00014 // ( 0,0 0,1 0,2 ) 00015 // ( 1,0 1,1 1,2 ) element i,j 00016 // ( 2,0 2,1 2,2 ) 00017 // 00018 // a) simple arrays double[3][3] [i][j] AMat33 00019 // b) clipper::Mat33<double> (i,j) CMat33 00020 // c) std::vector<double> M33(9) elements in order 00, 01, 02, 10 ... 00021 // M33(k) k = i*3 + j VMat33 00022 // 00023 // Routines here are:- 00024 // 1) converters Set[type] (other_type) 00025 // (not all combinations written yet!) 00026 // 2) tests of identity is_mat33_ident(Matrix) 00027 // 3) format reindex matrix FormatReindex_as_hkl(Matrix) 00028 // 4) vector utilities Modulus(Vec3<T>) 00029 // IntVec(Vec3<double>) 00030 // 5) LSQ general matrix to superimpose vector set y = [U] x 00031 // not necessarily orthonormal (no constraints) 00032 // 00033 00034 #ifndef MATVEC_UTILS 00035 #define MATVEC_UTILS 00036 00037 #include <cstring> 00038 00039 // Clipper 00040 #include <clipper/clipper.h> 00041 using clipper::Vec3; 00042 using clipper::Mat33; 00043 using clipper::RTop; 00044 typedef clipper::Vec3<double> DVect3; 00045 typedef clipper::Mat33<double> DMat33; 00046 00047 #include "csymlib.h" // CCP4 symmetry stuff 00048 00049 #include <assert.h> 00050 #define ASSERT assert 00051 00052 namespace MVutil { 00053 //-------------------------------------------------------------- 00054 // Make clipper::Mat33 from various things 00055 // 2D array 00056 // float -> double 00057 clipper::Mat33<double> SetCMat33(const float m[3][3]); 00058 clipper::Mat33<double> SetCMat33(const float m[9]); 00059 // double 00060 clipper::Mat33<double> SetCMat33(const double m[3][3]); 00061 // vector<double> (9) 00062 clipper::Mat33<double> SetCMat33(const std::vector<double>& m); 00063 // CCP4 symop 00064 clipper::Mat33<double> SetCMat33(CSym::ccp4_symop& op); 00065 // Diagonal matrix 00066 clipper::Mat33<double> SetDiagCMat33(const clipper::Vec3<double>& v); 00067 //-------------------------------------------------------------- 00068 // Convert Mat33 to vector type 00069 std::vector<double> SetVMat33(const clipper::Mat33<double>& R); 00070 std::vector<double> SetVMat33(const float m[3][3]); 00071 //-------------------------------------------------------------- 00072 //-------------------------------------------------------------- 00073 bool is_mat33_ident(const clipper::Mat33<float>& R); 00074 // Returns true if 3x3 matrix R is identity 00075 //-------------------------------------------------------------- 00076 bool is_mat33_ident(const clipper::Mat33<double>& R); 00077 // Returns true if 3x3 matrix R is identity 00078 //-------------------------------------------------------------- 00079 bool is_rtop_ident(const clipper::RTop<double>& R); 00080 // Returns true if 3x3 matrix R is identity 00081 //-------------------------------------------------------------- 00082 //-------------------------------------------------------------- 00083 // Format reciprocal space matrix (eg reindex or symmetry) 00084 // as h,k,l 00085 std::string FormatSymop_as_hkl(const clipper::Symop& op, 00086 const std::string& brackets = "[]"); 00087 std::string FormatReindex_as_hkl(const clipper::RTop<double>& op, 00088 const std::string& brackets = "[]"); 00089 //-------------------------------------------------------------- 00090 template<class T> inline T Modulus( const clipper::Vec3<T>& v) 00091 { return sqrt(v*v); } 00092 //-------------------------------------------------------------- 00093 clipper::Vec3<int> IntVec(const clipper::Vec3<double>& vector); 00094 //-------------------------------------------------------------- 00095 template<class T> int LargestVectorElement(const clipper::Vec3<T>& v) 00096 // Return index 0,1,2 of largest element in vector 00097 { 00098 if (v[0] > v[1]) { 00099 if (v[0] > v[2]) { 00100 return 0; // v[0] biggest 00101 } else { 00102 return 2; // v(2] biggest 00103 } 00104 } else if (v[1] > v[2]) { 00105 return 1; // v[1] biggest 00106 } 00107 return 2; // v(2] biggest 00108 } 00109 //-------------------------------------------------------------- 00110 // Generate LSQ general matrix transforming one set of 3-vectors into another 00111 // if the 3x3 matrix [U] transforms x -> y, with a possible error e 00112 // then e = [U] x - y 00113 // 00114 // Return [U] & mean residual E 00115 // 00116 // Minimising E = 1/2N Sum(e.e) ie least-squares gives 00117 // [U] = [R] [S]^-1 00118 // where 00119 // [S]ij = Sum(xi xj) i, j = 1,3 for 3-vector 00120 // [R]ij = Sum(yi xj) 00121 // On entry: 00122 // x, y corresponding vector lists 00123 // AvoidPlanar if true, generate a few orthogonal vectors to avoid the ambiguities 00124 // introduced in all x vectors are coplanar (not included in residual E) 00125 // but only if needed (only valid for orthogonal frames) 00126 DMat33 LsqTransform(const std::vector<DVect3>& x, const std::vector<DVect3>& y, 00127 double& E, const bool& AvoidPlanar=false); 00128 //-------------------------------------------------------------- 00129 // Get axis direction for rotation matrix R 00130 // Returns normalised axis direction 00131 clipper::Vec3<double> AxisDirection(const clipper::Mat33<double>& R); 00132 //-------------------------------------------------------------- 00133 } 00134 00135 #endif 00136
