Subversion Repositories gelsvn

Rev

Rev 325 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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