fixing core matrix calculations

TODO.md

@@ -0,0 +1,12 @@
+# core
+
+* rewrite matrix and associated code
+
+
+# scripting
+
+* refactor script into something like a runtime - script
+
+
+
+

src/CMakeLists.txt

@@ -1,14 +1,13 @@
+#add_subdirectory(deps)
 
-
-add_subdirectory(deps)
 add_subdirectory(core)
-add_subdirectory(scene)
-add_subdirectory(system)
-add_subdirectory(io)
+#add_subdirectory(scene)
+#add_subdirectory(system)
+#add_subdirectory(io)
 
 #add_subdirectory(ui)
-add_subdirectory(scripting)
-add_subdirectory(visual)
-add_subdirectory(geometry)
+#add_subdirectory(scripting)
+#add_subdirectory(visual)
+#add_subdirectory(geometry)
 
-add_subdirectory(engine)
+#add_subdirectory(engine)

src/core/CMakeLists.txt

@@ -1,35 +1,36 @@
 
 set(hdrs
     include/pw/core/debug.hpp
-    include/pw/core/axisangle.hpp
+#    include/pw/core/axisangle.hpp
     include/pw/core/core.hpp
     include/pw/core/math.hpp
     include/pw/core/matrixbase.hpp
-    include/pw/core/matrix.hpp
-    include/pw/core/vector.hpp
-    include/pw/core/quaternion.hpp
-    include/pw/core/serialize.hpp
+#    include/pw/core/matrix.hpp
+#    include/pw/core/vector.hpp
+#    include/pw/core/quaternion.hpp
+#    include/pw/core/serialize.hpp
     include/pw/core/image.hpp
     include/pw/core/point.hpp
     include/pw/core/rect.hpp
     include/pw/core/size.hpp
     include/pw/core/timer.hpp
-    include/pw/core/mesh.hpp
+#    include/pw/core/mesh.hpp
     include/pw/core/globals.hpp
     include/pw/core/image.hpp
     )
 
 set(srcs
-    src/buffer.cpp
+#    src/buffer.cpp
     src/image.cpp
     src/debug.cpp
-    src/mesh.cpp
+#    src/mesh.cpp
     src/core.cpp
-    src/serialize.cpp
+#    src/serialize.cpp
     src/timer.cpp
     src/image.cpp
     ${CMAKE_SOURCE_DIR}/README.md
     ${CMAKE_SOURCE_DIR}/LICENSE
+    ${CMAKE_SOURCE_DIR}/TODO.md
     )
 
 add_library(pwcore

src/core/include/pw/core/matrix.hpp

@@ -50,7 +50,13 @@ public:
 
     matrix_(std::initializer_list<T> il)
     {
-//        for (int i = 0;i < il.size();i++) (*this).at(i) = il[i];
+        this->set(il);
+    }
+
+    matrix_ set(std::initializer_list<T> il)
+    {
+        int i = 0;
+        for (T e : il) (*this).at(i++) = e;
     }
 
 	matrix_ transposed() const;
@@ -106,10 +112,6 @@ public:
 
 	matrix_<R,C,T>& copy_from_data(const T* src) { for (unsigned int i = 0; i < R*C; ++i) { (*this).at(i) = src[i]; } return *this; }
 
-	matrix_<R,C,T> operator * (const matrix_<R,C,T>& rhs) const {
-        return mul(rhs);
-	}
-
 	const matrix_<C,R,T> reshape() const {
 		matrix_<C,R,T> m;
 		for (unsigned int r = 0; r < R; ++r)
@@ -131,7 +133,7 @@ public:
 	void normalize();
 
     template <unsigned int bC>
-    matrix_<R,bC,T> mul(const matrix_<R,bC,T>& B) const
+    matrix_<R,bC,T> operator * (const matrix_<R,bC,T>& B) const
     {
         // aC == bR
         // set all null
@@ -235,7 +237,7 @@ inline matrix_<R,C,T> operator - (const matrix_<R,C,T>& a, const matrix_<R,C,T>&
 template <unsigned int R, unsigned int C,typename T>
 matrix_<R,C,T>& matrix_<R,C,T>::operator *= (const matrix_& rhs)
 {
-    *this = mul(rhs);
+    *this = *this * rhs;
 	return *this;
 }
 

src/core/include/pw/core/matrixbase.hpp

@@ -23,122 +23,327 @@
 #ifndef PW_CORE_MATRIXBASE_HPP
 #define PW_CORE_MATRIXBASE_HPP
 
-
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <numeric>
+#include <type_traits>
+#include <utility>
 #include <initializer_list>
 
 namespace pw {
 
-/**
- * \brief base class for all matrix and vector operations
- */
-template <typename T> class matrixbase {
-public:
+// http://coliru.stacked-crooked.com/a/34db679ee8c13a80
 
-    typedef T value_type;
-    typedef int offset_type;
-    typedef unsigned int size_type;
+template <typename T, typename Derived>
+struct matrixbase_ {
+    Derived& derived() { return static_cast<Derived&>(*this); }
+    const Derived& derived() const { return static_cast<const Derived&>(*this); }
 
-
-    matrixbase() {}
-
-
-	//! assignment constructor
-    explicit matrixbase(int rows, int cols,T* ptr,bool row_major,int data_offset = 0)
-        : _data(ptr)
-        , _data_offset(data_offset)
-        , _rows(rows)
-        , _cols(cols)
-        , _row_stride(row_major ? 1    : cols)
-        , _col_stride(row_major ? rows : 1)
-	{
+    std::size_t size() const {
+        return std::extent<decltype(Derived::data)>::value;
     }
 
-//    matrixbase(std::initializer_list<T> il)
-//    {
-//        for (int i = 0;i < il.size();i++) this->at(i) = il[i];
-//    }
+    Derived& fill(const T& v) {
+        std::fill(std::begin(derived().data), std::end(derived().data), T(v));
+        return derived();
+    }
 
-    inline const T& get_element(int e) const { return this->at(e); }
-    inline void set_element(int e,const T &v) { this->at(e) = v; }
+    Derived& zero() {
+        return derived().fill(0);
+    }
 
-    //! return number of rows
-    inline unsigned int rows() const { return _rows; }
-    //! return number of columns
-    inline unsigned int cols() const { return _cols; }
+    T trace() const {
+        return std::accumulate(std::begin(derived().data),std::end(derived().data),T());
+    }
 
-	//! get cell count
-    inline unsigned int cells() const { return _rows * _cols; }
+    using iterator = T*;
+    using const_iterator = const T*;
+    iterator begin() { return &derived().data[0]; }
+    iterator end() { return &derived().data[0] + size(); }
+    const_iterator begin() const { return &derived().data[0]; }
+    const_iterator end() const { return &derived().data[0] + size(); }
 
-    //! get data
-    inline T* data() { return _data + _data_offset; }
+    T& operator [] (std::size_t i) {
+        return derived().data[i];
+    }
 
-    //! get data
-    inline const T* data() const { return _data + _data_offset; }
+    const T& operator [] (std::size_t i) const {
+        return derived().data[i];
+    }
 
-    //! get item at index
-    inline T& at(unsigned int idx) { return data()[idx]; }
+    inline Derived& operator *= (const T& b) { for (auto & e : *this) e *= b; return derived(); }
+    inline Derived& operator /= (const T& b) { for (auto & e : *this) e *= b; return derived(); }
+    inline Derived& operator += (const T& b) { for (auto & e : *this) e *= b; return derived(); }
+    inline Derived& operator -= (const T& b) { for (auto & e : *this) e *= b; return derived(); }
+};
 
-    //! get item at index
-    inline const T& at(unsigned int idx) const { return data()[idx]; }
 
-    //! get item at position
-    inline T& at(unsigned int r,unsigned int c) { return data()[r * _row_stride + c * _col_stride]; }
+template <typename T, std::size_t R,std::size_t C,bool RowMajor = false>
+struct matrix_ : matrixbase_<T, matrix_<T, R, C>> {
 
-    //! get item at position
-    inline const T& at(unsigned int r,unsigned int c) const { return data()[r * _row_stride + c * _col_stride]; }
+    T data[R*C];
 
-    //! fill data
-    inline matrixbase& fill(const T& val) {
-        for (unsigned int r = 0; r < this->rows(); r++)
-            for (unsigned int c = 0; c < this->cols(); c++)
-                this->at(r,c) = val;
+    matrix_() = default;
+
+    matrix_(const matrix_<T,R,C,RowMajor>& other)
+    {
+        *this = other;
+    }
+
+    matrix_& operator = (const matrix_<T,R,C,RowMajor>& other)
+    {
+        for (size_t i = 0; i < other.size();i++) (*this)[i] = other[i];
         return *this;
     }
 
+    //! rows
+    inline std::size_t rows() const { return R; }
+
+    //! return number of columns
+    inline std::size_t cols() const { return C; }
+
+    //! get cell count
+    inline std::size_t coefficients() const { return R * C; }
+
+    inline size_t offset(int r,int c) const {
+        return (RowMajor) ? r * C + c : c * R + r;
+    }
+
+    inline T& operator () (std::size_t r, std::size_t c) {
+        return data[offset(r,c)];
+    }
+
+    inline const T& operator () (std::size_t r, std::size_t c) const {
+        return data[offset(r,c)];
+    }
+
     //! set identity
-    inline matrixbase& set_identity() {
+    inline matrix_& set_identity()
+    {
         for (unsigned int r = 0;r < rows(); r++)
             for (unsigned int c = 0; c < cols(); c++)
                 this->at(r,c) = (c == r) ? T(1) : T(0);
         return *this;
     }
 
-    //! set strides for row and column
-    inline void set_stride(int row_stride,int col_stride) {
-        _row_stride = row_stride; _col_stride = col_stride;
+    template<std::size_t Rs,std::size_t Cs,bool RowMajorSlice = RowMajor>
+    matrix_<T,Rs,Cs,RowMajorSlice> slice(std::size_t r,std::size_t c) const
+    {
+        matrix_<T,Rs,Cs,RowMajorSlice> s;
+        for (std::size_t ri = 0;ri < Rs;ri++)
+            for (std::size_t ci = 0;ci < Cs;ci++)
+                s(ri,ci) = (*this)(ri+r,ci+c);
+        return s;
     }
 
-    //! return row stride (elements between consequtive rows)
-    inline int row_stride() const { return _row_stride; }
+    template<std::size_t Rs,std::size_t Cs,bool RowMajorSlice = RowMajor>
+    matrix_& set_slice(const matrix_<T,Rs,Cs,RowMajorSlice>& s,std::size_t r,std::size_t c)
+    {
+        for (std::size_t ri = 0;ri < Rs;ri++)
+            for (std::size_t ci = 0;ci < Cs;ci++)
+                (*this)(ri+r,ci+c) = s(ri,ci);
+        return *this;
+    }
 
-    //! return column stride (elements between consequtive columns)
-    inline int col_stride() const { return _col_stride; }
+    T determinant() const {
+        T det(0);
+        for (size_t c = 0; c < C; c++) {
+            det += ((c % 2 == 0) ? (*this)(0,c) : -(*this)(0,c)) *  this->slice<R-1,C-1>(0,c).determinant();
+        }
+        return det;
+    }
 
-    //! get row major
-    inline bool is_row_major() const { return row_stride() > col_stride(); }
+    matrix_<T,C,R,RowMajor> transposed() const {
+        matrix_<T,C,R,RowMajor> res;
+        for (size_t r = this->rows();r-->0;)
+            for (size_t c = this->cols();c-->0;)
+                res(c,r) = (*this)(r,c);
+        return res;
+    }
 
-    //! check if squared
-    inline bool is_square() const { return (_rows == _cols); }
+    matrix_ inverse() const
+    {
+        matrix_ resMat;
 
-    //! flip row/column order
-    inline void flip_order() { set_stride(col_stride(),row_stride()); }
+        for ( unsigned int r = 0; r < C; ++r) {
+            for ( unsigned int j = 0; j < R; ++j) {
+                short sgn = ((r+j)%2) ? -1 : 1;
+                auto minor = this->slice<R-1,C-1>(r,j);
+                resMat(r,j) = minor.determinant() * sgn;
+            }
+        }
 
-    //! get the offset for data in the matrix
-    inline int get_data_offset(int r,int c) const { return (r * _row_stride + c * _col_stride) + _data_offset; }
+        resMat = resMat.transposed();
+        resMat *= T(1)/this->determinant();
+        return resMat;
+    }
 
-protected:
+    inline bool row_major() const {
+        return RowMajor;
+    }
 
-    T* _data;
-    int _data_offset;
-
-    int _cols;
-    int _rows;
-
-    int _row_stride;
-    int _col_stride;
+    inline bool square() const { return R == C; }
 
 };
 
+template <> inline
+float matrix_<float,1,1>::determinant() const
+{
+    return (*this)[0];
+}
+
+template <> inline
+double matrix_<double,1,1>::determinant() const
+{
+    return (*this)[0];
+}
+
+
+template <typename T, std::size_t N,bool RowMajor = false>
+struct vector_ : matrix_<T, N, 1, RowMajor>
+{
+    typedef matrix_<T, N, 1, RowMajor> derived_type;
+
+    vector_(const derived_type& rhs) : derived_type(rhs) {}
+
+};
+
+template <typename T, typename U, std::size_t N,bool RowMajor = false>
+auto operator * (const vector_<T, N, RowMajor>& a, const vector_<U, N, RowMajor>& b)
+-> vector_<decltype(a[0] * b[0]), N, RowMajor> {
+    vector_<decltype(a[0] * b[0]), N, RowMajor> result;
+    for (std::size_t i = 0; i < N; ++i) {
+        result[i] = a[i] * b[i];
+    }
+    return result;
+}
+
+template <typename T, typename U, std::size_t N,bool RowMajor = false>
+auto dot(const vector_<T, N, RowMajor>& a, const vector_<U, N, RowMajor>& b)
+-> decltype(a[0] * b[0]) {
+    auto product = a * b;
+    using V = decltype(product.x);
+    return std::accumulate(std::begin(product), std::end(product), V(0));
+}
+
+
+//
+//
+//
+
+template <typename T> using matrix2x2_ = matrix_<T, 2, 2>;
+using matrix2x2f = matrix_<float, 2, 2>;
+
+template <typename T> using vector2_ = vector_<T, 2>;
+using vector2f = vector2_<float>;
+
+
+///**
+// * \brief base class for all matrix and vector operations
+// */
+//template <typename T> class matrixbase {
+//public:
+
+//    typedef T value_type;
+//    typedef int offset_type;
+//    typedef unsigned int size_type;
+
+
+//    matrixbase() {}
+
+//	//! assignment constructor
+//    explicit matrixbase(int rows, int cols,T* ptr,bool row_major,int data_offset = 0)
+//        : _data(ptr)
+//        , _data_offset(data_offset)
+//        , _rows(rows)
+//        , _cols(cols)
+//        , _row_stride(row_major ? 1    : cols)
+//        , _col_stride(row_major ? rows : 1)
+//	{
+//	}
+
+//    inline const T& get_element(int e) const { return this->at(e); }
+//    inline void set_element(int e,const T &v) { this->at(e) = v; }
+
+//    //! return number of rows
+//    inline unsigned int rows() const { return _rows; }
+//    //! return number of columns
+//    inline unsigned int cols() const { return _cols; }
+
+//	//! get cell count
+//    inline unsigned int cells() const { return _rows * _cols; }
+
+//    //! get data
+//    inline T* data() { return _data + _data_offset; }
+
+//    //! get data
+//    inline const T* data() const { return _data + _data_offset; }
+
+//    //! get item at index
+//    inline T& at(unsigned int idx) { return data()[idx]; }
+
+//    //! get item at index
+//    inline const T& at(unsigned int idx) const { return data()[idx]; }
+
+//    //! get item at position
+//    inline T& at(unsigned int r,unsigned int c) { return data()[r * _row_stride + c * _col_stride]; }
+
+//    //! get item at position
+//    inline const T& at(unsigned int r,unsigned int c) const { return data()[r * _row_stride + c * _col_stride]; }
+
+//    //! fill data
+//    inline matrixbase& fill(const T& val) {
+//        for (unsigned int r = 0; r < this->rows(); r++)
+//            for (unsigned int c = 0; c < this->cols(); c++)
+//                this->at(r,c) = val;
+//        return *this;
+//    }
+
+//    //! set identity
+//    inline matrixbase& set_identity() {
+//        for (unsigned int r = 0;r < rows(); r++)
+//            for (unsigned int c = 0; c < cols(); c++)
+//                this->at(r,c) = (c == r) ? T(1) : T(0);
+//		return *this;
+//	}
+
+//    //! set strides for row and column
+//    inline void set_stride(int row_stride,int col_stride) {
+//        _row_stride = row_stride; _col_stride = col_stride;
+//    }
+
+//    //! return row stride (elements between consequtive rows)
+//    inline int row_stride() const { return _row_stride; }
+
+//    //! return column stride (elements between consequtive columns)
+//    inline int col_stride() const { return _col_stride; }
+
+//    //! get row major
+//    inline bool is_row_major() const { return row_stride() > col_stride(); }
+
+//    //! check if squared
+//    inline bool is_square() const { return (_rows == _cols); }
+
+//    //! flip row/column order
+//    inline void flip_order() { set_stride(col_stride(),row_stride()); }
+
+//    //! get the offset for data in the matrix
+//    inline int get_data_offset(int r,int c) const { return (r * _row_stride + c * _col_stride) + _data_offset; }
+
+//protected:
+
+//    T* _data;
+//    int _data_offset;
+
+//    int _cols;
+//    int _rows;
+
+//    int _row_stride;
+//    int _col_stride;
+
+//};
+
 }
 
 

src/core/include/pw/core/vector.hpp

@@ -57,6 +57,15 @@ public:
         return acos( nothr.dot(nself) );
     }
 
+
+    vector_<components+1,T> homogenous(const T& w) const {
+        return vector_<components+1,T>();
+    }
+
+    const vector_<2,T> project() const {
+        return vector_<2,T>(x()/z(),y()/z());
+    }
+
 };
 
 
@@ -98,13 +107,6 @@ public:
     T& z() { return (*this)(2); }
 
 
-    vector_<4,T> project(const T& w) const {
-        return vector_<4,T>({ x(),y(),z(),w } );
-    }
-
-    const vector_<2,T> unproject() const {
-        return vector_<2,T>(x()/z(),y()/z());
-    }
 
 	inline static vector3_<T> forward() { return vector3_<T> ( 0, 0,-1); }
 	inline static vector3_<T> backward() { return vector3_<T>( 0, 0, 1); }

src/core/src/core.cpp

@@ -6,16 +6,15 @@
 
 
 #include "pw/core/matrixbase.hpp"
-#include "pw/core/matrix.hpp"
+//#include "pw/core/matrix.hpp"
 
 
 void test_matrixbase() {
-
     using namespace pw;
 
-    matrix4x4f m;
+//    matrix4x4f m;
 
-    m.set_identity();
+//    m.set_identity();
 
 }
 

src/core/src/mesh.cpp

@@ -7,10 +7,10 @@ void mesh::apply(const matrix4x4 &m)
 {
     for (auto &v : _vertices)
     {
-//        v = m * v.project(1);
+//        v = vector4(m * v.project(1)).un_project();
 
                 auto vh = v.project(1);
-        m.mul(vh);
+        //        m.mul(vh);
     }
 }
 

src/core/tests/CMakeLists.txt

@@ -5,25 +5,25 @@ add_executable(pwcore_test_matrix
 target_link_libraries(pwcore_test_matrix
 	pwcore)
 
-add_executable(pwcore_test_vector
-	pwcore_test_vector.cpp
-	)
+#add_executable(pwcore_test_vector
+#	pwcore_test_vector.cpp
+#	)
 
-target_link_libraries(pwcore_test_vector
-	pwcore)
+#target_link_libraries(pwcore_test_vector
+#	pwcore)
 
 
-add_executable(pwcore_test_quaternion
-	pwcore_test_quaternion.cpp
-	)
+#add_executable(pwcore_test_quaternion
+#	pwcore_test_quaternion.cpp
+#	)
 
-target_link_libraries(pwcore_test_quaternion
-	pwcore)
+#target_link_libraries(pwcore_test_quaternion
+#	pwcore)
 
 
-add_executable(pwcore_test_axisangle
-	pwcore_test_axisangle.cpp
-	)
+#add_executable(pwcore_test_axisangle
+#	pwcore_test_axisangle.cpp
+#	)
 
-target_link_libraries(pwcore_test_axisangle
-	pwcore)
+#target_link_libraries(pwcore_test_axisangle
+#	pwcore)

src/core/tests/pwcore_test_matrix.cpp

@@ -1,11 +1,45 @@
 
-#include <pw/core/matrix.hpp>
-#include <pw/core/serialize.hpp>
+#include <pw/core/matrixbase.hpp>
+
+//#include <pw/core/matrix.hpp>
+//#include <pw/core/serialize.hpp>
+
+#include <pw/core/debug.hpp>
 
 #include <iostream>
 
 int main(int argc,char **argv) {
 
+
+    using namespace pw;
+
+    matrix2x2f m22; m22.zero();
+
+    m22(0,0) = 1; m22(0,1) = 2;
+    m22(1,0) = 3; m22(1,1) = 4;
+
+
+    debug::d() << "offset(0,1) col-major " << m22.offset(0,1);
+    debug::d() << "det " << m22.determinant();
+
+//    vector2f v2 = m22.slice<2,1>(0,0);
+//    m22.set_slice<2,1>(v2,0,0);
+
+//    v2 *= 2;
+//    v2 += 1;
+
+//    m22 *= 3;
+
+//    m22.determinant();
+
+//    debug::d() << sizeof(v2);
+
+//    auto m22_i = m22.inverse();
+
+
+
+#if 0
+
     pw::matrix4x4d m;
 
     m.set_identity();
@@ -48,4 +82,6 @@ int main(int argc,char **argv) {
 
 
     return 0;
+
+#endif
 }