Adding TMatrix<N,M> templated matrix class

radixmath/CMakeLists.txt

@@ -18,6 +18,7 @@ vector3d.cc
 SET(HEADERS
 constants.hh
 matrix.hh
+matrix.i.hh
 normal.hh
 point3d.hh
 ray.hh

radixmath/matrix.hh

 #ifndef RADIX_RADIXMATH_MATRIX_H
 #define RADIX_RADIXMATH_MATRIX_H
 #include <string>
+#include <array>
 #include "radixmath/constants.hh"
 namespace radix
 {
 
+typedef unsigned int tmatrix_index_t;
+/**
+ * @class TMatrix
+ * @brief TMatix<N rows, M column> Templated N by M matrix
+ */
+template<tmatrix_index_t N, tmatrix_index_t M>
+class TMatrix
+{
+public:
+    typedef tmatrix_index_t index_t;
+    typedef Real value_t;
+    typedef std::array<value_t, M> value_array_t;
+    typedef std::array<value_array_t, N> matrix_t;
+private:
+    matrix_t mData;
+public:
+
+    TMatrix();
+    bool equals(const TMatrix<N,M>& m, Real eps=kEpsilon) const;
+    value_array_t& operator[](size_t i) {return mData[i];}
+    const value_array_t& operator[](size_t i) const {return mData[i];}
+    template<index_t J>
+    TMatrix<N,J> operator*(const TMatrix<M,J>& m) const;
+    matrix_t& data() {return mData;}
+    const matrix_t& data() const {return mData;}
+}; // class TMatrix
+
     /**
      * @brief Matrix class for affine transformations
      */
@@ -83,4 +111,8 @@ namespace radix
 
     };
 } // namespace radix
+
+//
+// Include templated implementation
+#include "radixmath/matrix.i.hh"
 #endif

radixmath/matrix.i.hh

+#include <algorithm>
+#include <cmath>
+#include "radixmath/matrix.hh"
+namespace radix
+{
+
+template<tmatrix_index_t N, tmatrix_index_t M>
+TMatrix<N,M>::TMatrix()
+{
+    for(index_t i = 0; i < mData.size(); ++i)
+    {
+        std::fill(mData[i].begin(), mData[i].end(), 0.0);
+    }
+}
+
+template<tmatrix_index_t N, tmatrix_index_t M>
+bool TMatrix<N,M>::equals(const TMatrix<N,M>& m, Real eps) const
+{
+    //
+    // check data with eps
+    for(index_t i = 0; i < N; ++i)
+    {
+        for(index_t j = 0; j < M; ++j)
+        {
+            if(std::abs(mData[i][j]-m[i][j]) > eps)
+            {
+                return false;
+            }
+        }
+    }
+    return true;
+}
+template<tmatrix_index_t N, tmatrix_index_t M>
+template<tmatrix_index_t J>
+TMatrix<N,J> TMatrix<N,M>::operator*(const TMatrix<M,J>& m) const
+{
+    // result matrix
+    TMatrix<N,J> r;
+    value_t sum = 0.0;
+    for(index_t ni = 0; ni < N; ++ni)
+    {
+        for(index_t ji = 0; ji < J; ++ji)
+        {
+            for(index_t mi = 0; mi < M; ++mi)
+            {
+                sum += mData[ni][mi] * m[mi][ji];
+            }
+            r[ni][ji] = sum;
+            sum = 0.0;
+        }
+    }
+    return r;
+}
+
+} // namespace radix

radixmath/tests/tstMatrix.cc

@@ -14,4 +14,68 @@ TEST(radix, MatrixDivision){
         EXPECT_EQ(matrix.matrix[i][i], 0.5);
     }
 }
+TEST(radix, TMatrix)
+{
+    {
+        TMatrix<3,3> m1;
+        TMatrix<3,3> m2;
+        EXPECT_TRUE(m1.equals(m2));
+        // this doesn't compile since dimensions are different
+        //    TMatrix<4,3> m3;
+        //    EXPECT_TRUE(m3.equals(m1));
+    }
+
+    {
+        TMatrix<3,3> m1;
+        TMatrix<3,1> m2;
+        // multiplication should result in 3,1
+        TMatrix<3,1> m3 = m1*m2;
+        EXPECT_EQ(m3.data().size(), 3);
+        EXPECT_EQ(m3[0].size(), 1);
+    }
+
+    {
+        TMatrix<3,3> m1;
+        //
+        // initialize to identiy matrix
+        m1.data() = {{{1,0,0}
+                     ,{0,1,0}
+                     ,{0,0,1}
+                    }};
+        TMatrix<3,1> m2;
+        m2.data() = {{{0.964221}
+                     ,{1}
+                     ,{0.825213}
+                    }};
+        // multiplication should result in 3,1
+        TMatrix<3,1> m3 = m1*m2;
+        EXPECT_TRUE(m3.equals(m2));
+    }
+
+    {
+        TMatrix<3,3> m1;
+        //
+        // initialize to identiy matrix
+        m1.data() = {{{5,8,-4}
+                     ,{6,9,-5}
+                     ,{4,7,-2}
+                    }};
+        TMatrix<3,1> m2;
+        m2.data() = {{{2}
+                     ,{-3}
+                     ,{1}
+                    }};
+
+        TMatrix<3,1> solution;
+        solution.data() = {{{-18}
+                           ,{-20}
+                           ,{-15}
+                          }};
+        // multiplication should result in 3,1
+        TMatrix<3,1> m3 = m1*m2;
+        EXPECT_TRUE(m3.equals(solution));
+        EXPECT_TRUE(solution.equals(m3));
+    }
+
+}