https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet &lt;ugo.jardonnet(a)lrde.epita.fr&gt; Sandbox: ICP: Fix distance map version: 0.5s for 01.pbm over 02.pbm. * mln/algebra/vec.hh: Add check for dimension in to_point(). * sandbox/jardonnet/test/plotscript: Add gnuplot script. * sandbox/jardonnet/test/Makefile: Add g++ flags: Unexpected behavior disappear with -ffloat-store (float are not stored in register). * sandbox/jardonnet/test/check: Update test script. * sandbox/jardonnet/registration/quat7.hh: . * sandbox/jardonnet/registration/icp.hh: . * sandbox/jardonnet/registration/projection.hh: . * sandbox/jardonnet/registration/tools.hh: . mln/algebra/vec.hh | 4 +- sandbox/jardonnet/registration/icp.hh | 39 +++++++++++++-------------- sandbox/jardonnet/registration/projection.hh | 6 ++-- sandbox/jardonnet/registration/quat7.hh | 2 - sandbox/jardonnet/registration/tools.hh | 4 +- sandbox/jardonnet/test/Makefile | 9 +++++- sandbox/jardonnet/test/check | 6 ++++ sandbox/jardonnet/test/plotscript | 5 +++ 8 files changed, 47 insertions(+), 28 deletions(-) Index: mln/algebra/vec.hh --- mln/algebra/vec.hh (revision 1818) +++ mln/algebra/vec.hh (working copy) @@ -550,10 +550,10 @@ } - template <typename P, unsigned n> + template <typename P> inline P - to_point(const vec<n,float>& v) + to_point(const vec<P::dim,float>& v) { P tmp; for (unsigned i = 0; i < P::dim; ++i) Index: sandbox/jardonnet/test/plotscript --- sandbox/jardonnet/test/plotscript (revision 0) +++ sandbox/jardonnet/test/plotscript (revision 0) @@ -0,0 +1,5 @@ +set xlabel "k" +set ylabel "err" +set zlabel "dQk" +splot "log" +pause -1 \ No newline at end of file Index: sandbox/jardonnet/test/Makefile --- sandbox/jardonnet/test/Makefile (revision 1818) +++ sandbox/jardonnet/test/Makefile (working copy) @@ -17,8 +17,15 @@ g++ icp.cc $(FLAG) -o '+icp.exe' icp++: - g++ icp.cc -I../../.. -O3 -o '+icp.exe' + g++ icp.cc -ffloat-store -I../../.. -O3 -DNDEBUG -o '+icp.exe' + +icpsafe: + g++ icp.cc -fsignaling-nans -ffloat-store -I../../.. -O1 -o '+icp.exe' run: time ./+sub.exe . . ; time ./+gsub.exe . . +plot: icp + ./+icp.exe 01.pbm 02.pbm > log + gnuplot plotscript + Index: sandbox/jardonnet/test/check --- sandbox/jardonnet/test/check (revision 1818) +++ sandbox/jardonnet/test/check (working copy) @@ -3,11 +3,17 @@ touch log echo "### report : projection" >> log + echo "*** maps ***" >> log time ./+icp_maps.exe +01.pbm +02.pbm >> log + + echo "*** de_base ***" >> log time ./+icp_base.exe +01.pbm +02.pbm >> log + + echo "*** memo ***" >> log time ./+icp_memo.exe +01.pbm +02.pbm >> log + cat log Index: sandbox/jardonnet/registration/quat7.hh --- sandbox/jardonnet/registration/quat7.hh (revision 1818) +++ sandbox/jardonnet/registration/quat7.hh (working copy) @@ -99,7 +99,7 @@ // qR algebra::mat<P::dim,P::dim,float> Mk; - for (unsigned i = 0; i < C.npoints(); ++i) + for (size_t i = 0; i < C.npoints(); ++i) { algebra::vec<P::dim,float> Ci = C[i]; algebra::vec<P::dim,float> Xki = Xk[i]; Index: sandbox/jardonnet/registration/icp.hh --- sandbox/jardonnet/registration/icp.hh (revision 1818) +++ sandbox/jardonnet/registration/icp.hh (working copy) @@ -36,7 +36,7 @@ # include <mln/algebra/quat.hh> # include <mln/algebra/vec.hh> # include <mln/make/w_window.hh> -# include <mln/make/window3d.hh> +# include <mln/make/w_window3d.hh> # include "tools.hh" @@ -66,12 +66,12 @@ namespace impl { - template <typename P, typename T> + template <typename P, typename M> inline p_array<P> icp_(p_array<P>& C, - const p_array<P>& X, - lazy_map<T>& map) + const p_array<P>&, + M& map) { trace::entering("registration::impl::icp_"); @@ -88,23 +88,27 @@ //// step 1 k = 0; do { + std::cout << "step 2" << std::endl; //// step 2 - //err_bis = projection::fill_Xk(Ck, map, Xk); + projection::fill_Xk(Ck, map, Xk); //projection::de_base(Ck, X, Xk, err_bis); - projection::memo(Ck, X, Xk, map); + //projection::memo(Ck, X, Xk, map); mu_Xk = center(Xk); + //std::cout << "step 3" << std::endl; //// step 3 old_qk = qk; qk = match(C, mu_C, Xk, mu_Xk); + //std::cout << "step 4" << std::endl; //// step 4 qk.apply_on(C, Ck); // Ck+1 = qk(C) + //std::cout << "step err" << std::endl; //// err = d(Ck+1,Xk) err = rms(Ck, Xk); - std::cout << k << ' ' << err << std::endl; //plot file + std::cout << k << ' ' << err << ' ' << (qk - old_qk).sqr_norm() << std::endl; //plot file ++k; } while (k < 3 || (qk - old_qk).sqr_norm() > epsilon); @@ -132,17 +136,14 @@ I3d cloud = convert::to_image_3d(exact(cloud_)); const I3d surface = convert::to_image_3d(exact(surface_)); - /* + //create a pair (distance map, closest point) - bool w[27] = {true, true, true, true, true, true, true, true, true, - true, true, true, true, true, true, true, true, true, - true, true, true, true, true, true, true, true, true}; - window<mln_dpoint(I3d)> win3d = make::window3d(w); - fun::cham<point3d> fun; - w_window<mln_dpoint(I3d), float> chamfer = make::w_window(win3d, fun); - */ - //std::pair<mln_ch_value(I3d,float), mln_ch_value(I3d,mln_point(I3d)) > maps;// = - //dt::chamfer(surface, chamfer); + float w[27] = {1.4142, 1, 1.4142, 1.4142, 1, 1.4142, 0, 0, 0, + 1, 1, 1, 1, 1, 0, 0, 0, 0, + 1.4142, 1, 1.4142, 0, 0, 0, 0, 0, 0}; + w_window<mln_dpoint(I3d), float> chamfer = make::w_window3d(w); + std::pair<mln_ch_value(I3d,float), mln_ch_value(I3d,mln_point(I3d)) > map = + dt::chamfer(surface, chamfer); //build p_arrays. @@ -151,7 +152,7 @@ //build closest point map //lazy_map<I3d> map(enlarge(bigger(c.bbox(),x.bbox()),50)); - lazy_map<I3d> map(1000,1000,50); + //lazy_map<I3d> map(1000,1000,50); p_array<mln_point(I3d)> res = impl::icp_(c, x, map); @@ -162,7 +163,7 @@ { point2d p(res[i][0], res[i][1]); //FIXME: not necessary if output(res.bbox()) - //if (output.has(p)) + if (output.has(p)) output(p) = true; } Index: sandbox/jardonnet/registration/projection.hh --- sandbox/jardonnet/registration/projection.hh (revision 1818) +++ sandbox/jardonnet/registration/projection.hh (working copy) @@ -22,13 +22,13 @@ for (size_t i = 0; i < Ck.npoints(); ++i) { //FIXME: bof - //if (map.second.has(Ck[i])) - //{ + if (map.second.has(Ck[i])) + { //x[i] := Ck[i] closest point in X Xk.hook_()[i] = map.second(Ck[i]); //err := Distance between Ck[i] and its closest point err += map.first(Ck[i]); - //} + } } return err /= Ck.npoints(); } Index: sandbox/jardonnet/registration/tools.hh --- sandbox/jardonnet/registration/tools.hh (revision 1818) +++ sandbox/jardonnet/registration/tools.hh (working copy) @@ -133,8 +133,8 @@ image3d<T> to_image_3d(const image2d<T>& img) { - point3d pmin(img.domain().pmin()[0], img.domain().pmin()[1], -1); - point3d pmax(img.domain().pmax()[0], img.domain().pmax()[1], 1); + point3d pmin(img.domain().pmin()[0], img.domain().pmin()[1], 0); + point3d pmax(img.domain().pmax()[0], img.domain().pmax()[1], 0); image3d<T> img3d(box3d(pmin,pmax)); mln_piter(image3d<T>) p(img3d.domain());