[transform] use qvm instead of ublas

This commit is contained in:
barendgehrels 2017-01-04 11:55:56 +01:00
parent f486eee0f6
commit 138879adb9
5 changed files with 50 additions and 97 deletions

View File

@ -144,9 +144,7 @@ int main()
// Compose matrix for the two transformation
// Create transformer attached to the transformation matrix
ublas_transformer<double, 2, 2>
combined(boost::numeric::ublas::prod(rotate.matrix(), translate.matrix()));
//combined(rotate.matrix());
ublas_transformer<double, 2, 2> combined(rotate.matrix() * translate.matrix());
// Apply transformation to subject geometry point-by-point
model::polygon<point_2d> g4;

View File

@ -14,15 +14,8 @@
#ifndef BOOST_GEOMETRY_STRATEGIES_TRANSFORM_INVERSE_TRANSFORMER_HPP
#define BOOST_GEOMETRY_STRATEGIES_TRANSFORM_INVERSE_TRANSFORMER_HPP
// Remove the ublas checking, otherwise the inverse might fail
// (while nothing seems to be wrong)
#ifdef BOOST_UBLAS_TYPE_CHECK
#undef BOOST_UBLAS_TYPE_CHECK
#endif
#define BOOST_UBLAS_TYPE_CHECK 0
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/qvm/mat.hpp>
#include <boost/qvm/mat_operations.hpp>
#include <boost/geometry/strategies/transform/matrix_transformers.hpp>
@ -50,25 +43,7 @@ public :
template <typename Transformer>
inline inverse_transformer(Transformer const& input)
{
typedef boost::numeric::ublas::matrix<CalculationType> matrix_type;
// create a working copy of the input
matrix_type copy(input.matrix());
// create a permutation matrix for the LU-factorization
typedef boost::numeric::ublas::permutation_matrix<> permutation_matrix;
permutation_matrix pm(copy.size1());
// perform LU-factorization
int res = boost::numeric::ublas::lu_factorize<matrix_type>(copy, pm);
if( res == 0 )
{
// create identity matrix
this->m_matrix.assign(boost::numeric::ublas::identity_matrix<CalculationType>(copy.size1()));
// backsubstitute to get the inverse
boost::numeric::ublas::lu_substitute(copy, pm, this->m_matrix);
}
this->m_matrix = boost::qvm::inverse(input.matrix());
}
};

View File

