src/openbabel/matrix3x3.hpp

/**********************************************************************
matrix3x3.cpp - Handle 3D Rotation matrix.
 
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
 
This file is part of the Open Babel project.
For more information, see <http://openbabel.sourceforge.net/>
 
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
 
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
***********************************************************************/

#ifndef OB_MATRIX3x3_H
#define OB_MATRIX3x3_H

#include "oberror.hpp"

#if HAVE_IOSTREAM
#include <iostream>
#elif HAVE_IOSTREAM_H
#include <iostream.h>
#endif

#if HAVE_FSTREAM
#include <fstream>
#elif HAVE_FSTREAM_H
#include <fstream.h>
#endif

#include <math.h>

#include "obutil.hpp"
#include "vector3.hpp"

#ifndef PI
#define PI 3.1415926535897932384626433
#endif

#ifndef RAD_TO_DEG
#define RAD_TO_DEG 180.0/PI
#endif

#ifndef DEG_TO_RAD
#define DEG_TO_RAD PI/180.0
#endif

namespace OpenBabel
{

// class introduction in matrix3x3.cpp
class OBAPI matrix3x3
{
    //! Elements of the matrix
    /*! This array holds the matrix. The first index refers to the
    row, the second the column. */
    double ele[3][3];

public:
    //! constructs the zero-matrix
    matrix3x3(void)
    {
        ele[0][0] = 0.0;
        ele[0][1] = 0.0;
        ele[0][2] = 0.0;
        ele[1][0] = 0.0;
        ele[1][1] = 0.0;
        ele[1][2] = 0.0;
        ele[2][0] = 0.0;
        ele[2][1] = 0.0;
        ele[2][2] = 0.0;
    }

    //! constructs s times the unit matrix
    matrix3x3(double s)
    {
        ele[0][0] = s;
        ele[0][1] = 0.0;
        ele[0][2] = 0.0;
        ele[1][0] = 0.0;
        ele[1][1] = s;
        ele[1][2] = 0.0;
        ele[2][0] = 0.0;
        ele[2][1] = 0.0;
        ele[2][2] = s;
    }

    //! constructs a matrix from row vectors
    matrix3x3(vector3 row1,vector3 row2,vector3 row3)
    {
        ele[0][0] = row1.x();
        ele[0][1] = row1.y();
        ele[0][2] = row1.z();
        ele[1][0] = row2.x();
        ele[1][1] = row2.y();
        ele[1][2] = row2.z();
        ele[2][0] = row3.x();
        ele[2][1] = row3.y();
        ele[2][2] = row3.z();
    }

    //! constructs a matrix from a 3x3-array of doubles
    /*! constructs a matrix from a 3x3-array of doubles. The first
    index represents the row, the second index the column */
    matrix3x3(double d[3][3])
    {
        ele[0][0] = d[0][0];
        ele[0][1] = d[0][1];
        ele[0][2] = d[0][2];
        ele[1][0] = d[1][0];
        ele[1][1] = d[1][1];
        ele[1][2] = d[1][2];
        ele[2][0] = d[2][0];
        ele[2][1] = d[2][1];
        ele[2][2] = d[2][2];
    }

    //! access function
    /*! writes the matrix into the 1-dimensional array m, row by
    row. The array must be able to hold 9 doubles, otherwise your
    prgram will segfault. */
    void GetArray(double *m)
    {
        m[0] = ele[0][0];
        m[1] = ele[0][1];
        m[2] = ele[0][2];
        m[3] = ele[1][0];
        m[4] = ele[1][1];
        m[5] = ele[1][2];
        m[6] = ele[2][0];
        m[7] = ele[2][1];
        m[8] = ele[2][2];
    }

    //! Calculates the inverse of a matrix.
    matrix3x3 inverse(void) const throw(OBError);

    //! Calculates the transpose of a matrix.
    matrix3x3 transpose(void) const;

    //! generates a matrix for a random rotation
    void randomRotation(OBRandom &rnd);

    //! returns the determinant of the matrix
    double determinant() const;

    //! Checks if a matrix is symmetric
    bool isSymmetric(void) const;

