~ubuntu-branches/ubuntu/maverick/freecad/maverick

« back to all changes in this revision

Viewing changes to src/Mod/Mesh/App/WildMagic4/Wm4Matrix3.h

  • Committer: Bazaar Package Importer
  • Author(s): Teemu Ikonen
  • Date: 2009-07-16 18:37:41 UTC
  • Revision ID: james.westby@ubuntu.com-20090716183741-oww9kcxqrk991i1n
Tags: upstream-0.8.2237
ImportĀ upstreamĀ versionĀ 0.8.2237

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
// Wild Magic Source Code
 
2
// David Eberly
 
3
// http://www.geometrictools.com
 
4
// Copyright (c) 1998-2007
 
5
//
 
6
// This library is free software; you can redistribute it and/or modify it
 
7
// under the terms of the GNU Lesser General Public License as published by
 
8
// the Free Software Foundation; either version 2.1 of the License, or (at
 
9
// your option) any later version.  The license is available for reading at
 
10
// either of the locations:
 
11
//     http://www.gnu.org/copyleft/lgpl.html
 
12
//     http://www.geometrictools.com/License/WildMagicLicense.pdf
 
13
// The license applies to versions 0 through 4 of Wild Magic.
 
14
//
 
15
// Version: 4.0.2 (2006/09/05)
 
16
 
 
17
#ifndef WM4MATRIX3_H
 
18
#define WM4MATRIX3_H
 
19
 
 
20
// Matrix operations are applied on the left.  For example, given a matrix M
 
21
// and a vector V, matrix-times-vector is M*V.  That is, V is treated as a
 
22
// column vector.  Some graphics APIs use V*M where V is treated as a row
 
23
// vector.  In this context the "M" matrix is really a transpose of the M as
 
24
// represented in Wild Magic.  Similarly, to apply two matrix operations M0
 
25
// and M1, in that order, you compute M1*M0 so that the transform of a vector
 
26
// is (M1*M0)*V = M1*(M0*V).  Some graphics APIs use M0*M1, but again these
 
27
// matrices are the transpose of those as represented in Wild Magic.  You
 
28
// must therefore be careful about how you interface the transformation code
 
29
// with graphics APIS.
 
30
//
 
31
// For memory organization it might seem natural to chose Real[N][N] for the
 
32
// matrix storage, but this can be a problem on a platform/console that
 
33
// chooses to store the data in column-major rather than row-major format.
 
34
// To avoid potential portability problems, the matrix is stored as Real[N*N]
 
35
// and organized in row-major order.  That is, the entry of the matrix in row
 
36
// r (0 <= r < N) and column c (0 <= c < N) is stored at index i = c+N*r
 
37
// (0 <= i < N*N).
 
38
 
 
39
// The (x,y,z) coordinate system is assumed to be right-handed.  Coordinate
 
40
// axis rotation matrices are of the form
 
41
//   RX =    1       0       0
 
42
//           0     cos(t) -sin(t)
 
43
//           0     sin(t)  cos(t)
 
44
// where t > 0 indicates a counterclockwise rotation in the yz-plane
 
45
//   RY =  cos(t)    0     sin(t)
 
46
//           0       1       0
 
47
//        -sin(t)    0     cos(t)
 
48
// where t > 0 indicates a counterclockwise rotation in the zx-plane
 
49
//   RZ =  cos(t) -sin(t)    0
 
50
//         sin(t)  cos(t)    0
 
51
//           0       0       1
 
52
// where t > 0 indicates a counterclockwise rotation in the xy-plane.
 
53
 
 
54
#include "Wm4FoundationLIB.h"
 
55
#include "Wm4Vector3.h"
 
56
 
 
57
namespace Wm4
 
