Clipper
clipper_types.h
1 
4 //C Copyright (C) 2000-2006 Kevin Cowtan and University of York
5 //L
6 //L This library is free software and is distributed under the terms
7 //L and conditions of version 2.1 of the GNU Lesser General Public
8 //L Licence (LGPL) with the following additional clause:
9 //L
10 //L `You may also combine or link a "work that uses the Library" to
11 //L produce a work containing portions of the Library, and distribute
12 //L that work under terms of your choice, provided that you give
13 //L prominent notice with each copy of the work that the specified
14 //L version of the Library is used in it, and that you include or
15 //L provide public access to the complete corresponding
16 //L machine-readable source code for the Library including whatever
17 //L changes were used in the work. (i.e. If you make changes to the
18 //L Library you must distribute those, but you do not need to
19 //L distribute source or object code to those portions of the work
20 //L not covered by this licence.)'
21 //L
22 //L Note that this clause grants an additional right and does not impose
23 //L any additional restriction, and so does not affect compatibility
24 //L with the GNU General Public Licence (GPL). If you wish to negotiate
25 //L other terms, please contact the maintainer.
26 //L
27 //L You can redistribute it and/or modify the library under the terms of
28 //L the GNU Lesser General Public License as published by the Free Software
29 //L Foundation; either version 2.1 of the License, or (at your option) any
30 //L later version.
31 //L
32 //L This library is distributed in the hope that it will be useful, but
33 //L WITHOUT ANY WARRANTY; without even the implied warranty of
34 //L MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35 //L Lesser General Public License for more details.
36 //L
37 //L You should have received a copy of the CCP4 licence and/or GNU
38 //L Lesser General Public License along with this library; if not, write
39 //L to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
40 //L The GNU Lesser General Public can also be obtained by writing to the
41 //L Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
42 //L MA 02111-1307 USA
43 
44 
45 #ifndef CLIPPER_TYPES
46 #define CLIPPER_TYPES
47 
48 
49 #include "clipper_util.h"
50 #include "clipper_memory.h"
51 
52 
53 namespace clipper
54 {
55  // forward definitions
56  template<class T> class Mat33sym;
57 
58 
60 
64  class String : public std::string
65  {
66  public:
67  inline String() : std::string() {}
68  inline String( const std::string str ) : std::string( str ) {}
69  inline String( const char* str ) : std::string( str ) {}
70  String( const char* str, const int l );
71  String( const int i, const int w = 4 );
72  String( const long i, const int w = 4 );
73  String( const float f, const int w = 6, const int p = 6 );
76  String( const double f, const int w = 6, const int p = 6 );
77 
79  std::vector<String> split(const String sep) const;
81  String trim() const;
82 
84  String tail() const;
86  String head() const;
88  String nohead() const;
90  String notail() const;
91 
93  static String rational( const double f, const int b, const bool sign=false );
94 
95  int i() const;
96  long l() const;
97  ftype32 f32() const;
98  ftype64 f64() const;
99  ftype f() const;
100  ftype rational() const;
101  };
102 
103 
105  template<class T = ftype> class Vec3
106  {
107  public:
109  inline Vec3() {}
111  inline Vec3( const T& v0, const T& v1, const T& v2 )
112  { vec[0] = v0; vec[1] = v1; vec[2] = v2; }
114  template<class TT> explicit Vec3( const Vec3<TT>& v )
115  { vec[0] = TT(v[0]); vec[1] = TT(v[1]); vec[2] = TT(v[2]); }
117  bool equals( const Vec3<T>& v, const T& tol ) const;
119  inline const T& operator []( const int& i ) const { return vec[i]; }
121  inline T& operator []( const int& i ) { return vec[i]; }
123  inline Vec3<T> unit() const
124  { return (*this)*T(1.0/sqrt(ftype(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]))); }
126  inline static Vec3<T> zero() { return Vec3<T>( 0, 0, 0 ); }
128  inline static Vec3<T> null() { return Vec3<T>( T(Util::nan()), 0, 0 ); }
130  inline bool is_null() const { return Util::is_nan( vec[0] ); }
132  inline static T dot( const Vec3<T>& v1, const Vec3<T>& v2 )
133  { return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]); }
135  inline static Vec3<T> cross( const Vec3<T>& v1, const Vec3<T>& v2 )
136  { return Vec3<T>(v1[1]*v2[2]-v2[1]*v1[2],
137  v1[2]*v2[0]-v2[2]*v1[0],
138  v1[0]*v2[1]-v2[0]*v1[1]); }
140  String format() const
141  { return "("+String(vec[0],10,4)+","+String(vec[1],10,4)+","+String(vec[2],10,4)+")"; }
143  inline const Vec3<T>& operator += ( const Vec3<T>& v )
144  { vec[0] += v[0]; vec[1] += v[1]; vec[2] += v[2]; return (*this); }
146  inline const Vec3<T>& operator -= ( const Vec3<T>& v )
147  { vec[0] -= v[0]; vec[1] -= v[1]; vec[2] -= v[2]; return (*this); }
149  //-- friend int operator == ( const Vec3<T>& v1, const Vec3<T>& v2 );
151  //-- friend int operator != ( const Vec3<T>& v1, const Vec3<T>& v2 );
153  //-- friend Vec3<T> operator -( const Vec3<T>& v );
155  //-- friend Vec3<T> operator +( const Vec3<T>& v1, const Vec3<T> &v2 );
157  //-- friend Vec3<T> operator -( const Vec3<T>& v1, const Vec3<T>& v2 );
159  //-- friend Vec3<T> operator *( const T& s, const Vec3<T>& v1 );
161  //-- friend Vec3<T> operator *( const Vec3<T>& v1, const T& s );
163  //-- friend T operator *( const Vec3<T>& v1, const Vec3<T>& v2 );
164  private:
165  T vec[3];
166  };
167  template<class T> inline int operator == ( const Vec3<T>& v1, const Vec3<T>& v2 ) { return (v1[0]==v2[0] && v1[1]==v2[1] && v1[2]==v2[2]); }
168  template<class T> inline int operator != ( const Vec3<T>& v1, const Vec3<T>& v2 ) { return (v1[0]!=v2[0] || v1[1]!=v2[1] || v1[2]!=v2[2]); }
169  template<class T> inline Vec3<T> operator -( const Vec3<T>& v )
170  { return Vec3<T>( -v[0], -v[1], -v[2] ); }
171  template<class T> inline Vec3<T> operator +( const Vec3<T>& v1, const Vec3<T> &v2 ) { return Vec3<T>( v1[0]+v2[0], v1[1]+v2[1], v1[2]+v2[2]); }
172  template<class T> inline Vec3<T> operator -( const Vec3<T>& v1, const Vec3<T>& v2 ) { return Vec3<T>( v1[0]-v2[0], v1[1]-v2[1], v1[2]-v2[2]); }
173  template<class T> inline Vec3<T> operator *( const T& s, const Vec3<T>& v1 )
174  { return Vec3<T>( s*v1[0], s*v1[1], s*v1[2]); }
175  template<class T> inline Vec3<T> operator *( const Vec3<T>& v1, const T& s )
176  { return Vec3<T>( s*v1[0], s*v1[1], s*v1[2]); }
177  template<class T> inline T operator *( const Vec3<T>& v1, const Vec3<T>& v2 )
178  { return Vec3<T>::dot(v1,v2); }
179 
180 
182  template<class T = ftype> class Mat33
183  {
184  public:
186  inline Mat33() {}
188  inline Mat33( const T& m00, const T& m01, const T& m02, const T& m10, const T& m11, const T& m12, const T& m20, const T& m21, const T& m22 )
189  { mat[0][0] = m00; mat[0][1] = m01; mat[0][2] = m02; mat[1][0] = m10; mat[1][1] = m11; mat[1][2] = m12; mat[2][0] = m20; mat[2][1] = m21; mat[2][2] = m22; }
191  template<class TT> explicit Mat33( const Mat33<TT>& m );
193  template<class TT> explicit Mat33( const Mat33sym<TT>& m );
194  T det() const;
195  Mat33<T> inverse() const;
196  Mat33<T> transpose() const;
197  bool equals( const Mat33<T>& m, const T& tol ) const;
198  inline const T& operator ()( const int& i, const int& j ) const
199  { return mat[i][j]; }
200  inline T& operator ()( const int& i, const int& j )
201  { return mat[i][j]; }
202  String format() const
204  { return "|"+String(mat[0][0],10,4)+","+String(mat[0][1],10,4)+","+String(mat[0][2],10,4)+"|\n|"+String(mat[1][0],10,4)+","+String(mat[1][1],10,4)+","+String(mat[1][2],10,4)+"|\n|"+String(mat[2][0],10,4)+","+String(mat[2][1],10,4)+","+String(mat[2][2],10,4)+"|"; }
206  inline static Mat33<T> identity() { Mat33<T> m; m.mat[0][0] = m.mat[1][1] = m.mat[2][2] = 1; m.mat[0][1] = m.mat[0][2] = m.mat[1][0] = m.mat[1][2] = m.mat[2][0] = m.mat[2][1] = 0; return m; }
208  inline static Mat33<T> null() { Mat33<T> m; m.mat[0][0] = Util::nan(); return m; }
210  bool inline is_null() const { return Util::is_nan( mat[0][0] ); }
212 
213  //-- friend Vec3<T> operator *( const Mat33<T>& m, const Vec3<T>& v );
215 
217  //-- friend Vec3<T> operator *( const Vec3<T>& v, const Mat33<T>& m );
219  //-- friend Mat33<T> operator *(const Mat33<T>& m1, const Mat33<T>& m2);
221  //-- friend Mat33<T> operator +(const Mat33<T>& m1, const Mat33<T>& m2);
223  //-- friend Mat33<T> operator -(const Mat33<T>& m);
224  private:
225  T mat[3][3];
226  };
227  template<class T> inline Vec3<T> operator *( const Mat33<T>& m, const Vec3<T>& v )
228  { return Vec3<T>( m(0,0)*v[0]+m(0,1)*v[1]+m(0,2)*v[2],
229  m(1,0)*v[0]+m(1,1)*v[1]+m(1,2)*v[2],
230  m(2,0)*v[0]+m(2,1)*v[1]+m(2,2)*v[2] ); }
231  template<class T> inline Vec3<T> operator *( const Vec3<T>& v, const Mat33<T>& m )
232  { return Vec3<T>( v[0]*m(0,0)+v[1]*m(1,0)+v[2]*m(2,0),
233  v[0]*m(0,1)+v[1]*m(1,1)+v[2]*m(2,1),
234  v[0]*m(0,2)+v[1]*m(1,2)+v[2]*m(2,2) ); }
235  template<class T> inline Mat33<T> operator *( const Mat33<T>& m1, const Mat33<T>& m2 )
236  {
237  return Mat33<T> ( m1(0,0)*m2(0,0) + m1(0,1)*m2(1,0) + m1(0,2)*m2(2,0),
238  m1(0,0)*m2(0,1) + m1(0,1)*m2(1,1) + m1(0,2)*m2(2,1),
239  m1(0,0)*m2(0,2) + m1(0,1)*m2(1,2) + m1(0,2)*m2(2,2),
240  m1(1,0)*m2(0,0) + m1(1,1)*m2(1,0) + m1(1,2)*m2(2,0),
241  m1(1,0)*m2(0,1) + m1(1,1)*m2(1,1) + m1(1,2)*m2(2,1),
242  m1(1,0)*m2(0,2) + m1(1,1)*m2(1,2) + m1(1,2)*m2(2,2),
243  m1(2,0)*m2(0,0) + m1(2,1)*m2(1,0) + m1(2,2)*m2(2,0),
244  m1(2,0)*m2(0,1) + m1(2,1)*m2(1,1) + m1(2,2)*m2(2,1),
245  m1(2,0)*m2(0,2) + m1(2,1)*m2(1,2) + m1(2,2)*m2(2,2) );
246  }
247  template<class T> inline Mat33<T> operator +( const Mat33<T>& m1, const Mat33<T>& m2 )
248  { return Mat33<T>( m1(0,0)+m2(0,0), m1(0,1)+m2(0,1), m1(0,2)+m2(0,2),
249  m1(1,0)+m2(1,0), m1(1,1)+m2(1,1), m1(1,2)+m2(1,2),
250  m1(2,0)+m2(2,0), m1(2,1)+m2(2,1), m1(2,2)+m2(2,2) ); }
251  template<class T> inline Mat33<T> operator -( const Mat33<T>& m )
252  { return Mat33<T>( -m(0,0), -m(0,1), -m(0,2),
253  -m(1,0), -m(1,1), -m(1,2),
254  -m(2,0), -m(2,1), -m(2,2) ); }
255 
256 
258  template<class T = ftype> class Mat33sym
259  {
260  public:
262  inline Mat33sym() {}
264  template<class TT> explicit Mat33sym( const Mat33<TT>& m ) :
265  m00(m(0,0)), m11(m(1,1)), m22(m(2,2)),
266  m01(m(0,1)), m02(m(0,2)), m12(m(1,2)) {}
268  template<class TT> explicit Mat33sym( const Mat33sym<TT>& m ) :
269  m00(m.mat00()), m11(m.mat11()), m22(m.mat22()),
270  m01(m.mat01()), m02(m.mat02()), m12(m.mat12()) {}
272  inline Mat33sym( const T& c00, const T& c11, const T& c22,
273  const T& c01, const T& c02, const T& c12 ) :
274  m00(c00), m11(c11), m22(c22), m01(c01), m02(c02), m12(c12) {}
276  String format() const
277  { return "|"+String(m00,10,4)+","+String(m01,10,4)+","+String(m02,10,4)+"|\n|"+String(m01,10,4)+","+String(m11,10,4)+","+String(m12,10,4)+"|\n|"+String(m02,10,4)+","+String(m12,10,4)+","+String(m22,10,4)+"|"; }
279  inline static Mat33sym<T> identity()
280  { return Mat33sym<T>( 1, 1, 1, 0, 0, 0 ); }
282  inline static Mat33sym<T> null()
283  { return Mat33sym<T>(Util::nan(),0,0,0,0,0); }
285  inline bool is_null() const { return Util::is_nan( m00 ); }
287  T quad_form( const Vec3<T>& v ) const;
288  T det() const;
289  Mat33<T> sqrt() const;
290  Mat33sym<T> inverse() const;
291  inline const T& mat00() const { return m00; }
292  inline const T& mat11() const { return m11; }
293  inline const T& mat22() const { return m22; }
294  inline const T& mat01() const { return m01; }
295  inline const T& mat02() const { return m02; }
296  inline const T& mat12() const { return m12; }
297  const T& operator ()( const int& i, const int& j ) const;
300  //-- friend Vec3<T> operator *( const Mat33sym<T>& m, const Vec3<T>& v );
302  //-- friend Mat33sym<T> operator +( const Mat33sym<T>& m1, const Mat33sym<T>& m2 );
304  //-- friend Mat33sym<T> operator -( const Mat33sym<T>& m );
305  private:
306  T m00, m11, m22, m01, m02, m12;
307  };
308  template<class T> inline Vec3<T> operator *( const Mat33sym<T>& m, const Vec3<T>& v )
309  { return Vec3<T>( m.mat00()*v[0]+m.mat01()*v[1]+m.mat02()*v[2],
310  m.mat01()*v[0]+m.mat11()*v[1]+m.mat12()*v[2],
311  m.mat02()*v[0]+m.mat12()*v[1]+m.mat22()*v[2] ); }
312  template<class T> inline Mat33sym<T> operator +( const Mat33sym<T>& m1, const Mat33sym<T>& m2 )
313  { return Mat33sym<T>( m1.mat00()+m2.mat00(), m1.mat11()+m2.mat11(),
314  m1.mat22()+m2.mat22(), m1.mat01()+m2.mat01(),
315  m1.mat02()+m2.mat02(), m1.mat12()+m2.mat12() ); }
316  template<class T> inline Mat33sym<T> operator -( const Mat33sym<T>& m )
317  { return Mat33sym<T>( -m.mat00(), -m.mat11(), -m.mat22(),
318  -m.mat01(), -m.mat02(), -m.mat12() ); }
319 
320 
322  template<class T = ftype> class RTop
323  {
324  public:
326  inline RTop() {}
328  inline explicit RTop( const Mat33<T>& r ) : rot_( r ), trn_( Vec3<T>::zero() ) {}
330  inline RTop( const Mat33<T>& r, const Vec3<T>& t ) : rot_( r ), trn_( t ) {}
332  RTop<T> inverse() const
333  { Mat33<T> minv = rot().inverse(); return RTop<T>(minv, -(minv*trn())); }
335  inline bool equals( const RTop<T>& m, const T& tol ) const
336  { return ( rot().equals(m.rot(),tol) && trn().equals(m.trn(),tol) ); }
337  inline const Mat33<T>& rot() const { return rot_; }
338  inline const Vec3<T>& trn() const { return trn_; }
339  inline Mat33<T>& rot() { return rot_; }
340  inline Vec3<T>& trn() { return trn_; }
341  inline static RTop<T> identity()
343  { return RTop<T>( Mat33<T>::identity(), Vec3<T>::zero() ); }
345  inline static RTop<T> null()
346  { return RTop<T>( Mat33<T>::null(), Vec3<T>::null() ); }
348  inline bool is_null() const { return rot_.is_null() || trn_.is_null(); }
350  String format() const { return rot_.format() + "\n" + trn_.format(); }
352  //-- friend Vec3<T> operator *( const RTop<T>& r, const Vec3<T>& v );
354  //-- friend RTop<T> operator *( const RTop<T>& r1, const RTop<T>& r2 );
355  private:
356  Mat33<T> rot_; Vec3<T> trn_;
357  };
358  template<class T> inline Vec3<T> operator *( const RTop<T>& r, const Vec3<T>& v ) { return r.rot()*v + r.trn(); }
359  template<class T> inline RTop<T> operator *( const RTop<T>& r1, const RTop<T>& r2 ) { return RTop<T>( r1.rot()*r2.rot(), r1.rot()*r2.trn()+r1.trn() ); }
360 
361 
362 
364  template<class T = ftype> class Array2d
365  {
366  public:
368  inline Array2d() { d1_ = d2_ = 0; }
370  inline Array2d( const int& d1, const int& d2 ) { resize( d1, d2 ); }
372  inline Array2d( const int& d1, const int& d2, T val )
373  { resize( d1, d2, val ); }
375  void inline resize( const int& d1, const int& d2 )
376  { d1_ = d1; d2_ = d2; data.resize( size() ); }
378  void inline resize( const int& d1, const int& d2, const T& val )
379  { d1_ = d1; d2_ = d2; data.resize( size(), val ); }
380  inline int size() const { return d1_ * d2_; }
381  inline const int& rows() const { return d1_; }
382  inline const int& cols() const { return d2_; }
383  inline const T& operator () ( const int& i1, const int& i2 ) const
385  { return data[ i1*d2_ + i2 ]; }
387  inline T& operator () ( const int& i1, const int& i2 )
388  { return data[ i1*d2_ + i2 ]; }
389  protected:
390  std::vector<T> data;
391  int d1_, d2_;
392  };
393 
394 
396  template<class T = ftype> class Matrix : public Array2d<T>
397  {
398  public:
400  inline Matrix() {}
402  inline Matrix( const int& d1, const int& d2 )
403  { Array2d<T>::resize( d1, d2 ); }
405  inline Matrix( const int& d1, const int& d2, T val )
406  { Array2d<T>::resize( d1, d2, val ); }
408  std::vector<T> solve( const std::vector<T>& b ) const;
410  std::vector<T> eigen( const bool sort = true );
412 
413  //-- friend std::vector<T> operator *( const Matrix<T>& m, const std::vector<T>& v );
414  };
415  template<class T> std::vector<T> operator *( const Matrix<T>& m, const std::vector<T>& v )
416  {
417  if ( m.cols() != v.size() )
418  Message::message( Message_fatal( "Matrix*vector dimension mismatch" ) );
419  std::vector<T> r( m.rows() );
420  int i, j; T s;
421  for ( i = 0; i < m.rows(); i++ ) {
422  s = T(0);
423  for ( j = 0; j < m.cols(); j++ ) s += m(i,j) * v[j];
424  r[i] = s;
425  }
426  return r;
427  }
428 
429 
430 
431  // template implementations
432 
433  template<class T> bool Vec3<T>::equals( const Vec3<T>& v, const T& tol ) const
434  {
435  return ( pow(vec[0]-v[0],T(2)) + pow(vec[1]-v[1],T(2)) +
436  pow(vec[2]-v[2],T(2)) <= pow(tol,T(2)) );
437  }
438 
439  template<class T> template<class TT> Mat33<T>::Mat33( const Mat33<TT>& m )
440  {
441  mat[0][0]=T(m(0,0)); mat[0][1]=T(m(0,1)); mat[0][2]=T(m(0,2));
442  mat[1][0]=T(m(1,0)); mat[1][1]=T(m(1,1)); mat[1][2]=T(m(1,2));
443  mat[2][0]=T(m(2,0)); mat[2][1]=T(m(2,1)); mat[2][2]=T(m(2,2));
444  }
445 
446  template<class T> template<class TT> Mat33<T>::Mat33( const Mat33sym<TT>& m )
447  {
448  mat[0][0]=T(m.mat00());
449  mat[1][1]=T(m.mat11());
450  mat[2][2]=T(m.mat22());
451  mat[0][1]=mat[1][0]=T(m.mat01());
452  mat[0][2]=mat[2][0]=T(m.mat02());
453  mat[1][2]=mat[2][1]=T(m.mat12());
454  }
455 
456  template<class T> bool Mat33<T>::equals( const Mat33<T>& m, const T& tol ) const
457  {
458  return ( pow(mat[0][0]-m(0,0),T(2)) + pow(mat[0][1]-m(0,1),T(2)) +
459  pow(mat[0][2]-m(0,2),T(2)) + pow(mat[1][0]-m(1,0),T(2)) +
460  pow(mat[1][1]-m(1,1),T(2)) + pow(mat[1][2]-m(1,2),T(2)) +
461  pow(mat[2][0]-m(2,0),T(2)) + pow(mat[2][1]-m(2,1),T(2)) +
462  pow(mat[2][2]-m(2,2),T(2)) <= pow(tol,T(2)) );
463  }
464 
465  template<class T> T Mat33<T>::det() const
466  {
467  return ( mat[0][0]*(mat[1][1]*mat[2][2] - mat[1][2]*mat[2][1]) +
468  mat[0][1]*(mat[1][2]*mat[2][0] - mat[1][0]*mat[2][2]) +
469  mat[0][2]*(mat[1][0]*mat[2][1] - mat[1][1]*mat[2][0]) );
470  }
471 
472  template<class T> Mat33<T> Mat33<T>::inverse() const
473  {
474  T d = det();
475  Mat33<T> inv;
476  inv(0,0) = ( mat[1][1]*mat[2][2] - mat[1][2]*mat[2][1] ) / d;
477  inv(1,0) = ( mat[1][2]*mat[2][0] - mat[1][0]*mat[2][2] ) / d;
478  inv(2,0) = ( mat[1][0]*mat[2][1] - mat[1][1]*mat[2][0] ) / d;
479  inv(0,1) = ( mat[2][1]*mat[0][2] - mat[2][2]*mat[0][1] ) / d;
480  inv(1,1) = ( mat[2][2]*mat[0][0] - mat[2][0]*mat[0][2] ) / d;
481  inv(2,1) = ( mat[2][0]*mat[0][1] - mat[2][1]*mat[0][0] ) / d;
482  inv(0,2) = ( mat[0][1]*mat[1][2] - mat[0][2]*mat[1][1] ) / d;
483  inv(1,2) = ( mat[0][2]*mat[1][0] - mat[0][0]*mat[1][2] ) / d;
484  inv(2,2) = ( mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0] ) / d;
485  return inv;
486  }
487 
488  template<class T> Mat33<T> Mat33<T>::transpose() const
489  {
490  Mat33<T> t;
491  t(0,0) = mat[0][0]; t(0,1) = mat[1][0]; t(0,2) = mat[2][0];
492  t(1,0) = mat[0][1]; t(1,1) = mat[1][1]; t(1,2) = mat[2][1];
493  t(2,0) = mat[0][2]; t(2,1) = mat[1][2]; t(2,2) = mat[2][2];
494  return t;
495  }
496 
497  template<class T> T Mat33sym<T>::det() const
498  {
499  return ( m00*(m11*m22 - m12*m12) + m01*(m12*m02 - m01*m22) +
500  m02*(m01*m12 - m11*m02) );
501  }
502 
503  template<class T> Mat33<T> Mat33sym<T>::sqrt() const
504  {
505  Mat33<T> half( Mat33sym<T>( 0.5, 0.5, 0.5, 0.0, 0.0, 0.0 ) );
506  Mat33<T> target( *this );
507  Mat33<T> result( target );
508  result(1,0) = result(2,0) = result(2,1) = 0.0;
509  for ( int i = 0; i < 10; i++ )
510  result = half * ( result.inverse() * target + result );
511  return result;
512  }
513 
514  template<class T> Mat33sym<T> Mat33sym<T>::inverse() const
515  {
516  T d = det();
517  return Mat33sym<T> ( ( m11*m22 - m12*m12 ) / d,
518  ( m22*m00 - m02*m02 ) / d,
519  ( m00*m11 - m01*m01 ) / d,
520  ( m12*m02 - m22*m01 ) / d,
521  ( m01*m12 - m02*m11 ) / d,
522  ( m02*m01 - m00*m12 ) / d );
523  }
524 
525  template<class T> T Mat33sym<T>::quad_form( const Vec3<T>& v ) const
526  {
527  return ( v[0]*( v[0]*m00 + 2*(v[1]*m01+v[2]*m02) ) +
528  v[1]*( v[1]*m11 + 2*(v[2]*m12) ) + v[2]*v[2]*m22 );
529  }
530 
531  template<class T> const T& Mat33sym<T>::operator ()( const int& i, const int& j ) const
532  {
533  switch (i) {
534  case 0:
535  switch (j) {
536  case 0: return m00;
537  case 1: return m01;
538  default: return m02;
539  }
540  case 1:
541  switch (j) {
542  case 0: return m01;
543  case 1: return m11;
544  default: return m12;
545  }
546  default:
547  switch (j) {
548  case 0: return m02;
549  case 1: return m12;
550  default: return m22;
551  }
552  }
553  }
554 
555 
556 // complex matrix methods
557 
558 
561 template<class T> std::vector<T> Matrix<T>::solve( const std::vector<T>& b ) const
562 {
563  if ( Array2d<T>::rows() != Array2d<T>::cols() )
564  Message::message( Message_fatal("Matrix.solve() matrix not square") );
565  if ( Array2d<T>::rows() != b.size() )
566  Message::message( Message_fatal("Matrix.solve() matrix/vector mismatch") );
567  const int n = Array2d<T>::rows();
568 
569  // solve for X by Gaussian elimination
570  T s, pivot;
571  int i, j, k;
572 
573  Matrix<T> a = *this;
574  std::vector<T> x = b;
575  for ( i = 0; i < n; i++ ) {
576  // pick largest pivot
577  j = i;
578  for ( k = i+1; k < n; k++ )
579  if ( fabs(a(k,i)) > fabs(a(j,i)) ) j = k;
580  // swap rows
581  for ( k = 0; k < n; k++ )
582  Util::swap( a(i,k), a(j,k) );
583  Util::swap( x[i], x[j] );
584  // perform elimination
585  pivot = a(i,i);
586  for ( j = 0; j < n; j++ ) {
587  if ( j != i ) {
588  s = a(j,i) / pivot;
589  for ( k = i+1; k < n; k++ ) a(j,k) = a(j,k) - s*a(i,k);
590  x[j] = x[j] - s*x[i];
591  }
592  }
593  }
594  for ( i = 0; i < n; i++ ) x[i] /= a(i,i);
595  return x;
596 }
597 
598 
604 template<class T> std::vector<T> Matrix<T>::eigen( const bool sort )
605 {
606  if ( Array2d<T>::rows() != Array2d<T>::cols() )
607  Message::message( Message_fatal( "Matrix.eigen() matrix not square" ) );
608  const int n = Array2d<T>::rows();
609 
610  int cyc, j, q, p;
611  T spp, spq, t, s, c, theta, tau, h, ap, aq, a_pq;
612 
613  Matrix<T>& mat = *this;
614  Matrix evec( n, n, 0.0 );
615  std::vector<T> eval( n );
616  std::vector<T> b( n );
617  std::vector<T> z( n );
618 
619  // Set evec to identity, eval & b to diagonal, z to 0.
620  for( p = 0; p < n; p++ ) {
621  evec(p,p) = 1.0;
622  eval[p] = b[p] = mat(p,p);
623  }
624 
625  for ( cyc = 1; cyc <= 50; cyc++ ) {
626 
627  // calc sum of diagonal, off-diagonal
628  spp = spq = 0.0;
629  for ( p=0; p<n-1; p++ ) {
630  for ( q=p+1; q<n; q++ )
631  spq += fabs( mat(p,q) );
632  spp += fabs( mat(p,p) );
633  }
634  if ( spq <= 1.0e-12 * spp ) break;
635 
636  // zero z
637  for ( p = 0; p < n; p++ ) z[p] = 0.0;
638 
639  // now try and reduce each off-diagonal element in turn
640  for( p=0; p<n-1; p++ ) {
641  for( q=p+1; q<n; q++ ) {
642  a_pq = mat( p, q );
643  h = eval[q] - eval[p];
644  if ( fabs(a_pq) > 1.0e-12*fabs(h) ) {
645  theta = 0.5*h/a_pq;
646  t = 1.0/(fabs(theta) + sqrt(1.0 + theta*theta));
647  if ( theta < 0.0 ) t = -t;
648  } else {
649  t = a_pq/h;
650  }
651 
652  // calc trig properties
653  c = 1.0/sqrt(1.0+t*t);
654  s = t*c;
655  tau = s/(1.0+c);
656  h = t * a_pq;
657 
658  // update eigenvalues
659  z[p] -= h;
660  z[q] += h;
661  eval[p] -= h;
662  eval[q] += h;
663 
664  // rotate the upper diagonal of the matrix
665  mat( p, q ) = 0.0;
666  for ( j = 0; j < p; j++ ) {
667  ap = mat( j, p );
668  aq = mat( j, q );
669  mat( j, p ) = ap - s * ( aq + ap * tau );
670  mat( j, q ) = aq + s * ( ap - aq * tau );
671  }
672  for ( j = p+1; j < q; j++ ) {
673  ap = mat( p, j );
674  aq = mat( j, q );
675  mat( p, j ) = ap - s * ( aq + ap * tau );
676  mat( j, q ) = aq + s * ( ap - aq * tau );
677  }
678  for ( j = q+1; j < n; j++ ) {
679  ap = mat( p, j );
680  aq = mat( q, j );
681  mat( p, j ) = ap - s * ( aq + ap * tau );
682  mat( q, j ) = aq + s * ( ap - aq * tau );
683  }
684  // apply corresponding rotation to result
685  for ( j = 0; j < n; j++ ) {
686  ap = evec( j, p );
687  aq = evec( j, q );
688  evec( j, p ) = ap - s * ( aq + ap * tau );
689  evec( j, q ) = aq + s * ( ap - aq * tau );
690  }
691  }
692  }
693 
694  for ( p = 0; p < n; p++ ) {
695  b[p] += z[p];
696  eval[p] = b[p];
697  }
698  }
699 
700  // sort the eigenvalues
701  if ( sort ) {
702  for ( p = 0; p < n; p++ ) {
703  j = p; // set j to index of largest remaining eval
704  for ( q = p+1; q < n; q++ )
705  if ( eval[q] < eval[j] ) j = q;
706  Util::swap( eval[p], eval[j] ); // now swap evals, evecs
707  for ( q = 0; q < n; q++ )
708  Util::swap( evec( q, p ), evec( q, j ) );
709  }
710  }
711 
712  // return eigenvalues
713  mat = evec;
714  return eval;
715 }
716 
717 } // namespace clipper
718 
719 #endif
Array2d(const int &d1, const int &d2)
constructor
Definition: clipper_types.h:370
Mat33(const T &m00, const T &m01, const T &m02, const T &m10, const T &m11, const T &m12, const T &m20, const T &m21, const T &m22)
constructor
Definition: clipper_types.h:188
ftype rational() const
convert from rational to ftype
Definition: clipper_types.cpp:145
const T & mat01() const
element (0,1)
Definition: clipper_types.h:294
String trim() const
Return copy of string without leading and trailing blanks.
Definition: clipper_types.cpp:86
const T & mat00() const
element (0,0)
Definition: clipper_types.h:291
Mat33< T > sqrt() const
square root
Definition: clipper_types.h:503
static Vec3< T > zero()
return zero vector
Definition: clipper_types.h:126
Fatal message (level = 9)
Definition: clipper_message.h:129
T det() const
determinant
Definition: clipper_types.h:465
std::vector< T > solve(const std::vector< T > &b) const
equation solver (square matrices only)
Definition: clipper_types.h:561
int size() const
size
Definition: clipper_types.h:380
Simple 2-d array class.
Definition: clipper_types.h:364
static Vec3< T > cross(const Vec3< T > &v1, const Vec3< T > &v2)
Vector cross product.
Definition: clipper_types.h:135
static Mat33< T > identity()
return identity matrix
Definition: clipper_types.h:206
bool is_null() const
test for null operator
Definition: clipper_types.h:348
ftype f() const
convert to ftype
Definition: clipper_types.cpp:142
String(const char *str)
constructor: from char*
Definition: clipper_types.h:69
RTop(const Mat33< T > &r, const Vec3< T > &t)
constructor: from rotation and translation
Definition: clipper_types.h:330
Mat33sym(const Mat33sym< TT > &m)
constructor: from Mat33sym
Definition: clipper_types.h:268
String head() const
remove trailing path element
Definition: clipper_types.cpp:100
static RTop< T > identity()
return identity operator
Definition: clipper_types.h:342
static Vec3< T > null()
return null vector (only valid for floating point types)
Definition: clipper_types.h:128
Compressed form for 3x3 symmetric matrix class.
Definition: clipper_types.h:56
static Mat33sym< T > null()
return null matrix (only valid for floating point types)
Definition: clipper_types.h:282
bool equals(const Mat33< T > &m, const T &tol) const
test equality
Definition: clipper_types.h:456
bool is_null() const
test for null matrix (only valid for floating point types)
Definition: clipper_types.h:285
const T & mat22() const
element (2,2)
Definition: clipper_types.h:293
T quad_form(const Vec3< T > &v) const
return quadratic form with vector
Definition: clipper_types.h:525
Vec3< T > unit() const
return unit vector with same direction as this vector
Definition: clipper_types.h:123
const T & mat11() const
element (1,1)
Definition: clipper_types.h:292
const Vec3< T > & operator-=(const Vec3< T > &v)
subtract another vector from this one
Definition: clipper_types.h:146
const Vec3< T > & trn() const
get translation
Definition: clipper_types.h:338
static const ftype & nan()
fast Util::nan() value
Definition: clipper_util.h:67
void resize(const int &d1, const int &d2, const T &val)
resize
Definition: clipper_types.h:378
Vec3(const T &v0, const T &v1, const T &v2)
constructor: from individual values
Definition: clipper_types.h:111
String nohead() const
get leading path element
Definition: clipper_types.cpp:103
long l() const
convert to long
Definition: clipper_types.cpp:133
String format() const
return formatted String representation
Definition: clipper_types.h:276
String format() const
return formatted String representation
Definition: clipper_types.h:203
void resize(const int &d1, const int &d2)
resize
Definition: clipper_types.h:375
ftype64 ftype
ftype definition for floating point representation
Definition: clipper_precision.h:58
Mat33sym(const T &c00, const T &c11, const T &c22, const T &c01, const T &c02, const T &c12)
constructor: from coefficients
Definition: clipper_types.h:272
const T & operator()(const int &i1, const int &i2) const
read accessor
Definition: clipper_types.h:384
static Mat33sym< T > identity()
return identity matrix
Definition: clipper_types.h:279
Vec3(const Vec3< TT > &v)
constructor: copy/convert
Definition: clipper_types.h:114
static void swap(T &a, T &b)
swap the contents of two objects
Definition: clipper_util.h:148
Mat33< T > inverse() const
inverse
Definition: clipper_types.h:472
Mat33sym< T > inverse() const
inverse
Definition: clipper_types.h:514
static T dot(const Vec3< T > &v1, const Vec3< T > &v2)
Vector dot product (equivalent to *)
Definition: clipper_types.h:132
const Vec3< T > & operator+=(const Vec3< T > &v)
add another vector to this one
Definition: clipper_types.h:143
Matrix()
null constructor
Definition: clipper_types.h:400
bool is_null() const
test for null vector
Definition: clipper_types.h:130
String notail() const
remove leading path element
Definition: clipper_types.cpp:109
RTop< T > inverse() const
inverse
Definition: clipper_types.h:332
std::vector< String > split(const String sep) const
String splitter - a very simple parser component.
Definition: clipper_types.cpp:71
const T & mat02() const
element (0,2)
Definition: clipper_types.h:295
const int & cols() const
number of cols
Definition: clipper_types.h:382
String extension with simple parsing methods.
Definition: clipper_types.h:64
String tail() const
get trailing path element
Definition: clipper_types.cpp:97
ftype64 f64() const
convert to double
Definition: clipper_types.cpp:139
3-vector class
Definition: clipper_types.h:105
static bool is_nan(const ftype32 f)
fast Util::nan() test
Definition: clipper_util.h:82
Mat33()
null constructor
Definition: clipper_types.h:186
Mat33< T > transpose() const
transpose
Definition: clipper_types.h:488
bool is_null() const
test for null matrix (only valid for floating point types)
Definition: clipper_types.h:210
static Mat33< T > null()
return null matrix (only valid for floating point types)
Definition: clipper_types.h:208
Rotation-translation operator.
Definition: clipper_types.h:322
static void message(const T &message)
pass a message
Definition: clipper_message.h:93
const T & operator()(const int &i, const int &j) const
access elements as 3x3 matrix (inefficient)
Definition: clipper_types.h:531
std::vector< T > eigen(const bool sort=true)
eigenvalue calculation (square symmetric matrices only)
Definition: clipper_types.h:604
Matrix(const int &d1, const int &d2)
constructor
Definition: clipper_types.h:402
ftype32 f32() const
convert to float
Definition: clipper_types.cpp:136
T det() const
determinant
Definition: clipper_types.h:497
String()
null constructor
Definition: clipper_types.h:67
Mat33< T > & rot()
set rotation
Definition: clipper_types.h:339
const Mat33< T > & rot() const
get rotation
Definition: clipper_types.h:337
Matrix(const int &d1, const int &d2, T val)
constructor
Definition: clipper_types.h:405
RTop()
null constructor
Definition: clipper_types.h:326
Array2d(const int &d1, const int &d2, T val)
constructor
Definition: clipper_types.h:372
String format() const
return formatted String representation
Definition: clipper_types.h:140
RTop(const Mat33< T > &r)
constructor: from rotation
Definition: clipper_types.h:328
const T & mat12() const
element (1,2)
Definition: clipper_types.h:296
Mat33sym(const Mat33< TT > &m)
constructor: from Mat33 (does not check for symmetry)
Definition: clipper_types.h:264
bool equals(const RTop< T > &m, const T &tol) const
test equality with some tolerance
Definition: clipper_types.h:335
General matrix class: like Array2d but with numerical methods.
Definition: clipper_types.h:396
Array2d()
null constructor
Definition: clipper_types.h:368
Mat33sym()
null constructor
Definition: clipper_types.h:262
static RTop< T > null()
return identity operator
Definition: clipper_types.h:345
const T & operator[](const int &i) const
get element
Definition: clipper_types.h:119
String format() const
return formatted String representation
Definition: clipper_types.h:350
const int & rows() const
number of rows
Definition: clipper_types.h:381
Vec3< T > & trn()
set translation
Definition: clipper_types.h:340
int i() const
convert to int
Definition: clipper_types.cpp:130
String(const std::string str)
constructor: from string
Definition: clipper_types.h:68
bool equals(const Vec3< T > &v, const T &tol) const
test equality
Definition: clipper_types.h:433
Vec3()
null constructor
Definition: clipper_types.h:109
3x3-matrix class
Definition: clipper_types.h:182
const T & operator()(const int &i, const int &j) const
get element
Definition: clipper_types.h:198