| 1 |
/* ----------------------------------------------------------------------- *
|
1 |
/* ----------------------------------------------------------------------- *
|
| 2 |
* This file is part of GEL, http://www.imm.dtu.dk/GEL
|
2 |
* This file is part of GEL, http://www.imm.dtu.dk/GEL
|
| 3 |
* Copyright (C) the authors and DTU Informatics
|
3 |
* Copyright (C) the authors and DTU Informatics
|
| 4 |
* For license and list of authors, see ../../doc/intro.pdf
|
4 |
* For license and list of authors, see ../../doc/intro.pdf
|
| 5 |
* @brief Abstract 4x4 floating point matrix class
|
5 |
* @brief Abstract 4x4 floating point matrix class
|
| 6 |
* ----------------------------------------------------------------------- */
|
6 |
* ----------------------------------------------------------------------- */
|
| 7 |
|
7 |
|
| 8 |
/** @file ArithSqMat4x4Float.h
|
8 |
/** @file ArithSqMat4x4Float.h
|
| 9 |
* @brief Abstract 4x4 floating point matrix class
|
9 |
* @brief Abstract 4x4 floating point matrix class
|
| 10 |
*/
|
10 |
*/
|
| 11 |
|
11 |
|
| 12 |
#ifndef __CGLA_ARITHSQMAT4X4FLOAT_H
|
12 |
#ifndef __CGLA_ARITHSQMAT4X4FLOAT_H
|
| 13 |
#define __CGLA_ARITHSQMAT4X4FLOAT_H
|
13 |
#define __CGLA_ARITHSQMAT4X4FLOAT_H
|
| 14 |
|
14 |
|
| 15 |
#include "ExceptionStandard.h"
|
15 |
#include "ExceptionStandard.h"
|
| 16 |
#include "CGLA.h"
|
16 |
#include "CGLA.h"
|
| 17 |
#include "Vec3f.h"
|
17 |
#include "Vec3f.h"
|
| 18 |
#include "Vec3Hf.h"
|
18 |
#include "Vec3Hf.h"
|
| 19 |
#include "Vec4f.h"
|
19 |
#include "Vec4f.h"
|
| 20 |
#include "ArithSqMatFloat.h"
|
20 |
#include "ArithSqMatFloat.h"
|
| 21 |
|
21 |
|
| 22 |
|
22 |
|
| 23 |
namespace CGLA
|
23 |
namespace CGLA
|
| 24 |
{
|
24 |
{
|
| 25 |
CGLA_DERIVEEXCEPTION(Mat4x4fException)
|
25 |
CGLA_DERIVEEXCEPTION(Mat4x4fException)
|
| 26 |
CGLA_DERIVEEXCEPTION(Mat4x4fNotAffine)
|
26 |
CGLA_DERIVEEXCEPTION(Mat4x4fNotAffine)
|
| 27 |
CGLA_DERIVEEXCEPTION(Mat4x4fSingular)
|
27 |
CGLA_DERIVEEXCEPTION(Mat4x4fSingular)
|
| 28 |
|
28 |
|
| 29 |
/** \brief 4 by 4 float matrix template.
|
29 |
/** \brief 4 by 4 float matrix template.
|
| 30 |
|
30 |
|
| 31 |
this class template is useful for transformations such as perspective
|
31 |
this class template is useful for transformations such as perspective
|
| 32 |
projections or translation where 3x3 matrices do not suffice. */
|
32 |
projections or translation where 3x3 matrices do not suffice. */
|
| 33 |
template<class VT, class M>
|
33 |
template<class VT, class M>
|
| 34 |
class ArithSqMat4x4Float: public ArithSqMatFloat<VT, M, 4>
|
34 |
class ArithSqMat4x4Float: public ArithSqMatFloat<VT, M, 4>
|
| 35 |
{
|
35 |
{
|
| 36 |
public:
|
36 |
public:
|
| 37 |
|
37 |
|
| 38 |
/// Vector type
|
38 |
/// Vector type
|
| 39 |
typedef VT VectorType;
|
39 |
typedef VT VectorType;
|
| 40 |
|
40 |
|
| 41 |
/// The type of a matrix element
|
41 |
/// The type of a matrix element
|
| 42 |
typedef typename VT::ScalarType ScalarType;
|
42 |
typedef typename VT::ScalarType ScalarType;
|
| 43 |
|
43 |
|
| 44 |
public:
|
44 |
public:
|
| 45 |
|
45 |
|
| 46 |
/// Construct a Mat4x4f from four V vectors
|
46 |
/// Construct a Mat4x4f from four V vectors
|
| 47 |
ArithSqMat4x4Float(VT a, VT b, VT c, VT d):
|
47 |
ArithSqMat4x4Float(VT a, VT b, VT c, VT d):
|
| 48 |
ArithSqMatFloat<VT, M, 4> (a,b,c,d) {}
|
48 |
ArithSqMatFloat<VT, M, 4> (a,b,c,d) {}
|
| 49 |
|
49 |
|
| 50 |
/// Construct the NAN matrix
|
50 |
/// Construct the NAN matrix
|
| 51 |
ArithSqMat4x4Float() {}
|
51 |
ArithSqMat4x4Float() {}
|
| 52 |
|
52 |
|
| 53 |
/// Construct matrix where all values are equal to constructor argument.
|
53 |
/// Construct matrix where all values are equal to constructor argument.
|
| 54 |
explicit ArithSqMat4x4Float(ScalarType _a):
|
54 |
explicit ArithSqMat4x4Float(ScalarType _a):
|
| 55 |
ArithSqMatFloat<VT,M,4>(_a) {}
|
55 |
ArithSqMatFloat<VT,M,4>(_a) {}
|
| 56 |
|
56 |
|
| 57 |
/** Multiply vector onto matrix. Here the fourth coordinate
|
57 |
/** Multiply vector onto matrix. Here the fourth coordinate
|
| 58 |
is se to 0. This removes any translation from the matrix.
|
58 |
is se to 0. This removes any translation from the matrix.
|
| 59 |
Useful if one wants to transform a vector which does not
|
59 |
Useful if one wants to transform a vector which does not
|
| 60 |
represent a point but a direction. Note that this is not
|
60 |
represent a point but a direction. Note that this is not
|
| 61 |
correct for transforming normal vectors if the matric
|
61 |
correct for transforming normal vectors if the matric
|
| 62 |
contains anisotropic scaling. */
|
62 |
contains anisotropic scaling. */
|
| 63 |
template<class T, class VecT>
|
63 |
template<class T, class VecT>
|
| 64 |
const VecT mul_3D_vector(const ArithVec3Float<T,VecT>& v_in) const
|
64 |
const VecT mul_3D_vector(const ArithVec3Float<T,VecT>& v_in) const
|
| 65 |
{
|
65 |
{
|
| 66 |
VT v_out = (*this) * VT(v_in[0],v_in[1],v_in[2],0);
|
66 |
VT v_out = (*this) * VT(v_in[0],v_in[1],v_in[2],0);
|
| 67 |
return VecT(v_out[0],v_out[1],v_out[2]);
|
67 |
return VecT(v_out[0],v_out[1],v_out[2]);
|
| 68 |
}
|
68 |
}
|
| 69 |
|
69 |
|
| 70 |
/** Multiply 3D point onto matrix. Here the fourth coordinate
|
70 |
/** Multiply 3D point onto matrix. Here the fourth coordinate
|
| 71 |
becomes 1 to ensure that the point is translated. Note that
|
71 |
becomes 1 to ensure that the point is translated. Note that
|
| 72 |
the vector is converted back into a Vec3f without any
|
72 |
the vector is converted back into a Vec3f without any
|
| 73 |
division by w. This is deliberate: Typically, w=1 except
|
73 |
division by w. This is deliberate: Typically, w=1 except
|
| 74 |
for projections. If we are doing projection, we can use
|
74 |
for projections. If we are doing projection, we can use
|
| 75 |
project_3D_point instead */
|
75 |
project_3D_point instead */
|
| 76 |
template<class T, class VecT>
|
76 |
template<class T, class VecT>
|
| 77 |
const VecT mul_3D_point(const ArithVec3Float<T,VecT> & v_in) const
|
77 |
const VecT mul_3D_point(const ArithVec3Float<T,VecT> & v_in) const
|
| 78 |
{
|
78 |
{
|
| 79 |
VT v_out = (*this) * VT(v_in[0],v_in[1],v_in[2],1);
|
79 |
VT v_out = (*this) * VT(v_in[0],v_in[1],v_in[2],1);
|
| 80 |
return VecT(v_out[0],v_out[1],v_out[2]);
|
80 |
return VecT(v_out[0],v_out[1],v_out[2]);
|
| 81 |
}
|
81 |
}
|
| 82 |
|
82 |
|
| 83 |
/** Multiply 3D point onto matrix. We set w=1 before
|
83 |
/** Multiply 3D point onto matrix. We set w=1 before
|
| 84 |
multiplication and divide by w after multiplication. */
|
84 |
multiplication and divide by w after multiplication. */
|
| 85 |
template<class T, class VecT>
|
85 |
template<class T, class VecT>
|
| 86 |
const VecT project_3D_point(const ArithVec3Float<T,VecT>& v_in) const
|
86 |
const VecT project_3D_point(const ArithVec3Float<T,VecT>& v_in) const
|
| 87 |
{
|
87 |
{
|
| 88 |
VT v_out = (*this) * VT(v_in[0],v_in[1],v_in[2],1);
|
88 |
VT v_out = (*this) * VT(v_in[0],v_in[1],v_in[2],1);
|
| 89 |
v_out.de_homogenize();
|
89 |
v_out.de_homogenize();
|
| 90 |
return VecT(v_out[0],v_out[1],v_out[2]);
|
90 |
return VecT(v_out[0],v_out[1],v_out[2]);
|
| 91 |
}
|
91 |
}
|
| 92 |
|
92 |
|
| 93 |
};
|
93 |
};
|
| 94 |
|
94 |
|
| 95 |
/** Compute the adjoint of a matrix. This is the matrix where each
|
95 |
/** Compute the adjoint of a matrix. This is the matrix where each
|
| 96 |
entry is the subdeterminant of 'in' where the row and column of
|
96 |
entry is the subdeterminant of 'in' where the row and column of
|
| 97 |
the element is removed. Use mostly to compute the inverse */
|
97 |
the element is removed. Use mostly to compute the inverse */
|
| 98 |
template<class VT, class M>
|
98 |
template<class VT, class M>
|
| 99 |
M adjoint(const ArithSqMat4x4Float<VT,M>&);
|
99 |
M adjoint(const ArithSqMat4x4Float<VT,M>&);
|
| 100 |
|
100 |
|
| 101 |
/** Compute the determinant of a 4x4 matrix. The code below is what
|
101 |
/** Compute the determinant of a 4x4 matrix. The code below is what
|
| 102 |
I found to be most robust. The original implementation used direct
|
102 |
I found to be most robust. The original implementation used direct
|
| 103 |
computation of the 3x3 sub-determinants and I also tried a direct
|
103 |
computation of the 3x3 sub-determinants and I also tried a direct
|
| 104 |
computation based on enumerating all permutations. The code below
|
104 |
computation based on enumerating all permutations. The code below
|
| 105 |
is better. */
|
105 |
is better. */
|
| 106 |
template<class V, class M>
|
106 |
template<class V, class M>
|
| 107 |
inline double determinant(const ArithSqMat4x4Float<V,M>& m)
|
107 |
inline double determinant(const ArithSqMat4x4Float<V,M>& m)
|
| 108 |
{
|
108 |
{
|
| 109 |
/* return m[0][0] * (m[1][1]*(m[2][2]*m[3][3]-m[3][2]*m[2][3])- */
|
109 |
/* return m[0][0] * (m[1][1]*(m[2][2]*m[3][3]-m[3][2]*m[2][3])- */
|
| 110 |
/* m[2][1]*(m[1][2]*m[3][3]-m[3][2]*m[1][3])+ */
|
110 |
/* m[2][1]*(m[1][2]*m[3][3]-m[3][2]*m[1][3])+ */
|
| 111 |
/* m[3][1]*(m[1][2]*m[2][3]-m[2][2]*m[1][3])) */
|
111 |
/* m[3][1]*(m[1][2]*m[2][3]-m[2][2]*m[1][3])) */
|
| 112 |
/* - m[1][0] * (m[0][1]*(m[2][2]*m[3][3]-m[3][2]*m[2][3])- */
|
112 |
/* - m[1][0] * (m[0][1]*(m[2][2]*m[3][3]-m[3][2]*m[2][3])- */
|
| 113 |
/* m[2][1]*(m[0][2]*m[3][3]-m[3][2]*m[0][3])+ */
|
113 |
/* m[2][1]*(m[0][2]*m[3][3]-m[3][2]*m[0][3])+ */
|
| 114 |
/* m[3][1]*(m[0][2]*m[2][3]-m[2][2]*m[0][3])) */
|
114 |
/* m[3][1]*(m[0][2]*m[2][3]-m[2][2]*m[0][3])) */
|
| 115 |
/* + m[2][0] * (m[0][1]*(m[1][2]*m[3][3]-m[3][2]*m[1][3])- */
|
115 |
/* + m[2][0] * (m[0][1]*(m[1][2]*m[3][3]-m[3][2]*m[1][3])- */
|
| 116 |
/* m[1][1]*(m[0][2]*m[3][3]-m[3][2]*m[0][3])+ */
|
116 |
/* m[1][1]*(m[0][2]*m[3][3]-m[3][2]*m[0][3])+ */
|
| 117 |
/* m[3][1]*(m[0][2]*m[1][3]-m[1][2]*m[0][3])) */
|
117 |
/* m[3][1]*(m[0][2]*m[1][3]-m[1][2]*m[0][3])) */
|
| 118 |
/* - m[3][0] * (m[0][1]*(m[1][2]*m[2][3]-m[2][2]*m[1][3])- */
|
118 |
/* - m[3][0] * (m[0][1]*(m[1][2]*m[2][3]-m[2][2]*m[1][3])- */
|
| 119 |
/* m[1][1]*(m[0][2]*m[2][3]-m[2][2]*m[0][3])+ */
|
119 |
/* m[1][1]*(m[0][2]*m[2][3]-m[2][2]*m[0][3])+ */
|
| 120 |
/* m[2][1]*(m[0][2]*m[1][3]-m[1][2]*m[0][3])); */
|
120 |
/* m[2][1]*(m[0][2]*m[1][3]-m[1][2]*m[0][3])); */
|
| 121 |
typedef typename M::ScalarType ScalarType;
|
121 |
typedef typename M::ScalarType ScalarType;
|
| 122 |
ScalarType a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4;
|
122 |
ScalarType a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4;
|
| 123 |
|
123 |
|
| 124 |
/* assign to individual variable names to aid selecting */
|
124 |
/* assign to individual variable names to aid selecting */
|
| 125 |
/* correct elements */
|
125 |
/* correct elements */
|
| 126 |
|
126 |
|
| 127 |
a1 = m[0][0]; b1 = m[1][0]; c1 = m[2][0]; d1 = m[3][0];
|
127 |
a1 = m[0][0]; b1 = m[1][0]; c1 = m[2][0]; d1 = m[3][0];
|
| 128 |
a2 = m[0][1]; b2 = m[1][1]; c2 = m[2][1]; d2 = m[3][1];
|
128 |
a2 = m[0][1]; b2 = m[1][1]; c2 = m[2][1]; d2 = m[3][1];
|
| 129 |
a3 = m[0][2]; b3 = m[1][2]; c3 = m[2][2]; d3 = m[3][2];
|
129 |
a3 = m[0][2]; b3 = m[1][2]; c3 = m[2][2]; d3 = m[3][2];
|
| 130 |
a4 = m[0][3]; b4 = m[1][3]; c4 = m[2][3]; d4 = m[3][3];
|
130 |
a4 = m[0][3]; b4 = m[1][3]; c4 = m[2][3]; d4 = m[3][3];
|
| 131 |
|
131 |
|
| 132 |
return
|
132 |
return
|
| 133 |
a1 * (b2*(c3*d4-d3*c4)-c2*(b3*d4-d3*b4)+d2*(b3*c4-c3*b4))
|
133 |
a1 * (b2*(c3*d4-d3*c4)-c2*(b3*d4-d3*b4)+d2*(b3*c4-c3*b4))
|
| 134 |
- b1 * (a2*(c3*d4-d3*c4)-c2*(a3*d4-d3*a4)+d2*(a3*c4-c3*a4))
|
134 |
- b1 * (a2*(c3*d4-d3*c4)-c2*(a3*d4-d3*a4)+d2*(a3*c4-c3*a4))
|
| 135 |
+ c1 * (a2*(b3*d4-d3*b4)-b2*(a3*d4-d3*a4)+d2*(a3*b4-b3*a4))
|
135 |
+ c1 * (a2*(b3*d4-d3*b4)-b2*(a3*d4-d3*a4)+d2*(a3*b4-b3*a4))
|
| 136 |
- d1 * (a2*(b3*c4-c3*b4)-b2*(a3*c4-c3*a4)+c2*(a3*b4-b3*a4));
|
136 |
- d1 * (a2*(b3*c4-c3*b4)-b2*(a3*c4-c3*a4)+c2*(a3*b4-b3*a4));
|
| 137 |
|
137 |
|
| 138 |
}
|
138 |
}
|
| 139 |
|
139 |
|
| 140 |
/// Compute the inverse matrix of a Mat4x4f.
|
140 |
/// Compute the inverse matrix of a Mat4x4f.
|
| 141 |
template<class VT, class M>
|
141 |
template<class VT, class M>
|
| 142 |
M invert(const ArithSqMat4x4Float<VT,M>&);
|
142 |
M invert(const ArithSqMat4x4Float<VT,M>&);
|
| 143 |
|
143 |
|
| 144 |
/// Compute the inverse matrix of a Mat4x4f that is affine
|
144 |
/// Compute the inverse matrix of a Mat4x4f that is affine
|
| 145 |
template<class VT, class M>
|
145 |
template<class VT, class M>
|
| 146 |
M invert_affine(const ArithSqMat4x4Float<VT,M>&);
|
146 |
M invert_affine(const ArithSqMat4x4Float<VT,M>&);
|
| 147 |
|
147 |
|
| 148 |
}
|
148 |
}
|
| 149 |
#endif
|
149 |
#endif
|
| 150 |
|
150 |
|
| 151 |
|
151 |
|
| 152 |
|
152 |
|
| 153 |
|
153 |
|
| 154 |
|
154 |
|
| 155 |
|
155 |
|
| 156 |
|
156 |
|
| 157 |
|
157 |
|