put the new matrix implementation in place

This commit is contained in:
Hartmut Seichter 2019-01-18 21:45:23 +01:00
parent 55b7361717
commit 9dd862018b
6 changed files with 293 additions and 220 deletions

View file

@ -5,8 +5,8 @@ set(hdrs
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/matrix.hpp
include/pw/core/vector.hpp
# include/pw/core/quaternion.hpp
# include/pw/core/serialize.hpp
include/pw/core/image.hpp
@ -14,7 +14,7 @@ set(hdrs
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
)
@ -23,7 +23,7 @@ set(srcs
# src/buffer.cpp
src/image.cpp
src/debug.cpp
# src/mesh.cpp
src/mesh.cpp
src/core.cpp
# src/serialize.cpp
src/timer.cpp

View file

@ -32,6 +32,197 @@
namespace pw {
template <typename T, std::size_t R,std::size_t C,bool RowMajor = false>
struct matrix_ : matrixbase_<T, matrix_<T, R, C>> {
T data[R*C];
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;
}
template <typename... Arguments>
matrix_(Arguments ...values)
: data {values... }
{
static_assert(sizeof...(Arguments) == R*C,
"Incorrect number of arguments");
}
//! 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 this->size(); }
inline size_t offset(size_t r,size_t 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 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;
}
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;
}
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;
}
template<std::size_t Rs,std::size_t Cs,bool RowMajorSlice = RowMajor>
matrix_<T,Rs,Cs,RowMajorSlice> minor(std::size_t r0,std::size_t c0) const
{
matrix_<T,Rs,Cs,RowMajorSlice> m;
size_t r = 0;
for (size_t ri = 0; ri < R; ri++) {
size_t c = 0;
if (ri == r0)
continue;
for (size_t ci = 0; ci < C; ci++)
{
if (ci == c0)
continue;
m(r,c) = (*this)(ri,ci);
c++;
}
r++;
}
return m;
}
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->minor<R-1,C-1,RowMajor>(0,c).determinant();
return det;
}
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;
}
matrix_<T,C,R,RowMajor> inverse() const {
T invDet = T(1) / this->determinant();
matrix_<T,C,R,RowMajor> inv;
for (int j = 0; j < C; j++)
for (int i = 0; i < R; i++)
{
const T minorDet = this->minor<R-1,C-1,RowMajor>(j,i).determinant();
const T coFactor = ((i + j) % 2 == 1) ? -minorDet : minorDet;
inv(i, j) = invDet * coFactor;
}
return inv;
}
inline bool row_major() const {
return RowMajor;
}
inline bool square() const { return R == C; }
};
template <> inline
float matrix_<float,1,1>::determinant() const
{
return (*this)(0,0);
}
template <> inline
double matrix_<double,1,1>::determinant() const
{
return (*this)(0,0);
}
template <typename T, std::size_t R, std::size_t Ca,std::size_t Cb>
auto operator * (const matrix_<T, R, Ca>& A,
const matrix_<T, R, Cb>& B
)
{
matrix_<T,R,Cb> result; result.zero(); // zero the output
for (size_t r = 0; r < R; r++)
for (size_t c = 0; c < Cb; c++)
for (size_t iI = 0; iI < R; iI++)
result(r,c) += A(r,iI) * B(iI,c); // inner product
return result;
}
//
//
//
template <typename T> using matrix2x2_ = matrix_<T, 2, 2>;
template <typename T> using matrix3x3_ = matrix_<T, 3, 3>;
template <typename T> using matrix4x4_ = matrix_<T, 4, 4>;
using matrix2x2f = matrix_<float, 2, 2>;
using matrix2x2d = matrix_<double, 2, 2>;
using matrix2x2 = matrix_<real_t, 2, 2>;
using matrix3x3f = matrix_<float, 3, 3>;
using matrix3x3d = matrix_<double, 3, 3>;
using matrix3x3 = matrix_<real_t, 3, 3>;
using matrix4x4f = matrix_<float, 4, 4>;
using matrix4x4d = matrix_<double, 4, 4>;
using matrix4x4 = matrix_<real_t, 4, 4>;
}
#if __OLD
template <unsigned int R, unsigned int C, typename T>
class matrix_ : public matrixbase<T> {
@ -633,6 +824,7 @@ typedef matrix4x4_<float> matrix4x4f;
}
#endif
#endif

View file