@ -48,7 +48,8 @@ template
class map_transformer
: public ublas_transformer<CalculationType, Dimension1, Dimension2>
{
typedef boost::numeric::ublas::matrix<CalculationType> M;
typedef boost::qvm::mat<CalculationType, Dimension1 + 1, Dimension2 + 1> M;
typedef boost::qvm::mat<CalculationType, 3, 3> matrix33;
public :
template <typename B, typename D>
@ -76,26 +77,26 @@ private :
{
// Translate to a coordinate system centered on world coordinates (-wx, -wy)
M t1(3,3);
t1(0,0) = 1; t1(0,1) = 0; t1(0,2) = -wx;
t1(1,0) = 0; t1(1,1) = 1; t1(1,2) = -wy;
t1(2,0) = 0; t1(2,1) = 0; t1(2,2) = 1;
matrix33 t1;
qvm::A<0,0>(t1) = 1; qvm::A<0,1>(t1) = 0; qvm::A<0,2>(t1) = -wx;
qvm::A<1,0>(t1) = 0; qvm::A<1,1>(t1) = 1; qvm::A<1,2>(t1) = -wy;
qvm::A<2,0>(t1) = 0; qvm::A<2,1>(t1) = 0; qvm::A<2,2>(t1) = 1;
// Scale the map
M s(3,3);
s(0,0) = scalex; s(0,1) = 0; s(0,2) = 0;
s(1,0) = 0; s(1,1) = scaley; s(1,2) = 0;
s(2,0) = 0; s(2,1) = 0; s(2,2) = 1;
matrix33 s;
qvm::A<0,0>(s) = scalex; qvm::A<0,1>(s) = 0; qvm::A<0,2>(s) = 0;
qvm::A<1,0>(s) = 0; qvm::A<1,1>(s) = scaley; qvm::A<1,2>(s) = 0;
qvm::A<2,0>(s) = 0; qvm::A<2,1>(s) = 0; qvm::A<2,2>(s) = 1;
// Translate to a coordinate system centered on the specified pixels (+px, +py)
M t2(3, 3);
t2(0,0) = 1; t2(0,1) = 0; t2(0,2) = px;
t2(1,0) = 0; t2(1,1) = 1; t2(1,2) = py;
t2(2,0) = 0; t2(2,1) = 0; t2(2,2) = 1;
matrix33 t2;
qvm::A<0,0>(t2) = 1; qvm::A<0,1>(t2) = 0; qvm::A<0,2>(t2) = px;
qvm::A<1,0>(t2) = 0; qvm::A<1,1>(t2) = 1; qvm::A<1,2>(t2) = py;
qvm::A<2,0>(t2) = 0; qvm::A<2,1>(t2) = 0; qvm::A<2,2>(t2) = 1;
// Calculate combination matrix in two steps
this->m_matrix = boost::numeric::ublas::prod(s, t1);
this->m_matrix = boost::numeric::ublas::prod(t2, this->m_matrix);
this->m_matrix = s * t1;
this->m_matrix = t2 * this->m_matrix;
}
@ -140,20 +141,20 @@ private :
if (Mirror)
{
// Mirror in y-direction
M m(3,3);
m(0,0) = 1; m(0,1) = 0; m(0,2) = 0;
m(1,0) = 0; m(1,1) = -1; m(1,2) = 0;
m(2,0) = 0; m(2,1) = 0; m(2,2) = 1;
matrix33 m;
qvm::A<0,0>(m) = 1; qvm::A<0,1>(m) = 0; qvm::A<0,2>(m) = 0;
qvm::A<1,0>(m) = 0; qvm::A<1,1>(m) = -1; qvm::A<1,2>(m) = 0;
qvm::A<2,0>(m) = 0; qvm::A<2,1>(m) = 0; qvm::A<2,2>(m) = 1;
// Translate in y-direction such that it fits again
M y(3, 3);
y(0,0) = 1; y(0,1) = 0; y(0,2) = 0;
y(1,0) = 0; y(1,1) = 1; y(1,2) = height;
y(2,0) = 0; y(2,1) = 0; y(2,2) = 1;
matrix33 y;
qvm::A<0,0>(y) = 1; qvm::A<0,1>(y) = 0; qvm::A<0,2>(y) = 0;
qvm::A<1,0>(y) = 0; qvm::A<1,1>(y) = 1; qvm::A<1,2>(y) = height;
qvm::A<2,0>(y) = 0; qvm::A<2,1>(y) = 0; qvm::A<2,2>(y) = 1;
// Calculate combination matrix in two steps
this->m_matrix = boost::numeric::ublas::prod(m, this->m_matrix);
this->m_matrix = boost::numeric::ublas::prod(y, this->m_matrix);
this->m_matrix = m * this->m_matrix;
this->m_matrix = y * this->m_matrix;
}
}
};

View File

