[Olena-patches] cleanup-2008 2475: Start to integrate registration in Milena.

https://svn.lrde.epita.fr/svn/oln/branches/cleanup-2008/milena Index: ChangeLog from Ugo Jardonnet &lt;ugo.jardonnet(a)lrde.epita.fr&gt; Start to integrate registration in Milena. * mln/binarization/binarization.hh: Remove whitesapce. * mln/registration: New. * mln/registration/multiscale.hh: Multiscale registration of tow images. * mln/registration/registration.hh: Simple registration of tow images. * mln/registration/icp.hh: ICP registration of two p_array. icp.hh | 159 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ multiscale.hh | 114 ++++++++++++++++++++++++++++++++++++++++ registration.hh | 119 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 392 insertions(+) Index: mln/binarization/binarization.hh Index: mln/registration/multiscale.hh --- mln/registration/multiscale.hh (revision 0) +++ mln/registration/multiscale.hh (revision 0) @@ -0,0 +1,114 @@ +// Copyright (C) 2008 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_REGISTRATION_MULTISCALE_HH +# define MLN_REGISTRATION_MULTISCALE_HH + +# include <mln/registration/icp.hh> + +namespace mln +{ + + namespace registration + { + + template <typename I, typename J> + inline + composed< rotation<I::psite::dim, float>, translation<I::Psite::dim, float> > + multiscale(const Image<I>& cloud, + const Image<J>& surface, + const float q, + const unsigned nb_it); + + +# ifndef MLN_INCLUDE_ONLY + + namespace impl + { + + template <typename I, typename J> + inline + composed< rotation<P::dim, float>, translation<P::dim, float> > + multiscale_(const Image<I>& cloud, + const Image<J>& surface, + const float q, + const unsigned nb_it) + { + // Shuffle cloud + shuffle(cloud); + + //working box + const box<point2d> working_box = + enlarge(bigger(geom::bbox(c), geom::bbox(x)), 100); + + //make a lazy_image map via function closest_point + closest_point<mln_psite_(image2db)> fun(x, working_box); //FIXME: to port + lazy_image<image2d<bool>, closest_point<mln_psite_(image2db)>, box2d > + map(fun, fun.domain()); + + //init rigid transform qk + composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > qk; + + //run registration + for (int e = nb_it-1; e >= 0; e--) + { + unsigned l = cloud.nsites() / std::pow(q, e); + l = (l<1) ? 1 : l; + impl::registration_(cloud, map, qk, l, 1e-3); + } + return qk; + } + + } + + template <typename I, typename J> + inline + composed< rotation<P::dim, float>, translation<P::dim, float> > + multiscale(const Image<I>& cloud, + const Image<J>& surface, + const float q, + const unsigned nb_it) + { + trace::entering("registration::registration"); + + mln_precondition(P::dim == 3 || P::dim == 2); + + impl::multiscale_(cloud, surface, q, nb_it); + + trace::exiting("registration::registration"); + + return qk; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::registration + +} // end of namespace mln + + +#endif // ! MLN_REGISTRATION_MULTISCALE_HH Index: mln/registration/registration.hh --- mln/registration/registration.hh (revision 0) +++ mln/registration/registration.hh (revision 0) @@ -0,0 +1,119 @@ +// Copyright (C) 2008 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_REGISTRATION_REGISTRATION_HH +# define MLN_REGISTRATION_REGISTRATION_HH + +/*! \file mln/registration/registration.hh + * + * \brief image registration + */ + +# include <mln/registration/icp.hh> +# include <mln/core/image/lazy_image.hh> + +namespace mln +{ + + namespace registration + { + + + using namespace fun::x2x::geom; + + /*! Register an image \p cloud over the image \p surface. + */ + template <typename I, typename J> + inline + composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > + registration(const Image<I>& cloud, + const Image<J>& surface); + + +# ifndef MLN_INCLUDE_ONLY + + namespace impl + { + + template <typename I, typename J> + inline + composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > + registration_(const Image<I>& cloud, + const Image<J>& surface) + { + p_array<mln_psite(I)> c = convert::to< p_array<mln_psite(I)> >(cloud); + p_array<mln_psite(J)> x = convert::to< p_array<mln_psite(I)> >(surface); + + //init rigid transform qk + composed< rotation<I::P::dim, float>, translation<I::P::dim, float> > qk; + + //working box + const box<point2d> working_box = + enlarge(bigger(geom::bbox(c), geom::bbox(x)), 100); + + //make a lazy_image map via function closest_point + closest_point<mln_psite(I)> fun(x, working_box); //FIXME: to port + lazy_image<mln_ch_value(I,bool), closest_point<mln_psite(I)>, box2d > + map(fun, fun.domain()); + + //run registration + registration::icp(c, map, qk, l, 1e-3); + + } + + } + + template <typename I, typename J> + inline + composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > + registration(const Image<I>& cloud, + const Image<J>& surface) + { + trace::entering("registration::registration"); + + mln_precondition(I::psite::dim == J::psite::dim); + mln_precondition(I::psite::dim == 3 || J::psite::dim == 2); + + composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > + qk = impl::registration_(cloud, surface); + + trace::exiting("registration::registration"); + + return qk; + } + + +# endif // ! MLN_INCLUDE_ONLY + + + } // end of namespace mln::registration + + +} // end of namespace mln + + +#endif // ! MLN_REGISTRATION_REGISTRATION_HH Index: mln/registration/icp.hh --- mln/registration/icp.hh (revision 0) +++ mln/registration/icp.hh (revision 0) @@ -0,0 +1,159 @@ +// Copyright (C) 2008 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_REGISTRATION_ICP_HH +# define MLN_REGISTRATION_ICP_HH + +/*! \file mln/registration/icp.hh + * + * \brief Register an image over an another using the ICP algorithm. + */ + +#include <mln/fun/x2x/all.hh> + +namespace mln +{ + + namespace registration + { + + using namespace fun::x2x::geom; + + /*! Register point in \p c using a map of closest points \p map + * + * \param[in] c The cloud of points. + * \param[in] map The map of closest points. + * \param[in] epsilon ICP stops if sqr_norm(qk - qk-1) / + * sqr_norm(qk) > epsilon + * \param[out] qk The rigid transformation obtained. + * + * \pre \p ima has to be initialized. + */ + template <typename P, typename M> + inline + composed<rotation<2u, float>, translation<2u, float> > + icp(const p_array<P>& c, const M& map, + const float epsilon = 1e-3); + + /*! + * \fixme + */ + template <typename P, typename M, typename T> + inline + void + icp_subset(const p_array<P>& c, + const M& map, + T& qk, + const unsigned c_length, + const float epsilon = 1e-3); + + +# ifndef MLN_INCLUDE_ONLY + + namespace impl + { + + template <typename P, typename M, typename T> + inline + void + icp_(const p_array<P>& c, + const unsigned c_length, + const M& map, + T& qk, + const float epsilon = 1e-3) + { + trace::entering("registration::impl::icp_"); + + buffer<4,T> buf_qk; + buffer<3,float> buf_dk; + + float d_k = 10000; + p_array<P> Ck(c); + + algebra::vec<P::dim,float> mu_C = center(c, c_length), mu_Xk; + + buf_qk.store(qk); + + qk.apply_on(c, Ck, c_length); + + unsigned int k = 0; + do { + buf_dk.store(d_k); + + //compute qk + qk = match(c, mu_C, Ck, map, c_length); + buf_qk.store(qk); + + //Ck+1 = qk(c) + qk.apply_on(c, Ck, c_length); + + // d_k = d(Yk, Pk+1) + // = d(closest(qk-1(P)), qk(P)) + d_k = rms(c, map, c_length, buf_qk[1], qk); + + k++; + } while ((qk - buf_qk[1]).sqr_norm() / qk.sqr_norm() > epsilon); + + trace::exiting("registration::impl::icp_"); + } + + } // end of namespace mln::registration::impl + + + template <typename P, typename M, typename T> + inline + T + icp(const p_array<P>& c, + const M& map, + const float epsilon = 1e-3) + { + T qk; + impl::icp_(c, map, qk, c.length(), epsilon); + return qk; + } + + template <typename P, typename M, typename T> + inline + void + icp_subset(const p_array<P>& c, + const unsigned c_length, + const M& map, + T& qk, + const float epsilon = 1e-3) + { + impl::icp_(c, map, qk, c_length, epsilon); + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::registration + +} // end of namespace mln + + + +#endif // ! MLN_REGISTRATION_ICP_HH