@ -30,6 +30,7 @@
#include <type_traits>
#include <utility>
#include <initializer_list>
#include <cmath>
#include <iostream>
@ -54,9 +55,26 @@ struct matrixbase_ {
}
T trace() const {
return std::accumulate(std::begin(derived().data),std::end(derived().data),T());
return std::accumulate(std::begin(derived().data),std::end(derived().data),T(0));
}
inline T squared_norm() const {
return std::accumulate(std::begin(derived().data),std::end(derived().data), T(0),
[&](const T& a,const T& b){
return a + b * b;
});
}
inline T norm() const {
return std::sqrt(squared_norm());
}
inline Derived& normalize() {
(*this) /= this->norm();
return derived();
}
using iterator = T*;
using const_iterator = const T*;
iterator begin() { return &derived().data[0]; }
@ -73,226 +91,18 @@ struct matrixbase_ {
}
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(); }
};
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(); }
template <typename T, std::size_t R,std::size_t C,bool RowMajor = false>
struct matrix_ : matrixbase_<T, matrix_<T, R, C>> {
T data[R*C];
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;
}
template <typename... Arguments>
matrix_(Arguments ...values)
: data {values... }
{
static_assert(sizeof...(Arguments) == R*C,
"Incorrect number of arguments");
}
//! 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(size_t r,size_t 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 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;
}
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;
}
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;
}
template<std::size_t Rs,std::size_t Cs,bool RowMajorSlice = RowMajor>
matrix_<T,Rs,Cs,RowMajorSlice> minor(std::size_t r0,std::size_t c0) const
{
matrix_<T,Rs,Cs,RowMajorSlice> m;
size_t r = 0;
for (size_t ri = 0; ri < R; ri++) {
size_t c = 0;
if (ri == r0)
continue;
for (size_t ci = 0; ci < C; ci++)
{
if (ci == c0)
continue;
m(r,c) = (*this)(ri,ci);
c++;
}
r++;
}
return m;
}
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->minor<R-1,C-1,RowMajor>(0,c).determinant();
return det;
}
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;
}
matrix_<T,C,R,RowMajor> inverse() const {
T invDet = T(1) / this->determinant();
matrix_<T,C,R,RowMajor> inv;
for (int j = 0; j < C; j++)
for (int i = 0; i < R; i++)
{
const T minorDet = this->minor<R-1,C-1,RowMajor>(j,i).determinant();
const T coFactor = ((i + j) % 2 == 1) ? -minorDet : minorDet;
inv(i, j) = invDet * coFactor;
}
return inv;
}
inline bool row_major() const {
return RowMajor;
}
inline bool square() const { return R == C; }
};
template <> inline
float matrix_<float,1,1>::determinant() const
{
return (*this)(0,0);
}
template <> inline
double matrix_<double,1,1>::determinant() const
{
return (*this)(0,0);
}
template <typename T, std::size_t R, std::size_t Ca,std::size_t Cb>
auto operator * (const matrix_<T, R, Ca>& A,
const matrix_<T, R, Cb>& B
)
{
matrix_<T,R,Cb> result; result.zero(); // zero the output
for (size_t r = 0; r < R; r++)
for (size_t c = 0; c < Cb; c++)
for (size_t iI = 0; iI < R; iI++)
result(r,c) += A(r,iI) * B(iI,c); // inner product
return result;
}
//
//
//
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_() = default;
vector_(const derived_type& rhs) : derived_type(rhs) {}
template <typename... Arguments>
vector_(Arguments ...values)
: derived_type( { values...} )
{
static_assert(sizeof...(Arguments) == N,
"Incorrect number of arguments");
}
};
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>;
template <typename T> using matrix3x3_ = matrix_<T, 3, 3>;
template <typename T> using matrix4x4_ = matrix_<T, 4, 4>;
using matrix2x2f = matrix_<float, 2, 2>;
template <typename T> using vector2_ = vector_<T, 2>;
template <typename T> using vector3_ = vector_<T, 3>;
template <typename T> using vector4_ = vector_<T, 4>;
using vector2f = vector2_<float>;
///**

View file

@ -27,6 +27,70 @@
namespace pw {
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_() = default;
vector_(const derived_type& rhs) : derived_type(rhs) {}
template <typename... Arguments>
vector_(Arguments ...values)
: derived_type( { values...} )
{
static_assert(sizeof...(Arguments) == N,
"Incorrect number of arguments");
}
};
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 vector2_ = vector_<T, 2>;
template <typename T> using vector3_ = vector_<T, 3>;
template <typename T> using vector4_ = vector_<T, 4>;
using vector2f = vector2_<float>;
using vector2d = vector2_<double>;
using vector2 = vector2_<real_t>;
using vector3f = vector3_<float>;
using vector3d = vector3_<double>;
using vector3 = vector3_<real_t>;
using vector4f = vector4_<float>;
using vector4d = vector4_<double>;
using vector4 = vector4_<real_t>;
}
#if ___OLDSTUFF
template <unsigned int components,typename T>
class vector_ : public matrix_<components,1,T> {
public:
@ -213,6 +277,6 @@ typedef vector4_<unsigned int> vector4ui;
#endif
#endif

View file

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

View file

@ -1,7 +1,8 @@
#include <pw/core/matrixbase.hpp>
//#include <pw/core/matrix.hpp>
#include <pw/core/matrix.hpp>
#include <pw/core/vector.hpp>
//#include <pw/core/serialize.hpp>
#include <pw/core/debug.hpp>
@ -62,6 +63,12 @@ int main(int argc,char **argv) {
std::cout << "v2_f " << v2_f << std::endl;
std::cout << "v2_b " << v2_b << std::endl;
std::cout << "v2_b.norm " << v2_b.norm() << std::endl;
v2_b.normalize();
std::cout << "v2_b.normalized " << v2_b << std::endl;