@ -22,27 +22,8 @@
#include <cstddef>
// Remove the ublas checking, otherwise the inverse might fail
// (while nothing seems to be wrong)
#ifdef BOOST_UBLAS_TYPE_CHECK
#undef BOOST_UBLAS_TYPE_CHECK
#endif
#define BOOST_UBLAS_TYPE_CHECK 0
#include <boost/numeric/conversion/cast.hpp>
#if defined(__clang__)
// Avoid warning about unused UBLAS function: boost_numeric_ublas_abs
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wunused-function"
#endif
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#if defined(__clang__)
#pragma clang diagnostic pop
#endif
#include <boost/qvm/mat.hpp>
#include <boost/qvm/mat_access.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_dimension.hpp>
@ -85,7 +66,7 @@ class ublas_transformer<CalculationType, 2, 2>
{
protected :
typedef CalculationType ct;
typedef boost::numeric::ublas::matrix<ct> matrix_type;
typedef boost::qvm::mat<ct, 3, 3> matrix_type;
matrix_type m_matrix;
public :
@ -94,11 +75,10 @@ public :
ct const& m_0_0, ct const& m_0_1, ct const& m_0_2,
ct const& m_1_0, ct const& m_1_1, ct const& m_1_2,
ct const& m_2_0, ct const& m_2_1, ct const& m_2_2)
: m_matrix(3, 3)
{
m_matrix(0,0) = m_0_0; m_matrix(0,1) = m_0_1; m_matrix(0,2) = m_0_2;
m_matrix(1,0) = m_1_0; m_matrix(1,1) = m_1_1; m_matrix(1,2) = m_1_2;
m_matrix(2,0) = m_2_0; m_matrix(2,1) = m_2_1; m_matrix(2,2) = m_2_2;
qvm::A<0,0>(m_matrix) = m_0_0; qvm::A<0,1>(m_matrix) = m_0_1; qvm::A<0,2>(m_matrix) = m_0_2;
qvm::A<1,0>(m_matrix) = m_1_0; qvm::A<1,1>(m_matrix) = m_1_1; qvm::A<1,2>(m_matrix) = m_1_2;
qvm::A<2,0>(m_matrix) = m_2_0; qvm::A<2,1>(m_matrix) = m_2_1; qvm::A<2,2>(m_matrix) = m_2_2;
}
inline ublas_transformer(matrix_type const& matrix)
@ -106,7 +86,7 @@ public :
{}
inline ublas_transformer() : m_matrix(3, 3) {}
inline ublas_transformer() {}
template <typename P1, typename P2>
inline bool apply(P1 const& p1, P2& p2) const
@ -117,8 +97,8 @@ public :
ct const& c1 = get<0>(p1);
ct const& c2 = get<1>(p1);
ct p2x = c1 * m_matrix(0,0) + c2 * m_matrix(0,1) + m_matrix(0,2);
ct p2y = c1 * m_matrix(1,0) + c2 * m_matrix(1,1) + m_matrix(1,2);
ct p2x = c1 * qvm::A<0,0>(m_matrix) + c2 * qvm::A<0,1>(m_matrix) + qvm::A<0,2>(m_matrix);
ct p2y = c1 * qvm::A<1,0>(m_matrix) + c2 * qvm::A<1,1>(m_matrix) + qvm::A<1,2>(m_matrix);
typedef typename geometry::coordinate_type<P2>::type ct2;
set<0>(p2, boost::numeric_cast<ct2>(p2x));
@ -159,7 +139,7 @@ class ublas_transformer<CalculationType, 3, 3>
{
protected :
typedef CalculationType ct;
typedef boost::numeric::ublas::matrix<ct> matrix_type;
typedef boost::qvm::mat<ct, 4, 4> matrix_type;
matrix_type m_matrix;
public :
@ -169,15 +149,14 @@ public :
ct const& m_2_0, ct const& m_2_1, ct const& m_2_2, ct const& m_2_3,
ct const& m_3_0, ct const& m_3_1, ct const& m_3_2, ct const& m_3_3
)
: m_matrix(4, 4)
{
m_matrix(0,0) = m_0_0; m_matrix(0,1) = m_0_1; m_matrix(0,2) = m_0_2; m_matrix(0,3) = m_0_3;
m_matrix(1,0) = m_1_0; m_matrix(1,1) = m_1_1; m_matrix(1,2) = m_1_2; m_matrix(1,3) = m_1_3;
m_matrix(2,0) = m_2_0; m_matrix(2,1) = m_2_1; m_matrix(2,2) = m_2_2; m_matrix(2,3) = m_2_3;
m_matrix(3,0) = m_3_0; m_matrix(3,1) = m_3_1; m_matrix(3,2) = m_3_2; m_matrix(3,3) = m_3_3;
qvm::A<0,0>(m_matrix) = m_0_0; qvm::A<0,1>(m_matrix) = m_0_1; qvm::A<0,2>(m_matrix) = m_0_2; qvm::A<0,3>(m_matrix) = m_0_3;
qvm::A<1,0>(m_matrix) = m_1_0; qvm::A<1,1>(m_matrix) = m_1_1; qvm::A<1,2>(m_matrix) = m_1_2; qvm::A<1,3>(m_matrix) = m_1_3;
qvm::A<2,0>(m_matrix) = m_2_0; qvm::A<2,1>(m_matrix) = m_2_1; qvm::A<2,2>(m_matrix) = m_2_2; qvm::A<2,3>(m_matrix) = m_2_3;
qvm::A<3,0>(m_matrix) = m_3_0; qvm::A<3,1>(m_matrix) = m_3_1; qvm::A<3,2>(m_matrix) = m_3_2; qvm::A<3,3>(m_matrix) = m_3_3;
}
inline ublas_transformer() : m_matrix(4, 4) {}
inline ublas_transformer() {}
template <typename P1, typename P2>
inline bool apply(P1 const& p1, P2& p2) const

View File

@ -25,11 +25,6 @@
#include <boost/test/included/unit_test.hpp>
#include <boost/assign/list_of.hpp>
#include <boost/core/ignore_unused.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/tuple/tuple.hpp>
#include <boost/geometry/core/point_type.hpp>
#include <boost/geometry/core/tags.hpp>
@ -46,6 +41,11 @@
#include <boost/geometry/io/wkt/wkt.hpp>
#include <boost/geometry/io/dsv/write.hpp>
#include <boost/assign/list_of.hpp>
#include <boost/core/ignore_unused.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/tuple/tuple.hpp>
namespace bg = ::boost::geometry;
namespace ba = ::boost::assign;