1 /**************************************************************************/
2 /*  Copyright 2009 Tim Day                                                */
3 /*                                                                        */
4 /*  This file is part of Fracplanet                                       */
5 /*                                                                        */
6 /*  Fracplanet is free software: you can redistribute it and/or modify    */
7 /*  it under the terms of the GNU General Public License as published by  */
8 /*  the Free Software Foundation, either version 3 of the License, or     */
9 /*  (at your option) any later version.                                   */
10 /*                                                                        */
11 /*  Fracplanet is distributed in the hope that it will be useful,         */
12 /*  but WITHOUT ANY WARRANTY; without even the implied warranty of        */
13 /*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         */
14 /*  GNU General Public License for more details.                          */
15 /*                                                                        */
16 /*  You should have received a copy of the GNU General Public License     */
17 /*  along with Fracplanet.  If not, see <http://www.gnu.org/licenses/>.   */
18 /**************************************************************************/
19 
20 /*! \file
21   \brief Interface for class Matrix33.
22 */
23 
24 #ifndef _matrix33_h_
25 #define _matrix33_h_
26 
27 #include "common.h"
28 #include "xyz.h"
29 
30 //! Class to hold 3x3 matrics
31 class Matrix33
32 {
33  public:
34 
35   //! Column vectors of matrix
36   XYZ basis[3];
37 
38   //! Null constructor
39   /*! NB Doesn't set up identity or anything useful
40    */
Matrix33()41   Matrix33()
42     {}
43 
44   //! Construct from column vectors
Matrix33(const XYZ & rx,const XYZ & ry,const XYZ & rz)45   Matrix33(const XYZ& rx,const XYZ& ry,const XYZ& rz)
46     {
47       basis[0]=rx;
48       basis[1]=ry;
49       basis[2]=rz;
50     }
51 
52   //! Destructor.
~Matrix33()53   ~Matrix33()
54     {}
55 
assign(const Matrix33 & m)56   void assign(const Matrix33& m)
57     {
58       basis[0]=m.basis[0];
59       basis[1]=m.basis[1];
60       basis[2]=m.basis[2];
61     }
62 
63   //! Access a given element
element(uint row,uint col)64   const float& element(uint row,uint col) const
65     {
66       return basis[col].element(row);
67     }
68 
69   //! Access a given element
element(uint row,uint col)70   float& element(uint row,uint col)
71     {
72       return basis[col].element(row);
73     }
74 
75   //! Extract copy of the first row of matrix
row0()76   const XYZ row0() const
77     {
78       return XYZ(basis[0].x,basis[1].x,basis[2].x);
79     }
80 
81   //! Extract copy of the second row of matrix
row1()82   const XYZ row1() const
83     {
84       return XYZ(basis[0].y,basis[1].y,basis[2].y);
85     }
86 
87   //! Extract copy of the third row of matrix
row2()88   const XYZ row2() const
89     {
90       return XYZ(basis[0].z,basis[1].z,basis[2].z);
91     }
92 
93   //! Cofactor of an element
94   float cofactor(uint row,uint col) const;
95 
96   //! Determinant of matrix
97   float determinant() const;
98 
99   //! Return inverse of matrix
100   const Matrix33 inverted() const;
101 };
102 
103 //! Multiplication by scalar
104 inline const Matrix33 operator*(float k,const Matrix33& m)
105 {
106   return Matrix33(k*m.basis[0],k*m.basis[1],k*m.basis[2]);
107 }
108 
109 //! Multiplication by scalar
110 inline const Matrix33 operator*(const Matrix33& m,float k)
111 {
112   return Matrix33(m.basis[0]*k,m.basis[1]*k,m.basis[2]*k);
113 }
114 
115 //! Division by scalar
116 inline const Matrix33 operator/(const Matrix33& m,float k)
117 {
118   return Matrix33(m.basis[0]/k,m.basis[1]/k,m.basis[2]/k);
119 }
120 
121 //! Apply matrix to vector
122 inline const XYZ operator*(const Matrix33& m,const XYZ& v)
123 {
124   return XYZ
125     (
126      m.row0()%v,
127      m.row1()%v,
128      m.row2()%v
129      );
130 }
131 
132 //! Matrix multiplication
133 inline const Matrix33 operator*(const Matrix33& a,const Matrix33& b)
134 {
135   return Matrix33
136     (
137      a*b.basis[0],
138      a*b.basis[1],
139      a*b.basis[2]
140      );
141 }
142 
143 //! 3x3 Identity matrix
144 class Matrix33Identity : public Matrix33
145 {
146  public:
147 
Matrix33Identity()148   Matrix33Identity()
149     {
150       basis[0]=XYZ(1.0f,0.0f,0.0f);
151       basis[1]=XYZ(0.0f,1.0f,0.0f);
152       basis[2]=XYZ(0.0f,0.0f,1.0f);
153     }
154 };
155 
156 class Matrix33RotateAboutX : public Matrix33
157 {
158  public:
159 
Matrix33RotateAboutX(float angle)160   Matrix33RotateAboutX(float angle)
161     {
162       const float ca=cos(angle);
163       const float sa=sin(angle);
164       basis[0]=XYZ(1.0f,0.0f,0.0f);
165       basis[1]=XYZ(0.0f, ca,  sa);
166       basis[2]=XYZ(0.0f,-sa,  ca);
167     }
168 };
169 
170 class Matrix33RotateAboutY : public Matrix33
171 {
172  public:
173 
Matrix33RotateAboutY(float angle)174   Matrix33RotateAboutY(float angle)
175     {
176       const float ca=cos(angle);
177       const float sa=sin(angle);
178       basis[0]=XYZ(ca,  0.0f,-sa);
179       basis[1]=XYZ(0.0f,1.0f,0.0f);
180       basis[2]=XYZ(sa  ,0.0f,ca);
181     }
182 };
183 
184 class Matrix33RotateAboutZ : public Matrix33
185 {
186  public:
187 
Matrix33RotateAboutZ(float angle)188   Matrix33RotateAboutZ(float angle)
189     {
190       const float ca=cos(angle);
191       const float sa=sin(angle);
192       basis[0]=XYZ( ca,sa,0.0f);
193       basis[1]=XYZ(-sa,ca,0.0f);
194       basis[2]=XYZ(0.0f,0.0f,1.0f);
195     }
196 };
197 
198 class Matrix33RotateAboutAxis : public Matrix33
199 {
200  public:
201 
202   //! NB Axis must be normalized.
Matrix33RotateAboutAxis(const XYZ & axis,float angle)203   Matrix33RotateAboutAxis(const XYZ& axis,float angle)
204     {
205       // Want 2 vectors perpendicular to axis TODO: Check for degenerate cases
206       const XYZ axis_ortho0((axis*XYZ(1.0,0.0,0.0)).normalised());
207       const XYZ axis_ortho1(axis*axis_ortho0);
208 
209       // The matrix which rotates identity basis to axis&orthos.  z axis goes to passed in axis
210       const Matrix33 xyz_to_axis
211 	(
212 	 axis_ortho0,
213 	 axis_ortho1,
214 	 axis
215 	 );
216 
217       const Matrix33 axis_to_xyz(xyz_to_axis.inverted());
218 
219       assign(xyz_to_axis*Matrix33RotateAboutZ(angle)*axis_to_xyz);
220     }
221 };
222 
223 #endif
224