[Olena-patches] 3225: Fix stuff for registration.

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet &lt;ugo.jardonnet(a)lrde.epita.fr&gt; Fix stuff for registration. * mln/core/image/lazy_image.hh: Fix box type. * mln/fun/x2x/rotation.hh: Fix matrix creation. * mln/registration/get_rot.hh: Assert 3d. * mln/registration/registration.hh: p_array convertion. * mln/registration/icp.hh: Remove debug print. * sandbox/lazzara/igr.cc: Add useless things. mln/core/image/lazy_image.hh | 4 ++-- mln/fun/x2x/rotation.hh | 2 +- mln/registration/get_rot.hh | 28 ++++++++++++++++++++-------- mln/registration/icp.hh | 14 ++------------ mln/registration/registration.hh | 6 +++--- sandbox/lazzara/igr.cc | 34 ++++++++++++++++++++++++++++++++-- 6 files changed, 60 insertions(+), 28 deletions(-) Index: mln/core/image/lazy_image.hh --- mln/core/image/lazy_image.hh (revision 3224) +++ mln/core/image/lazy_image.hh (working copy) @@ -120,7 +120,7 @@ void init_(const F& fun, const B& box); /// Return domain of lazyd_image. - const box2d& domain() const; + const box<mln_psite(I)>& domain() const; /// Test if a pixel value is accessible at \p p. bool has(const mln_psite(I)&) const; @@ -223,7 +223,7 @@ template <typename I, typename F, typename B> inline - const box2d& + const box<mln_psite(I)>& lazy_image<I,F,B>::domain() const { return this->data_->bb_; Index: mln/fun/x2x/rotation.hh --- mln/fun/x2x/rotation.hh (revision 3224) +++ mln/fun/x2x/rotation.hh (working copy) @@ -192,7 +192,7 @@ xy + zw, 1.f - x2 - z2, yz - xw, 0, xz - yw, yz + xw, 1.f - x2 - y2, 0, 0, 0, 0, 1}; - this->m_(make::mat<4,4>(t)); + this->m_ = make::mat<4,4>(t); } template <unsigned n, typename C> Index: mln/registration/get_rot.hh --- mln/registration/get_rot.hh (revision 3224) +++ mln/registration/get_rot.hh (working copy) @@ -61,12 +61,22 @@ const M& map, const algebra::vec<2u,float>& mu_xk) { - (void) c; - (void) mu_c; - (void) ck; - (void) map; - (void) mu_xk; - //FIXME: Write 2d version of rotation computation betwenn two p_array + assert(0 && "TODO"); + + //////SLOW VERSION + + //// V1 + // M1 := c covariance + // V1 := greatest eigen vector of M1 + + //// V2 + // M2 := c covariance + // V2 := greatest eigen vector of M2 + + //// alpha + // cos(alpha) = (V1.V2) / (|V1|.|V2|) + + //FIXME: Write 2d version of rotation computation between two p_arrays return fun::x2x::rotation<2u, float>(); } @@ -85,8 +95,8 @@ algebra::mat<3u,3u,float> Mk(literal::zero); for (unsigned i = 0; i < c.nsites(); ++i) { - algebra::vec<3u,float> ci = c[i]; - algebra::vec<3u,float> xki = map(ck[i]); + algebra::vec<3u,float> ci = convert::to< algebra::vec<3u,float> >(c[i]); + algebra::vec<3u,float> xki = convert::to< algebra::vec<3u,float> >(map(ck[i])); Mk += make::mat(ci - mu_c) * trans(make::mat(xki - mu_xk)); } Mk /= c.nsites(); @@ -121,6 +131,8 @@ algebra::quat qR(literal::zero); qR = math::jacobi(Qk); + std::cout << qR << std::endl; + return fun::x2x::rotation<3u, float>(qR); } Index: mln/registration/registration.hh --- mln/registration/registration.hh (revision 3224) +++ mln/registration/registration.hh (working copy) @@ -66,8 +66,8 @@ const J& surface) { //FIXME: Use convert::to< p_array<mln_psite(I)> >() - p_array<mln_psite(I)> c = convert::to_p_array(cloud); - p_array<mln_psite(J)> x = convert::to_p_array(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(J)> >(surface); //init rigid transform qk composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > qk; @@ -78,7 +78,7 @@ //make a lazy_image map via function closest_point fun::x2p::closest_point<mln_psite(I)> fun(x, working_box); - lazy_image<I, fun::x2p::closest_point<mln_psite(I)>, box2d> + lazy_image< I, fun::x2p::closest_point<mln_psite(I)>, box<mln_psite(I)> > map(fun, fun.domain()); //run registration Index: mln/registration/icp.hh --- mln/registration/icp.hh (revision 3224) +++ mln/registration/icp.hh (working copy) @@ -169,18 +169,10 @@ //work on a reduced version of c_ p_array<P> c = c_; - - std::cout << "c_length" << c_length << std::endl; - std::cout << "-nsites" << c_.nsites() << std::endl; - - std::cout << "bad point" << c[9100] << std::endl; - //FIXME: Use c_ and specify c_length to function using c_ c.hook_std_vector_().resize(c_length); - std::cout << "+nsites" << c.nsites() << std::endl; - std::cout << "bad point" << c[9100] << std::endl; - + //FIXME: avoid use of ck p_array<P> ck(c); algebra::vec<P::dim,float> mu_c = center(c); @@ -188,7 +180,6 @@ apply(qk, c, ck); - unsigned int k = 0; do { std::cout << "*------" << std::endl; buf_dk.store(d_k); @@ -204,13 +195,12 @@ // = d(closest(qk-1(P)), qk(P)) d_k = internal::rms(c, map, buf_qk[1], qk); - k++; //FIXME: Add matrix norm //} while ((qk - buf_qk[1]).sqr_norm() / qk.sqr_norm() > epsilon); std::cout << d_k << std::endl; - } while (d_k > 9.2); + } while (d_k > 9.2 + epsilon); trace::exiting("registration::impl::icp_"); } Index: sandbox/lazzara/igr.cc --- sandbox/lazzara/igr.cc (revision 3224) +++ sandbox/lazzara/igr.cc (working copy) @@ -1,4 +1,5 @@ #include <mln/essential/2d.hh> +#include <mln/core/image/image3d.hh> #include <mln/binarization/binarization.hh> #include <mln/labeling/background.hh> @@ -21,7 +22,38 @@ }; +template <typename T> +inline +mln::image3d<T> +to_image3d(const mln::image2d<T>& img) +{ + mln::point3d pmin(img.domain().pmin()[0], img.domain().pmin()[1], 0); + mln::point3d pmax(img.domain().pmax()[0], img.domain().pmax()[1], 0); + mln::image3d<T> img3d(mln::box3d(pmin,pmax)); + mln_piter(mln::image3d<T>) p(img3d.domain()); + for_all(p) + img3d(p) = exact(img)(mln::point2d(p[0],p[1])); + + return img3d; +} + +template <typename T> +inline +mln::point<mln::grid::square, T> +to_point2d(const mln::point<mln::grid::cube, T>& p) +{ + return mln::point<mln::grid::square, T>(p[0], p[1]); +} + +//point2d -> point3d +template <typename T> +inline +mln::point<mln::grid::cube, T> +to_point3d(const mln::point<mln::grid::square, T>& p) +{ + return mln::point<mln::grid::cube, T>(p[0], p[1], 0); +} template <typename I> mln_concrete(I) @@ -132,7 +164,5 @@ for_all(p) res(qk(p.to_site().to_vec())) = true; - io::pbm::save(slice(res,0), "04_registered.pbm"); - }