Subversion Repositories gelsvn

Rev

Rev 113 | Rev 327 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 113 Rev 324
Line 1... Line 1...
1 
#include "ArithSqMat4x4Float.h"
 1 
#include "ArithSqMat4x4Float.h"
2 
#include "Mat4x4f.h"
 2 
#include "Mat4x4f.h"
3 
#include "Mat4x4d.h"
 3 
#include "Mat4x4d.h"
4 
 
 4 
 
- 
 
 5 
namespace CGLA {
- 
 
 6
5 
namespace
 7 
namespace
6 
{
 8 
{
7
9
8 
	/* Aux func. computes 3x3 determinants. */
 10 
	/* Aux func. computes 3x3 determinants. */
9 
	template<class T>
 11 
	template<class T>
Line 18... Line 20...
18 
			a0*b2*c1 -
 20 
			a0*b2*c1 -
19 
			a1*b0*c2 ;
 21 
			a1*b0*c2 ;
20 
	}
 22 
	}
21 
}
 23 
}
22 
 
 24 
 
23 
 
 - 
 
- 
 
 25 
		/** Compute the determinant of a 4x4 matrix. The code below is what
- 
 
 26 
				I found to be most robust. The original implementation used direct
- 
 
 27 
				computation of the 3x3 sub-determinants and I also tried a direct
- 
 
 28 
				computation based on enumerating all permutations. The code below
24 
namespace CGLA {
 29 
				is better. */
25
- 
 
26 
 
 - 
 
27
- 
 
28 
	template<class V, class M>
 30 
	template<class V, class M>
29 
	double
31 
	double
30 
	determinant(const ArithSqMat4x4Float<V,M>& m)
 32 
	determinant(const ArithSqMat4x4Float<V,M>& m)
31 
	{
 33 
	{
32 
		typedef typename M::ScalarType ScalarType;
 34 
		typedef typename M::ScalarType ScalarType;
Line 34... Line 36...
34 
    ScalarType a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4;
 36 
    ScalarType a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4;
35
37
36 
    /* assign to individual variable names to aid selecting */
 38 
    /* assign to individual variable names to aid selecting */
37 
		/*  correct elements */
 39 
		/*  correct elements */
38
40
39 
		a1 = m[0][0]; b1 = m[0][1];
41 
		a1 = m[0][0]; b1 = m[1][0]; c1 = m[2][0]; d1 = m[3][0];
40 
		c1 = m[0][2]; d1 = m[0][3];
 - 
 
41
- 
 
42 
		a2 = m[1][0]; b2 = m[1][1];
42 
		a2 = m[0][1]; b2 = m[1][1]; c2 = m[2][1]; d2 = m[3][1];
43 
		c2 = m[1][2]; d2 = m[1][3];
 - 
 
44
- 
 
45 
		a3 = m[2][0]; b3 = m[2][1];
43 
		a3 = m[0][2]; b3 = m[1][2]; c3 = m[2][2]; d3 = m[3][2];
46 
		c3 = m[2][2]; d3 = m[2][3];
 - 
 
47
- 
 
48 
		a4 = m[3][0]; b4 = m[3][1];
- 
 
49 
		c4 = m[3][2]; d4 = m[3][3];
 44 
		a4 = m[0][3]; b4 = m[1][3]; c4 = m[2][3]; d4 = m[3][3];
50 
 
 45 
 
51 
		ans =
46 
		ans =
52 
			a1 * d3x3f( b2, b3, b4, c2, c3, c4, d2, d3, d4)
 47 
				a1 * (b2*(c3*d4-d3*c4)-c2*(b3*d4-d3*b4)+d2*(b3*c4-c3*b4))
53 
			- b1 * d3x3f( a2, a3, a4, c2, c3, c4, d2, d3, d4)
 48 
				- b1 * (a2*(c3*d4-d3*c4)-c2*(a3*d4-d3*a4)+d2*(a3*c4-c3*a4))
54 
			+ c1 * d3x3f( a2, a3, a4, b2, b3, b4, d2, d3, d4)
 49 
				+ c1 * (a2*(b3*d4-d3*b4)-b2*(a3*d4-d3*a4)+d2*(a3*b4-b3*a4))
55 
			- d1 * d3x3f( a2, a3, a4, b2, b3, b4, c2, c3, c4);
 50 
				- d1 * (a2*(b3*c4-c3*b4)-b2*(a3*c4-c3*a4)+c2*(a3*b4-b3*a4));
56 
		return ans;
 51 
		return ans;
57 
	}
 52 
	}
58 
 
 53 
 
59 
 
 54 
 
60 
	template<class V, class M>
 55 
	template<class V, class M>