https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Ugo Jardonnet <jardonnet@lrde.epita.fr> Sandbox: ICP: clean-up before milena update. Prepare patch for math and geom modules. * jardonnet/registration/power_it.hh: Add comment . * jardonnet/registration/cov.hh: New: clean math::covariance. * jardonnet/registration/variance.hh: New: clean math::variance. * jardonnet/registration/exp_val.hh: New: clean math::expected value. * jardonnet/registration/cross_cov.hh: Update: clean cross covariance. * jardonnet/registration/center.hh: New: clean geom::center. Update rendering for seminar. * jardonnet/registration/save.hh: Update rendering. * jardonnet/test/img/c5.pbm, * jardonnet/test/img/x5.pbm: Update images . registration/center.hh | 73 ++++++++++++++++++++++++++++++++++++++++ registration/cov.hh | 82 ++++++++++++++++++++++++++++++++++++++++++++++ registration/cross_cov.hh | 6 +-- registration/exp_val.hh | 77 +++++++++++++++++++++++++++++++++++++++++++ registration/power_it.hh | 3 + registration/save.hh | 52 ++++++++++++++--------------- registration/variance.hh | 27 +++++++++++++++ test/icp_ref.cc | 2 - 8 files changed, 292 insertions(+), 30 deletions(-) Index: jardonnet/test/icp_ref.cc --- jardonnet/test/icp_ref.cc (revision 2044) +++ jardonnet/test/icp_ref.cc (working copy) @@ -60,7 +60,7 @@ } //working box - const box_< point_<grid::cube, float> > working_box = enlarge(bigger(c.bbox(),x.bbox()),100); + const box_< point_<grid::cube, float> > working_box = enlarge(bigger(c.bbox(),x.bbox()),50); // FIXME : TODO : map : vec<3,float> -> point closest_point< point_<grid::cube, float> > map(x, working_box); Index: jardonnet/test/img/c5.pbm Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: jardonnet/test/img/x5.pbm Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: jardonnet/registration/cov.hh --- jardonnet/registration/cov.hh (revision 0) +++ jardonnet/registration/cov.hh (revision 0) @@ -0,0 +1,82 @@ +// 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_MATH_COV_HH +# define MLN_MATH_COV_HH + +/*! \file mln/math/cov.hh + * + * \brief Define the covariance (cov) routine. + */ + +# include "exp_val.hh" + +namespace mln +{ + + namespace math + { + + template <typename P> + algebra::mat<P::dim,P::dim,float> + cov(const p_array<P>& a1, + const p_array<P>& a2); + +# ifndef MLN_INCLUDE_ONLY + + template <typename P> + algebra::mat<P::dim,P::dim,float> + cov(const p_array<P>& a1, + const p_array<P>& a2) + { + mln_precondition(a1.npoints() == a2.npoints()); + + //centers of mass + algebra::vec<P::dim,float> mu_a1 = math::exp_value(a1); + algebra::vec<P::dim,float> mu_a2 = math::exp_value(a2); + + //covariance matrix + algebra::mat<P::dim,P::dim,float> Mk(literal::zero); + for (unsigned i = 0; i < a1.npoints(); ++i) + { + // FIXME: ugly cast. + algebra::vec<P::dim,float> a1i = a1[i]; + algebra::vec<P::dim,float> a2i = a2[i]; + Mk += make::mat(a1i - mu_a1) * trans(make::mat(a2i - mu_a2)); + } + + return Mk / a1.npoints(); + } + + # endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::math + +} // end of namespace mln + + +#endif // ! MLN_MATH_COV_HH Index: jardonnet/registration/power_it.hh --- jardonnet/registration/power_it.hh (revision 2044) +++ jardonnet/registration/power_it.hh (working copy) @@ -12,6 +12,9 @@ namespace mln { + /** + * Return the biggest eigen vector. + */ template <uint n> algebra::vec<n,float> power_it(algebra::mat<n,n,float>& A) { Index: jardonnet/registration/variance.hh --- jardonnet/registration/variance.hh (revision 0) +++ jardonnet/registration/variance.hh (revision 0) @@ -0,0 +1,27 @@ +// 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. + Index: jardonnet/registration/save.hh --- jardonnet/registration/save.hh (revision 2044) +++ jardonnet/registration/save.hh (working copy) @@ -5,11 +5,16 @@ # include <mln/io/ppm/save.hh> # include <mln/io/pbm/save.hh> # include <mln/draw/all.hh> +# include <mln/morpho/erosion.hh> +# include <mln/make/window2d.hh> # include <string> # include "quat7.hh" # include "tools.hh" # include "power_it.hh" +# include "center.hh" +# include "cov.hh" + namespace mln { @@ -55,17 +60,11 @@ level::fill(out, literal::white); //plot mu_Ck - algebra::vec<P::dim,float> mu_Ck = center(ck, ck.npoints()); + algebra::vec<P::dim,float> mu_Ck = geom::center(ck); draw::plot(out, point2d(mu_Ck[0], mu_Ck[1]), literal::green); //Ck orientation - algebra::mat<P::dim,P::dim,float> Mk(literal::zero); - for (unsigned i = 0; i < ck.npoints(); ++i) - { - algebra::vec<P::dim,float> Cki = ck[i]; - Mk += make::mat(Cki - mu_Ck) * trans(make::mat(Cki - mu_Ck)); - } - Mk /= c.npoints(); + algebra::mat<P::dim,P::dim,float> Mk = math::cov(ck,ck); algebra::vec<3,float> vck = power_it(Mk); draw::line(out, point2d(mu_Ck[0], mu_Ck[1]), point2d(mu_Ck[0]+vck[0]*10, mu_Ck[1]+vck[1]*10), @@ -77,31 +76,16 @@ xk.append(map(ck[i])); //plot mu_Xk - algebra::vec<P::dim,float> mu_Xk = center(xk, xk.npoints()); + algebra::vec<P::dim,float> mu_Xk = geom::center(xk); draw::plot(out, point2d(mu_Xk[0], mu_Xk[1]), literal::blue); //Xk orientation - algebra::mat<P::dim,P::dim,float> Mxk(literal::zero); - for (unsigned i = 0; i < xk.npoints(); ++i) - { - algebra::vec<P::dim,float> Xki = xk[i]; - Mxk += make::mat(Xki - mu_Xk) * trans(make::mat(Xki - mu_Xk)); - } - Mxk /= c.npoints(); + algebra::mat<P::dim,P::dim,float> Mxk = math::cov(xk,xk); algebra::vec<3,float> vxk = power_it(Mxk); draw::line(out, point2d(mu_Xk[0], mu_Xk[1]), point2d(mu_Xk[0]+vxk[0]*10, mu_Xk[1]+vxk[1]*10), literal::red); - - //ck in green - for (unsigned i = 0; i < ck.npoints(); i++) - { - point2d p(ck[i][0], ck[i][1]); - if (out.has(p)) - out(p) = literal::green; - } - //x in black for (unsigned i = 0; i < x.npoints(); i++) { @@ -110,7 +94,7 @@ out(p) = literal::black; } - //xk in blue + //xk in red for (unsigned i = 0; i < xk.npoints(); i++) { point2d p(xk[i][0], xk[i][1]); @@ -118,6 +102,22 @@ out(p) = literal::red; } + //ck in green + for (unsigned i = 0; i < ck.npoints(); i++) + { + point2d p(ck[i][0], ck[i][1]); + if (out.has(p)) + out(p) = literal::green; + } + + /* + FIXME: + bool vals[] = {1, 1, 1, + 1, 1, 1, + 1, 1, 1}; + morpho::erosion(out, make::window2d(vals)); + */ + //save std::stringstream oss; oss << "step_" << id++ << ".ppm"; Index: jardonnet/registration/exp_val.hh --- jardonnet/registration/exp_val.hh (revision 0) +++ jardonnet/registration/exp_val.hh (revision 0) @@ -0,0 +1,77 @@ +// 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_MATH_EXP_VAL_HH +# define MLN_MATH_EXP_VAL_HH + +/*! \file mln/math/exp_val.hh + * + * \brief Define the expected value (exp_val) routine. + */ + +# include <mln/algebra/vec.hh> + +namespace mln +{ + + namespace math + { + + template <typename P> + algebra::vec<P::dim,float> + exp_value(const p_array<P>& a); + +# ifndef MLN_INCLUDE_ONLY + + template <typename P> + inline + algebra::vec<P::dim,float> + exp_value(const p_array<P>& a) + { + if (a.npoints() == 0) + return P(); + + algebra::vec<P::dim,float> c(literal::zero); + for (unsigned i = 0; i < a.npoints(); ++i) + { + // FIXME : Ugly. + algebra::vec<P::dim,float> ai = a[i]; + c += ai; + } + + return c / a.npoints(); + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::math + +} // end of namespace mln + + +#endif // ! MLN_MATH_ABS_HH Index: jardonnet/registration/cross_cov.hh --- jardonnet/registration/cross_cov.hh (revision 2044) +++ jardonnet/registration/cross_cov.hh (working copy) @@ -1,5 +1,5 @@ -#ifndef MLN_ACCU_CROSS_COV_HH_HH -# define MLN_ACCU_CROSS_COV_HH_HH +#ifndef MLN_ACCU_CROSS_COV_HH +# define MLN_ACCU_CROSS_COV_HH namespace mln @@ -15,4 +15,4 @@ } -#endif // ! MLN_ACCU_CROSS_COV_HH_HH +#endif // ! MLN_ACCU_CROSS_COV_HH Index: jardonnet/registration/center.hh --- jardonnet/registration/center.hh (revision 0) +++ jardonnet/registration/center.hh (revision 0) @@ -0,0 +1,73 @@ +// 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_GEOM_CENTER_HH +# define MLN_GEOM_CENTER_HH + +/*! \file mln/math/exp_val.hh + * + * \brief Define the center routine. + */ + +namespace mln +{ + + namespace geom + { + + template <typename P> + P center(const p_array<P>& a); + +# ifndef MLN_INCLUDE_ONLY + + template <typename P> + inline + P exp_value(const p_array<P>& a) + { + if (a.npoints() == 0) + return P(); + + algebra::vec<P::dim,float> c(literal::zero); + for (unsigned i = 0; i < a.npoints(); ++i) + { + // FIXME : Ugly. + algebra::vec<P::dim,float> ai = a[i]; + c += ai; + } + + return algebra::to_point<P>(c / a.npoints()); + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::geom + +} // end of namespace mln + + +#endif // ! MLN_MATH_ABS_HH