[Olena-patches] 1784: Icp legacy. algebra:: now contains types vec, mat and quat.

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet &lt;ugo.jardonnet(a)lrde.epita.fr&gt; Icp legacy. algebra:: now contains types vec, mat and quat. * tests/metal_mat.cc: . * tests/metal_pow.cc: . * tests/all.cc: . * tests/mat.cc: . * tests/metal_vec.cc: . * tests/core/h_vec.cc: . * tests/core/tr_image.cc: . * tests/stack.cc: . * tests/dpoint1d.cc: . * tests/dpoint2d.cc: . * tests/dpoint3d.cc: . * tests/fun/x2x/translation.cc: . * tests/fun/x2x/rotation.cc: . * tests/fun/x2x/composed.cc: . * tests/fun/vv2v/min.cc: . * tests/fun/vv2v/max.cc: . * tests/fun/v2v/norm.cc: . * tests/vec.cc: . * tests/point1d.cc: . * tests/point2d.cc: . * tests/norm/l1.cc: . * tests/norm/l2.cc: . * tests/norm/linfty.cc: . * tests/interpolated.cc: . * mln/trait/ch_value.hh: . * mln/trait/value_.hh: . * mln/core/point.hh: . * mln/core/line_graph_image.hh: . * mln/core/h_vec.hh: . * mln/core/tr_image.hh: . * mln/core/interpolated.hh: . * mln/core/bgraph_image.hh: . * mln/core/trait/op_mult.hh: . * mln/core/graph_image.hh: . * mln/core/h_mat.hh: . * mln/core/dpoint.hh: . * mln/core/plain.hh: . * mln/metal/mat.hh: . * mln/metal/math/pow.hh: . * mln/metal/math/max.hh: . * mln/metal/math/all.hh: . * mln/metal/math/sqrt.hh: . * mln/metal/all.hh: . * mln/metal/vec.hh: . * mln/linear/sobel.hh: . * mln/linear/gaussian.hh: . * mln/value/graylevel.hh: . * mln/value/graylevel_f.hh: . * mln/value/float01_.hh: . * mln/value/quat.hh: Remove. * mln/value/int_s.hh: . * mln/value/int_u.hh: . * mln/value/internal/gray_.hh: . * mln/value/internal/gray_f.hh: . * mln/value/int_u_sat.hh: . * mln/value/stack.hh: . * mln/value/rgb.hh: . * mln/value/label.hh: . * mln/make/vec.hh: . * mln/make/mat.hh: . * mln/make/voronoi.hh: . * mln/geom/seeds2tiling_roundness.hh: . * mln/geom/seeds2tiling.hh: . * mln/fun/x2x/composed.hh: . * mln/fun/x2x/translation.hh: . * mln/fun/x2x/rotation.hh: . * mln/fun/internal/selector.hh: . * mln/algebra/mat.hh: New. * mln/algebra/quat.hh: New. * mln/norm/linfty.hh: . * mln/norm/l1.hh: . * mln/norm/l2.hh: . * sandbox/duhamel/slow_seed2tiling.cc: . * sandbox/duhamel/labeling_algo.hh: . * sandbox/duhamel/mesh_image.hh: . * sandbox/jardonnet/test/icp.cc: New. * sandbox/jardonnet/test/gaussian_subsampling.cc: . * sandbox/jardonnet/test/Makefile: . * sandbox/jardonnet/subsampling/gaussian_subsampling.hh: . * sandbox/jardonnet/subsampling/sub_sampled_image.hh: . * sandbox/jardonnet/TODO: . * sandbox/jardonnet/registration/quat7.hh: New. * sandbox/jardonnet/registration/cloud.hh: New. * sandbox/jardonnet/registration/jacobi.hh: New. * sandbox/jardonnet/registration/icp.hh: . * sandbox/jardonnet/registration/projection.hh: New. * sandbox/jardonnet/registration/quat: New. * sandbox/jardonnet/registration/quat/all.hh: New. * sandbox/jardonnet/registration/quat/misc.hh: New. * sandbox/jardonnet/registration/quat/interpol.hh: New. * sandbox/jardonnet/registration/quat/rotation.hh: New. * sandbox/vigouroux/color/my_yuv.hh: . * sandbox/vigouroux/color/my_hsi.hh: . * sandbox/vigouroux/color/my_hsl.hh: . * sandbox/vigouroux/color/my_xyz.hh: . * sandbox/vigouroux/color/my_hsv.hh: . * sandbox/vigouroux/color/my_yiq.hh: . * sandbox/folio/dmap.cc: . * sandbox/garrigues/image_identity/interpolated.hh: . * sandbox/garrigues/image_identity/interpolated.cc: . mln/algebra/mat.hh | 487 +++++++++++++ mln/algebra/quat.hh | 648 ++++++++++++++++++ mln/core/bgraph_image.hh | 2 mln/core/dpoint.hh | 14 mln/core/graph_image.hh | 2 mln/core/h_mat.hh | 12 mln/core/h_vec.hh | 22 mln/core/interpolated.hh | 10 mln/core/line_graph_image.hh | 2 mln/core/plain.hh | 2 mln/core/point.hh | 14 mln/core/tr_image.hh | 20 mln/core/trait/op_mult.hh | 28 mln/fun/internal/selector.hh | 4 mln/fun/x2x/composed.hh | 2 mln/fun/x2x/rotation.hh | 22 mln/fun/x2x/translation.hh | 22 mln/geom/seeds2tiling.hh | 2 mln/geom/seeds2tiling_roundness.hh | 2 mln/linear/gaussian.hh | 4 mln/linear/sobel.hh | 2 mln/make/mat.hh | 23 mln/make/vec.hh | 36 - mln/make/voronoi.hh | 2 mln/metal/all.hh | 6 mln/metal/math/all.hh | 10 mln/metal/math/max.hh | 8 mln/metal/math/pow.hh | 8 mln/metal/math/sqrt.hh | 6 mln/norm/l1.hh | 12 mln/norm/l2.hh | 12 mln/norm/linfty.hh | 14 mln/trait/ch_value.hh | 2 mln/trait/value_.hh | 4 mln/value/float01_.hh | 2 mln/value/graylevel.hh | 10 mln/value/graylevel_f.hh | 2 mln/value/int_s.hh | 6 mln/value/int_u.hh | 2 mln/value/int_u_sat.hh | 4 mln/value/internal/gray_.hh | 6 mln/value/internal/gray_f.hh | 6 mln/value/label.hh | 2 mln/value/rgb.hh | 28 mln/value/stack.hh | 30 sandbox/duhamel/labeling_algo.hh | 6 sandbox/duhamel/mesh_image.hh | 2 sandbox/duhamel/slow_seed2tiling.cc | 2 sandbox/folio/dmap.cc | 2 sandbox/garrigues/image_identity/interpolated.cc | 6 sandbox/garrigues/image_identity/interpolated.hh | 10 sandbox/jardonnet/TODO | 6 sandbox/jardonnet/registration/cloud.hh | 49 + sandbox/jardonnet/registration/icp.hh | 53 - sandbox/jardonnet/registration/jacobi.hh | 108 +++ sandbox/jardonnet/registration/projection.hh | 52 + sandbox/jardonnet/registration/quat/all.hh | 9 sandbox/jardonnet/registration/quat/interpol.hh | 100 ++ sandbox/jardonnet/registration/quat/misc.hh | 19 sandbox/jardonnet/registration/quat/rotation.hh | 60 + sandbox/jardonnet/registration/quat7.hh | 148 ++++ sandbox/jardonnet/subsampling/gaussian_subsampling.hh | 4 sandbox/jardonnet/subsampling/sub_sampled_image.hh | 4 sandbox/jardonnet/test/Makefile | 7 sandbox/jardonnet/test/gaussian_subsampling.cc | 2 sandbox/jardonnet/test/icp.cc | 16 sandbox/vigouroux/color/my_hsi.hh | 2 sandbox/vigouroux/color/my_hsl.hh | 2 sandbox/vigouroux/color/my_hsv.hh | 2 sandbox/vigouroux/color/my_xyz.hh | 2 sandbox/vigouroux/color/my_yiq.hh | 2 sandbox/vigouroux/color/my_yuv.hh | 2 tests/all.cc | 2 tests/core/h_vec.cc | 6 tests/core/tr_image.cc | 2 tests/dpoint1d.cc | 2 tests/dpoint2d.cc | 2 tests/dpoint3d.cc | 2 tests/fun/v2v/norm.cc | 4 tests/fun/vv2v/max.cc | 4 tests/fun/vv2v/min.cc | 4 tests/fun/x2x/composed.cc | 2 tests/fun/x2x/rotation.cc | 2 tests/fun/x2x/translation.cc | 2 tests/interpolated.cc | 6 tests/mat.cc | 10 tests/metal_mat.cc | 16 tests/metal_pow.cc | 8 tests/metal_vec.cc | 8 tests/norm/l1.cc | 6 tests/norm/l2.cc | 6 tests/norm/linfty.cc | 4 tests/point1d.cc | 2 tests/point2d.cc | 2 tests/stack.cc | 2 tests/vec.cc | 12 96 files changed, 2041 insertions(+), 316 deletions(-) Index: tests/metal_mat.cc --- tests/metal_mat.cc (revision 1783) +++ tests/metal_mat.cc (working copy) @@ -27,10 +27,10 @@ /*! \file tests/metal_mat.cc * - * \brief Test on mln::metal::mat. + * \brief Test on mln::algebra::mat. */ -#include <mln/metal/mat.hh> +#include <mln/algebra/mat.hh> #include <mln/value/int_u8.hh> @@ -47,18 +47,18 @@ // tab2[6] = {2, 5, 1, 0, 7, 2}, // tab3[6] = {3, 1, 6, 2, 1, 0}; - // metal::mat<3,6,int> mat36 = make::mat<3,6,18>(tab1); - // metal::mat<2,3,int> mat23_1 = make::mat<2,3,6>(tab2); - // metal::mat<2,3,int> mat23_2 = make::mat<2,3,6>(tab3); + // algebra::mat<3,6,int> mat36 = make::mat<3,6,18>(tab1); + // algebra::mat<2,3,int> mat23_1 = make::mat<2,3,6>(tab2); + // algebra::mat<2,3,int> mat23_2 = make::mat<2,3,6>(tab3); - // metal::mat<2,3,float> mat23_3 = mat23_1 - mat23_2; + // algebra::mat<2,3,float> mat23_3 = mat23_1 - mat23_2; // std::cout << mat23_3 << std::endl << mat23_3 * mat36 << std::endl; - using metal::vec; + using algebra::vec; vec<2,int> v = make::vec(5,1); - using metal::mat; + using algebra::mat; mat<2,2, vec<2,int> > mv; mv.set_all(v); // std::cout << mv << std::endl; Index: tests/metal_pow.cc --- tests/metal_pow.cc (revision 1783) +++ tests/metal_pow.cc (working copy) @@ -27,12 +27,12 @@ /*! \file tests/metal_pow.cc * - * \brief Test on mln::metal::math. + * \brief Test on mln::algebra::math. */ #include <iostream> #include <mln/core/contract.hh> -#include <mln/metal/math/pow.hh> +#include <mln/algebra/math/pow.hh> int main() @@ -40,8 +40,8 @@ using namespace mln; using namespace mln::metal; - int res = metal::math::pow_int<2,3>::value; + int res = algebra::math::pow_int<2,3>::value; mln_assertion(res == 8); - std::cout << metal::math::pow< int_<2>, int_<3> >::ret().name() << std::endl; + std::cout << algebra::math::pow< int_<2>, int_<3> >::ret().name() << std::endl; } Index: tests/all.cc --- tests/all.cc (revision 1783) +++ tests/all.cc (working copy) @@ -48,7 +48,7 @@ #include <mln/math/all.hh> //#include <mln/set/all.hh> #include <mln/draw/all.hh> -#include <mln/metal/math/all.hh> +#include <mln/algebra/math/all.hh> #include <mln/metal/all.hh> #include <mln/morpho/all.hh> #include <mln/io/pfm/all.hh> Index: tests/mat.cc --- tests/mat.cc (revision 1783) +++ tests/mat.cc (working copy) @@ -27,13 +27,13 @@ /*! \file tests/mat.cc * - * \brief Tests on mln::metal::mat. + * \brief Tests on mln::algebra::mat. */ #include <iostream> #include <mln/fun/i2v/all_to.hh> -#include <mln/metal/mat.hh> +#include <mln/algebra/mat.hh> #include <mln/core/h_mat.hh> @@ -42,14 +42,14 @@ { using namespace mln; - metal::mat<1,3,float> m1(all_to(4.)); - metal::mat<2,2,float> m2 = metal::mat<2,2,float>::Id; + algebra::mat<1,3,float> m1(all_to(4.)); + algebra::mat<2,2,float> m2 = metal::mat<2,2,float>::Id; h_mat<1,float> hm1(m2); h_mat<2,float> hm2; h_mat<3,float> hm3(all_to(1.5)); - metal::mat<4,4,float> m4 = hm3; + algebra::mat<4,4,float> m4 = hm3; std::cout << "m1 = " << m1 << ";" << std::endl; std::cout << "m2 = " << m2 << ";" << std::endl; Index: tests/metal_vec.cc --- tests/metal_vec.cc (revision 1783) +++ tests/metal_vec.cc (working copy) @@ -27,10 +27,10 @@ /*! \file tests/metal_vec.cc * - * \brief Test on mln::metal::vec. + * \brief Test on mln::algebra::vec. */ -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/value/int_u8.hh> @@ -40,8 +40,8 @@ { using namespace mln; - metal::vec<3,int> v_int = make::vec(3, 6, 7); - metal::vec<3,float> v_f = make::vec(2.6, 1.9, 5.2); + algebra::vec<3,int> v_int = make::vec(3, 6, 7); + algebra::vec<3,float> v_f = make::vec(2.6, 1.9, 5.2); mln_assertion((v_int + v_f) == ((v_f + v_int))); mln_assertion((v_f / 3) == ((3 * v_f) / 9)); Index: tests/core/h_vec.cc --- tests/core/h_vec.cc (revision 1783) +++ tests/core/h_vec.cc (working copy) @@ -36,7 +36,7 @@ using namespace mln; -void run_in_3d(const metal::vec<3, int>&) +void run_in_3d(const algebra::vec<3, int>&) { } @@ -53,7 +53,7 @@ int main() { - metal::vec<3, int> x; + algebra::vec<3, int> x; h_vec<3, int> w = x.to_h_vec(); typedef h_vec<3, int> p3d; @@ -65,7 +65,7 @@ run_in_3d_h(k.to_h_vec()); { - metal::vec<3,float> v; + algebra::vec<3,float> v; h_vec<3,float> w(v.to_h_vec()); w = v.to_h_vec(); foo(v.to_h_vec()); Index: tests/core/tr_image.cc --- tests/core/tr_image.cc (revision 1783) +++ tests/core/tr_image.cc (working copy) @@ -58,7 +58,7 @@ for_all(p) { - metal::vec<3,int> vec = (image3d<int_u8>::point)p; + algebra::vec<3,int> vec = (image3d<int_u8>::point)p; if (inter.has(vec)) out(p) = inter(vec); else Index: tests/stack.cc --- tests/stack.cc (revision 1783) +++ tests/stack.cc (working copy) @@ -44,7 +44,7 @@ image2d<int> ima5(b), ima1(b); point2d p = make::point2d(0, 0); - metal::vec<2, int> v = make::vec(5, 1); + algebra::vec<2, int> v = make::vec(5, 1); value::stack(ima5, ima1)(p) = v; mln_assertion(value::stack(ima5, ima1)(p) == v); Index: tests/dpoint1d.cc --- tests/dpoint1d.cc (revision 1783) +++ tests/dpoint1d.cc (working copy) @@ -48,6 +48,6 @@ mln_assertion(dp == q - p); mln_assertion(q == p + dp); - metal::vec<1, float> v = dp; + algebra::vec<1, float> v = dp; mln_assertion(v[0] == 3); } Index: tests/dpoint2d.cc --- tests/dpoint2d.cc (revision 1783) +++ tests/dpoint2d.cc (working copy) @@ -48,7 +48,7 @@ mln_assertion(dp == q - p); mln_assertion(q == p + dp); - metal::vec<2, float> v = dp; + algebra::vec<2, float> v = dp; mln_assertion(v[0] / 2 == 1.5); p += dp; Index: tests/dpoint3d.cc --- tests/dpoint3d.cc (revision 1783) +++ tests/dpoint3d.cc (working copy) @@ -48,6 +48,6 @@ mln_assertion(dp == q - p); mln_assertion(q == p + dp); - metal::vec<3, float> v = dp; + algebra::vec<3, float> v = dp; mln_assertion(v[1] == 6); } Index: tests/fun/x2x/translation.cc --- tests/fun/x2x/translation.cc (revision 1783) +++ tests/fun/x2x/translation.cc (working copy) @@ -46,7 +46,7 @@ b = 0, c = 2.9; - metal::vec<3,float> vec1 = make::vec(a, b, c); + algebra::vec<3,float> vec1 = make::vec(a, b, c); fun::x2x::translation<3,float> tr1(all_to(1.6)); std::cout << vec1 << std::endl; Index: tests/fun/x2x/rotation.cc --- tests/fun/x2x/rotation.cc (revision 1783) +++ tests/fun/x2x/rotation.cc (working copy) @@ -59,7 +59,7 @@ for_all(p) { - metal::vec<2,float> v = rot1.inv()((point2d::vec_t)(point2d)p); + algebra::vec<2,float> v = rot1.inv()((point2d::vec_t)(point2d)p); if (inter.owns_(v)) out(p) = inter(v); else Index: tests/fun/x2x/composed.cc --- tests/fun/x2x/composed.cc (revision 1783) +++ tests/fun/x2x/composed.cc (working copy) @@ -48,7 +48,7 @@ b = 0, c = 2.9; - metal::vec<3,float> vec1 = make::vec(a, b, c); + algebra::vec<3,float> vec1 = make::vec(a, b, c); fun::x2x::translation<3,float> tr(all_to(1.6)); fun::x2x::rotation<3,float> rot(0.3, 1); Index: tests/fun/vv2v/min.cc --- tests/fun/vv2v/min.cc (revision 1783) +++ tests/fun/vv2v/min.cc (working copy) @@ -30,7 +30,7 @@ #include <cassert> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/fun/vv2v/min.hh> int main() @@ -43,7 +43,7 @@ // work. #if 0 // Vectors. - typedef mln::metal::vec<3, int> vec_t; + typedef mln::algebra::vec<3, int> vec_t; mln::fun::vv2v::min<vec_t> min_vec_t; vec_t t, u; t.set (1, -2, 3); Index: tests/fun/vv2v/max.cc --- tests/fun/vv2v/max.cc (revision 1783) +++ tests/fun/vv2v/max.cc (working copy) @@ -30,7 +30,7 @@ #include <cassert> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/fun/vv2v/max.hh> int main() @@ -43,7 +43,7 @@ // work. #if 0 // Vectors. - typedef mln::metal::vec<3, int> vec_t; + typedef mln::algebra::vec<3, int> vec_t; mln::fun::vv2v::max<vec_t> max_vec_t; vec_t t, u; t.set (1, -2, 3); Index: tests/fun/v2v/norm.cc --- tests/fun/v2v/norm.cc (revision 1783) +++ tests/fun/v2v/norm.cc (working copy) @@ -30,7 +30,7 @@ * \brief Test the norm functors. */ -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/fun/v2v/norm.hh> #include <tests/norm/common.hh> @@ -38,7 +38,7 @@ int main() { - typedef mln::metal::vec<3, int> vec_t; + typedef mln::algebra::vec<3, int> vec_t; // L1-norm. { Index: tests/vec.cc --- tests/vec.cc (revision 1783) +++ tests/vec.cc (working copy) @@ -27,13 +27,13 @@ /*! \file tests/vec.cc * - * \brief Tests on mln::metal::vec. + * \brief Tests on mln::algebra::vec. */ #include <iostream> #include <mln/fun/i2v/all_to.hh> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/core/h_vec.hh> @@ -42,8 +42,8 @@ { using namespace mln; - metal::vec<1,float> v1(all_to(4.)); - metal::vec<2,float> v2 = make::vec(6., 2.8); + algebra::vec<1,float> v1(all_to(4.)); + algebra::vec<2,float> v2 = make::vec(6., 2.8); h_vec<1,float> hv1; h_vec<2,float> hv2 = v2.to_h_vec(); // Immersion into homogeneous. @@ -52,8 +52,8 @@ hv3 += make::vec(0., 0., 0., 0.5); - metal::vec<3,float> v3 = hv3.to_vec(); // Back from homogeneous. - metal::vec<4,float> v4 = hv3; + algebra::vec<3,float> v3 = hv3.to_vec(); // Back from homogeneous. + algebra::vec<4,float> v4 = hv3; std::cout << "v1 = " << v1 << ";" << std::endl; std::cout << "v2 = " << v2 << ";" << std::endl; Index: tests/point1d.cc --- tests/point1d.cc (revision 1783) +++ tests/point1d.cc (working copy) @@ -43,7 +43,7 @@ // assignment p[0] = 4; - metal::vec<1,int> v = p; + algebra::vec<1,int> v = p; std::cout << v << std::endl; p.ind() += 1; Index: tests/point2d.cc --- tests/point2d.cc (revision 1783) +++ tests/point2d.cc (working copy) @@ -74,5 +74,5 @@ for (unsigned i = 0; i < p.dim; ++i) mln_assertion(q[i] == 0); - std::cout << 3.4 * metal::vec<2, int>(p) << std::endl; + std::cout << 3.4 * algebra::vec<2, int>(p) << std::endl; } Index: tests/norm/l1.cc --- tests/norm/l1.cc (revision 1783) +++ tests/norm/l1.cc (working copy) @@ -32,7 +32,7 @@ #include <tests/norm/common.hh> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/norm/l1.hh> @@ -58,12 +58,12 @@ int main() { - typedef mln::metal::vec<3, int> vec_t; + typedef mln::algebra::vec<3, int> vec_t; // Reference value. int d = (5 - 1) + (1 + 2) + 3; - // Tests using mln::metal::vec. + // Tests using mln::algebra::vec. vec_t t, u; t.set(1, -2, 3); u.set(5, 1, 0); Index: tests/norm/l2.cc --- tests/norm/l2.cc (revision 1783) +++ tests/norm/l2.cc (working copy) @@ -35,7 +35,7 @@ #include <tests/norm/common.hh> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/norm/l2.hh> @@ -61,14 +61,14 @@ int main() { - typedef mln::metal::vec<3, int> vec_t; + typedef mln::algebra::vec<3, int> vec_t; // Reference value. float d = std::sqrt((4 - 2) * (4 - 2) + (1 + 2) * (1 + 2) + (0 - 3) * (0 - 3)); - // Tests using mln::metal::vec. + // Tests using mln::algebra::vec. vec_t t, u; t.set(2, -2, 3); u.set(4, 1, 0); Index: tests/norm/linfty.cc --- tests/norm/linfty.cc (revision 1783) +++ tests/norm/linfty.cc (working copy) @@ -32,7 +32,7 @@ #include <tests/norm/common.hh> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/math/abs.hh> #include <mln/norm/linfty.hh> @@ -59,7 +59,7 @@ int main() { - typedef mln::metal::vec<3, int> vec_t; + typedef mln::algebra::vec<3, int> vec_t; // Reference value. float d = std::max(std::abs(4 - 2), Index: tests/interpolated.cc --- tests/interpolated.cc (revision 1783) +++ tests/interpolated.cc (working copy) @@ -35,7 +35,7 @@ #include <mln/core/image2d.hh> #include <mln/core/interpolated.hh> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/level/fill.hh> @@ -60,8 +60,8 @@ interpolated< image2d<float> > inter(f); - metal::vec<2, float> v1 = make::vec(2.3, 0.6); - metal::vec<2, float> v2 = make::vec(3.2, 1.8); + algebra::vec<2, float> v1 = make::vec(2.3, 0.6); + algebra::vec<2, float> v2 = make::vec(3.2, 1.8); debug::println(f); Index: mln/trait/ch_value.hh --- mln/trait/ch_value.hh (revision 1783) +++ mln/trait/ch_value.hh (working copy) @@ -98,7 +98,7 @@ { /* FIXME: The code used to read - typedef metal::vec<n, V> value; + typedef algebra::vec<n, V> value; typedef mln_ch_value(I, value) ret; here. But this is wrong IMHO (Roland). Changing the value Index: mln/trait/value_.hh --- mln/trait/value_.hh (revision 1783) +++ mln/trait/value_.hh (working copy) @@ -37,11 +37,11 @@ # include <string> # include <mln/metal/int.hh> -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/pow.hh> # include <mln/metal/if.hh> # include <mln/trait/value/all.hh> -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/pow.hh> # define mln_trait_value_nature(V) typename mln::trait::value_< V >::nature Index: mln/core/point.hh --- mln/core/point.hh (revision 1783) +++ mln/core/point.hh (working copy) @@ -38,7 +38,7 @@ # include <mln/fun/i2v/all_to.hh> # include <mln/metal/bool.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/metal/converts_to.hh> # include <mln/core/h_vec.hh> @@ -62,7 +62,7 @@ template <typename M, typename C> struct point_to_ { - typedef metal::vec<M::dim, C> metal_vec; + typedef algebra::vec<M::dim, C> metal_vec; typedef h_vec<M::dim, C> h_vec; }; @@ -143,17 +143,17 @@ point_<M,C>& operator-=(const dpoint& dp); /// Type of the array of coordinates. - typedef metal::vec<M::dim, C> vec_t; + typedef algebra::vec<M::dim, C> vec_t; /// Hook to coordinates. operator typename internal::point_to_<M, C>::metal_vec () const; - operator metal::vec<M::dim, float> () const; + operator algebra::vec<M::dim, float> () const; /// Transform to point in homogene coordinate system. h_vec<M::dim, C> to_h_vec() const; protected: - metal::vec<M::dim, C> coord_; + algebra::vec<M::dim, C> coord_; }; @@ -312,9 +312,9 @@ template <typename M, typename C> inline - point_<M,C>::operator metal::vec<M::dim, float> () const + point_<M,C>::operator algebra::vec<M::dim, float> () const { - metal::vec<dim, float> tmp; + algebra::vec<dim, float> tmp; for (unsigned i = 0; i < dim; ++i) tmp[i] = coord_[i]; return tmp; Index: mln/core/line_graph_image.hh --- mln/core/line_graph_image.hh (revision 1783) +++ mln/core/line_graph_image.hh (working copy) @@ -34,7 +34,7 @@ # include <mln/trait/images.hh> # include <mln/core/internal/image_primary.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/core/p_line_graph.hh> # include <mln/core/line_graph_psite.hh> # include <mln/value/set.hh> Index: mln/core/h_vec.hh --- mln/core/h_vec.hh (revision 1783) +++ mln/core/h_vec.hh (working copy) @@ -33,7 +33,7 @@ * \brief Definition of a vector with homogeneous coordinates. */ -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/literal/one.hh> namespace mln @@ -86,7 +86,7 @@ * */ template <unsigned d, typename C> - struct h_vec : public metal::vec<d + 1, C> + struct h_vec : public algebra::vec<d + 1, C> { /// Dimension is the 'natural' one (3 for 3D), not the one of the vector (dim + 1). enum { dim = d }; @@ -94,12 +94,12 @@ /// Constructor without argument. h_vec(); /// Constructor with the underlying vector. - h_vec(const metal::vec<d+1, C>& other); + h_vec(const algebra::vec<d+1, C>& other); - h_vec& operator=(const metal::vec<d+1, C>& rhs); + h_vec& operator=(const algebra::vec<d+1, C>& rhs); /// Back to the natural (non-homogeneous) space. - metal::vec<d,C> to_vec() const; + algebra::vec<d,C> to_vec() const; }; @@ -119,18 +119,18 @@ template <unsigned d, typename C> inline - h_vec<d,C>::h_vec(const metal::vec<d+1, C>& other) - : metal::vec<d+1, C>(other) + h_vec<d,C>::h_vec(const algebra::vec<d+1, C>& other) + : algebra::vec<d+1, C>(other) { } template <unsigned d, typename C> inline - h_vec<d,C>& h_vec<d,C>::operator=(const metal::vec<d+1, C>& rhs) + h_vec<d,C>& h_vec<d,C>::operator=(const algebra::vec<d+1, C>& rhs) { if (& rhs == this) return *this; - this->metal::vec<d+1, C>::operator=(rhs); + this->algebra::vec<d+1, C>::operator=(rhs); return *this; } @@ -154,12 +154,12 @@ template <unsigned d, typename C> inline - metal::vec<d,C> h_vec<d,C>::to_vec() const + algebra::vec<d,C> h_vec<d,C>::to_vec() const { const C w = this->data_[d]; mln_assertion(w != 0); - metal::vec<d,C> tmp; + algebra::vec<d,C> tmp; for (unsigned i = 0; i < d; ++i) tmp[i] = this->data_[i] / w; return tmp; Index: mln/core/tr_image.hh --- mln/core/tr_image.hh (revision 1783) +++ mln/core/tr_image.hh (working copy) @@ -38,7 +38,7 @@ # include <cmath> # include <mln/core/internal/image_identity.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/value/set.hh> @@ -114,18 +114,18 @@ using super_::owns_; /// Test if a pixel value is accessible at \p v. - bool owns_(const mln::metal::vec<I::point::dim, float>& v) const; + bool owns_(const mln::algebra::vec<I::point::dim, float>& v) const; using super_::has; /// Test if a pixel value is belonging to image at \p v. - bool has(const mln::metal::vec<I::point::dim, float>& v) const; + bool has(const mln::algebra::vec<I::point::dim, float>& v) const; /// Read-only access of pixel value at point site \p p. /// Mutable access is only OK for reading (not writing). using super_::operator(); - mln_value(I) operator()(const mln::metal::vec<I::point::dim, float>& v) const; + mln_value(I) operator()(const mln::algebra::vec<I::point::dim, float>& v) const; void set_tr(T& tr); }; @@ -173,10 +173,10 @@ template <typename T, typename I> inline - bool tr_image<T,I>::owns_(const metal::vec<I::point::dim, float>& v) const + bool tr_image<T,I>::owns_(const algebra::vec<I::point::dim, float>& v) const { mln_point(I) p; - metal::vec<I::point::dim, float> v2 = this->data_->tr_.inv()(v); + algebra::vec<I::point::dim, float> v2 = this->data_->tr_.inv()(v); for (unsigned i = 0; i < I::point::dim; ++i) p[i] = static_cast<int>(round(v2[i])); return this->delegatee_().owns_(p); @@ -184,10 +184,10 @@ template <typename T, typename I> inline - bool tr_image<T,I>::has(const metal::vec<I::point::dim, float>& v) const + bool tr_image<T,I>::has(const algebra::vec<I::point::dim, float>& v) const { mln_point(I) p; - metal::vec<I::point::dim, float> v2 = this->data_->tr_.inv()(v); + algebra::vec<I::point::dim, float> v2 = this->data_->tr_.inv()(v); for (unsigned i = 0; i < I::point::dim; ++i) p[i] = static_cast<int>(round(v2[i])); return this->delegatee_()->domain().has(p); @@ -196,10 +196,10 @@ template <typename T, typename I> inline mln_value(I) - tr_image<T,I>::operator()(const metal::vec<I::point::dim, float>& v) const + tr_image<T,I>::operator()(const algebra::vec<I::point::dim, float>& v) const { mln_point(I) p; - metal::vec<I::point::dim, float> v2 = this->data_->tr_.inv()(v); + algebra::vec<I::point::dim, float> v2 = this->data_->tr_.inv()(v); for (unsigned i = 0; i < I::point::dim; ++i) p[i] = static_cast<int>(round(v2[i])); mln_assertion(this->delegatee_()->owns_(p)); Index: mln/core/interpolated.hh --- mln/core/interpolated.hh (revision 1783) +++ mln/core/interpolated.hh (working copy) @@ -37,7 +37,7 @@ # include <cmath> # include <mln/core/internal/image_identity.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/value/set.hh> namespace mln @@ -101,13 +101,13 @@ using super_::owns_; /// Test if a pixel value is accessible at \p v. - bool owns_(const mln::metal::vec<I::point::dim, float>& v) const; + bool owns_(const mln::algebra::vec<I::point::dim, float>& v) const; /// Read-only access of pixel value at point site \p p. /// Mutable access is only OK for reading (not writing). using super_::operator(); - mln_value(I) operator()(const mln::metal::vec<I::point::dim, float>& v) const; + mln_value(I) operator()(const mln::algebra::vec<I::point::dim, float>& v) const; /// Give the set of values of the image. @@ -156,7 +156,7 @@ template <typename I> inline - bool interpolated<I>::owns_(const mln::metal::vec<I::point::dim, float>& v) const + bool interpolated<I>::owns_(const mln::algebra::vec<I::point::dim, float>& v) const { mln_point(I) p; for (unsigned i = 0; i < I::point::dim; ++i) @@ -167,7 +167,7 @@ template <typename I> inline mln_value(I) - interpolated<I>::operator()(const mln::metal::vec<I::point::dim, float>& v) const + interpolated<I>::operator()(const mln::algebra::vec<I::point::dim, float>& v) const { mln_point(I) p; for (unsigned i = 0; i < I::point::dim; ++i) Index: mln/core/bgraph_image.hh --- mln/core/bgraph_image.hh (revision 1783) +++ mln/core/bgraph_image.hh (working copy) @@ -33,7 +33,7 @@ # include <mln/trait/images.hh> # include <mln/core/internal/image_primary.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/core/p_bgraph.hh> # include <mln/core/bgraph_psite.hh> Index: mln/core/trait/op_mult.hh --- mln/core/trait/op_mult.hh (revision 1783) +++ mln/core/trait/op_mult.hh (working copy) @@ -101,43 +101,43 @@ }; template <unsigned n, typename T, typename U> - struct op_mult<metal::vec<n, T>, U> + struct op_mult<algebra::vec<n, T>, U> { - typedef metal::vec<n, mln_op_mult(T, U)> ret; + typedef algebra::vec<n, mln_op_mult(T, U)> ret; }; template <typename U, unsigned n, typename T> - struct op_mult<U, metal::vec<n, T> > + struct op_mult<U, algebra::vec<n, T> > { - typedef metal::vec<n, mln_op_mult(T, U)> ret; + typedef algebra::vec<n, mln_op_mult(T, U)> ret; }; template <unsigned n, unsigned m, typename T, typename U> - struct op_mult<metal::mat<n, m, T>, U> + struct op_mult<algebra::mat<n, m, T>, U> { - typedef metal::mat<n, m, mln_op_mult(T, U)> ret; + typedef algebra::mat<n, m, mln_op_mult(T, U)> ret; }; template <typename U, unsigned n, unsigned m, typename T> - struct op_mult<U, metal::mat<n, m, T> > + struct op_mult<U, algebra::mat<n, m, T> > { - typedef metal::mat<n, m, mln_op_mult(T, U)> ret; + typedef algebra::mat<n, m, mln_op_mult(T, U)> ret; }; template <unsigned n, unsigned o, typename T, unsigned m, typename U> - struct op_mult<metal::mat<n, o, T>, metal::mat<o, m, U> > + struct op_mult<algebra::mat<n, o, T>, algebra::mat<o, m, U> > { - typedef metal::mat<n, m, mln_op_mult(T, U)> ret; + typedef algebra::mat<n, m, mln_op_mult(T, U)> ret; }; template <unsigned m, unsigned n, typename T, typename U> - struct op_mult<metal::mat<m, n, T>, metal::vec<n, U> > + struct op_mult<algebra::mat<m, n, T>, algebra::vec<n, U> > { - typedef metal::mat<m, 1, mln_op_mult(T, U)> ret; + typedef algebra::mat<m, 1, mln_op_mult(T, U)> ret; }; template <unsigned n, typename U, unsigned m, typename T> - struct op_mult< metal::vec<n, U>, metal::mat<n, m, T> > + struct op_mult< algebra::vec<n, U>, algebra::mat<n, m, T> > { - typedef metal::mat<1, m, mln_op_mult(T, U)> ret; + typedef algebra::mat<1, m, mln_op_mult(T, U)> ret; }; Index: mln/core/graph_image.hh --- mln/core/graph_image.hh (revision 1783) +++ mln/core/graph_image.hh (working copy) @@ -34,7 +34,7 @@ # include <mln/trait/images.hh> # include <mln/core/internal/image_primary.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/core/p_graph.hh> # include <mln/core/graph_psite.hh> # include <mln/value/set.hh> Index: mln/core/h_mat.hh --- mln/core/h_mat.hh (revision 1783) +++ mln/core/h_mat.hh (working copy) @@ -33,7 +33,7 @@ * \brief Definition of a matrix with homogeneous coordinates. */ -# include <mln/metal/mat.hh> +# include <mln/algebra/mat.hh> namespace mln @@ -43,7 +43,7 @@ * */ template <unsigned d, typename T> - struct h_mat : public metal::mat<d+1, d+1, T> + struct h_mat : public algebra::mat<d+1, d+1, T> { /// Dimension is the 'natural' one (3 for 3D), not the one of the vector (dim + 1) enum { N = d, @@ -53,7 +53,7 @@ /// Constructor without argument. h_mat(); /// Constructor with the underlying matrix. - h_mat(const metal::mat<d+1, d+1, T>& x); + h_mat(const algebra::mat<d+1, d+1, T>& x); }; @@ -62,14 +62,14 @@ template <unsigned d, typename T> inline h_mat<d,T>::h_mat() - : metal::mat<d+1, d+1, T>(metal::mat<d+1, d+1, T>::Id) + : algebra::mat<d+1, d+1, T>(algebra::mat<d+1, d+1, T>::Id) { } template <unsigned d, typename T> inline - h_mat<d,T>::h_mat(const metal::mat<d+1, d+1, T>& x) - : metal::mat<d+1, d+1, T>(x) + h_mat<d,T>::h_mat(const algebra::mat<d+1, d+1, T>& x) + : algebra::mat<d+1, d+1, T>(x) { } Index: mln/core/dpoint.hh --- mln/core/dpoint.hh (revision 1783) +++ mln/core/dpoint.hh (working copy) @@ -36,7 +36,7 @@ # include <mln/core/concept/dpoint.hh> # include <mln/core/internal/coord_impl.hh> # include <mln/fun/i2v/all.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/metal/converts_to.hh> @@ -113,15 +113,15 @@ /// Set all coordinates to the value \p c. void set_all(C c); - /// Conversion towards a metal::vec. + /// Conversion towards a algebra::vec. template <typename Q> - operator metal::vec<M::dim, Q>() const; + operator algebra::vec<M::dim, Q>() const; /// Explicit conversion. - metal::vec<M::dim, C> to_vec() const; + algebra::vec<M::dim, C> to_vec() const; protected: - metal::vec<M::dim, C> coord_; + algebra::vec<M::dim, C> coord_; }; @@ -232,14 +232,14 @@ template <typename M, typename C> template <typename Q> inline - dpoint_<M,C>::operator metal::vec<M::dim, Q> () const + dpoint_<M,C>::operator algebra::vec<M::dim, Q> () const { return coord_; } template <typename M, typename C> inline - metal::vec<M::dim, C> + algebra::vec<M::dim, C> dpoint_<M,C>::to_vec() const { return coord_; Index: mln/core/plain.hh --- mln/core/plain.hh (revision 1783) +++ mln/core/plain.hh (working copy) @@ -39,7 +39,7 @@ # include <mln/core/internal/image_identity.hh> # include <mln/core/clone.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> namespace mln Index: mln/metal/math/pow.hh --- mln/metal/math/pow.hh (revision 1783) +++ mln/metal/math/pow.hh (working copy) @@ -28,7 +28,7 @@ #ifndef MLN_METAL_MATH_POW_HH # define MLN_METAL_MATH_POW_HH -/*! \file mln/metal/math/pow.hh +/*! \file mln/algebra/math/pow.hh * * \brief Definition of the 'power' static function. */ @@ -37,8 +37,8 @@ # include <mln/metal/int.hh> -# define mlc_pow(X, N) typename mln::metal::math::pow< X, N >::ret -# define mlc_pow_int(x, n) mln::metal::math::pow_int< x, n >::value +# define mlc_pow(X, N) typename mln::algebra::math::pow< X, N >::ret +# define mlc_pow_int(x, n) mln::algebra::math::pow_int< x, n >::value @@ -105,7 +105,7 @@ }; - } // end of namespace mln::metal::math + } // end of namespace mln::algebra::math } // end of namespace mln::metal Index: mln/metal/math/max.hh --- mln/metal/math/max.hh (revision 1783) +++ mln/metal/math/max.hh (working copy) @@ -28,7 +28,7 @@ #ifndef MLN_METAL_MATH_MAX_HH # define MLN_METAL_MATH_MAX_HH -/*! \file mln/metal/math/max.hh +/*! \file mln/algebra/math/max.hh * * \brief Definition of the 'max' static function. */ @@ -36,8 +36,8 @@ # include <mln/metal/bool.hh> # include <mln/metal/int.hh> -# define mlc_max(X, Y) typename mln::metal::math::max< X, Y >::ret -# define mlc_max_int(x, y) mln::metal::math::max_int< x, y >::value +# define mlc_max(X, Y) typename mln::algebra::math::max< X, Y >::ret +# define mlc_max_int(x, y) mln::algebra::math::max_int< x, y >::value namespace mln { @@ -69,7 +69,7 @@ }; - } // end of namespace mln::metal::math + } // end of namespace mln::algebra::math } // end of namespace mln::metal Index: mln/metal/math/all.hh --- mln/metal/math/all.hh (revision 1783) +++ mln/metal/math/all.hh (working copy) @@ -28,7 +28,7 @@ #ifndef MLN_METAL_MATH_ALL_HH # define MLN_METAL_MATH_ALL_HH -/*! \file mln/metal/math/all.hh +/*! \file mln/algebra/math/all.hh * * \brief Include all static mathematical functions. */ @@ -43,7 +43,7 @@ /// Namespace of static mathematical functions. namespace math { - /// \internal Implementation namespace of metal::math namespace. + /// \internal Implementation namespace of algebra::math namespace. namespace impl {} } @@ -53,9 +53,9 @@ -# include <mln/metal/math/pow.hh> -# include <mln/metal/math/sqrt.hh> -# include <mln/metal/math/max.hh> +# include <mln/algebra/math/pow.hh> +# include <mln/algebra/math/sqrt.hh> +# include <mln/algebra/math/max.hh> // ... Index: mln/metal/math/sqrt.hh --- mln/metal/math/sqrt.hh (revision 1783) +++ mln/metal/math/sqrt.hh (working copy) @@ -28,7 +28,7 @@ #ifndef MLN_METAL_MATH_SQRT_HH # define MLN_METAL_MATH_SQRT_HH -/*! \file mln/metal/math/sqrt.hh +/*! \file mln/algebra/math/sqrt.hh * * \brief Definition of the 'sqrt' static function. */ @@ -82,7 +82,7 @@ { }; - } // end of namespace mln::metal::math::impl + } // end of namespace mln::algebra::math::impl template <int n> struct sqrt_int : impl::sqrt_int_if_< n, (n >= 0) > @@ -102,7 +102,7 @@ }; - } // end of namespace mln::metal::math + } // end of namespace mln::algebra::math } // end of namespace mln::metal Index: mln/metal/all.hh --- mln/metal/all.hh (revision 1783) +++ mln/metal/all.hh (working copy) @@ -74,10 +74,10 @@ # include <mln/metal/unqualif.hh> # include <mln/metal/is_unqualif.hh> -# include <mln/metal/vec.hh> -# include <mln/metal/mat.hh> +# include <mln/algebra/vec.hh> +# include <mln/algebra/mat.hh> -# include <mln/metal/math/all.hh> +# include <mln/algebra/math/all.hh> // FIXME: Remove the following includes below! # include <mln/metal/same_coord.hh> Index: mln/linear/sobel.hh --- mln/linear/sobel.hh (revision 1783) +++ mln/linear/sobel.hh (working copy) @@ -243,7 +243,7 @@ // The type of a component of a vector from the gradient. typedef mln_sum(mln_value(I)) gradient_val_t; // The type of a vector from the gradient. - typedef mln::metal::vec<I::point::dim, gradient_val_t> gradient_vec_t; + typedef mln::algebra::vec<I::point::dim, gradient_val_t> gradient_vec_t; return sobel_norm(input, fun::v2v::l1_norm<gradient_vec_t, gradient_val_t>()); } Index: mln/linear/gaussian.hh --- mln/linear/gaussian.hh (revision 1783) +++ mln/linear/gaussian.hh (working copy) @@ -327,7 +327,7 @@ inline void gaussian(const Image<I>& input, float sigma, - Image<O>& out) + Image<O>& output) { mln_precondition(exact(input).has_data()); mln_precondition(exact(output).has_data()); @@ -339,7 +339,7 @@ 0.6318f, 1.997f, sigma); impl::gaussian_common_(mln_trait_value_nature(mln_value(I))(), - input, coef, sigma, out); + input, coef, sigma, output); } # endif // ! MLN_INCLUDE_ONLY Index: mln/value/graylevel.hh --- mln/value/graylevel.hh (revision 1783) +++ mln/value/graylevel.hh (working copy) @@ -38,8 +38,8 @@ # include <mln/value/ops.hh> # include <mln/core/contract.hh> -# include <mln/metal/math/pow.hh> -# include <mln/metal/math/max.hh> +# include <mln/algebra/math/pow.hh> +# include <mln/algebra/math/max.hh> # include <mln/metal/bexpr.hh> # include <mln/literal/ops.hh> @@ -513,7 +513,7 @@ inline graylevel<n>::graylevel(const literal::medium_gray_t&) { - this->v_ = metal::math::pow_int<2, n - 1>::value; + this->v_ = algebra::math::pow_int<2, n - 1>::value; } template <unsigned n> @@ -521,7 +521,7 @@ graylevel<n>& graylevel<n>::operator=(const literal::medium_gray_t&) { - this->v_ = metal::math::pow_int<2, n - 1>::value; + this->v_ = algebra::math::pow_int<2, n - 1>::value; return *this; } @@ -555,7 +555,7 @@ float graylevel<n>::to_float() const { - static const float denom = float(metal::math::pow_int<2, n>::value) - 1.f; + static const float denom = float(algebra::math::pow_int<2, n>::value) - 1.f; return float(this->v_) / denom; } Index: mln/value/graylevel_f.hh --- mln/value/graylevel_f.hh (revision 1783) +++ mln/value/graylevel_f.hh (working copy) @@ -38,7 +38,7 @@ # include <mln/value/ops.hh> # include <mln/core/contract.hh> -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/pow.hh> # include <mln/metal/bexpr.hh> # include <mln/literal/ops.hh> Index: mln/value/float01_.hh --- mln/value/float01_.hh (revision 1783) +++ mln/value/float01_.hh (working copy) @@ -35,7 +35,7 @@ # include <iostream> # include <mln/core/contract.hh> -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/pow.hh> # include <mln/metal/bexpr.hh> # include <mln/value/int_u.hh> Index: mln/value/int_s.hh --- mln/value/int_s.hh (revision 1783) +++ mln/value/int_s.hh (working copy) @@ -35,7 +35,7 @@ # include <mln/value/ops.hh> -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/pow.hh> # include <mln/value/internal/value_like.hh> # include <mln/value/concept/integer.hh> # include <mln/value/internal/encoding.hh> @@ -180,7 +180,7 @@ inline int_s<n>::int_s(int i) { - static const int max = metal::math::pow_int<2, n-1>::value - 1; + static const int max = algebra::math::pow_int<2, n-1>::value - 1; static const int min = - max; mln_precondition(i >= min); mln_precondition(i <= max); @@ -192,7 +192,7 @@ int_s<n>& int_s<n>::operator=(int i) { - static const int max = metal::math::pow_int<2, n-1>::value - 1; + static const int max = algebra::math::pow_int<2, n-1>::value - 1; static const int min = - max; mln_precondition(i >= min); mln_precondition(i <= max); Index: mln/value/int_u.hh --- mln/value/int_u.hh (revision 1783) +++ mln/value/int_u.hh (working copy) @@ -35,7 +35,7 @@ # include <mln/value/ops.hh> -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/pow.hh> # include <mln/value/internal/value_like.hh> # include <mln/value/internal/encoding.hh> # include <mln/value/concept/integer.hh> Index: mln/value/internal/gray_.hh --- mln/value/internal/gray_.hh (revision 1783) +++ mln/value/internal/gray_.hh (working copy) @@ -37,8 +37,8 @@ # include <iostream> # include <cmath> -# include <mln/metal/math/max.hh> -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/max.hh> +# include <mln/algebra/math/pow.hh> # include <mln/value/concept/integer.hh> @@ -312,7 +312,7 @@ inline gray_<n>::operator graylevel_f() const { - static const float denom = float(metal::math::pow_int<2, n>::value) - 1.f; + static const float denom = float(algebra::math::pow_int<2, n>::value) - 1.f; return graylevel_f(float(this->v_) / denom); } Index: mln/value/internal/gray_f.hh --- mln/value/internal/gray_f.hh (revision 1783) +++ mln/value/internal/gray_f.hh (working copy) @@ -38,7 +38,7 @@ # include <mln/value/ops.hh> # include <mln/core/contract.hh> -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/pow.hh> # include <mln/metal/bexpr.hh> # include <mln/literal/ops.hh> @@ -225,7 +225,7 @@ template <unsigned n> gray_f::gray_f(const gray_<n>& rhs) { - static const float denom = float(metal::math::pow_int<2, n>::value) - 1.f; + static const float denom = float(algebra::math::pow_int<2, n>::value) - 1.f; this->v_ = float(rhs.value()) / denom; } @@ -233,7 +233,7 @@ gray_f& gray_f::operator=(const gray_<n>& rhs) { - static const float denom = float(metal::math::pow_int<2, n>::value) - 1.f; + static const float denom = float(algebra::math::pow_int<2, n>::value) - 1.f; this->v_ = float(rhs.value()) / denom; return *this; } Index: mln/value/int_u_sat.hh --- mln/value/int_u_sat.hh (revision 1783) +++ mln/value/int_u_sat.hh (working copy) @@ -34,7 +34,7 @@ * behavior. */ -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/pow.hh> # include <mln/value/internal/value_like.hh> # include <mln/value/concept/integer.hh> # include <mln/value/internal/encoding.hh> @@ -59,7 +59,7 @@ struct value_< mln::value::int_u_sat<n> > { // FIXME: Overhaul these traits (see other value traits). - static const std::size_t card = metal::math::pow_int<2, n>::value; + static const std::size_t card = algebra::math::pow_int<2, n>::value; static const mln::value::int_u_sat<n> min() { return 0; } static const mln::value::int_u_sat<n> max() { return card - 1; } static const unsigned nbits = n; Index: mln/value/stack.hh --- mln/value/stack.hh (revision 1783) +++ mln/value/stack.hh (working copy) @@ -35,7 +35,7 @@ # include <mln/core/internal/image_value_morpher.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/value/set.hh> # include <mln/value/proxy.hh> @@ -57,8 +57,8 @@ struct data_< value::stack_image<n, I> > { public: - data_(const metal::vec<n,I>& imas); - metal::vec<n,I> imas_; + data_(const algebra::vec<n,I>& imas); + algebra::vec<n,I> imas_; I& ima_; }; @@ -83,7 +83,7 @@ template <unsigned n, typename I> struct helper_stack_image_lvalue_< n, const I > { - typedef metal::vec<n, mln_value(I)> ret; + typedef algebra::vec<n, mln_value(I)> ret; static ret make(stack_image<n, const I>& ima, const mln_psite(I)& p) { return ima.read_(p); @@ -102,7 +102,7 @@ template <unsigned n, typename I> struct image_< mln::value::stack_image<n, I> > : default_image_morpher_< I, - metal::vec<n, mln_value(I)>, + algebra::vec<n, mln_value(I)>, mln::value::stack_image<n, I> > { // FIXME: We shall carefully define the missing required traits @@ -155,7 +155,7 @@ typedef mln_pset(I) pset; /// Value associated type. - typedef metal::vec<n, mln_value(I)> value; + typedef algebra::vec<n, mln_value(I)> value; /// Return type of read-only access. typedef value rvalue; @@ -173,7 +173,7 @@ /// Constructors. /// \{ - stack_image(const metal::vec<n,I>& imas); + stack_image(const algebra::vec<n,I>& imas); stack_image(); /// \} @@ -215,7 +215,7 @@ template <unsigned n, typename I> inline - data_< value::stack_image<n,I> >::data_(const metal::vec<n,I>& imas) + data_< value::stack_image<n,I> >::data_(const algebra::vec<n,I>& imas) : imas_(imas), ima_(imas_[0]) { @@ -236,7 +236,7 @@ template <unsigned n, typename I> inline - stack_image<n,I>::stack_image(const metal::vec<n,I>& imas) + stack_image<n,I>::stack_image(const algebra::vec<n,I>& imas) { this->data_ = new mln::internal::data_< stack_image<n, I> >(imas); for (unsigned i = 0; i < n; ++i) @@ -256,11 +256,11 @@ template <unsigned n, typename I> inline - metal::vec<n, mln_value(I)> + algebra::vec<n, mln_value(I)> stack_image<n,I>::read_(const psite& p) const { mln_precondition(this->owns_(p)); - metal::vec<n, mln_value(I)> tmp; + algebra::vec<n, mln_value(I)> tmp; for (unsigned i = 0; i < n; ++i) tmp[i] = this->data_->imas_[i].operator()(p); return tmp; @@ -268,7 +268,7 @@ template <unsigned n, typename I> inline - metal::vec<n, mln_value(I)> + algebra::vec<n, mln_value(I)> stack_image<n,I>::operator()(const psite& p) const { return read_(p); @@ -295,7 +295,7 @@ template <unsigned n, typename I> inline - const mln::value::set< metal::vec<n, mln_value(I)> >& + const mln::value::set< algebra::vec<n, mln_value(I)> >& stack_image<n,I>::values() const { return vset::the(); @@ -309,7 +309,7 @@ stack(const Image<I>& ima1, const Image<I>& ima2) { mln_precondition(exact(ima1).domain() == exact(ima2).domain()); - metal::vec<2, const I> imas; + algebra::vec<2, const I> imas; imas[0] = exact(ima1); imas[1] = exact(ima2); return imas; @@ -321,7 +321,7 @@ stack(Image<I>& ima1, Image<I>& ima2) { mln_precondition(exact(ima1).domain() == exact(ima2).domain()); - metal::vec<2, I> imas; + algebra::vec<2, I> imas; imas[0] = exact(ima1); imas[1] = exact(ima2); return imas; Index: mln/value/rgb.hh --- mln/value/rgb.hh (revision 1783) +++ mln/value/rgb.hh (working copy) @@ -38,7 +38,7 @@ # include <mln/value/concept/vectorial.hh> # include <mln/value/int_u.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> namespace mln @@ -137,7 +137,7 @@ typedef trait::value::kind::color kind; typedef mln_value_quant_from_(card) quant; - typedef metal::vec<3, float> sum; + typedef algebra::vec<3, float> sum; }; } // end of namespace trait @@ -155,9 +155,9 @@ : public Vectorial< rgb<n> > , - public internal::value_like_< metal::vec< 3, int_u<n> >, // Equivalent. - metal::vec< 3, int_u<n> >, // Encoding. - metal::vec< 3, int >, // Interoperation. + public internal::value_like_< algebra::vec< 3, int_u<n> >, // Equivalent. + algebra::vec< 3, int_u<n> >, // Encoding. + algebra::vec< 3, int >, // Interoperation. rgb<n> > // Exact. { public: @@ -179,15 +179,15 @@ /// Constructor from component values. rgb<n>(int r, int g, int b); - /// Constructor from a metal::vec. - rgb<n>(const metal::vec<3, int>& rhs); - rgb<n>(const metal::vec<3, unsigned>& rhs); - rgb<n>(const metal::vec<3, int_u<n> >& rhs); + /// Constructor from a algebra::vec. + rgb<n>(const algebra::vec<3, int>& rhs); + rgb<n>(const algebra::vec<3, unsigned>& rhs); + rgb<n>(const algebra::vec<3, int_u<n> >& rhs); // Conversion to the interoperation type. - operator metal::vec<3, int>() const { return this->v_; } + operator algebra::vec<3, int>() const { return this->v_; } // Conversion to the sum type. - operator metal::vec<3, float>() const { return this->v_; } + operator algebra::vec<3, float>() const { return this->v_; } /// \{ Constructors with literals. rgb<n>(const literal::white_t&); @@ -293,21 +293,21 @@ template <unsigned n> inline - rgb<n>::rgb(const metal::vec<3, int>& v) + rgb<n>::rgb(const algebra::vec<3, int>& v) { this->v_ = v; } template <unsigned n> inline - rgb<n>::rgb(const metal::vec<3, unsigned>& v) + rgb<n>::rgb(const algebra::vec<3, unsigned>& v) { this->v_ = v; } template <unsigned n> inline - rgb<n>::rgb(const metal::vec<3, int_u<n> >& v) + rgb<n>::rgb(const algebra::vec<3, int_u<n> >& v) { this->v_ = v; } Index: mln/value/label.hh --- mln/value/label.hh (revision 1783) +++ mln/value/label.hh (working copy) @@ -33,7 +33,7 @@ * \brief Define a generic class for labels. */ -# include <mln/metal/math/pow.hh> +# include <mln/algebra/math/pow.hh> # include <mln/value/concept/symbolic.hh> # include <mln/value/internal/value_like.hh> # include <mln/value/internal/convert.hh> Index: mln/make/vec.hh --- mln/make/vec.hh (revision 1783) +++ mln/make/vec.hh (working copy) @@ -30,10 +30,10 @@ /*! \file mln/make/vec.hh * - * \brief Routine to construct an mln::metal::vec. + * \brief Routine to construct an mln::algebra::vec. */ -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/core/concept/function.hh> namespace mln @@ -43,16 +43,16 @@ { - /*! \brief Create an mln::metal::vec<n,T>. + /*! \brief Create an mln::algebra::vec<n,T>. * * \param[in] v_0 First coordinate. * * \return A 1D vector. */ template <typename T> - metal::vec<1, T> vec(const T& v_0); + algebra::vec<1, T> vec(const T& v_0); - /*! \brief Create an mln::metal::vec<2,T>. + /*! \brief Create an mln::algebra::vec<2,T>. * * \param[in] v_0 First coordinate. * \param[in] v_1 Second coordinate. @@ -60,9 +60,9 @@ * \return A 2D vector. */ template <typename T> - metal::vec<2, T> vec(const T& v_0, const T& v_1); + algebra::vec<2, T> vec(const T& v_0, const T& v_1); - /*! \brief Create an mln::metal::vec<3,T>. + /*! \brief Create an mln::algebra::vec<3,T>. * * \param[in] v_0 First coordinate. * \param[in] v_1 Second coordinate. @@ -71,9 +71,9 @@ * \return A 3D vector. */ template <typename T> - metal::vec<3, T> vec(const T& v_0, const T& v_1, const T& v_2); + algebra::vec<3, T> vec(const T& v_0, const T& v_1, const T& v_2); - /*! \brief Create an mln::metal::vec<4,T>. + /*! \brief Create an mln::algebra::vec<4,T>. * * \param[in] v_0 First coordinate. * \param[in] v_1 Second coordinate. @@ -83,25 +83,25 @@ * \return A 4D vector. */ template <typename T> - metal::vec<4, T> vec(const T& v_0, const T& v_1, const T& v_2, const T& v_3); + algebra::vec<4, T> vec(const T& v_0, const T& v_1, const T& v_2, const T& v_3); # ifndef MLN_INCLUDE_ONLY template <typename T> inline - metal::vec<1, T> vec(const T& v_0) + algebra::vec<1, T> vec(const T& v_0) { - metal::vec<1, T> tmp; + algebra::vec<1, T> tmp; tmp[0] = v_0; return tmp; } template <typename T> inline - metal::vec<2, T> vec(const T& v_0, const T& v_1) + algebra::vec<2, T> vec(const T& v_0, const T& v_1) { - metal::vec<2, T> tmp; + algebra::vec<2, T> tmp; tmp[0] = v_0; tmp[1] = v_1; return tmp; @@ -109,9 +109,9 @@ template <typename T> inline - metal::vec<3, T> vec(const T& v_0, const T& v_1, const T& v_2) + algebra::vec<3, T> vec(const T& v_0, const T& v_1, const T& v_2) { - metal::vec<3, T> tmp; + algebra::vec<3, T> tmp; tmp[0] = v_0; tmp[1] = v_1; tmp[2] = v_2; @@ -120,9 +120,9 @@ template <typename T> inline - metal::vec<4, T> vec(const T& v_0, const T& v_1, const T& v_2, const T& v_3) + algebra::vec<4, T> vec(const T& v_0, const T& v_1, const T& v_2, const T& v_3) { - metal::vec<4, T> tmp; + algebra::vec<4, T> tmp; tmp[0] = v_0; tmp[1] = v_1; tmp[2] = v_2; Index: mln/make/mat.hh --- mln/make/mat.hh (revision 1783) +++ mln/make/mat.hh (working copy) @@ -30,17 +30,17 @@ /*! \file mln/make/mat.hh * - * \brief Routine to construct an mln::metal::mat. + * \brief Routine to construct an mln::algebra::mat. */ -# include <mln/metal/mat.hh> +# include <mln/algebra/mat.hh> namespace mln { namespace make { - /*! \brief Create an mln::metal::mat<n,m,T>. + /*! \brief Create an mln::algebra::mat<n,m,T>. * * \param[in] tab Tab of value. * @@ -48,21 +48,32 @@ * with n and m, the dimensions oh the matrix. */ template <unsigned n, unsigned m, unsigned N, typename T> - metal::mat<n,m,T> mat(const T tab[N]); + algebra::mat<n,m,T> mat(const T tab[N]); + + template <unsigned n, unsigned m, typename T> + algebra::mat<n,m,T> mat(algebra::vec<n,T> v); # ifndef MLN_INCLUDE_ONLY template <unsigned n, unsigned m, unsigned N, typename T> inline - metal::mat<n,m,T> mat(const T tab[N]) + algebra::mat<n,m,T> mat(const T tab[N]) { mln_precondition(n * m == N); - metal::mat<n,m,T> tmp; + algebra::mat<n,m,T> tmp; for (unsigned i = 0; i < N; ++i) tmp(i / m, i % m) = tab[i]; return tmp; } + template <unsigned n, typename T> + algebra::mat<n,1,T> mat(algebra::vec<n,T> v) + { + algebra::mat<n,1,T> tmp; + for (unsigned i = 0; i < n; i++) + tmp(i,1) = v[i]; + return tmp; + } # endif // ! MLN_INCLUDE_ONLY Index: mln/make/voronoi.hh --- mln/make/voronoi.hh (revision 1783) +++ mln/make/voronoi.hh (working copy) @@ -77,7 +77,7 @@ Image<I>& orig_, const Neighborhood<N>& nbh) { - typedef metal::vec<2,float> X; + typedef algebra::vec<2,float> X; typedef mln_value(I) V; typedef mln_psite(I) P; Index: mln/geom/seeds2tiling_roundness.hh --- mln/geom/seeds2tiling_roundness.hh (revision 1783) +++ mln/geom/seeds2tiling_roundness.hh (working copy) @@ -41,7 +41,7 @@ # include <mln/core/clone.hh> # include <mln/accu/mean.hh> # include <mln/estim/min_max.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/geom/chamfer.hh> Index: mln/geom/seeds2tiling.hh --- mln/geom/seeds2tiling.hh (revision 1783) +++ mln/geom/seeds2tiling.hh (working copy) @@ -40,7 +40,7 @@ # include <mln/core/clone.hh> # include <mln/accu/mean.hh> # include <mln/estim/min_max.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> namespace mln Index: mln/fun/x2x/composed.hh --- mln/fun/x2x/composed.hh (revision 1783) +++ mln/fun/x2x/composed.hh (working copy) @@ -35,7 +35,7 @@ # include <mln/core/concept/function.hh> # include <mln/fun/internal/x2x_linear_impl.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/metal/is.hh> # include <mln/metal/bexpr.hh> # include <mln/core/h_mat.hh> Index: mln/fun/x2x/translation.hh --- mln/fun/x2x/translation.hh (revision 1783) +++ mln/fun/x2x/translation.hh (working copy) @@ -35,7 +35,7 @@ # include <mln/core/concept/function.hh> # include <mln/fun/internal/x2x_linear_impl.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/core/h_mat.hh> # include <mln/fun/i2v/all.hh> @@ -54,11 +54,11 @@ template <unsigned n, typename C> struct translation - : internal::x2x_linear_impl_< metal::vec<n,C>, translation<n,C> > + : internal::x2x_linear_impl_< algebra::vec<n,C>, translation<n,C> > , public Bijection_x2x< translation<n,C> > { - typedef fun::internal::x2x_linear_impl_< metal::vec<n,C>, translation<n,C> > super_; + typedef fun::internal::x2x_linear_impl_< algebra::vec<n,C>, translation<n,C> > super_; /// Type of the inverse function. typedef translation<n,C> invert; @@ -68,19 +68,19 @@ /// Constructor without argument. translation(); /// Constructor with the translation vector. - translation(const metal::vec<n,C>& t); + translation(const algebra::vec<n,C>& t); using super_::operator(); /// Perform the translation of the given vector - metal::vec<n,C> operator()(const metal::vec<n,C>& v) const; + algebra::vec<n,C> operator()(const algebra::vec<n,C>& v) const; /// Set a net translation vector. - void set_t(const metal::vec<n,C>& t); + void set_t(const algebra::vec<n,C>& t); protected: void update(); - metal::vec<n,C> t_; + algebra::vec<n,C> t_; }; @@ -94,7 +94,7 @@ template <unsigned n, typename C> inline - translation<n,C>::translation(const metal::vec<n,C>& t) + translation<n,C>::translation(const algebra::vec<n,C>& t) :t_(t) { this->m_ = h_mat<n,C>::Id; @@ -103,8 +103,8 @@ template <unsigned n, typename C> inline - metal::vec<n,C> - translation<n,C>::operator()(const metal::vec<n,C>& v) const + algebra::vec<n,C> + translation<n,C>::operator()(const algebra::vec<n,C>& v) const { return v + t_; } @@ -122,7 +122,7 @@ template <unsigned n, typename C> inline void - translation<n,C>::set_t(const metal::vec<n,C>& t) + translation<n,C>::set_t(const algebra::vec<n,C>& t) { this->t_ = t; this->update(); Index: mln/fun/x2x/rotation.hh --- mln/fun/x2x/rotation.hh (revision 1783) +++ mln/fun/x2x/rotation.hh (working copy) @@ -35,8 +35,8 @@ # include <mln/core/concept/function.hh> # include <mln/fun/internal/x2x_linear_impl.hh> -# include <mln/metal/vec.hh> -# include <mln/metal/mat.hh> +# include <mln/algebra/vec.hh> +# include <mln/algebra/mat.hh> # include <cmath> namespace mln @@ -53,10 +53,10 @@ */ template <unsigned n, typename C> struct rotation - : internal::x2x_linear_impl_< metal::vec<n,C>, rotation<n,C> > + : internal::x2x_linear_impl_< algebra::vec<n,C>, rotation<n,C> > , public Bijection_x2x< rotation<n,C> > { - typedef fun::internal::x2x_linear_impl_< metal::vec<n,C>, rotation<n,C> > super_; + typedef fun::internal::x2x_linear_impl_< algebra::vec<n,C>, rotation<n,C> > super_; /// Type of the inverse function. typedef rotation<n,C> invert; @@ -70,7 +70,7 @@ using super_::operator(); /// Perform the rotation of the given vector. - metal::vec<n,C> operator()(const metal::vec<n,C>& v) const; + algebra::vec<n,C> operator()(const algebra::vec<n,C>& v) const; /// Set a new grade alpha. void set_alpha(float alpha); @@ -106,12 +106,12 @@ template <unsigned n, typename C> inline - metal::vec<n,C> - rotation<n,C>::operator()(const metal::vec<n,C>& v) const + algebra::vec<n,C> + rotation<n,C>::operator()(const algebra::vec<n,C>& v) const { - metal::mat<n+1,1,C> hmg; - metal::mat<n+1,1,C> tmp; - metal::vec<n,C> res; + algebra::mat<n+1,1,C> hmg; + algebra::mat<n+1,1,C> tmp; + algebra::vec<n,C> res; for (unsigned i = 0; i < n; ++i) hmg(i,0) = v[i]; @@ -158,7 +158,7 @@ { const float cos_a = cos(alpha_); const float sin_a = sin(alpha_); - const metal::vec<4,float> vec = make::vec(cos_a, -sin_a, sin_a, cos_a); + const algebra::vec<4,float> vec = make::vec(cos_a, -sin_a, sin_a, cos_a); unsigned k = 0; for (unsigned i = 0; i < n; ++i) Index: mln/fun/internal/selector.hh --- mln/fun/internal/selector.hh (revision 1783) +++ mln/fun/internal/selector.hh (working copy) @@ -38,7 +38,7 @@ # include <mln/metal/unqualif.hh> # include <mln/metal/if.hh> # include <mln/metal/is_a.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> namespace mln @@ -152,7 +152,7 @@ }; template <unsigned n, typename T> - struct tag_< metal::vec<n,T> > + struct tag_< algebra::vec<n,T> > { enum { value = x_ }; }; Index: mln/algebra/mat.hh --- mln/algebra/mat.hh (revision 0) +++ mln/algebra/mat.hh (revision 0) @@ -0,0 +1,487 @@ +// Copyright (C) 2006 EPITA Research and Development Laboratory +// +// 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. + +#ifndef MLN_ALGEBRA_MAT_HH +# define MLN_ALGEBRA_MAT_HH + +/*! + * \file mln/algebra/mat.hh + * + * \brief Definition of a generic matrix class. + */ + +# include <iostream> + +# include <mln/core/concept/object.hh> +# include <mln/core/concept/function.hh> +# include <mln/core/contract.hh> +# include <mln/trait/all.hh> +# include <mln/trait/value_.hh> +# include <mln/algebra/vec.hh> + + +// FIXME: Document. + + + +namespace mln +{ + + + // Fwd decl. + namespace algebra { + template <unsigned n, unsigned m, typename T> class mat; + } + + + namespace trait + { + + template <unsigned n, unsigned m, typename T> + struct value_< algebra::mat<n,m,T> > + { + typedef trait::value::nature::matrix nature; + typedef trait::value::kind::data kind; + + enum { + nbits = n * m * mln_nbits(T), + card = n * m * mln_card(T) + }; + typedef mln_value_quant_from_(card) quant; + + typedef algebra::mat<n, m, mln_sum(T)> sum; + }; + + } // end of namespace mln::trait + + + + namespace algebra + { + + template <unsigned n, unsigned m, typename T> + class mat : public Object< mat<n,m,T> > + { + public: + + typedef T coord; + enum { N = n, + M = m, + dim = n * m }; + + static const mat<n,m,T> Id; + + mat() + { + } + + template <typename U> + mat(const mat<n,m,U>& rhs); + + /// Constructor; coordinates are set by function \p f. + template <typename F> + mat(const Function_i2v<F>& f); + + template <typename U> + mat& operator=(const mat<n,m,U>& rhs); + + const T& operator()(unsigned i, unsigned j) const; + + T& operator()(unsigned i, unsigned j); + + void set_all(const T& val); + + unsigned size() const; + + static mat identity(); + + private: + T data_[n][m]; + }; + + } + + + namespace trait + { + + // Unarys. + + template < template<class> class Name, + unsigned n, unsigned m, typename T > + struct set_precise_unary_< Name, algebra::mat<n,m,T> > + { + typedef algebra::mat<n, m, mln_trait_unary(Name, T)> ret; + }; + + // Default for binarys; works for (+), (-), comparisons, and promote. + + template < template<class, class> class Name, + unsigned n, unsigned m, typename T, typename U> + struct set_precise_binary_< Name, algebra::mat<n,m,T>, algebra::mat<n,m,U> > + { + typedef algebra::mat<n, m, mln_trait_binary(Name, T, U)> ret; + }; + + // mat * mat + + template < unsigned n, unsigned o, typename T, + unsigned m, typename U > + struct set_precise_binary_< op::times, algebra::mat<n,o,T>, algebra::mat<o,m,U> > + { + typedef algebra::mat<n, m, mln_sum_x(T, U)> ret; + }; + + template < unsigned n, typename T, typename U > + struct set_precise_binary_< op::times, algebra::mat<n,n,T>, algebra::mat<n,n,U> > + { // Disambiguate between both previous defs. + typedef algebra::mat<n, n, mln_sum_x(T, U)> ret; + }; + + // mat * vec + + template < unsigned n, unsigned m, typename T, + typename U > + struct set_precise_binary_< op::times, algebra::mat<n,m,T>, algebra::vec<m,U> > + { + typedef algebra::vec<n, mln_sum_x(T, U)> ret; + }; + + // mat * s + + template < template<class, class> class Name, + unsigned n, unsigned m, typename T, + typename S > + struct set_precise_binary_< Name, algebra::mat<n,m,T>, mln::value::scalar_<S> > + { + typedef algebra::mat<n, m, mln_trait_binary(Name, T, S)> ret; + }; + + template < template<class, class> class Name, + unsigned n, unsigned m, typename T, + typename S > + struct set_binary_< Name, + mln::Object, algebra::mat<n,m,T>, + mln::value::Scalar, S > + { + typedef algebra::mat<n, m, mln_trait_binary(Name, T, S)> ret; + }; + + } // end of namespace mln::trait + + + + namespace algebra + { + + // == + + template <unsigned n, unsigned m, typename T, typename U> + bool + operator==(mat<n,m,T>& lhs, const mat<n,m,U>& rhs); + + // + + + template <unsigned n, unsigned m, typename T, typename U> + inline + mat<n, m, mln_trait_op_plus(T,U)> + operator+(const mat<n,m,T>& lhs, const mat<n,m,U>& rhs); + + // - + + template <unsigned n, unsigned m, typename T, typename U> + mat<n, m, mln_trait_op_minus(T,U)> + operator-(const mat<n,m,T>& lhs, const mat<n,m,U>& rhs); + + // - (unary) + + template <unsigned n, unsigned m, typename T> + mat<n, m, mln_trait_op_uminus(T)> + operator-(const mat<n,m,T>& lhs); + + // mat * mat + + template <unsigned n, unsigned o, typename T, + unsigned m, typename U> + mat<n, m, mln_sum_x(T,U)> + operator*(const mat<n,o,T>& lhs, const mat<o,m,U>& rhs); + + // mat * vec + + template <unsigned n, unsigned m, typename T, + typename U> + vec<n, mln_sum_x(T,U)> + operator*(const mat<n,m,T>& lhs, const vec<m,U>& rhs); + + // mat * s + + template <unsigned n, unsigned m, typename T, + typename S> + mat<n, m, mln_trait_op_times(T,S)> + operator*(const mat<n,m,T>& lhs, const value::scalar_<S>& s); + + // mat / s + + template <unsigned n, unsigned m, typename T, typename S> + mat<n, m, mln_trait_op_div(T,S)> + operator/(const mat<n,m,T>& lhs, const value::scalar_<S>& s); + + // << + + template <unsigned n, unsigned m, typename T> + std::ostream& + operator<<(std::ostream& ostr, const mat<n,m,T>& v); + + + +# ifndef MLN_INCLUDE_ONLY + + template <unsigned n, unsigned m, typename T> + 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() + { + static mat<n,m,T> id_; + static bool flower = true; + if (flower) + { + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + id_(i, j) = (i == j); + flower = false; + } + return id_; + } + + template <unsigned n, unsigned m, typename T> + template <typename U> + inline + mat<n,m,T>::mat(const mat<n,m,U>& rhs) + { + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + data_[i][j] = rhs(i, j); + } + + template <unsigned n, unsigned m, typename T> + template <typename F> + inline + mat<n,m,T>::mat(const Function_i2v<F>& f_) + { + mlc_converts_to(mln_result(F), T)::check(); + const F& f = exact(f_); + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + data_[i][j] = f(i * n + j); + } + + template <unsigned n, unsigned m, typename T> + template <typename U> + inline + mat<n,m,T>& + mat<n,m,T>::operator=(const mat<n,m,U>& rhs) + { + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + data_[i][j] = rhs(i, j); + return *this; + } + + template <unsigned n, unsigned m, typename T> + inline + const T& + mat<n,m,T>::operator()(unsigned i, unsigned j) const + { + mln_precondition(i < n && j < m); + return data_[i][j]; + } + + template <unsigned n, unsigned m, typename T> + inline + T& + mat<n,m,T>::operator()(unsigned i, unsigned j) + { + mln_precondition(i < n && j < m); + return data_[i][j]; + } + + template <unsigned n, unsigned m, typename T> + inline + void mat<n,m,T>::set_all(const T& val) + { + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + data_[i][j] = val; + } + + template <unsigned n, unsigned m, typename T> + inline + unsigned mat<n,m,T>::size() const + { + return n * m; + } + + + // Operators. + + + template <unsigned n, unsigned m, typename T, typename U> + inline + bool + operator==(const mat<n,m,T>& lhs, const mat<n,m,U>& rhs) + { + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + if (lhs(i, j) != rhs(i, j)) + return false; + return true; + } + + template <unsigned n, unsigned m, typename T, typename U> + inline + mat<n, m, mln_trait_op_plus(T,U)> + operator+(const mat<n,m,T>& lhs, const mat<n,m,U>& rhs) + { + mat<n, m, mln_trait_op_plus(T,U)> tmp; + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + tmp(i, j) = lhs(i, j) + rhs(i, j); + return tmp; + } + + template <unsigned n, unsigned m, typename T, typename U> + inline + mat<n,m, mln_trait_op_minus(T,U)> + operator-(const mat<n,m,T>& lhs, const mat<n,m,U>& rhs) + { + mat<n,m, mln_trait_op_minus(T,U)> tmp; + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + tmp(i, j) = lhs(i, j) - rhs(i, j); + return tmp; + } + + template <unsigned n, unsigned m, typename T> + inline + mat<n,m, mln_trait_op_uminus(T)> + operator-(const mat<n,m,T>& rhs) + { + mat<n,m, mln_trait_op_uminus(T)> tmp; + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; i < m; ++i) + tmp(i, j) = - rhs(i, j); + return tmp; + } + + template <unsigned n, unsigned o, typename T, + unsigned m, typename U> + inline + mat<n, m, mln_sum_x(T,U)> + operator*(const mat<n,o,T>& lhs, const mat<o,m,U>& rhs) + { + mat<n,m, mln_sum_x(T,U)> tmp; + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + { + tmp(i, j) = literal::zero; + for (unsigned k = 0; k < o; ++k) + tmp(i, j) += lhs(i, k) * rhs(k, j); + } + return tmp; + } + + template <unsigned n, unsigned m, typename T, + typename U> + inline + vec<n, mln_sum_x(T,U)> + operator*(const mat<n,m,T>& lhs, const vec<m,U>& rhs) + { + vec<n, mln_sum_x(T,U)> tmp; + for (unsigned i = 0; i < n; ++i) + { + mln_sum_x(T,U) sum(literal::zero); + for (unsigned j = 0; j < m; ++j) + sum += lhs(i, j) * rhs[j]; + tmp[i] = sum; + } + return tmp; + } + + template <unsigned n, unsigned m, typename T, typename S> + inline + mat<n, m, mln_trait_op_times(T,S)> + operator*(const mat<n,m,T>& lhs, const value::scalar_<S>& s_) + { + S s = s_.to_equiv(); + mat<n, m, mln_trait_op_times(T,S)> tmp; + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + tmp(i, j) = lhs(i, j) * s; + return tmp; + } + + template <unsigned n, unsigned m, typename T, typename S> + inline + mat<n,m, mln_trait_op_div(T,S)> + operator/(const mat<n,m,T>& lhs, const value::scalar_<S>& s_) + { + S s = s_.to_equiv(); + mat<n,m, mln_trait_op_times(T,S)> tmp; + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + tmp(i,j) = lhs(i, j) / s; + return tmp; + } + + // << + + template <unsigned n, unsigned m, typename T> + inline + std::ostream& + operator<<(std::ostream& ostr, const mat<n,m,T>& v) + { + for (unsigned i = 0; i < n; ++i) + { + ostr << '['; + for (unsigned j = 0; j < m; ++j) + ostr << debug::format(v(i, j)) << (j == m - 1 ? "]" : ", "); + ostr << std::endl; + } + return ostr; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::algebra + +} // end of namespace mln + +# include <mln/make/mat.hh> + +#endif // ! MLN_ALGEBRA_MAT_HH Index: mln/algebra/quat.hh --- mln/algebra/quat.hh (revision 0) +++ mln/algebra/quat.hh (revision 0) @@ -0,0 +1,648 @@ +// Copyright (C) 2007 EPITA Research and Development Laboratory +// +// 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. + +#ifndef MLN_ALGEBRA_QUAT_HH +# define MLN_ALGEBRA_QUAT_HH + +/*! \file mln/algebra/quat.hh + * + * \brief Define a class for quaternion values. + */ + +# include <cmath> + +# include <mln/value/ops.hh> + +# include <mln/value/concept/vectorial.hh> +# include <mln/value/internal/value_like.hh> +# include <mln/trait/value_.hh> + +# include <mln/algebra/vec.hh> +# include <mln/norm/l2.hh> + + +namespace mln +{ + + // Fwd decls. + namespace value { class quat; } + namespace literal { struct zero_t; struct one_t; } + + + namespace trait + { + + // quat OP quat + + template <> + struct set_precise_binary_< op::plus, mln::algebra::quat, mln::algebra::quat > + { + typedef mln::algebra::quat ret; + }; + + template <> + struct set_precise_binary_< op::minus, mln::algebra::quat, mln::algebra::quat > + { + typedef mln::algebra::quat ret; + }; + + template <> + struct set_precise_binary_< op::times, mln::algebra::quat, mln::algebra::quat > + { + typedef mln::algebra::quat ret; + }; + + // quat OP scalar + + template < typename S > + struct set_precise_binary_< op::times, mln::algebra::quat, mln::value::scalar_<S> > + { + typedef mln::algebra::quat ret; + }; + + template < typename S > + struct set_precise_binary_< op::div, mln::algebra::quat, mln::value::scalar_<S> > + { + typedef mln::algebra::quat ret; + }; + + + // 'quat' as a value. + + + template <> + struct value_< mln::algebra::quat > + { + typedef trait::value::nature::vectorial nature; + typedef trait::value::kind::data kind; + typedef trait::value::quant::high quant; + + enum { + nbits = 4 * sizeof(float), + card = 0 + }; + + typedef mln::algebra::quat sum; + }; + + + } // end of namespace mln::trait + + + + namespace value + { + + // FIXME doesn't compile + + class quat + : + public Vectorial< quat > + , + public internal::value_like_< algebra::vec<4, float>, // Equivalent. + algebra::vec<4, float>, // Encoding. + algebra::vec<4, float>, // Interoperation. + quat > // Exact. + { + public: + + /// Constructor without argument. + quat(); + + /// Constructor with components. + quat(float s, float x, float y, float z); + + /// Constructor from a scalar and a 3D vector. + quat(float s, const algebra::vec<3,float>& v); + + + /// Constructor from a 4D vector. + quat(const algebra::vec<4,float>& v); + + /// Assignment from a 4D vector. + quat& operator=(const algebra::vec<4,float>& v); + + + /// \{ Constructors/assignments with literals zero and one. + quat(const literal::zero_t&); + quat& operator=(const literal::zero_t&); + quat(const literal::one_t&); + quat& operator=(const literal::one_t&); + /// \} + + + /// Explicit conversion to a 4D algebra::vec. + const algebra::vec<4,float>& to_vec() const; + + + /// Give the scalar part. + float s() const; + + /// Access to the scalar part. + float& s(); + + const algebra::vec<3,float>& v() const; + algebra::vec<3,float>& v(); + + void set_v(float x, float y, float z); + + /// Scalar product. + float sprod(const quat& rhs) const; + + /// Test if it is a unit quaternion. + bool is_unit() const; + + /// Test if the quaternion is null. + bool is_null() const; + + /// Test if it is a pure quaternion. + bool is_pure() const; + + /// Give the conjugate. + quat conj() const; + + /// Give the invert. + quat inv() const; // FIXME: Rename inv as invert. + + /// Transform into unit quaternion. + quat& set_unit(); + + /// Transform into unit quaternion. + template <typename T> + void set_unit(float theta, const algebra::vec<3,T>& uv); + + // only for unit quaternions described by theta and uv such as: + // q = ( cos(theta), sin(theta) * uv ) + + quat(unsigned one, float theta, const algebra::vec<3,float>& uv); + + float theta() const; + void set_theta(float theta); + + algebra::vec<3,float> uv() const; + void set_uv(const algebra::vec<3,float>& uv); + }; + + + // Operators. + + std::ostream& operator<<(std::ostream& ostr, const quat& q); + + quat operator+(const quat& lhs, const quat& rhs); + quat operator-(const quat& lhs, const quat& rhs); + quat operator*(const quat& lhs, const quat& rhs); + template <typename S> quat operator*(const quat& lhs, const value::scalar_<S>& rhs); + template <typename S> quat operator/(const quat& lhs, const value::scalar_<S>& rhs); + + // overloaded math procs + + quat log(const quat& q); + quat exp(const quat& q); + quat pow(const quat& q, double t); + template <typename T> + bool about_equal(const T& f, const T& q); + bool about_equal(const quat& p, const quat& q); + + + // Misc. + + bool interpol_ok(const quat& p, const quat& q, float h); + + + // Linear Quaternion Interpolation. + + quat lerp(const quat& p, const quat& q, float h); + + + // Spherical Linear Quaternion Interpolation. + + quat slerp(const quat& p, const quat& q, float h); + + quat slerp_2(const quat& p, const quat& q, float h); + + quat slerp_3(const quat& p, const quat& q, float h); + + quat slerp_4(const quat& p, const quat& q, float h); + + quat slerp_5(const quat& p, const quat& q, float h); + + + +# ifndef MLN_INCLUDE_ONLY + + // Constructors. + + inline + quat::quat() + { + } + + inline + quat::quat(float s, float x, float y, float z) + { + v_[0] = s; + set_v(x, y, z); + } + + inline + quat::quat(float s, const algebra::vec<3,float>& v) + { + v_[0] = s; + this->v() = v; + } + + inline + quat::quat(const algebra::vec<4,float>& v) + { + this->v_ = v; + } + + inline + quat& + quat::operator=(const algebra::vec<4,float>& v) + { + this->v_ = v; + return *this; + } + + + // With literals. + + inline + quat::quat(const literal::zero_t&) + { + v_.set_all(0); + } + + inline + quat& + quat::operator=(const literal::zero_t&) + { + v_.set_all(0); + return *this; + } + + inline + quat::quat(const literal::one_t&) + { + s() = 1; + v().set_all(0); + } + + inline + quat& + quat::operator=(const literal::one_t&) + { + s() = 1; + v().set_all(0); + return *this; + } + + + inline + const algebra::vec<4,float>& + quat::to_vec() const + { + return this->v_; + } + + inline + float + quat::s() const + { + return this->v_[0]; + } + + inline + float& + quat::s() + { + return this->v_[0]; + } + + inline + const algebra::vec<3, float>& + quat::v() const + { + return *(const algebra::vec<3, float>*)(const void*)(& this->v_[1]); + // return make::vec(this->v_[1], this->v_[2], this->v_[3]); + } + + inline + algebra::vec<3, float>& + quat::v() + { + return *(algebra::vec<3, float>*)(void*)(& this->v_[1]); + } + + inline + void quat::set_v(float x, float y, float z) + { + this->v_[1] = x; + this->v_[2] = y; + this->v_[3] = z; + } + + inline + float + quat::sprod(const quat& rhs) const + { + return v_ * rhs.to_vec(); + } + + inline + bool quat::is_unit() const + { + return about_equal(norm::l2(v_), 1.f); + } + + inline + bool quat::is_null() const + { + return about_equal(norm::l2(v_), 0.f); + } + + inline + bool quat::is_pure() const + { + return about_equal(v_[0], 0.f); + } + + inline + quat quat::conj() const + { + return quat(s(), - v()); + } + + inline + quat quat::inv() const + { + assert(! is_null()); + float f = norm::l2(v_); + return conj().to_vec() / (f * f); + } + + inline + quat& quat::set_unit() + { + v_.normalize(); + return *this; + } + + template <typename T> + inline + void quat::set_unit(float theta, const algebra::vec<3,T>& uv) + { + static const float pi = 3.14159265358979323846; + + mln_precondition(theta > - pi - mln_epsilon(float) + && theta < pi + mln_epsilon(float)); + mln_precondition(about_equal(norm::l2(uv), 1.f)); + + this->v_[0] = cos(theta); + float sint = sin(theta); + this->v_[1] = uv[0] * sint; + this->v_[2] = uv[1] * sint; + this->v_[3] = uv[2] * sint; + } + + // only for unit quaternions described by theta and uv such as: + // q = ( cos(theta), sin(theta) * uv ) + + inline + quat::quat(unsigned one, float theta, const algebra::vec<3,float>& uv) + { + mln_precondition(one == 1); + set_unit(theta, uv); + } + + inline + float quat::theta() const + { + mln_precondition(is_unit()); + return acos(s()); + } + + inline + void quat::set_theta(float theta) + { + mln_precondition(is_unit()); + set_unit(theta, uv()); + } + + inline + algebra::vec<3, float> quat::uv() const + { + mln_precondition(is_unit()); + algebra::vec<3, float> w = v(); + return w.normalize(); + } + + inline + void quat::set_uv(const algebra::vec<3,float>& uv) + { + mln_precondition(is_unit()); + set_unit(theta(), uv); + } + + + // Operators. + + inline + std::ostream& operator<<(std::ostream& ostr, const quat& q) + { + return ostr << q.to_vec(); + } + + inline + quat operator+(const quat& lhs, const quat& rhs) + { + quat tmp(lhs.to_vec() + rhs.to_vec()); + return tmp; + } + + inline + quat operator-(const quat& lhs, const quat& rhs) + { + quat tmp(lhs.to_vec() - rhs.to_vec()); + return tmp; + } + + inline + quat operator*(const quat& lhs, const quat& rhs) + { + quat tmp(lhs.s() * rhs.s() - lhs.v() * rhs.v(), + algebra::vprod(lhs.v(), rhs.v()) + lhs.s() * rhs.v() + rhs.s() * lhs.v()); + return tmp; + } + + template <typename S> + inline + quat operator*(const quat& lhs, const value::scalar_<S>& rhs) + { + mlc_converts_to(S, float)::check(); + quat tmp(lhs.to_vec() * rhs.to_equiv()); + return tmp; + } + + template <typename S> + inline + quat operator/(const quat& lhs, const value::scalar_<S>& rhs_) + { + mlc_converts_to(S, float)::check(); + float rhs = rhs_.to_equiv(); + mln_precondition(rhs != 0.f); + quat tmp(lhs.to_vec() / rhs); + return tmp; + } + + + // overloaded math procs + + inline + quat log(const quat& q) + { + mln_precondition(q.is_unit()); + return quat(0.f, q.theta() * q.uv()); + } + + + inline + quat exp(const quat& q) + { + mln_precondition(about_equal(q.s(), 0.f)); + algebra::vec<3, float> v = q.v(); + float theta = norm::l2(v); + mln_precondition(!about_equal(theta, 0.f)); + algebra::vec<3, float> uv = v / theta; + return quat(cos(theta), sin(theta) * uv); + } + + + inline + quat pow(const quat& q, double t) + { + mln_precondition(q.is_unit()); + return exp(t * log(q)); + } + + template <typename T> + inline + bool about_equal(const T& f, const T& q) + { + // FIXME: Use abs! + if (f > q) + return (f - q ) < mln_epsilon(T); + return (q - f) < mln_epsilon(T); + } + + inline + bool about_equal(const quat& p, const quat& q) + { + return about_equal<float>(norm::l2(p.to_vec() - q.to_vec()), 0); + } + + // Misc. + + inline + bool interpol_ok(const quat& p, const quat& q, float h) + { + return + p.is_unit() && + q.is_unit() && + h >= 0 && + h <= 1; + } + + + // Linear Quaternion Interpolation. + + inline + quat lerp(const quat& p, const quat& q, float h) + { + assert(interpol_ok(p, q, h)); + return (1 - h) * p + h * q; + } + + + // Spherical Linear Quaternion Interpolation. + + inline + quat slerp(const quat& p, const quat& q, float h) + { + assert(interpol_ok(p, q, h)); + float omega = acos(p.sprod(q)); + return + about_equal(omega, 0.f) ? + lerp(p, q, h) : + quat((sin((1-h)*omega) * p + sin(h*omega) * q) / sin(omega)); + } + + inline + quat slerp_2(const quat& p, const quat& q, float h) + { + assert(interpol_ok(p, q, h)); + quat tmp = p * value::pow(p.conj() * q, h); + assert(about_equal(tmp, slerp(p, q, h))); + return tmp; + } + + inline + quat slerp_3(const quat& p, const quat& q, float h) + { + assert(interpol_ok(p, q, h)); + quat tmp = value::pow(p * q.conj(), 1 - h) * q; + assert(about_equal(tmp, slerp(p, q, h))); + return tmp; + } + + inline + quat slerp_4(const quat& p, const quat& q, float h) + { + assert(interpol_ok(p, q, h)); + quat tmp = value::pow(q * p.conj(), h) * p; + assert(about_equal(tmp, slerp(p, q, h))); + return tmp; + } + + inline + quat slerp_5(const quat& p, const quat& q, float h) + { + assert(interpol_ok(p, q, h)); + quat tmp = q * value::pow(q.conj() * p, 1 - h); + assert(about_equal(tmp, slerp(p, q, h))); + return tmp; + } + + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::value + +} // end of namespace mln + +#endif // ! MLN_ALGEBRA_QUAT_HH Index: mln/norm/linfty.hh --- mln/norm/linfty.hh (revision 1783) +++ mln/norm/linfty.hh (working copy) @@ -36,7 +36,7 @@ */ # include <mln/math/abs.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> namespace mln @@ -51,7 +51,7 @@ C linfty(const C (&vec)[n]); template <unsigned n, typename C> - C linfty(const metal::vec<n,C>& vec); + C linfty(const algebra::vec<n,C>& vec); /// \} @@ -61,8 +61,8 @@ C linfty_distance(const C (&vec1)[n], const C (&vec2)[n]); template <unsigned n, typename C> - C linfty_distance(const metal::vec<n,C>& vec1, - const metal::vec<n,C>& vec2); + C linfty_distance(const algebra::vec<n,C>& vec1, + const algebra::vec<n,C>& vec2); /// \} @@ -119,7 +119,7 @@ template <unsigned n, typename C> inline - C linfty(const metal::vec<n,C>& vec) + C linfty(const algebra::vec<n,C>& vec) { return impl::linfty_<n, C>(vec); } @@ -133,8 +133,8 @@ template <unsigned n, typename C> inline - C linfty_distance(const metal::vec<n,C>& vec1, - const metal::vec<n,C>& vec2) + C linfty_distance(const algebra::vec<n,C>& vec1, + const algebra::vec<n,C>& vec2) { return impl::linfty_distance_<n, C>(vec1, vec2); } Index: mln/norm/l1.hh --- mln/norm/l1.hh (revision 1783) +++ mln/norm/l1.hh (working copy) @@ -35,7 +35,7 @@ */ # include <mln/math/abs.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> // FIXME: Use mln_sum_x (to be renamed as mln_sum_product) instead of // mln_sum. @@ -52,7 +52,7 @@ mln_sum(C) l1(const C (&vec)[n]); template <unsigned n, typename C> - mln_sum(C) l1(const metal::vec<n,C>& vec); + mln_sum(C) l1(const algebra::vec<n,C>& vec); /// \} /// L1-norm distance between vectors \a vec1 and \a vec2. @@ -61,8 +61,8 @@ mln_sum(C) l1_distance(const C (&vec1)[n], const C (&vec2)[n]); template <unsigned n, typename C> - mln_sum(C) l1_distance(const metal::vec<n,C>& vec1, - const metal::vec<n,C>& vec2); + mln_sum(C) l1_distance(const algebra::vec<n,C>& vec1, + const algebra::vec<n,C>& vec2); /// \} @@ -110,7 +110,7 @@ template <unsigned n, typename C> inline mln_sum(C) - l1(const metal::vec<n,C>& vec) + l1(const algebra::vec<n,C>& vec) { return impl::l1_<n, C>(vec); } @@ -126,7 +126,7 @@ template <unsigned n, typename C> inline mln_sum(C) - l1_distance(const metal::vec<n,C>& vec1, const metal::vec<n,C>& vec2) + l1_distance(const algebra::vec<n,C>& vec1, const algebra::vec<n,C>& vec2) { return impl::l1_distance_<n, C>(vec1, vec2); } Index: mln/norm/l2.hh --- mln/norm/l2.hh (revision 1783) +++ mln/norm/l2.hh (working copy) @@ -36,7 +36,7 @@ # include <mln/math/sqr.hh> # include <mln/math/sqrt.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> // FIXME: Use mln_sum_x (to be renamed as mln_sum_product) instead of @@ -54,7 +54,7 @@ mln_sum(C) l2(const C (&vec)[n]); template <unsigned n, typename C> - mln_sum(C) l2(const metal::vec<n,C>& vec); + mln_sum(C) l2(const algebra::vec<n,C>& vec); /// \} /// L2-norm distance between vectors \a vec1 and \p vec2. @@ -63,8 +63,8 @@ mln_sum(C) l2_distance(const C (&vec1)[n], const C (&vec2)[n]); template <unsigned n, typename C> - mln_sum(C) l2_distance(const metal::vec<n,C>& vec1, - const metal::vec<n,C>& vec2); + mln_sum(C) l2_distance(const algebra::vec<n,C>& vec1, + const algebra::vec<n,C>& vec2); /// \} @@ -113,7 +113,7 @@ template <unsigned n, typename C> inline mln_sum(C) - l2(const metal::vec<n,C>& vec) + l2(const algebra::vec<n,C>& vec) { return impl::l2_<n, C>(vec); } @@ -129,7 +129,7 @@ template <unsigned n, typename C> inline mln_sum(C) - l2_distance(const metal::vec<n,C>& vec1, const metal::vec<n,C>& vec2) + l2_distance(const algebra::vec<n,C>& vec1, const algebra::vec<n,C>& vec2) { return impl::l2_distance_<n, C>(vec1, vec2); } Index: sandbox/duhamel/slow_seed2tiling.cc --- sandbox/duhamel/slow_seed2tiling.cc (revision 1783) +++ sandbox/duhamel/slow_seed2tiling.cc (working copy) @@ -60,7 +60,7 @@ # include <mln/core/clone.hh> # include <mln/accu/mean.hh> # include <mln/estim/min_max.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/norm/l2.hh> # include <mln/norm/l1.hh> Index: sandbox/duhamel/labeling_algo.hh --- sandbox/duhamel/labeling_algo.hh (revision 1783) +++ sandbox/duhamel/labeling_algo.hh (working copy) @@ -17,7 +17,7 @@ #include <algorithm> #include <utility> #include <map> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> namespace mln { @@ -64,7 +64,7 @@ graph_with_no_border (Image<I>& ima_, const Neighborhood<N>& nbh) { - typedef metal::vec<2,float> X; + typedef algebra::vec<2,float> X; typedef mln_value(I) V; typedef mln_psite(I) P; @@ -111,7 +111,7 @@ Image<I>& orig_, const Neighborhood<N>& nbh) { - typedef metal::vec<2,float> X; + typedef algebra::vec<2,float> X; typedef mln_value(I) V; typedef mln_psite(I) P; Index: sandbox/duhamel/mesh_image.hh --- sandbox/duhamel/mesh_image.hh (revision 1783) +++ sandbox/duhamel/mesh_image.hh (working copy) @@ -36,7 +36,7 @@ # include <cmath> # include <mln/core/internal/image_identity.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include "mesh_p.hh" # include "mesh_psite.hh" # include <vector> Index: sandbox/jardonnet/test/icp.cc --- sandbox/jardonnet/test/icp.cc (revision 0) +++ sandbox/jardonnet/test/icp.cc (revision 0) @@ -0,0 +1,16 @@ +#include <mln/core/image2d.hh> + +#include <mln/io/ppm/load.hh> +#include <mln/io/ppm/save.hh> + +#include <sandbox/jardonnet/registration/icp.hh> + +int main(int, char*) +{ + using namespace mln; + + image2d< value::rgb8 > img; + + registration::icp(img,img); +} + Index: sandbox/jardonnet/test/gaussian_subsampling.cc --- sandbox/jardonnet/test/gaussian_subsampling.cc (revision 1783) +++ sandbox/jardonnet/test/gaussian_subsampling.cc (working copy) @@ -15,7 +15,7 @@ io::pgm::load(img, "../../../img/lena.pgm"); - image2d<int_u8> output = subsampling::gaussian_subsampling(img, make::dpoint2d(1,1), argc); + image2d<int_u8> output = subsampling::gaussian_subsampling(img, 0.1, make::dpoint2d(1,1), argc); io::pgm::save(output, "gsub.pgm"); } Index: sandbox/jardonnet/test/Makefile --- sandbox/jardonnet/test/Makefile (revision 1783) +++ sandbox/jardonnet/test/Makefile (working copy) @@ -1,7 +1,7 @@ SOURCE=test.cc subsampling.cc -FLAG=-Wall -W -I../../.. -g +FLAG=-Wall -W -pedantic -I../../.. -g -all: sub gsub +all: sub gsub gau icp sub: @@ -13,5 +13,9 @@ gau: g++ gaussian.cc $(FLAG) -o '+gau.exe' +icp: + g++ icp.cc $(FLAG) -o '+icp.exe' + run: time ./+sub.exe . . ; time ./+gsub.exe . . \ No newline at end of file + Index: sandbox/jardonnet/subsampling/gaussian_subsampling.hh --- sandbox/jardonnet/subsampling/gaussian_subsampling.hh (revision 1783) +++ sandbox/jardonnet/subsampling/gaussian_subsampling.hh (working copy) @@ -55,7 +55,7 @@ template <typename I> inline mln_concrete(I) - gaussian_subsampling(const Image<I>& input, float sigma + gaussian_subsampling(const Image<I>& input, float sigma, const mln_dpoint(I)& first_p, const mln_coord(I)& gap); @@ -76,7 +76,7 @@ mln_concrete(I) output(geom::nrows(input) / gap, geom::ncols(input) / gap); //FIXME : only for image2d. - linear::gaussian(input, 0.1, temp); + linear::gaussian(input, sigma, temp); output = impl::subsampling_(exact(temp), first_p, gap); trace::exiting("subsampling::gaussian_subsampling"); Index: sandbox/jardonnet/subsampling/sub_sampled_image.hh --- sandbox/jardonnet/subsampling/sub_sampled_image.hh (revision 1783) +++ sandbox/jardonnet/subsampling/sub_sampled_image.hh (working copy) @@ -34,7 +34,7 @@ # include <mln/core/internal/image_morpher.hh> # include <mln/convert/to_dpoint.hh> -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> # include <mln/value/set.hh> @@ -214,7 +214,7 @@ sub_sampled_image<I>::translate_coords_(const mln_point(I)& p) const { - return mln_point(I)(metal::vec<2, int>(p) * gap + metal::vec<2, int>(first_p)); + return mln_point(I)(algebra::vec<2, int>(p) * gap + algebra::vec<2, int>(first_p)); } Index: sandbox/jardonnet/TODO --- sandbox/jardonnet/TODO (revision 1783) +++ sandbox/jardonnet/TODO (working copy) @@ -1,10 +1,8 @@ image2d< value::rgb8 > img == out - -- - - - - - - - - - +- gaussian.cc: In function 'int main(int, char*)': gaussian.cc:22: error: no match for 'operator==' in 'img == out' +- - - -*/*/*/*/*/*/*/* \ No newline at end of file Index: sandbox/jardonnet/registration/quat7.hh --- sandbox/jardonnet/registration/quat7.hh (revision 0) +++ sandbox/jardonnet/registration/quat7.hh (revision 0) @@ -0,0 +1,148 @@ + +#ifndef QUAT7_HH +# define QUAT7_HH + +# include <assert.h> +# include <algorithm> + +# include <mln/algebra/mat.hh> + +# include "quat/all.hh" +# include "jacobi.hh" + +# include <mln/norm/l2.hh> + +# include <mln/trait/all.hh> + + +// FIXME : move elsewhere +namespace mln +{ + namespace algebra + { + + template<unsigned n, unsigned m, typename T> + mat<m,n,T> + trans(const mat<n,m,T>& matrice) + { + mat<m,n,T> tmp; + for (unsigned i = 0; i < n; ++i) + for (unsigned j = 0; j < m; ++j) + tmp(j,i) = matrice(i,j); + return tmp; + } + + // trace + + template<unsigned n, typename T> inline + float tr(const mat<n,n,T>& m) + { + float f = 0.f; + for (unsigned i = 0; i < n; ++i) + f += m(i,i); + return f; + } + + } +} + +namespace mln +{ + + struct quat7 + { + value::quat _qR; + vec3f _qT; + + quat7() + { + } + + quat7(const value::quat& qR, const vec3f& qT) : + _qR(qR), + _qT(qT) + { + } + + float sqr_norm() const + { + return norm::l2(_qR.to_vec()) + norm::l2(_qT); + } + + quat7 operator-(const quat7& rhs) const + { + return quat7(_qR - rhs._qR, _qT - rhs._qT); + } + + // quat7 is an object-function + + vec3f operator()(const vec3f& v) const + { + return rotate(_qR, v) + _qT; + } + + void apply_on(const std::vector< vec3f >& input, std::vector< vec3f >& output) const + { + assert(input.size() == output.size()); + assert(_qR.is_unit()); + + std::transform(input.begin(), input.end(), + output.begin(), + *this); + } + }; + + + // very usefull routine + + template <unsigned p, unsigned q, unsigned n, unsigned m> + void put(const algebra::mat<p,q,float>& in, // a matrix to put into... + unsigned row, unsigned col, // top-left location + algebra::mat<n,m,float>& inout) // ...a matrix to modify + { + assert(row + p <= n && col + q <= m); + for (unsigned i = 0; i < p; ++i) + for (unsigned j = 0; j < q; ++j) + inout(row + i, col + j) = in(i,j); + } + + + quat7 match(const vecs_t& P, + const vec3f& mu_P, + const vecs_t& Xk, + const vec3f& mu_Xk) + { + assert(P.size() == Xk.size()); + + // qR + + algebra::mat<3,3,float> Ck; + for (unsigned i = 0; i < P.size(); ++i) + Ck += make::mat(P[i] - mu_P) * trans(make::mat(Xk[i] - mu_Xk)); + Ck /= P.size(); + + const algebra::mat<3,3,float> Ak = Ck - trans(Ck); + + const float v[3] = {Ak(1,2), Ak(2,0), Ak(0,1)}; + const algebra::mat<3,1,float> D = make::mat<3,1,3,float>(v); // FIXME why <...> + + algebra::mat<4,4,float> Qk; + Qk(0,0) = tr(Ck); + put(trans(D), 0,1, Qk); + put(D, 1,0, Qk); + + put(Ck + trans(Ck) - algebra::mat<3,3,float>::identity() * tr(Ck), 1,1, Qk); + + value::quat qR; + jacobi(Qk, qR); + + // qT + + const vec3f qT = mu_Xk - rotate(qR, mu_P); + + return quat7(qR, qT); + } + +} + +#endif // ndef QUAT7_HH Index: sandbox/jardonnet/registration/cloud.hh --- sandbox/jardonnet/registration/cloud.hh (revision 0) +++ sandbox/jardonnet/registration/cloud.hh (revision 0) @@ -0,0 +1,49 @@ +#ifndef CLOUD_HH +# define CLOUD_HH + +# include <assert.h> +# include <string> +# include <vector> +# include <fstream> +# include <sstream> + +# include <mln/algebra/vec.hh> + +namespace mln +{ + + namespace registration + { + + typedef algebra::vec<3, float> vec3f; + + + vec3f center(const std::vector< vec3f >& vecs) + { + vec3f c; + for (size_t i = 0; i < vecs.size(); ++i) + c += vecs[i]; + return c / vecs.size(); + } + + + // FIXME : move + float sqr_norm(const vec3f& v) + { + return v[1] * v[1] + v[2] * v[2] + v[3] * v[3]; + } + + float rms(const std::vector< vec3f >& vecs1, + const std::vector< vec3f >& vecs2) + { + assert(vecs1.size() == vecs2.size()); + float f = 0.f; + for (size_t i = 0; i < vecs1.size(); ++i) + f += sqr_norm(vecs1[i] - vecs2[i]); + return f / vecs1.size(); + } + + } +} + +#endif // ndef CLOUD_HH Index: sandbox/jardonnet/registration/jacobi.hh --- sandbox/jardonnet/registration/jacobi.hh (revision 0) +++ sandbox/jardonnet/registration/jacobi.hh (revision 0) @@ -0,0 +1,108 @@ + +#ifndef JACOBI_HH +# define JACOBI_HH + + +#include <mln/algebra/mat.hh> + +// from num. rec. in C + + +#define rotateJacobi(a,i,j,k,l) g=a(i,j);h=a(k,l);a(i,j)=g-s*(h+g*tau); \ + a(k,l)=h+s*(g-h*tau); + +void jacobi(mln::algebra::mat<4,4,float> a, mln::value::quat& q) +{ + float dd, d[4]; + mln::algebra::mat<4,4,float> v; + int j,iq,ip,i = 0; + float tresh,theta,tau,t,sm,s,h,g,c,b[4],z[4]; + for (ip=0;ip<4;ip++) { + for (iq=0;iq<4;iq++) v(ip,iq)=0.0; + v(ip,ip)=1.0; + } + for (ip=0;ip<4;ip++) { + b[ip]=d[ip]=a(ip,ip); + z[ip]=0.0; + } + while (1) { + sm=0.0; + for (ip=0;ip<3;ip++) { + for (iq=ip+1;iq<4;iq++) + sm += fabs(a(ip,iq)); + } + if (sm < 1e-12) { + dd = d[0]; + iq = 0; + for (ip=1;ip<4;ip++) + if (d[ip]>dd) { + iq = ip; + dd = d[ip]; + } + q = mln::value::quat(v(0,iq), + v(1,iq), + v(2,iq), + v(3,iq)); + q.set_unit(); + return; + } + if (i < 4) { + i++; + tresh=0.0125*sm; + } + else + tresh=0.0; + for (ip=0;ip<3;ip++) { + for (iq=ip+1;iq<4;iq++) { + + /* unusefull */ + g=100.0*fabs(a(ip,iq)); + if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) + && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) + a(ip,iq)=0.0; + /* unusefull */ + + else if (fabs(a(ip,iq)) > tresh) { + h=d[iq]-d[ip]; + if ((float)(fabs(h)+g) == (float)fabs(h)) + t=(a(ip,iq))/h; + else { + theta=0.5*h/(a(ip,iq)); + t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); + if (theta < 0.0) t = -t; + } + c=1.0/sqrt(1+t*t); + s=t*c; + tau=s/(1.0+c); + h=t*a(ip,iq); + z[ip] -= h; + z[iq] += h; + d[ip] -= h; + d[iq] += h; + a(ip,iq)=0.0; + for (j=0;j<=ip-1;j++) { + rotateJacobi(a,j,ip,j,iq) + } + for (j=ip+1;j<=iq-1;j++) { + rotateJacobi(a,ip,j,j,iq) + } + for (j=iq+1;j<4;j++) { + rotateJacobi(a,ip,j,iq,j) + } + for (j=0;j<4;j++) { + rotateJacobi(v,j,ip,j,iq) + } + } + } + } + for (ip=0;ip<4;ip++) { + b[ip] += z[ip]; + d[ip]=b[ip]; + z[ip]=0.0; + } + } +} + + + +#endif // ndef JACOBI_HH Index: sandbox/jardonnet/registration/icp.hh --- sandbox/jardonnet/registration/icp.hh (revision 1783) +++ sandbox/jardonnet/registration/icp.hh (working copy) @@ -25,17 +25,25 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_REGISTRATION_REGISTRATION_HH -# define MLN_REGISTRATION_REGISTRATION_HH +#ifndef MLN_REGISTRATION_ICP_HH +# define MLN_REGISTRATION_ICP_HH -/*! \file mln/binarization/threshold.hh +/*! \file mln/registration/icp.hh * - * \brief Produce a subsampled image + * \brief image registration */ # include <mln/value/quat.hh> # include <mln/algebra/vec.hh> + +typedef mln::algebra::vec<3, float> vec3f; +typedef std::vector< vec3f > vecs_t; + +#include "cloud.hh" +#include "quat7.hh" +#include "projection.hh" + namespace mln { @@ -46,7 +54,7 @@ * * */ - template <typename I, template J> + template <typename I, typename J> inline void icp(const Image<I>& cloud, @@ -57,34 +65,43 @@ namespace impl { - template <typename I, typename J> inline void - icp_(const I& P, - const J& X) + icp_(const vecs_t& P, + const vecs_t& X) { trace::entering("registration::impl::icp_"); - mln_concrete(I) Pk(cloud.domain()); + unsigned int k; + quat7 old_qk, qk; + float err, err_bis; + + vecs_t Pk(P.size()), Xk(Pk.size()); + vec3f mu_P = center(P), mu_Xk; const float epsilon = 1e-3; - float err, err_bis; - quat old_qk, qk; - unsigned int k; + //step 1 k = 0; Pk = P; + do { - //projection + //step 2 FIXME : etienne + projection::de_base(Pk, X, Xk, mu_Xk, err_bis); + // step 3 old_qk = qk; - //qk = match(P, mu_P, Xk, mu_Xk); + qk = match(P, mu_P, Xk, mu_Xk); + + // step 4 + qk.apply_on(P, Pk); // Pk+1 = qk(P) + + // err = d(Pk+1,Xk) + err = rms(Pk, Xk); - // error = ++k; } while (k < 3 || (qk - old_qk).sqr_norm() > epsilon); - trace::exiting("registration::impl::icp_"); } @@ -102,9 +119,9 @@ mln_precondition(exact(cloud).has_data()); mln_precondition(exact(surface).has_data()); + vecs_t a,b; // FIXME : to built. - - output = impl::icp_(exact(cloud), exact(surface)); + impl::icp_(a, b); trace::exiting("registration::icp"); } Index: sandbox/jardonnet/registration/projection.hh --- sandbox/jardonnet/registration/projection.hh (revision 0) +++ sandbox/jardonnet/registration/projection.hh (revision 0) @@ -0,0 +1,52 @@ + +#ifndef PROJECTION_HH +# define PROJECTION_HH + +# include <assert.h> + +namespace mln +{ + + namespace projection + { + + template <class Pk_t, class X_t, class Xk_t> + void de_base(// input + const Pk_t& Pk, + const X_t& X, + // inout + Xk_t& Xk, + // output + algebra::vec<3, float>& mu_Xk, + float& err) + { + assert(Pk.size() == Xk.size()); + + err = 0.f; + mu_Xk = make::vec(0,0,0); + + for (size_t i = 0; i < Pk.size(); ++i) + { + algebra::vec<3,float> best_x = X[0]; + float best_d = norm::l2(Pk[i] - best_x); + for (size_t j = 1; j < X.size(); ++j) + { + float d = norm::l2(Pk[i] - X[j]); + if (d < best_d) + { + best_d = d; + best_x = X[j]; + } + } + Xk[i] = best_x; + mu_Xk += Xk[i]; + err += best_d; + } + mu_Xk /= Pk.size(); + err /= Pk.size(); + } + + } +} + +#endif // ndef PROJECTION_HH Index: sandbox/jardonnet/registration/quat/all.hh --- sandbox/jardonnet/registration/quat/all.hh (revision 0) +++ sandbox/jardonnet/registration/quat/all.hh (revision 0) @@ -0,0 +1,9 @@ +#ifndef QUAT_ALL_HH +# define QUAT_ALL_HH + + +# include "interpol.hh" +# include "rotation.hh" + + +#endif // ndef QUAT_ALL_HH Index: sandbox/jardonnet/registration/quat/misc.hh --- sandbox/jardonnet/registration/quat/misc.hh (revision 0) +++ sandbox/jardonnet/registration/quat/misc.hh (revision 0) @@ -0,0 +1,19 @@ +#ifndef MISC_HH +# define MISC_HH + + +inline +float epsilon() +{ + static const float e = 1e-5; + return e; +} + +inline +bool about_equal(float val1, float val2) +{ + return fabs(val1 - val2) < epsilon(); +} + + +#endif // ndef MISC_HH Index: sandbox/jardonnet/registration/quat/interpol.hh --- sandbox/jardonnet/registration/quat/interpol.hh (revision 0) +++ sandbox/jardonnet/registration/quat/interpol.hh (revision 0) @@ -0,0 +1,100 @@ +#ifndef INTERPOL_HH +# define INTERPOL_HH + +# include <vector> +# include "misc.hh" + +# include <mln/value/quat.hh> + +// misc + +namespace mln +{ + + inline + value::quat slerp_2(const value::quat& p, const value::quat& q, float h) + { + assert(interpol_ok(p, q, h)); + value::quat tmp = p * pow(p.conj() * q, h); + assert(about_equal(tmp, slerp(p, q, h))); + return tmp; + } + + inline + value::quat slerp_3(const value::quat& p, const value::quat& q, float h) + { + assert(interpol_ok(p, q, h)); + value::quat tmp = pow(p * q.conj(), 1 - h) * q; + assert(about_equal(tmp, slerp(p, q, h))); + return tmp; + } + + inline + value::quat slerp_4(const value::quat& p, const value::quat& q, float h) + { + assert(interpol_ok(p, q, h)); + value::quat tmp = pow(q * p.conj(), h) * p; + assert(about_equal(tmp, slerp(p, q, h))); + return tmp; + } + + inline + value::quat slerp_5(const value::quat& p, const value::quat& q, float h) + { + assert(interpol_ok(p, q, h)); + value::quat tmp = q * pow(q.conj() * p, 1 - h); + assert(about_equal(tmp, slerp(p, q, h))); + return tmp; + } + + + // Spherical Spline Value::Quaternion Interpolation + + inline + value::quat squad(const std::vector<value::quat>& q, + const std::vector<value::quat>& s, + unsigned i, + float h) + { + assert(i < q.size() && s.size() == q.size()); + return slerp(slerp(q[i], q[i+1], h), + slerp(s[i], s[i+1], h), + 2 * h * (1 - h)); + } + + + inline + value::quat s_from_q(const std::vector<value::quat>& q, unsigned i) + { + assert(i < q.size()); + assert(q[i].is_unit()); + if (i == 0 || i == q.size() - 1) + return q[i]; + assert(q[i-1].is_unit() && q[i+1].is_unit()); + return q[i] * exp(-(log(q[i].inv() * q[i+1]) + log(q[i].inv() * q[i-1])) + / 4); + } + + + inline + std::vector<value::quat> do_squad(const std::vector<value::quat>& q, + unsigned n) + { + assert(n > 1); + std::vector<value::quat> s(q.size()), res; + + for (unsigned i = 0; i < q.size(); ++i) + s[i] = s_from_q(q, i); + + res.push_back(q[0]); + for (unsigned i = 0; i < q.size() - 1; ++i) + { + for (unsigned j = 1; j < n; ++j) + res.push_back(squad(q, s, i, float(j) / n)); + res.push_back(q[i + 1]); + } + return res; + } +} + +#endif // ndef INTERPOL_HH Index: sandbox/jardonnet/registration/quat/rotation.hh --- sandbox/jardonnet/registration/quat/rotation.hh (revision 0) +++ sandbox/jardonnet/registration/quat/rotation.hh (revision 0) @@ -0,0 +1,60 @@ +#ifndef ROTATION_HH +# define ROTATION_HH + +# include <stdlib.h> +//# include "quat.hh" + +//# include "mat.hh" + +# include <mln/algebra/mat.hh> +# include <mln/algebra/vec.hh> +# include <mln/make/vec.hh> +# include <mln/make/mat.hh> +# include <mln/value/quat.hh> + + +// FIXME: rotation should be an abstract class +// and derived classes encapsulate either a quaternion or a algebra::matrix +namespace mln +{ + + vec3f rotate(const value::quat& q, const vec3f& p) + { + + return (q * value::quat(0. ,p) * q.inv()).v(); + } + + + bool check_rotation(const algebra::mat<3,3,float>& mat, + const value::quat& q) + { + assert(q.is_unit()); + vec3f + tmp = make::vec(rand(), rand(), rand()), + p = tmp / norm::l2(tmp), + p_rot_1 = rotate(q, p), + p_rot_2 = mat * p; + return about_equal(norm::l2(p_rot_1 - p_rot_2), 0.f); + } + + + algebra::mat<3,3,float> quat2mat(const value::quat& q) + { + assert(q.is_unit()); + float + w = q.to_vec()[0], + x = q.to_vec()[1], x2 = 2*x*x, xw = 2*x*w, + y = q.to_vec()[2], y2 = 2*y*y, xy = 2*x*y, yw = 2*y*w, + z = q.to_vec()[3], z2 = 2*z*z, xz = 2*x*z, yz = 2*y*z, zw = 2*z*w; + float t[9] = {1.f - y2 - z2, xy - zw, xz + yw, + xy + zw, 1.f - x2 - z2, yz - xw, + xz - yw, yz + xw, 1.f - x2 - y2}; + algebra::mat<3,3,float> tmp = make::mat<3,3,9,float>(t); + // postcondition + assert(check_rotation(tmp, q)); + return tmp; + } +} + + +#endif // ndef ROTATION_HH Index: sandbox/vigouroux/color/my_yuv.hh --- sandbox/vigouroux/color/my_yuv.hh (revision 1783) +++ sandbox/vigouroux/color/my_yuv.hh (working copy) @@ -2,7 +2,7 @@ // #include <mln/value/concept/vectorial.hh> // #include <mln/value/int_u.hh> -// #include <mln/metal/vec.hh> +// #include <mln/algebra/vec.hh> #ifndef MLN_VALUE_YUV_HH # define MLN_VALUE_YUV_HH Index: sandbox/vigouroux/color/my_hsi.hh --- sandbox/vigouroux/color/my_hsi.hh (revision 1783) +++ sandbox/vigouroux/color/my_hsi.hh (working copy) @@ -2,7 +2,7 @@ #include <mln/value/concept/vectorial.hh> #include <mln/value/int_u.hh> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/value/float01_8.hh> Index: sandbox/vigouroux/color/my_hsl.hh --- sandbox/vigouroux/color/my_hsl.hh (revision 1783) +++ sandbox/vigouroux/color/my_hsl.hh (working copy) @@ -2,7 +2,7 @@ #include <mln/value/concept/vectorial.hh> #include <mln/value/int_u.hh> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #ifndef MLN_VALUE_HSL_HH # define MLN_VALUE_HSL_HH Index: sandbox/vigouroux/color/my_xyz.hh --- sandbox/vigouroux/color/my_xyz.hh (revision 1783) +++ sandbox/vigouroux/color/my_xyz.hh (working copy) @@ -2,7 +2,7 @@ // #include <mln/value/concept/vectorial.hh> // #include <mln/value/int_u.hh> -// #include <mln/metal/vec.hh> +// #include <mln/algebra/vec.hh> #ifndef MLN_VALUE_XYZ_HH # define MLN_VALUE_XYZ_HH Index: sandbox/vigouroux/color/my_hsv.hh --- sandbox/vigouroux/color/my_hsv.hh (revision 1783) +++ sandbox/vigouroux/color/my_hsv.hh (working copy) @@ -2,7 +2,7 @@ #include <mln/value/concept/vectorial.hh> #include <mln/value/int_u.hh> -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #ifndef MLN_VALUE_HSV_HH # define MLN_VALUE_HSV_HH Index: sandbox/vigouroux/color/my_yiq.hh --- sandbox/vigouroux/color/my_yiq.hh (revision 1783) +++ sandbox/vigouroux/color/my_yiq.hh (working copy) @@ -2,7 +2,7 @@ // #include <mln/value/concept/vectorial.hh> // #include <mln/value/int_u.hh> -// #include <mln/metal/vec.hh> +// #include <mln/algebra/vec.hh> #ifndef MLN_VALUE_YIQ_HH # define MLN_VALUE_YIQ_HH Index: sandbox/folio/dmap.cc --- sandbox/folio/dmap.cc (revision 1783) +++ sandbox/folio/dmap.cc (working copy) @@ -34,7 +34,7 @@ for_all(q) if (input(q) == true) // q is in the object. { - metal::vec<2,int> vp = p.to_point(), vq = q.to_point(); + algebra::vec<2,int> vp = p.to_point(), vq = q.to_point(); min.take(norm::l2_distance(vp, vq)); } output(p) = min; Index: sandbox/garrigues/image_identity/interpolated.hh --- sandbox/garrigues/image_identity/interpolated.hh (revision 1783) +++ sandbox/garrigues/image_identity/interpolated.hh (working copy) @@ -36,7 +36,7 @@ # include <cmath> # include "image_identity.hh" -# include <mln/metal/vec.hh> +# include <mln/algebra/vec.hh> namespace mln @@ -99,13 +99,13 @@ using super_::owns_; /// Test if a pixel value is accessible at \p v. - bool owns_(const mln::metal::vec<I::point::dim, float>& v) const; + bool owns_(const mln::algebra::vec<I::point::dim, float>& v) const; /// Read-only access of pixel value at point site \p p. /// Mutable access is only OK for reading (not writing). using super_::operator(); - mln_value(I) operator()(const mln::metal::vec<I::point::dim, float>& v) const; + mln_value(I) operator()(const mln::algebra::vec<I::point::dim, float>& v) const; /// Give the set of values of the image. @@ -149,7 +149,7 @@ } template <typename I> - bool interpolated<I>::owns_(const mln::metal::vec<I::point::dim, float>& v) const + bool interpolated<I>::owns_(const mln::algebra::vec<I::point::dim, float>& v) const { mln_point(I) p; for (unsigned i = 0; i < I::point::dim; ++i) @@ -159,7 +159,7 @@ template <typename I> mln_value(I) - interpolated<I>::operator()(const mln::metal::vec<I::point::dim, float>& v) const + interpolated<I>::operator()(const mln::algebra::vec<I::point::dim, float>& v) const { mln_point(I) p; for (unsigned i = 0; i < I::point::dim; ++i) Index: sandbox/garrigues/image_identity/interpolated.cc --- sandbox/garrigues/image_identity/interpolated.cc (revision 1783) +++ sandbox/garrigues/image_identity/interpolated.cc (working copy) @@ -35,7 +35,7 @@ #include <mln/core/image2d_b.hh> #include "interpolated.hh" -#include <mln/metal/vec.hh> +#include <mln/algebra/vec.hh> #include <mln/level/fill.hh> @@ -60,8 +60,8 @@ interpolated< image2d_b<float> > inter(f); - metal::vec<2, float> v1 = make::vec(2.3, 0.6); - metal::vec<2, float> v2 = make::vec(3.2, 1.8); + algebra::vec<2, float> v1 = make::vec(2.3, 0.6); + algebra::vec<2, float> v2 = make::vec(3.2, 1.8); debug::println(f);