    //! Checks if a matrix is orthogonal
    /*! This method checks if a matrix describes an orthogonal
    transformation, i.e. if all column vectors are normalized and
    are mutually orthogonal. An orthogonal transformation is a
    transformation the preserves length and angle. 

    The check is performed using the method isUnitMatrix() to
    check if
    \code
    *this * transpose()
    \endcode
    is a unit matrix. The criterion is therefore numerically quite
    tight. */
    bool isOrthogonal(void) const
    {
        return (*this * transpose()).isUnitMatrix();
    };

    //! Checks if a matrix is diagonal
    bool isDiagonal(void) const;

    //! Checks if a matrix is the unit matrix
    bool isUnitMatrix(void) const;

    //! access function
    /*! \warning row or column are not in the range 0..2, random
    results are returned, and your program may even
    segfault. (Stefan Kebekus)

    \todo Replace this method with a more fool-proof version.
    */
    double Get(int row,int column) const
    {
      if (row >= 0 && row <= 2 && column >= 0 && column <= 2)
        return(ele[row][column]);
      else
        return 0.0f;
    }

    //! access function
    /*! \warning if row or column are not in the range 0..2, random
    variables are overwritten, and your program may
    segfault. (Stefan Kebekus)

    \todo Replace this method with a more fool-proof version.
    */
    void Set(int row,int column, double v)
    {
      if (row >= 0 && row <= 2 && column >= 0 && column <= 2)
        ele[row][column]= v;
    }

    //! access function
    /*! \warning If column is not in the range 0..2, the vector
    remains unchanged and an exception is thrown. */
    void SetColumn(int column, const vector3 &v) throw(OBError);

    //! access function
    /*! \warning If column is not in the range 0..2, the vector
    remains unchanged and an exception is thrown. */
    void SetRow(int row, const vector3 &v) throw(OBError);

    //! access function
    /*! \warning If col is not in the range 0..2, an exception is
    thrown. */
    vector3 GetColumn(unsigned int col) const throw(OBError);

    //! access function
    /*! \warning If row is not in the range 0..2, an exception is
    thrown. */
    vector3 GetRow(unsigned int row) const throw(OBError);


    //! divides all entries of the matrix by a scalar c
    matrix3x3 &operator/=(const double &c);

    void SetupRotMat(double,double,double);

    //! calculates a matrix that represents reflection on a plane
    void PlaneReflection(const vector3 &norm);

    //! calculates a rotation matrix
    void RotAboutAxisByAngle(const vector3 &axis, const double angle);

    void FillOrth(double,double,double,double,double,double);

    //! find the eigenvalues and -vectors of a symmetric matrix
    matrix3x3 findEigenvectorsIfSymmetric(vector3 &eigenvals) const throw(OBError);

    //! matrix-vector multiplication
    friend OBAPI vector3 operator *(const matrix3x3 &,const vector3 &);

    //! matrix-matrix multiplication
    friend OBAPI matrix3x3 operator *(const matrix3x3 &,const matrix3x3 &);

    friend OBAPI std::ostream& operator<< ( std::ostream&, const matrix3x3 & ) ;

    //! eigenvalue calculation
    static void jacobi(unsigned int n, double *a, double *d, double *v);
};

OBAPI vector3 center_coords(double*,int);
}

#endif // OB_MATRIX3x3_H