58
{
 
59
 
 
60
template <class Real>
 
61
class Matrix3
 
62
{
 
63
public:
 
64
    // If bZero is true, create the zero matrix.  Otherwise, create the
 
65
    // identity matrix.
 
66
    Matrix3 (bool bZero = true);
 
67
 
 
68
    // copy constructor
 
69
    Matrix3 (const Matrix3& rkM);
 
70
 
 
71
    // input Mrc is in row r, column c.
 
72
    Matrix3 (Real fM00, Real fM01, Real fM02,
 
73
             Real fM10, Real fM11, Real fM12,
 
74
             Real fM20, Real fM21, Real fM22);
 
75
 
 
76
    // Create a matrix from an array of numbers.  The input array is
 
77
    // interpreted based on the Boolean input as
 
78
    //   true:  entry[0..8]={m00,m01,m02,m10,m11,m12,m20,m21,m22} [row major]
 
79
    //   false: entry[0..8]={m00,m10,m20,m01,m11,m21,m02,m12,m22} [col major]
 
80
    Matrix3 (const Real afEntry[9], bool bRowMajor);
 
81
 
 
82
    // Create matrices based on vector input.  The Boolean is interpreted as
 
83
    //   true: vectors are columns of the matrix
 
84
    //   false: vectors are rows of the matrix
 
85
    Matrix3 (const Vector3<Real>& rkU, const Vector3<Real>& rkV,
 
86
        const Vector3<Real>& rkW, bool bColumns);
 
87
    Matrix3 (const Vector3<Real>* akV, bool bColumns);
 
88
 
 
89
    // create a diagonal matrix
 
90
    Matrix3 (Real fM00, Real fM11, Real fM22);
 
91
 
 
92
    // Create rotation matrices (positive angle - counterclockwise).  The
 
93
    // angle must be in radians, not degrees.
 
94
    Matrix3 (const Vector3<Real>& rkAxis, Real fAngle);
 
95
 
 
96
    // create a tensor product U*V^T
 
97
    Matrix3 (const Vector3<Real>& rkU, const Vector3<Real>& rkV);
 
98
 
 
99
    // create various matrices
 
100
    Matrix3& MakeZero ();
 
101
    Matrix3& MakeIdentity ();
 
102
    Matrix3& MakeDiagonal (Real fM00, Real fM11, Real fM22);
 
103
    Matrix3& FromAxisAngle (const Vector3<Real>& rkAxis, Real fAngle);
 
104
    Matrix3& MakeTensorProduct (const Vector3<Real>& rkU,
 
105
        const Vector3<Real>& rkV);
 
106
 
 
107
    // member access
 
108
    operator const Real* () const;
 
109
    operator Real* ();
 
110
    const Real* operator[] (int iRow) const;
 
111
    Real* operator[] (int iRow);
 
112
    Real operator() (int iRow, int iCol) const;
 
113
    Real& operator() (int iRow, int iCol);
 
114
    void SetRow (int iRow, const Vector3<Real>& rkV);
 
115
    Vector3<Real> GetRow (int iRow) const;
 
116
    void SetColumn (int iCol, const Vector3<Real>& rkV);
 
117
    Vector3<Real> GetColumn (int iCol) const;
 
118
    void GetColumnMajor (Real* afCMajor) const;
 
119
 
 
120
    // assignment
 
121
    Matrix3& operator= (const Matrix3& rkM);
 
122
 
 
123
    // comparison
 
124
    bool operator== (const Matrix3& rkM) const;
 
125
    bool operator!= (const Matrix3& rkM) const;
 
126
    bool operator<  (const Matrix3& rkM) const;
 
127
    bool operator<= (const Matrix3& rkM) const;
 
128
    bool operator>  (const Matrix3& rkM) const;
 
129
    bool operator>= (const Matrix3& rkM) const;
 
130
 
 
131
    // arithmetic operations
 
132
    Matrix3 operator+ (const Matrix3& rkM) const;
 
133
    Matrix3 operator- (const Matrix3& rkM) const;
 
134
    Matrix3 operator* (const Matrix3& rkM) const;
 
135
    Matrix3 operator* (Real fScalar) const;
 
136
    Matrix3 operator/ (Real fScalar) const;
 
137
    Matrix3 operator- () const;
 
138
 
 
139
    // arithmetic updates
 
140
    Matrix3& operator+= (const Matrix3& rkM);
 
141
    Matrix3& operator-= (const Matrix3& rkM);
 
142
    Matrix3& operator*= (Real fScalar);
 
143
    Matrix3& operator/= (Real fScalar);
 
144
 
 
145
    // matrix times vector
 
146
    Vector3<Real> operator* (const Vector3<Real>& rkV) const;  // M * v
 
147
 
 
148
    // other operations
 
149
    Matrix3 Transpose () const;  // M^T
 
150
    Matrix3 TransposeTimes (const Matrix3& rkM) const;  // this^T * M
 
151
    Matrix3 TimesTranspose (const Matrix3& rkM) const;  // this * M^T
 
152
    Matrix3 Inverse () const;
 
153
    Matrix3 Adjoint () const;
 
154
    Real Determinant () const;
 
155
    Real QForm (const Vector3<Real>& rkU,
 
156
        const Vector3<Real>& rkV) const;  // u^T*M*v
 
157
    Matrix3 TimesDiagonal (const Vector3<Real>& rkDiag) const;  // M*D
 
158
    Matrix3 DiagonalTimes (const Vector3<Real>& rkDiag) const;  // D*M
 
159
 
 
160
    // The matrix must be a rotation for these functions to be valid.  The
 
161
    // last function uses Gram-Schmidt orthonormalization applied to the
 
162
    // columns of the rotation matrix.  The angle must be in radians, not
 
163
    // degrees.
 
164
    void ToAxisAngle (Vector3<Real>& rkAxis, Real& rfAngle) const;
 
165
    void Orthonormalize ();
 
166
 
 
167
    // The matrix must be symmetric.  Factor M = R * D * R^T where
 
168
    // R = [u0|u1|u2] is a rotation matrix with columns u0, u1, and u2 and
 
169
    // D = diag(d0,d1,d2) is a diagonal matrix whose diagonal entries are d0,
 
170
    // d1, and d2.  The eigenvector u[i] corresponds to eigenvector d[i].
 
171
    // The eigenvalues are ordered as d0 <= d1 <= d2.
 
172
    void EigenDecomposition (Matrix3& rkRot, Matrix3& rkDiag) const;
 
173
 
 
174
    // Create rotation matrices from Euler angles.
 
175
    Matrix3& FromEulerAnglesXYZ (Real fXAngle, Real fYAngle, Real fZAngle);
 
176
    Matrix3& FromEulerAnglesXZY (Real fXAngle, Real fZAngle, Real fYAngle);
 
177
    Matrix3& FromEulerAnglesYXZ (Real fYAngle, Real fXAngle, Real fZAngle);
 
178
    Matrix3& FromEulerAnglesYZX (Real fYAngle, Real fZAngle, Real fXAngle);
 
179
    Matrix3& FromEulerAnglesZXY (Real fZAngle, Real fXAngle, Real fYAngle);
 
180
    Matrix3& FromEulerAnglesZYX (Real fZAngle, Real fYAngle, Real fXAngle);
 
181
 
 
182
    // Extract Euler angles from rotation matrices.  The return value is
 
183
    // 'true' iff the factorization is unique relative to certain angle
 
184
    // ranges.  That is, if (U,V,W) is some permutation of (X,Y,Z), the angle
 
185
    // ranges for the outputs from ToEulerAnglesUVW(uAngle,vAngle,wAngle) are
 
186
    // uAngle in [-pi,pi], vAngle in [-pi/2,pi/2], and wAngle in [-pi,pi].  If
 
187
    // the function returns 'false', wAngle is 0 and vAngle is either pi/2 or
 
188
    // -pi/2.
 
189
    bool ToEulerAnglesXYZ (Real& rfXAngle, Real& rfYAngle, Real& rfZAngle)
 
190
        const;
 
191
    bool ToEulerAnglesXZY (Real& rfXAngle, Real& rfZAngle, Real& rfYAngle)
 
192
        const;
 
193
    bool ToEulerAnglesYXZ (Real& rfYAngle, Real& rfXAngle, Real& rfZAngle)
 
194
        const;
 
195
    bool ToEulerAnglesYZX (Real& rfYAngle, Real& rfZAngle, Real& rfXAngle)
 
196
        const;
 
197
    bool ToEulerAnglesZXY (Real& rfZAngle, Real& rfXAngle, Real& rfYAngle)
 
198
        const;
 
199
    bool ToEulerAnglesZYX (Real& rfZAngle, Real& rfYAngle, Real& rfXAngle)
 
200
        const;
 
201
 
 
202
    // SLERP (spherical linear interpolation) without quaternions.  Computes
 
203
    // R(t) = R0*(Transpose(R0)*R1)^t.  If Q is a rotation matrix with
 
204
    // unit-length axis U and angle A, then Q^t is a rotation matrix with
 
205
    // unit-length axis U and rotation angle t*A.
 
206
    Matrix3& Slerp (Real fT, const Matrix3& rkR0, const Matrix3& rkR1);
 
207
 
 
208
    // Singular value decomposition, M = L*S*R, where L and R are orthogonal
 
209
    // and S is a diagonal matrix whose diagonal entries are nonnegative.
 
210
    void SingularValueDecomposition (Matrix3& rkL, Matrix3& rkS,
 
211
        Matrix3& rkR) const;
 
212
    void SingularValueComposition (const Matrix3& rkL, const Matrix3& rkS,
 
213
        const Matrix3& rkR);
 
214
 
 
215
    // factor M = Q*D*U with orthogonal Q, diagonal D, upper triangular U
 
216
    void QDUDecomposition (Matrix3& rkQ, Matrix3& rkD, Matrix3& rkU) const;
 
217
 
 
218
    // special matrices
 
219
    WM4_FOUNDATION_ITEM static const Matrix3 ZERO;
 
220
    WM4_FOUNDATION_ITEM static const Matrix3 IDENTITY;
 
221
 
 
222
private:
 
223
    // Support for eigendecomposition.  The Tridiagonalize function applies
 
224
    // a Householder transformation to the matrix.  If that transformation
 
225
    // is the identity (the matrix is already tridiagonal), then the return
 
226
    // value is 'false'.  Otherwise, the transformation is a reflection and
 
227
    // the return value is 'true'.  The QLAlgorithm returns 'true' iff the
 
228
    // QL iteration scheme converged.
 
229
    bool Tridiagonalize (Real afDiag[3], Real afSubd[2]);
 
230
    bool QLAlgorithm (Real afDiag[3], Real afSubd[2]);
 
231
 
 
232
    // support for singular value decomposition
 
233
    static void Bidiagonalize (Matrix3& rkA, Matrix3& rkL, Matrix3& rkR);
 
234
    static void GolubKahanStep (Matrix3& rkA, Matrix3& rkL, Matrix3& rkR);
 
235
 
 
236
    // support for comparisons
 
237
    int CompareArrays (const Matrix3& rkM) const;
 
238
 
 
239
    Real m_afEntry[9];
 
240
};
 
241
 
 
242
// c * M
 
243
template <class Real>
 
244
Matrix3<Real> operator* (Real fScalar, const Matrix3<Real>& rkM);
 
245
 
 
246
// v^T * M
 
247
template <class Real>
 
248
Vector3<Real> operator* (const Vector3<Real>& rkV, const Matrix3<Real>& rkM);
 
249
 
 
250
#include "Wm4Matrix3.inl"
 
251
 
 
252
typedef Matrix3<float> Matrix3f;
 
253
typedef Matrix3<double> Matrix3d;
 
254
 
 
255
}
 
256
 
 
257
#endif