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);

    //! 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);

    //! 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);

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

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


    //! 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);

    //! 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:	778
Committed:	Fri Dec 2 20:11:25 2005 UTC (19 years, 11 months ago) by gezelter
File size:	7231 byte(s)
Log Message:	Errors are no longer thrown with consts (fixes compilation bug on the Mac).
#	User	Rev	Content
1	tim	741	/**********************************************************************
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	gezelter	778	matrix3x3 inverse(void);
143	tim	741
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	gezelter	778	void SetColumn(int column, const vector3 &v);
212	tim	741
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	gezelter	778	void SetRow(int row, const vector3 &v);
217	tim	741
218			//! access function
219			/*! \warning If col is not in the range 0..2, an exception is
220			thrown. */
221	gezelter	778	vector3 GetColumn(unsigned int col);
222	tim	741
223			//! access function
224			/*! \warning If row is not in the range 0..2, an exception is
225			thrown. */
226	gezelter	778	vector3 GetRow(unsigned int row);
227	tim	741
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	gezelter	778	matrix3x3 findEigenvectorsIfSymmetric(vector3 &eigenvals);
244	tim	741
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.