ndarray_cg/d2/mat4x4/
general.rs

1use ndarray::Dimension;
2
3use crate::*;
4
5fn minor
6< 
7  E : MatEl + nd::NdFloat, 
8  Descriptor : mat::Descriptor 
9>
10( 
11  from : &Mat4< E, Descriptor >, 
12  to : &mut Mat3< E, Descriptor >, 
13  i : usize, 
14  j : usize 
15)
16where 
17Mat4< E, Descriptor > : RawSliceMut< Scalar = E > + IndexingRef< Scalar = E, Index = Ix2 >,
18Mat3< E, Descriptor > : RawSliceMut< Scalar = E >
19{
20  for( id, ( _, v ) ) in from
21  .iter_indexed_unstable()
22  .filter( 
23    | ( id, _ ) |
24    { 
25      let ( r, c ) = id.into_pattern();
26      r != i && c != j
27    } 
28  ).enumerate()
29  {
30    to.raw_slice_mut()[ id ] = *v;
31  }
32}
33
34fn cofactor
35< 
36  E : MatEl + nd::NdFloat, 
37  Descriptor : mat::Descriptor 
38>
39( 
40  from : &Mat4< E, Descriptor >, 
41  to : &mut Mat3< E, Descriptor >,  
42  i : usize, 
43  j : usize 
44) -> E
45where 
46Mat4< E, Descriptor > : 
47  RawSliceMut< Scalar = E > + 
48  IndexingRef< Scalar = E, Index = Ix2 >,
49Mat3< E, Descriptor > : 
50  RawSliceMut< Scalar = E > + 
51  ScalarRef< Scalar = E, Index = Ix2 > + 
52  ConstLayout< Index = Ix2 > + 
53  IndexingMut< Scalar = E, Index = Ix2 >
54{
55  let k = E::from( ( -1i32 ).pow( ( i + j ) as u32 ) ).unwrap();
56  minor( from, to, i, j );
57  k * to.determinant()
58}
59
60impl< E, Descriptor > Mat< 4, 4, E, Descriptor > 
61where 
62E : MatEl + nd::NdFloat,
63Descriptor : mat::Descriptor,
64Self : ScalarMut< Scalar = E, Index = Ix2 > +
65       RawSliceMut< Scalar = E > + 
66       ConstLayout< Index = Ix2 > + 
67       IndexingMut< Scalar = E, Index = Ix2 >
68{
69  /// Converts the matrix to an array
70  pub fn to_array( &self ) -> [ E; 16 ]
71  {
72    self.raw_slice().try_into().unwrap()
73  }
74
75
76  /// Computes the determinant of the matrix
77  pub fn determinant( &self ) -> E
78  where 
79    Mat< 3, 3, E, Descriptor > : 
80      RawSliceMut< Scalar = E > +
81      ScalarMut< Scalar = E, Index = Ix2 > + 
82      ConstLayout< Index = Ix2 > + 
83      IndexingMut< Scalar = E, Index = Ix2 >
84  {
85    let _a11 = *self.scalar_ref( Ix2( 0, 0 ) );
86    let _a12 = *self.scalar_ref( Ix2( 0, 1 ) );
87    let _a13 = *self.scalar_ref( Ix2( 0, 2 ) );
88    let _a14 = *self.scalar_ref( Ix2( 0, 3 ) );
89
90    let mut m = Mat3::< E, Descriptor >::default();
91
92    minor( self, &mut m, 0, 0 );
93    let _det11 = m.determinant();
94    minor( self, &mut m, 0, 1 );
95    let _det12 = m.determinant();
96    minor( self, &mut m, 0, 2 );
97    let _det13 = m.determinant();
98    minor( self, &mut m, 0, 3 );
99    let _det14 = m.determinant();
100
101    _a11 * _det11 - _a12 * _det12 + _a13 * _det13 - _a14 * _det14
102  }
103
104  /// Computes the inverse of the matrix.
105  /// If the determinant is zero - return `None`
106  pub fn inverse( &self ) -> Option< Self >
107  where 
108    Mat< 3, 3, E, Descriptor > : 
109      RawSliceMut< Scalar = E > +
110      ScalarMut< Scalar = E, Index = Ix2 > + 
111      ConstLayout< Index = Ix2 > + 
112      IndexingMut< Scalar = E, Index = Ix2 >
113  {
114    let det = self.determinant();
115
116    if det == E::zero() { return None; }
117
118    let mut cfm = Mat3::default();
119    let mut cf = | i, j |
120    {
121      cofactor( self, &mut cfm, i, j )
122    };
123
124    let adj = Self::from_column_major
125    ([
126      cf( 0, 0 ), cf( 0, 1 ), cf( 0, 2 ), cf( 0, 3 ),
127      cf( 1, 0 ), cf( 1, 1 ), cf( 1, 2 ), cf( 1, 3 ),
128      cf( 2, 0 ), cf( 2, 1 ), cf( 2, 2 ), cf( 2, 3 ),
129      cf( 3, 0 ), cf( 3, 1 ), cf( 3, 2 ), cf( 3, 3 ),
130    ]);
131
132    Some( adj / det )
133  }
134}