[Olena-patches] 3691: Improve vector and matrix interoperability.

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; Improve vector and matrix interoperability. * mln/algebra/mat.hh (t): New method to transpose; replace... (trans): Remove procedure. (operator*): New overloads to improve interoperatibility between vectors and matrices. * mln/registration/get_rot.hh: Update. * mln/registration/icp.hh: Update. * mln/algebra/vec.hh (mat<n, 1, U>): New conversion op, ctor, and assignment; declarations only, definitions are... * mln/algebra/mat.hh (mat<n, 1, U>): ...new. * tests/algebra/op_times.cc: New. * mln/labeling/compute.hh: Disable check to run on int_u. * mln/accu/sum.hh (take_as_init): New. * mln/accu/image/take_as_init.hh: New overload for image as arg. * tests/accu/image/take_as_init.cc: Augment. mln/accu/image/take_as_init.hh | 110 ++++++++++++++++++++++++++- mln/accu/sum.hh | 12 ++ mln/algebra/mat.hh | 158 ++++++++++++++++++++++++++++++++++----- mln/algebra/vec.hh | 18 ++++ mln/labeling/compute.hh | 2 mln/registration/get_rot.hh | 5 - mln/registration/icp.hh | 4 tests/accu/image/take_as_init.cc | 3 tests/algebra/op_times.cc | 72 +++++++++++++++++ 9 files changed, 357 insertions(+), 27 deletions(-) Index: mln/accu/image/take_as_init.hh --- mln/accu/image/take_as_init.hh (revision 3690) +++ mln/accu/image/take_as_init.hh (working copy) @@ -50,6 +50,11 @@ take_as_init(Image<I>& input, const mln_deduce(I, value, argument)& v); + template <typename I, typename J> + void + take_as_init(Image<I>& input, + const Image<J>& values); + # ifndef MLN_INCLUDE_ONLY @@ -57,7 +62,7 @@ namespace impl { - // Generic version. + // Generic versions. namespace generic { @@ -79,10 +84,29 @@ trace::exiting("accu::impl::image::generic::take_as_init"); } + template <typename I, typename J> + void + take_as_init(Image<I>& input_, + const Image<J>& values_) + { + trace::entering("accu::impl::image::generic::take_as_init"); + + I& input = exact(input_); + const I& values = exact(values_); + mln_precondition(input.is_valid()); + mln_precondition(values.is_valid()); + + mln_piter(I) p(input.domain()); + for_all(p) + input(p).take_as_init(values(p)); + + trace::exiting("accu::impl::image::generic::take_as_init"); + } + } // end of namespace mln::accu::image::impl::generic - // Fastest version. + // Fastest versions. template <typename I> void @@ -101,6 +125,26 @@ trace::exiting("accu::impl::image::take_as_init_fastest"); } + template <typename I, typename J> + void + take_as_init_fastest(Image<I>& input_, + const Image<J>& values_) + { + trace::entering("accu::impl::image::take_as_init_fastest"); + + I& input = exact(input_); + const I& values = exact(values_); + mln_precondition(input.is_valid()); + mln_precondition(values.is_valid()); + + mln_pixter(I) p_in(input); + mln_pixter(const J) p_v(values); + for_all_2(p_in, p_v) + p_in.val().take_as_init(p_v.val()); + + trace::exiting("accu::impl::image::take_as_init_fastest"); + } + } // end of namespace mln::accu::image::impl @@ -110,7 +154,10 @@ namespace internal { + // With a single value. + template <typename I, typename V> + inline void take_as_init_dispatch(trait::image::speed::any, Image<I>& input, const V& v) @@ -119,6 +166,7 @@ } template <typename I, typename V> + inline void take_as_init_dispatch(trait::image::speed::fastest, Image<I>& input, const V& v) @@ -127,6 +175,7 @@ } template <typename I, typename V> + inline void take_as_init_dispatch(Image<I>& input, const V& v) { @@ -134,12 +183,46 @@ input, v); } + // With an image of values. + + template <typename I, typename J> + inline + void + take_as_init_dispatch(trait::image::speed::any, + trait::image::speed::any, + Image<I>& input, const Image<J>& values) + { + impl::generic::take_as_init(input, values); + } + + template <typename I, typename J> + inline + void + take_as_init_dispatch(trait::image::speed::fastest, + trait::image::speed::fastest, + Image<I>& input, const Image<J>& values) + { + impl::take_as_init_fastest(input, values); + } + + template <typename I, typename J> + inline + void + take_as_init_dispatch(Image<I>& input, + const Image<J>& values) + { + take_as_init_dispatch(mln_trait_image_speed(I)(), + mln_trait_image_speed(J)(), + input, values); + } + } // end of namespace mln::accu::image::internal - // Facade. + // Facades. template <typename I> + inline void take_as_init(Image<I>& input, const mln_deduce(I, value, argument)& v) @@ -154,6 +237,27 @@ trace::exiting("accu::image::take_as_init"); } + template <typename I, typename J> + inline + void + take_as_init(Image<I>& input, + const Image<J>& values) + { + trace::entering("accu::image::take_as_init"); + + typedef mln_value(I) A; + mlc_is_a(A, Accumulator)::check(); + mlc_converts_to(mln_value(J), mln_argument(A))::check(); + + mln_precondition(exact(input).is_valid()); + mln_precondition(exact(values).is_valid()); + // mln_precondition(exact(values).domain() == exact(input).domain()); + + internal::take_as_init_dispatch(input, values); + + trace::exiting("accu::image::take_as_init"); + } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::accu::image Index: mln/accu/sum.hh --- mln/accu/sum.hh (revision 3690) +++ mln/accu/sum.hh (working copy) @@ -1,5 +1,5 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory -// (LRDE) +// Copyright (C) 2007, 2008, 2009 EPITA Research and Development +// Laboratory (LRDE) // // This file is part of the Olena Library. This library is free // software; you can redistribute it and/or modify it under the terms @@ -67,6 +67,7 @@ /// \{ void init(); void take(const argument& t); + void take_as_init(const argument& t); void take(const sum<T,S>& other); /// \} @@ -128,6 +129,13 @@ template <typename T, typename S> inline + void sum<T,S>::take_as_init(const argument& t) + { + s_ = static_cast<S>(t); + } + + template <typename T, typename S> + inline void sum<T,S>::take(const sum<T,S>& other) { Index: mln/algebra/mat.hh --- mln/algebra/mat.hh (revision 3690) +++ mln/algebra/mat.hh (working copy) @@ -82,6 +82,7 @@ namespace algebra { + template <unsigned n, unsigned m, typename T> class mat : public Object< mat<n,m,T> > { @@ -118,6 +119,9 @@ static mat identity(); + /// Return the transpose of the matrix. + mat<m,n,T> t() const; + private: T data_[n][m]; }; @@ -229,6 +233,11 @@ mat<n, m, mln_sum_product(T,U)> operator*(const mat<n,o,T>& lhs, const mat<o,m,U>& rhs); + template <unsigned o, typename T, + typename U> + mln_sum_product(T,U) + operator*(const mat<1,o,T>& lhs, const mat<o,1,U>& rhs); + // mat * vec template <unsigned n, unsigned m, typename T, @@ -236,6 +245,17 @@ vec<n, mln_sum_product(T,U)> operator*(const mat<n,m,T>& lhs, const vec<m,U>& rhs); + template <unsigned m, typename T, + typename U> + mln_sum_product(T,U) // scalar + operator*(const mat<1,m,T>& lhs, const vec<m,U>& rhs); + + // vec * mat + + template <unsigned m, typename T, typename U> + mat<m, m, mln_trait_op_times(T,U)> + operator*(const vec<m,T>& lhs, const mat<1,m,U>& rhs); + // mat * s template <unsigned n, unsigned m, typename T, @@ -255,11 +275,6 @@ std::ostream& operator<<(std::ostream& ostr, const mat<n,m,T>& v); - // transpose - - template<unsigned n, unsigned m, typename T> - mat<m,n,T> - trans(const mat<n,m,T>& matrix); // trace @@ -267,14 +282,61 @@ float tr(const mat<n,n,T>& m); + # ifndef MLN_INCLUDE_ONLY + + // vec -> mat + + template <unsigned n, typename T> + template <typename U> + inline + vec<n,T>::operator mat<n,1,U>() const + { + mlc_converts_to(T, U)::check(); + mat<n,1,U> tmp; + for (unsigned i = 0; i < n; ++i) + tmp(i, 0) = data_[i]; + return tmp; + } + + + // mat -> vec + + template <unsigned n, typename T> + template <typename U> + inline + vec<n,T>::vec(const mat<n, 1, U>& rhs) + { + mlc_converts_to(T, U)::check(); + for (unsigned i = 0; i < n; ++i) + data_[i] = rhs(i, 0); + } + + template <unsigned n, typename T> + template <typename U> + inline + vec<n,T>& + vec<n,T>::operator=(const mat<n, 1, U>& rhs) + { + mlc_converts_to(T, U)::check(); + for (unsigned i = 0; i < n; ++i) + data_[i] = rhs(i, 0); + return *this; + } + + + + // Id + template <unsigned n, unsigned m, typename T> - const mat<n,m,T> mat<n,m,T>::Id = mat<n,m,T>::identity(); + const mat<n,m,T> + mat<n,m,T>::Id = mat<n,m,T>::identity(); template <unsigned n, unsigned m, typename T> inline - mat<n,m,T> mat<n,m,T>::identity() + mat<n,m,T> + mat<n,m,T>::identity() { static mat<n,m,T> id_; static bool flower = true; @@ -369,6 +431,18 @@ return n * m; } + template <unsigned n, unsigned m, typename T> + inline + mat<m,n,T> + mat<n,m,T>::t() const + { + mat<m,n,T> tmp; + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + tmp(j,i) = data_[i][j]; + return tmp; + } + // Operators. @@ -421,6 +495,8 @@ return tmp; } + // mat * mat + template <unsigned n, unsigned o, typename T, unsigned m, typename U> inline @@ -438,6 +514,20 @@ return tmp; } + template <unsigned o, typename T, + typename U> + inline + mln_sum_product(T,U) + operator*(const mat<1,o,T>& lhs, const mat<o,1,U>& rhs) + { + mln_sum_product(T,U) tmp(literal::zero); + for (unsigned k = 0; k < o; ++k) + tmp += lhs(0, k) * rhs(k, 0); + return tmp; + } + + // mat * vec + template <unsigned n, unsigned m, typename T, typename U> inline @@ -455,6 +545,35 @@ return tmp; } + template <unsigned m, typename T, + typename U> + inline + mln_sum_product(T,U) // scalar + operator*(const mat<1,m,T>& lhs, const vec<m,U>& rhs) + { + mln_sum_product(T,U) tmp(literal::zero); + for (unsigned j = 0; j < m; ++j) + tmp += lhs(0, j) * rhs[j]; + return tmp; + } + + // vec * mat + + template <unsigned m, typename T, + typename U> + inline + mat<m, m, mln_trait_op_times(T,U)> + operator*(const vec<m,T>& lhs, const mat<1,m,U>& rhs) + { + mat<m, m, mln_trait_op_times(T,U)> tmp; + for (unsigned i = 0; i < m; ++i) + for (unsigned j = 0; j < m; ++j) + tmp(i, j) = lhs[i] * rhs(0, j); + return tmp; + } + + // mat * s + template <unsigned n, unsigned m, typename T, typename S> inline mat<n, m, mln_trait_op_times(T,S)> @@ -499,16 +618,7 @@ } - template<unsigned n, unsigned m, typename T> - mat<m,n,T> - trans(const mat<n,m,T>& matrix) - { - mat<m,n,T> tmp; - for (unsigned i = 0; i < n; ++i) - for (unsigned j = 0; j < m; ++j) - tmp(j,i) = matrix(i,j); - return tmp; - } + // Trace. template<unsigned n, typename T> inline float tr(const mat<n,n,T>& m) @@ -519,6 +629,20 @@ return f; } + + // vec methods. + + template <unsigned n, typename T> + inline + mat<1,n,T> + vec<n,T>::t() const + { + mat<1,n,T> tmp; + for (unsigned i = 0; i < n; ++i) + tmp(0,i) = data_[i]; + return tmp; + } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::algebra Index: mln/algebra/vec.hh --- mln/algebra/vec.hh (revision 3690) +++ mln/algebra/vec.hh (working copy) @@ -59,6 +59,7 @@ namespace algebra { template <unsigned n, typename T> class vec; template <unsigned d, typename C> class h_vec; + template <unsigned n, unsigned m, typename T> class mat; } namespace literal { @@ -235,6 +236,9 @@ lower (sic) to zero. */ const vec<n, T>& normalize(); + /// Transposition. + mat<1, n, T> t() const; + /// Constructor; coordinates are set by function \p f. template <typename F> vec(const Function_i2v<F>& f); @@ -244,6 +248,19 @@ /// Origin value. static const vec<n, T> origin; + + + /// Conversion to a matrix. + template <typename U> + operator mat<n, 1, U>() const; + + /// Construction from a matrix. + template <typename U> + vec(const mat<n, 1, U>& rhs); + + /// Assignment from a matrix. + template <typename U> + vec& operator=(const mat<n, 1, U>& rhs); }; } // end of namespace mln::algebra @@ -632,6 +649,7 @@ # include <mln/make/vec.hh> +# include <mln/algebra/mat.hh> #endif // ! MLN_ALGEBRA_VEC_HH Index: mln/registration/get_rot.hh --- mln/registration/get_rot.hh (revision 3690) +++ mln/registration/get_rot.hh (working copy) @@ -1,4 +1,4 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory +// Copyright (C) 2008, 2009 EPITA Research and Development Laboratory // (LRDE) // // This file is part of the Olena Library. This library is free @@ -35,6 +35,7 @@ # include <mln/algebra/vec.hh> # include <mln/math/jacobi.hh> + namespace mln { @@ -97,7 +98,7 @@ { algebra::vec<3u,float> ci = convert::to< algebra::vec<3u,float> >(c[i]); algebra::vec<3u,float> xki = convert::to< algebra::vec<3u,float> >(map(ck[i])); - Mk += make::mat(ci - mu_c) * trans(make::mat(xki - mu_xk)); + Mk += (ci - mu_c) * (xki - mu_xk).t(); } Mk /= c.nsites(); Index: mln/registration/icp.hh --- mln/registration/icp.hh (revision 3690) +++ mln/registration/icp.hh (working copy) @@ -78,10 +78,10 @@ # include <mln/io/cloud/save.hh> + namespace mln { - namespace registration { @@ -530,7 +530,7 @@ vec3d_f P_i = p; vec3d_f Pk_i = qR.rotate(P_i) + qT; vec3d_f Yk_i = closest_point(Pk_i); - Spx += make::mat(P_i - mu_P) * trans(make::mat(Yk_i - mu_Yk)); + Spx += (P_i - mu_P) * (Yk_i - mu_Yk).t(); } Spx /= P_.nsites(); Index: mln/labeling/compute.hh --- mln/labeling/compute.hh (revision 3690) +++ mln/labeling/compute.hh (working copy) @@ -145,7 +145,7 @@ { mln_precondition(exact(input).is_valid()); mln_precondition(exact(label).is_valid()); - mlc_is_a(mln_value(L), mln::value::Symbolic)::check(); + // mlc_is_a(mln_value(L), mln::value::Symbolic)::check(); (void) a; (void) input; (void) label; Index: tests/accu/image/take_as_init.cc --- tests/accu/image/take_as_init.cc (revision 3690) +++ tests/accu/image/take_as_init.cc (working copy) @@ -45,4 +45,7 @@ accu::image::take_as_init(ima, 3); debug::println(ima); + + accu::image::take_as_init(ima, ima); + debug::println(ima); } Index: tests/algebra/op_times.cc --- tests/algebra/op_times.cc (revision 0) +++ tests/algebra/op_times.cc (revision 0) @@ -0,0 +1,72 @@ +// Copyright (C) 2009 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +/// \file tests/algebra/op_times.cc +/// +/// Tests on mln::algebra::operator *. + +#include <iostream> +#include <mln/algebra/mat.hh> + + + +int main() +{ + using namespace mln; + + using algebra::vec; + using algebra::mat; + + // Debase tests. + { + algebra::vec<3, float> v3; + algebra::mat<2,3,float> m23; + algebra::mat<3,2,float> m32; + + algebra::mat<2,2,float> m22 = m23 * m32; + algebra::vec<2, float> v2 = m23 * v3; + + m22 = v2 * v2.t(); + } + + // Tests with horizontal and/or vertical matrices and/or vectors. + { + float f; + algebra::vec<3, float> v3; + algebra::mat<1,3,float> m13; + algebra::mat<3,1,float> m31; + + f = m13 * m31; + f = m31.t() * m13.t(); + f = m13 * v3; + f = v3.t() * v3; + + v3 = m31; + m31 = v3; + } + +}