//! \file matrix3x3.h
//! \brief Handle 3D Rotation matrix.
Revision:	2440
Committed:	Wed Nov 16 19:42:11 2005 UTC (19 years, 11 months ago) by tim
File size:	7345 byte(s)
Log Message:	adding openbabel
#	Content
1	/**********************************************************************
2	matrix3x3.cpp - Handle 3D Rotation matrix.
3
4	Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
5	Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
6
7	This file is part of the Open Babel project.
8	For more information, see <http://openbabel.sourceforge.net/>
9
10	This program is free software; you can redistribute it and/or modify
11	it under the terms of the GNU General Public License as published by
12	the Free Software Foundation version 2 of the License.
13
14	This program is distributed in the hope that it will be useful,
15	but WITHOUT ANY WARRANTY; without even the implied warranty of
16	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17	GNU General Public License for more details.
18	***********************************************************************/
19
20	#ifndef OB_MATRIX3x3_H
21	#define OB_MATRIX3x3_H
22
23	#include "oberror.hpp"
24
25	#if HAVE_IOSTREAM
26	#include <iostream>
27	#elif HAVE_IOSTREAM_H
28	#include <iostream.h>
29	#endif
30
31	#if HAVE_FSTREAM
32	#include <fstream>
33	#elif HAVE_FSTREAM_H
34	#include <fstream.h>
35	#endif
36
37	#include <math.h>
38
39	#include "obutil.hpp"
40	#include "vector3.hpp"
41
42	#ifndef PI
43	#define PI 3.1415926535897932384626433
44	#endif
45
46	#ifndef RAD_TO_DEG
47	#define RAD_TO_DEG 180.0/PI
48	#endif
49
50	#ifndef DEG_TO_RAD
51	#define DEG_TO_RAD PI/180.0
52	#endif
53
54	namespace OpenBabel
55	{
56
57	// class introduction in matrix3x3.cpp
58	class OBAPI matrix3x3
59	{
60	//! Elements of the matrix
61	/*! This array holds the matrix. The first index refers to the
62	row, the second the column. */
63	double ele[3][3];
64
65	public:
66	//! constructs the zero-matrix
67	matrix3x3(void)
68	{
69	ele[0][0] = 0.0;
70	ele[0][1] = 0.0;
71	ele[0][2] = 0.0;
72	ele[1][0] = 0.0;
73	ele[1][1] = 0.0;
74	ele[1][2] = 0.0;
75	ele[2][0] = 0.0;
76	ele[2][1] = 0.0;
77	ele[2][2] = 0.0;
78	}
79
80	//! constructs s times the unit matrix
81	matrix3x3(double s)
82	{
83	ele[0][0] = s;
84	ele[0][1] = 0.0;
85	ele[0][2] = 0.0;
86	ele[1][0] = 0.0;
87	ele[1][1] = s;
88	ele[1][2] = 0.0;
89	ele[2][0] = 0.0;
90	ele[2][1] = 0.0;
91	ele[2][2] = s;
92	}
93
94	//! constructs a matrix from row vectors
95	matrix3x3(vector3 row1,vector3 row2,vector3 row3)
96	{
97	ele[0][0] = row1.x();
98	ele[0][1] = row1.y();
99	ele[0][2] = row1.z();
100	ele[1][0] = row2.x();
101	ele[1][1] = row2.y();
102	ele[1][2] = row2.z();
103	ele[2][0] = row3.x();
104	ele[2][1] = row3.y();
105	ele[2][2] = row3.z();
106	}
107
108	//! constructs a matrix from a 3x3-array of doubles
109	/*! constructs a matrix from a 3x3-array of doubles. The first
110	index represents the row, the second index the column */
111	matrix3x3(double d[3][3])
112	{
113	ele[0][0] = d[0][0];
114	ele[0][1] = d[0][1];
115	ele[0][2] = d[0][2];
116	ele[1][0] = d[1][0];
117	ele[1][1] = d[1][1];
118	ele[1][2] = d[1][2];
119	ele[2][0] = d[2][0];
120	ele[2][1] = d[2][1];
121	ele[2][2] = d[2][2];
122	}
123
124	//! access function
125	/*! writes the matrix into the 1-dimensional array m, row by
126	row. The array must be able to hold 9 doubles, otherwise your
127	prgram will segfault. */
128	void GetArray(double *m)
129	{
130	m[0] = ele[0][0];
131	m[1] = ele[0][1];
132	m[2] = ele[0][2];
133	m[3] = ele[1][0];
134	m[4] = ele[1][1];
135	m[5] = ele[1][2];
136	m[6] = ele[2][0];
137	m[7] = ele[2][1];
138	m[8] = ele[2][2];
139	}
140
141	//! Calculates the inverse of a matrix.
142	matrix3x3 inverse(void) const throw(OBError);
143
144	//! Calculates the transpose of a matrix.
145	matrix3x3 transpose(void) const;
146
147	//! generates a matrix for a random rotation
148	void randomRotation(OBRandom &rnd);
149
150	//! returns the determinant of the matrix
151	double determinant() const;
152
153	//! Checks if a matrix is symmetric
154	bool isSymmetric(void) const;
155
156	//! Checks if a matrix is orthogonal
157	/*! This method checks if a matrix describes an orthogonal
158	transformation, i.e. if all column vectors are normalized and
159	are mutually orthogonal. An orthogonal transformation is a
160	transformation the preserves length and angle.
161
162	The check is performed using the method isUnitMatrix() to
163	check if
164	\code
165	this transpose()
166	\endcode
167	is a unit matrix. The criterion is therefore numerically quite
168	tight. */
169	bool isOrthogonal(void) const
170	{
171	return (this transpose()).isUnitMatrix();
172	};
173
174	//! Checks if a matrix is diagonal
175	bool isDiagonal(void) const;
176
177	//! Checks if a matrix is the unit matrix
178	bool isUnitMatrix(void) const;
179
180	//! access function
181	/*! \warning row or column are not in the range 0..2, random
182	results are returned, and your program may even
183	segfault. (Stefan Kebekus)
184
185	\todo Replace this method with a more fool-proof version.
186	*/
187	double Get(int row,int column) const
188	{
189	if (row >= 0 && row <= 2 && column >= 0 && column <= 2)
190	return(ele[row][column]);
191	else
192	return 0.0f;
193	}
194
195	//! access function
196	/*! \warning if row or column are not in the range 0..2, random
197	variables are overwritten, and your program may
198	segfault. (Stefan Kebekus)
199
200	\todo Replace this method with a more fool-proof version.
201	*/
202	void Set(int row,int column, double v)
203	{
204	if (row >= 0 && row <= 2 && column >= 0 && column <= 2)
205	ele[row][column]= v;
206	}
207
208	//! access function
209	/*! \warning If column is not in the range 0..2, the vector
210	remains unchanged and an exception is thrown. */
211	void SetColumn(int column, const vector3 &v) throw(OBError);
212
213	//! access function
214	/*! \warning If column is not in the range 0..2, the vector
215	remains unchanged and an exception is thrown. */
216	void SetRow(int row, const vector3 &v) throw(OBError);
217
218	//! access function
219	/*! \warning If col is not in the range 0..2, an exception is
220	thrown. */
221	vector3 GetColumn(unsigned int col) const throw(OBError);
222
223	//! access function
224	/*! \warning If row is not in the range 0..2, an exception is
225	thrown. */
226	vector3 GetRow(unsigned int row) const throw(OBError);
227
228
229	//! divides all entries of the matrix by a scalar c
230	matrix3x3 &operator/=(const double &c);
231
232	void SetupRotMat(double,double,double);
233
234	//! calculates a matrix that represents reflection on a plane
235	void PlaneReflection(const vector3 &norm);
236
237	//! calculates a rotation matrix
238	void RotAboutAxisByAngle(const vector3 &axis, const double angle);
239
240	void FillOrth(double,double,double,double,double,double);
241
242	//! find the eigenvalues and -vectors of a symmetric matrix
243	matrix3x3 findEigenvectorsIfSymmetric(vector3 &eigenvals) const throw(OBError);
244
245	//! matrix-vector multiplication
246	friend OBAPI vector3 operator *(const matrix3x3 &,const vector3 &);
247
248	//! matrix-matrix multiplication
249	friend OBAPI matrix3x3 operator *(const matrix3x3 &,const matrix3x3 &);
250
251	friend OBAPI std::ostream& operator<< ( std::ostream&, const matrix3x3 & ) ;
252
253	//! eigenvalue calculation
254	static void jacobi(unsigned int n, double a, double d, double *v);
255	};
256
257	OBAPI vector3 center_coords(double*,int);
258	}
259
260	#endif // OB_MATRIX3x3_H
261
262	//! \file matrix3x3.h
263	//! \brief Handle 3D Rotation matrix.