- Olena-patches

milena r1844: PSN algorithm
by Etienne FOLIO 07 Apr '08

07 Apr '08

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2008-04-07 Etienne FOLIO <folio(a)lrde.epita.fr> PSN algorithm. * sandbox/folio/chamfer.cc: Debug using an image to compare with PSN. * sandbox/folio/psn.cc: Propagation using a single neighborhood. --- chamfer.cc | 47 +++++++++---- psn.cc | 211 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 243 insertions(+), 15 deletions(-) Index: trunk/milena/sandbox/folio/psn.cc =================================================================== --- trunk/milena/sandbox/folio/psn.cc (revision 0) +++ trunk/milena/sandbox/folio/psn.cc (revision 1844) @@ -0,0 +1,211 @@ +// 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_DT_CHAMFER_HH +# define MLN_DT_CHAMFER_HH + +# include <queue> +# include <map> +# include <cmath> + +# include <mln/core/concept/image.hh> +# include <mln/core/concept/neighborhood.hh> +# include <mln/literal/zero.hh> + +namespace mln +{ + + namespace dt + { + + /*! Propagation using a single neighborhood (PSN). + * + * \param[in] input_ The input image. + * \param[in] chamfer The chamfer window to use for distance calcul. + * \return A pair <distance map, closest point map>. + * + * \pre \p img has to be initialized. + */ + template<typename I, typename N> + mln_ch_value(I, unsigned) + psn(const Image& input_, const N& nbh); + + +# ifndef MLN_INCLUDE_ONLY + + + namespace impl + { + + template <typename BP> + class compare + { + public: + bool + operator()(const BP& lhs, const BP& rhs) const + { + return lhs.second > rhs.second; + } + }; + + template <typename D> + unsigned + sq(const D& dp) + { + unsigned res = 0; + + for (unsigned i = 0; i < D::dim; ++i) + res += std::abs(dp[i]); // FIXME: dp[i] * dp[i]; + + return res; + } + + + } // end of namespace mln::dt::impl + + + + // Facade. + + template<typename I, typename N> + inline + mln_ch_value(I, unsigned) + psn(const Image& input_, const N& nbh) + { + const I& input = exact(input_); + mln_precondition(input.has_data()); + + mln_ch_value(I, unsigned) D; + initialize(D, input); + + static const unsigned M = 666; // FIXME + + // Initialization. + typedef mln_point(I) point; + typedef mln_dpoint(I) dpoint; + typedef std::pair<point, dpoint> BP; + + std::map<unsigned, std::queue<BP> > bucket; + unsigned bucket_size = 0; + + mln_fwd_piter(I) p(input.domain()); + for_all(p) + if (input(p)) + { + D(p) = literal::zero; + bucket[0].push(BP(p, literal::zero)); + ++bucket_size; + } + else + D(p) = M; + + unsigned d = 0; + + while (bucket_size != 0) + { + std::queue<BP>& bucket_d = bucket[d]; + + while (! bucket_d.empty()) + { + point p = bucket_d.front().first; + dpoint dp = bucket_d.front().second; + + bucket_d.pop(); + --bucket_size; + + if (D(p) == d) + { + mln_niter(N) n_(nbh, p); + + for_all(n_) + if (D.has(n_)) + { + dpoint n = n_ - p; + unsigned newD = impl::sq(dp + n); + + if (newD < D(p + n)) // p + n = p + n_ - p = n_ + { + D(n_) = newD; + bucket[newD].push(BP(p + n, dp + n)); + ++bucket_size; + } + } + } + } + ++d; + } + + return D; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::dt + +} // end of namespace mln + +#endif // ! MLN_DT_CHAMFER_HH + +#include <iostream> +#include <mln/core/image2d.hh> +#include <mln/debug/println.hh> +#include <mln/make/win_chamfer.hh> +#include <mln/level/fill.hh> +#include <mln/core/neighb2d.hh> + +#include <mln/io/pbm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/level/stretch.hh> +#include <mln/value/int_u8.hh> + +int main() +{ + using namespace mln; + +// image2d<bool> ima(5,5); +// bool vals[] = { 1, 1, 1, 0, 0, +// 0, 0, 1, 0, 0, +// 0, 0, 1, 0, 0, +// 0, 0, 0, 0, 0, +// 0, 0, 0, 0, 0 }; + +// level::fill(ima, vals); +// debug::println(ima); + + image2d<bool> ima = io::pbm::load("../../img/c01.pbm"); + + image2d<unsigned> out; + out = dt::psn(ima, c4()); + + image2d<value::int_u8> out2(out.domain()); + level::stretch(out, out2); + + io::pgm::save(out2, "out.pgm"); + +// std::cerr << "Distance:" << std::endl; +// debug::println(out); +} Index: trunk/milena/sandbox/folio/chamfer.cc =================================================================== --- trunk/milena/sandbox/folio/chamfer.cc (revision 1843) +++ trunk/milena/sandbox/folio/chamfer.cc (revision 1844) @@ -156,34 +156,51 @@ #endif // ! MLN_DT_CHAMFER_HH #include <iostream> -#include <mln/debug/println.hh> #include <mln/core/image2d.hh> +#include <mln/debug/println.hh> #include <mln/make/win_chamfer.hh> #include <mln/level/fill.hh> +#include <mln/io/pbm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/level/stretch.hh> +#include <mln/value/int_u8.hh> + int main() { using namespace mln; - w_window2d_int chamfer = make::mk_chamfer_3x3_int<3, 4>(); + w_window2d_int chamfer = make::mk_chamfer_3x3_int<1, 2>(); { - image2d<bool> ima(5,5); - bool vals[] = { 1, 1, 1, 0, 0, - 0, 0, 1, 0, 0, - 0, 0, 1, 0, 0, - 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0 }; +// image2d<bool> ima(5,5); +// bool vals[] = { 1, 1, 1, 0, 0, +// 0, 0, 1, 0, 0, +// 0, 0, 1, 0, 0, +// 0, 0, 0, 0, 0, +// 0, 0, 0, 0, 0 }; + +// level::fill(ima, vals); +// debug::println(ima); + +// std::pair<image2d<int>, image2d<mln_point_(image2d<bool>)> > out; +// for (int i = 0; i < 999; ++i) +// out = dt::chamfer(ima, chamfer); + +// std::cerr << "Distance:" << std::endl; +// debug::println(out.first); +// std::cerr << "PPP:" << std::endl; +// debug::println(out.second); - level::fill(ima, vals); - debug::println(ima); + image2d<bool> ima = io::pbm::load("../../img/c01.pbm"); - std::pair<image2d<int>, image2d<mln_point_(image2d<bool>)> > + std::pair<image2d<int>, image2d<mln_point_(image2d<bool>)> > out; out = dt::chamfer(ima, chamfer); - std::cerr << "Distance:" << std::endl; - debug::println(out.first); - std::cerr << "PPP:" << std::endl; - debug::println(out.second); + image2d<value::int_u8> out2(out.first.domain()); + level::stretch(out.first, out2); + + io::pgm::save(out2, "out.pgm"); + } }

1 0

1843: Sandbox: ICP: Remove intermediate p_array Xk.
by Ugo Jardonnet 05 Apr '08

05 Apr '08

Index: ChangeLog =================================================================== --- ChangeLog (revision 1842) +++ ChangeLog (working copy) @@ -1,3 +1,14 @@ +2008-04-05 Ugo Jardonnet <jardonnet(a)lrde.epita.fr> + + Sandbox: ICP: Remove intermediate p_array Xk. + + * sandbox/jardonnet/test/icp.cc: Move a big part of icp.hh here. + * sandbox/jardonnet/test/Makefile (icp++): Add rule. + * sandbox/jardonnet/registration/quat7.hh: Add direct access to the map. + * sandbox/jardonnet/registration/icp.hh: Remove Xk. + * sandbox/jardonnet/registration/tools.hh: Improve lazy_image interface + and checks. + 2008-04-05 Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: Add file for refactoring. Index: sandbox/jardonnet/test/icp.cc =================================================================== --- sandbox/jardonnet/test/icp.cc (revision 1842) +++ sandbox/jardonnet/test/icp.cc (working copy) @@ -4,17 +4,64 @@ #include <mln/io/pbm/save.hh> #include <sandbox/jardonnet/registration/icp.hh> +#include <sandbox/jardonnet/registration/tools.hh> -int main(int, char* argv[]) +void usage(char *argv[]) { + std::cout << "usage : " << argv[0] + << " cloud.pbm surface.pbm" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 3) + usage(argv); + using namespace mln; + typedef image3d< bool > I3d; image2d< bool > img1; image2d< bool > img2; - + + //load images io::pbm::load(img1, argv[1]); io::pbm::load(img2, argv[2]); + + //convert to image3d + I3d cloud = convert::to_image3d(img1); + I3d surface = convert::to_image3d(img2); - io::pbm::save(registration::icp(img1,img2), "./+registred.pbm"); + //build p_arrays. + p_array<mln_point_(I3d)> c = convert::to_p_array(cloud); + p_array<mln_point_(I3d)> x = convert::to_p_array(surface); + + //Make a lazy_image map via function closest_point + closest_point<mln_point_(I3d)> fun(x, box_<point3d>(1000,1000,1)); + lazy_image< closest_point<mln_point_(I3d)> > map(fun); + + registration::icp(c, map); + +#ifndef NDEBUG + std::cout << "closest_point(Ck[i]) = " << fun.i << std::endl; + std::cout << "Pts processed = " << registration::pts << std::endl; +#endif + + // //init output image + // //FIXME: Should be + // //mln_concrete(I) output(res.bbox()) = convert::to_image(res)? + // const box_<point2d> box2d(1000,1000,0); + // image2d<bool> output(box2d, 1); + + // //to 2d : projection (FIXME:if 3d) + // for (size_t i = 0; i < res.npoints(); i++) + // { + // point2d p(res[i][0], res[i][1]); + // if (output.has(p)) + // output(p) = true; + // } + + //io::pbm::save("./+registred.pbm"); } Index: sandbox/jardonnet/test/Makefile =================================================================== --- sandbox/jardonnet/test/Makefile (revision 1842) +++ sandbox/jardonnet/test/Makefile (working copy) @@ -16,8 +16,11 @@ run: time ./+sub.exe . . ; time ./+gsub.exe . . + icp: g++ icp.cc -I../../.. -g -o '+icp.exe' +icp++: + g++ icp.cc -I../../.. -O3 -DNDEBUG -o '+icp.exe' icp_max++: g++ icp_max.cc -I../../.. -O3 -DNDEBUG -pg -o '+icp_map.exe' Index: sandbox/jardonnet/registration/quat7.hh =================================================================== --- sandbox/jardonnet/registration/quat7.hh (revision 1842) +++ sandbox/jardonnet/registration/quat7.hh (working copy) @@ -59,10 +59,7 @@ template <typename P> void apply_on(const p_array& input, p_array& output) const { - assert(input.npoints() == output.npoints()); - - std::cout << _qR << std::endl; - + assert(input.npoints() == output.npoints()); assert(_qR.is_unit()); for (size_t i = 0; i < input.npoints(); i++) @@ -92,21 +89,28 @@ } - template <typename P> + template <typename P, typename M> quat7<P::dim> match(const p_array& C, const algebra::vec<P::dim,float>& mu_C, - const p_array& Xk, - const algebra::vec<P::dim,float>& mu_Xk) + const p_array& Ck, + M& map) { - assert(C.npoints() == Xk.npoints()); + //mu_Xk = center map(Ck) + algebra::vec<P::dim,float> mu_Xk(literal::zero); + for (size_t i = 0; i < Ck.npoints(); ++i) + { + algebra::vec<P::dim,float> xki = map(Ck[i]); + mu_Xk += xki; + } + mu_Xk /= Ck.npoints(); + // qR - algebra::mat<P::dim,P::dim,float> Mk(literal::zero); 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]; + algebra::vec<P::dim,float> Ci = C[i]; + algebra::vec<P::dim,float> Xki = map(Ck[i]); Mk += make::mat(Ci - mu_C) * trans(make::mat(Xki - mu_Xk)); } Mk /= C.npoints(); Index: sandbox/jardonnet/registration/icp_subsampled.hh =================================================================== --- sandbox/jardonnet/registration/icp_subsampled.hh (revision 1842) +++ sandbox/jardonnet/registration/icp_subsampled.hh (working copy) @@ -149,7 +149,7 @@ p_array<mln_point(I3d)> c = convert::to_p_array(cloud); p_array<mln_point(I3d)> x = convert::to_p_array(surface); - //init qk + //init rigid transform qk quat7<3> qk; //create subsampled p_array Index: sandbox/jardonnet/registration/icp.hh =================================================================== --- sandbox/jardonnet/registration/icp.hh (revision 1842) +++ sandbox/jardonnet/registration/icp.hh (working copy) @@ -51,6 +51,10 @@ namespace registration { +#ifndef NDEBUG + static unsigned pts = 0; +#endif + /*! Registration FIXME : doxy * * @@ -69,45 +73,52 @@ template <typename P, typename M> inline p_array - icp_(p_array& C, - const p_array&, - M& map) + icp_(const p_array& C, + M& map, + quat7<P::dim>& qk) { trace::entering("registration::impl::icp_"); +#ifndef NDEBUG + //display registred points + std::cout << "Register " + << C.npoints() << " points" << std::endl; + std::cout << "k\terror\tdqk" << std::endl; +#endif + unsigned int k; - quat7<P::dim> old_qk, qk; + quat7<P::dim> old_qk; float err; //float err_bis; p_array Ck(C), Xk(C); //FIXME: Xk copy C algebra::vec<P::dim,float> mu_C = center(C), mu_Xk; - + const float epsilon = 1;//1e-3; //// step 1 k = 0; do { - - //// step 2 - //projection::fill_Xk(Ck, map, Xk); - //projection::de_base(Ck, X, Xk, err_bis); - projection::memo(Ck, Xk, map); - - mu_Xk = center(Xk); - - //// step 3 + //compute qk old_qk = qk; - qk = match(C, mu_C, Xk, mu_Xk); + qk = match(C, mu_C, Ck, map); - //// step 4 - qk.apply_on(C, Ck); // Ck+1 = qk(C) + //Ck+1 = qk(C) + qk.apply_on(C, Ck); - //// err = d(Ck+1,Xk) + //err = d(Ck+1,Xk) err = rms(Ck, Xk); - std::cout << k << ' ' << err << ' ' << (qk - old_qk).sqr_norm() << std::endl; //plot file +#ifndef NDEBUG + //plot file + std::cout << k << '\t' << err << '\t' + << (qk - old_qk).sqr_norm() << '\t' + << std::endl; + //count the number of points processed + pts += Ck.npoints(); +#endif + ++k; } while (k < 3 || (qk - old_qk).sqr_norm() > epsilon); @@ -117,50 +128,28 @@ } // end of namespace mln::registration::impl - - //Only for 2d and 3d image - template <typename I, typename J> + + // Only for 3d images + template <typename P, typename M> inline - mln_concrete(I) //FIXME: should return something else ? qk ? - icp(const Image& cloud_, - const Image<J>& surface_) + quat7<P::dim> + icp(const p_array& cloud, + M& map) { trace::entering("registration::icp"); - mln_precondition(exact(cloud_).has_data()); - mln_precondition(exact(surface_).has_data()); + + mln_precondition(cloud.npoints() != 0); + mln_precondition(P::dim == 3); - //convert to image: time consuming - typedef image3d<mln_value(I)> I3d; - I3d cloud = convert::to_image_3d(exact(cloud_)); - const I3d surface = convert::to_image_3d(exact(surface_)); - - //build p_arrays. - p_array<mln_point(I3d)> c = convert::to_p_array(cloud); - p_array<mln_point(I3d)> x = convert::to_p_array(surface); - - //build closest point map - //lazy_map<I3d> map(enlarge(bigger(c.bbox(),x.bbox()),50)); - //lazy_map<I3d> map(1000,1000,50); - - c_point<mln_point(I3d)> fun(x); - //Make via function - lazy_image< c_point<mln_point(I3d)> > map(fun); + //init rigid transform qk + quat7<P::dim> qk; - p_array<mln_point(I3d)> res = impl::icp_(c, x, map); + //run icp + p_array res = impl::icp_(cloud, map, qk); - //to 2d : projection (FIXME:if 3d) - //mln_concrete(I) output = convert::to_image(res)? - mln_concrete(I) output(exact(cloud_).domain()); - for (size_t i = 0; i < res.npoints(); i++) - { - point2d p(res[i][0], res[i][1]); - //FIXME: not necessary if output(res.bbox()) - if (output.has(p)) - output(p) = true; - } - trace::exiting("registration::icp"); - return output; + + return qk; } # endif // ! MLN_INCLUDE_ONLY Index: sandbox/jardonnet/registration/tools.hh =================================================================== --- sandbox/jardonnet/registration/tools.hh (revision 1842) +++ sandbox/jardonnet/registration/tools.hh (working copy) @@ -11,19 +11,24 @@ { template <typename P> - struct c_point + struct closest_point { typedef P input; typedef P result; - c_point(const p_array& X) - : X(X) + closest_point(const p_array& X, const box_& box) + : X(X), box(box), i(0) { } result //inline operator () (const P& Ck) { + +#ifndef NDEBUG + ++i; +#endif + algebra::vec<P::dim,float> Cki = Ck; algebra::vec<P::dim,float> best_x = X[0]; float best_d = norm::l2(Cki - best_x); @@ -40,12 +45,14 @@ return algebra::to_point(best_x); } - const box_ domain() + const box_& domain() { - return X.bbox(); + return box; } const p_array& X; + const box_ box; + mutable unsigned i; }; @@ -56,29 +63,32 @@ { // Fun is potentially an image. lazy_image(F& fun) - : value(fun.domain()), is_known(fun.domain()), functor(fun) + : value(fun.domain()), is_known(fun.domain()), fun(fun) { } - // FIXME: remove this constructor - lazy_image(const F& fun, int nrows, int ncols) - : value(nrows, ncols), is_known(nrows,ncols), functor(fun) + // FIXME: hack, remove this constructor + lazy_image(F& fun, int nrows, int ncols, int nslis) + : value(nrows, ncols,1), is_known(nrows,ncols,1), fun(fun) { } const mln_result(F) //inline operator() (const typename F::input& p) { + mln_precondition(fun.domain().has(p)); + //FIXME: What about domain? if (is_known(p)) return value(p); - value(p) = functor(p); + value(p) = fun(p); + is_known(p) = true; return value(p); } //FIXME: 3d -> //mln_dim(ml_input(input)) image3d<mln_result(F)> value; image3d<bool> is_known; - F& functor; + F& fun; }; @@ -165,26 +175,27 @@ } } */ - + + //const return a const template <typename T> inline const image3d<T>& - to_image_3d(const image3d<T>& img) + to_image3d(const image3d<T>& img) { return img; } - + //non const return non const template <typename T> image3d<T>& - to_image_3d(image3d<T>& img) + to_image3d(image3d<T>& img) { return img; } - + template <typename T> inline image3d<T> - to_image_3d(const image2d<T>& img) + to_image3d(const image2d<T>& img) { point3d pmin(img.domain().pmin()[0], img.domain().pmin()[1], 0); point3d pmax(img.domain().pmax()[0], img.domain().pmax()[1], 0);

1 0

1842: Sandbox: ICP: Add file for refactoring.
by Ugo Jardonnet 05 Apr '08

05 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: Add file for refactoring. * sandbox/jardonnet/registration/icp.hh: New file for the refactored version. icp.hh | 173 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 173 insertions(+) Index: sandbox/jardonnet/registration/icp.hh --- sandbox/jardonnet/registration/icp.hh (revision 0) +++ sandbox/jardonnet/registration/icp.hh (revision 0) @@ -0,0 +1,173 @@ +// 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 image registration + */ + +# include <mln/algebra/quat.hh> +# include <mln/algebra/vec.hh> +# include <mln/make/w_window.hh> +# include <mln/make/w_window3d.hh> + +# include "tools.hh" + +# include "cloud.hh" +# include "quat7.hh" +# include "projection.hh" +# include "chamfer.hh" + +namespace mln +{ + + namespace registration + { + + /*! Registration FIXME : doxy + * + * + */ + template <typename I, typename J> + inline + mln_concrete(I) + icp(const Image& cloud, + const Image<J>& surface); + +# ifndef MLN_INCLUDE_ONLY + + namespace impl + { + + template <typename P, typename M> + inline + p_array + icp_(p_array& C, + const p_array&, + M& map) + { + trace::entering("registration::impl::icp_"); + + unsigned int k; + quat7<P::dim> old_qk, qk; + float err; + //float err_bis; + p_array Ck(C), Xk(C); //FIXME: Xk copy C + + algebra::vec<P::dim,float> mu_C = center(C), mu_Xk; + + const float epsilon = 1;//1e-3; + + //// step 1 + k = 0; + + do { + + //// step 2 + //projection::fill_Xk(Ck, map, Xk); + //projection::de_base(Ck, X, Xk, err_bis); + projection::memo(Ck, Xk, map); + + mu_Xk = center(Xk); + + //// step 3 + old_qk = qk; + qk = match(C, mu_C, Xk, mu_Xk); + + //// step 4 + qk.apply_on(C, Ck); // Ck+1 = qk(C) + + //// err = d(Ck+1,Xk) + err = rms(Ck, Xk); + std::cout << k << ' ' << err << ' ' << (qk - old_qk).sqr_norm() << std::endl; //plot file + + ++k; + } while (k < 3 || (qk - old_qk).sqr_norm() > epsilon); + + trace::exiting("registration::impl::icp_"); + return Ck; + } + + } // end of namespace mln::registration::impl + + + //Only for 2d and 3d image + template <typename I, typename J> + inline + mln_concrete(I) //FIXME: should return something else ? qk ? + icp(const Image& cloud_, + const Image<J>& surface_) + { + trace::entering("registration::icp"); + mln_precondition(exact(cloud_).has_data()); + mln_precondition(exact(surface_).has_data()); + + //convert to image: time consuming + typedef image3d<mln_value(I)> I3d; + I3d cloud = convert::to_image_3d(exact(cloud_)); + const I3d surface = convert::to_image_3d(exact(surface_)); + + //build p_arrays. + p_array<mln_point(I3d)> c = convert::to_p_array(cloud); + p_array<mln_point(I3d)> x = convert::to_p_array(surface); + + //build closest point map + //lazy_map<I3d> map(enlarge(bigger(c.bbox(),x.bbox()),50)); + //lazy_map<I3d> map(1000,1000,50); + + c_point<mln_point(I3d)> fun(x); + //Make via function + lazy_image< c_point<mln_point(I3d)> > map(fun); + + p_array<mln_point(I3d)> res = impl::icp_(c, x, map); + + //to 2d : projection (FIXME:if 3d) + //mln_concrete(I) output = convert::to_image(res)? + mln_concrete(I) output(exact(cloud_).domain()); + for (size_t i = 0; i < res.npoints(); i++) + { + point2d p(res[i][0], res[i][1]); + //FIXME: not necessary if output(res.bbox()) + if (output.has(p)) + output(p) = true; + } + + trace::exiting("registration::icp"); + return output; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::registration + +} // end of namespace mln + + +#endif // ! MLN_REGISTRATION_ICP_HH

1 0

1841: Sandbox: ICP: Improve lazy_image implementation.
by Ugo Jardonnet 04 Apr '08

04 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: Improve lazy_image implementation. * sandbox/jardonnet/test/Makefile: Add rules for the lazy version. * sandbox/jardonnet/registration/icp_map.hh: New dedicated version to map. * sandbox/jardonnet/registration/icp_lazy.hh: Rename lazy_map as lazy_image. Improve implementation. * sandbox/jardonnet/registration/icp.hh (icp_): Refactoring. * sandbox/jardonnet/registration/projection.hh (memo): Make a projection using the new lazy_image interface. * sandbox/jardonnet/registration/tools.hh (c_point): Add functor. registration/icp_lazy.hh | 2 registration/projection.hh | 28 ------------ registration/tools.hh | 103 ++++++++++++++++++++++++++++++++++----------- test/Makefile | 7 +-- test/icp_lazy.cc | 27 +++++++++++ 5 files changed, 110 insertions(+), 57 deletions(-) Index: sandbox/jardonnet/test/icp_lazy.cc --- sandbox/jardonnet/test/icp_lazy.cc (revision 1840) +++ sandbox/jardonnet/test/icp_lazy.cc (working copy) @@ -5,6 +5,33 @@ #include <sandbox/jardonnet/registration/icp_lazy.hh> +template <typename I> +class fun +{ + typedef mln_point(I) result; + + const operator () (mln_point(I)& p) + { + algebra::vec<P::dim,float> Cki = Ck[i]; + algebra::vec<P::dim,float> best_x = X[0]; + best_d = norm::l2(Cki - best_x); + for (size_t j = 1; j < X.npoints(); ++j) + { + algebra::vec<P::dim,float> Xj = X[j]; + float d = norm::l2(Cki - Xj); + if (d < best_d) + { + best_d = d; + best_x = Xj; + } + } + Xk.hook_()[i] = algebra::to_point(best_x); + map.map(Ck[i]) = Xk[i]; + map.known(Ck[i]) = true; + } +} + + int main(int, char* argv[]) { using namespace mln; Index: sandbox/jardonnet/test/Makefile --- sandbox/jardonnet/test/Makefile (revision 1840) +++ sandbox/jardonnet/test/Makefile (working copy) @@ -16,12 +16,11 @@ run: time ./+sub.exe . . ; time ./+gsub.exe . . - icp: - g++ icp.cc -I../../.. -g -o '+icp_map.exe' + g++ icp.cc -I../../.. -g -o '+icp.exe' -icp++: - g++ icp.cc -I../../.. -O3 -DNDEBUG -pg -o '+icp_map.exe' +icp_max++: + g++ icp_max.cc -I../../.. -O3 -DNDEBUG -pg -o '+icp_map.exe' icp_lazy++: g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -pg -o '+icp_lazy.exe' Index: sandbox/jardonnet/registration/icp_lazy.hh --- sandbox/jardonnet/registration/icp_lazy.hh (revision 1840) +++ sandbox/jardonnet/registration/icp_lazy.hh (working copy) @@ -89,8 +89,6 @@ k = 0; do { - //std::cout << "step 2" << std::endl; - //// step 2 projection::memo(Ck, X, Xk, map); Index: sandbox/jardonnet/registration/projection.hh --- sandbox/jardonnet/registration/projection.hh (revision 1840) +++ sandbox/jardonnet/registration/projection.hh (working copy) @@ -69,37 +69,13 @@ template <typename P, typename T> void memo(const p_array& Ck, - const p_array& X, p_array& Xk, - lazy_map<T>& map) // first: closest points, second: is_computed + lazy_image<T>& map) // first: closest points, second: is_computed { assert(Ck.npoints() == Xk.npoints()); for (size_t i = 0; i < Ck.npoints(); ++i) - { - float best_d; - if (map.known(Ck[i]) == false) - { - algebra::vec<P::dim,float> Cki = Ck[i]; - algebra::vec<P::dim,float> best_x = X[0]; - best_d = norm::l2(Cki - best_x); - for (size_t j = 1; j < X.npoints(); ++j) - { - algebra::vec<P::dim,float> Xj = X[j]; - float d = norm::l2(Cki - Xj); - if (d < best_d) - { - best_d = d; - best_x = Xj; - } - } - Xk.hook_()[i] = algebra::to_point(best_x); - map.map(Ck[i]) = Xk[i]; - map.known(Ck[i]) = true; - } - else - Xk.hook_()[i] = map.map(Ck[i]); - } + Xk.hook_()[i] = map(Ck[i]); } } // end of namespace projection Index: sandbox/jardonnet/registration/tools.hh --- sandbox/jardonnet/registration/tools.hh (revision 1840) +++ sandbox/jardonnet/registration/tools.hh (working copy) @@ -10,43 +10,96 @@ namespace mln { - - //FIXME: Shall we use something *really* lazy - template <typename I> - struct lazy_map - { template <typename P> - lazy_map(const box_& domain) - : map(domain), known(domain) - { } + struct c_point + { + typedef P input; + typedef P result; - lazy_map(int nrows, int ncols, int bdr = 3) - : map(nrows, ncols, bdr), known(nrows,ncols,bdr) + c_point(const p_array& X) + : X(X) { } - mln_ch_value(I, mln_point(I)) map; - mln_ch_value(I, bool) known; - }; + result + //inline + operator () (const P& Ck) + { + algebra::vec<P::dim,float> Cki = Ck; + algebra::vec<P::dim,float> best_x = X[0]; + float best_d = norm::l2(Cki - best_x); + for (size_t j = 1; j < X.npoints(); ++j) + { + algebra::vec<P::dim,float> Xj = X[j]; + float d = norm::l2(Cki - Xj); + if (d < best_d) + { + best_d = d; + best_x = Xj; + } + } + return algebra::to_point(best_x); + } + const box_ domain() + { + return X.bbox(); + } - // Point + const p_array& X; + }; - template <typename P> - P min(const P& a, const P& b) + + // FIXME: Should be a morpher ? + // we could acces domain of a lazy map, iterator etc... + template < typename F> + struct lazy_image + { + // Fun is potentially an image. + lazy_image(F& fun) + : value(fun.domain()), is_known(fun.domain()), functor(fun) { - if (a > b) - return b; - return a; } - template <typename P> - P max(const P& a, const P& b) - { - if (a < b) - return b; - return a; + // FIXME: remove this constructor + lazy_image(const F& fun, int nrows, int ncols) + : value(nrows, ncols), is_known(nrows,ncols), functor(fun) + { } + + const mln_result(F) + //inline + operator() (const typename F::input& p) + { + if (is_known(p)) + return value(p); + value(p) = functor(p); + return value(p); } + //FIXME: 3d -> //mln_dim(ml_input(input)) + image3d<mln_result(F)> value; + image3d<bool> is_known; + F& functor; + }; + + +// // Point + +// template <typename P> +// P min(const P& a, const P& b) +// { +// if (a > b) +// return b; +// return a; +// } + +// template <typename P> +// P max(const P& a, const P& b) +// { +// if (a < b) +// return b; +// return a; +// } + // Box

1 0

1840: sandbox: Continue documentation about image types and properties.
by Nicolas Ballas 04 Apr '08

04 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Nicolas Ballas <ballas(a)lrde.epita.fr> sandbox: Continue documentation about image types and properties. * sandbox/ballas/doc/image_tours.txt: update documentation. * sandbox/ballas/doc/draft.txt: New file about the images properties. draft.txt | 260 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ image_tours.txt | 160 ++++++++++++---------------------- 2 files changed, 319 insertions(+), 101 deletions(-) Index: sandbox/ballas/doc/image_tours.txt --- sandbox/ballas/doc/image_tours.txt (revision 1839) +++ sandbox/ballas/doc/image_tours.txt (working copy) @@ -7,36 +7,36 @@ ** Two types of properties - -*** Declaratif mode: (find an exemple). - +*** Declarative mode: (find an example). Here, the properties must be define by the image type programmer. They cannot be deduce from the image associated type. So we can't find the properties by -introspecting the image type (with introspection). +introspecting the image type. declaration definition ..... I check ------- ---------- | | / \ | Props | ----> | Image Type | -| | \ ___________/ +| | \ __________ / ------- -*** Automatic mode: (find an exemple). +*** Automatic mode: (find an example). Here, we get the image properties from the image associated type (site, pset...). - ----------- + ---------- / \ - | Image Type | automatic - and associated -----> Props - type - \ ___________/ + | | + | Image Type | + | and | automatic -----> Props + | associated| + | type | + \ __________ / @@ -53,36 +53,29 @@ (ex: dilatation/erosion). TODO: detail the example. -*** Implementation inherentence +*** Implementation inheritance Milena image types are property driven. It is possible to get a special behavior (not describe in the Image Concept) depends on the property define in the property associated for the image type. So, it is possible to recover some piece of interface depending on the image type properties. -For instance an Image2d has the bidimensional, and random value access. +For instance an Image2d has the bi-dimensional, and random value access. properties. Hence, Image2d automatically provides a row/column access: V at(Point2d<T>::coord x, Point2d<T>::coord y). - - *** Difference between a property and its values -*** Hierarchie of property (scalar/integer) +*** Hierarchies of property (scalar/integer) ** Implementation of properties in Milena -* Algorithm and Properties - -** Specialization - -** Static checking - -* Image Properties: +* List of the image properties: ** global Properties: + category: primary, { domain_morpher, value_morpher, identity_morpher } < morpher @@ -176,20 +169,17 @@ ** Image nD -Image on regular grid where grides nodes are points. -To gain efficienty, Images nD have a virtual borders. +Image on regular grid where grid nodes are points. +To gain efficiency, Images nD have a virtual borders. All these image are templated by T, the image value type. -*** Image1d<T> + +*** Common Associated type//properties between imageND **** Associated types: value = T -site = point1d -psite = point1d == point_<tick,int> -pset = box1d - **** Properties: --> global properties @@ -211,99 +201,56 @@ --> properties related to the pset access = random -space = one_d -size = regular +size = regular // this a global property support = aligned -**** Specific interface: - FIXME.... +*** Image1d<T> + +**** Associated types: +site = point1d +psite = point1d == point_<tick,int> +pset = box1d +**** Properties: +space = one_d +**** Specific interface: + FIXME.... *** Image2d<T> **** Associated types: -value = T - site = point2d psite = point2d == point_<square,int> pset = box2d **** Properties: ---> global properties - -category = primary -border = stored -neighb: none -data = raw -io = read_write -speed = fastest - ---> properties related to values - -// Depend on T FIXME, where the property are defined?? -kind = -quant = -value = - ---> properties related to the pset - -access = random space = two_d -size = regular -support = aligned + **** Specific interface: FIXME.... - - *** Image3d<T> **** Associated types: value = T -site = point3d -psite = point3d == point_<cube,int> -pset = box3d - **** Properties: ---> global properties - -category = primary -border = stored -neighb: none -data = raw -io = read_write -speed = fastest - ---> properties related to values - -// Depend on T FIXME, where the property are defined?? -kind = -quant = -value = - ---> properties related to the pset - -access = random space = three_d -size = regular -support = aligned **** Specific interface: FIXME.... - ** Image based on function The images based on function are templated by F and P, where @@ -355,32 +302,44 @@ **** Specific interface: FIXME.... - - - ** Run based Encoded image -All the run based encoded images have their definition domains encoded by ... -runs. These images types are not morpher types since these do not lie on -another image type. -A run is a "continuous" (note this is wrong, cf graph index) set of sites. -It is define by a site start, and a length. -All the run based encoded image are templated by P, T where P is a site type -and T a value type. - - +All the run based encoded images have their definition domains encoded by runs. +These images types are not morpher types since these do not lie on another +image type.A run is a "continuous" (note this is wrong, cf graph index) set of +sites.It is define by a site start, and a length. All the run based encoded +image are templated by P, T where P is a site type and T a value type. Note: type one shot (const)!!! Do we store the the null value in the rle encoding. (Do we compress all the images values, or only the image objects...) Is there any restriction on P and T? -Can we compress all the image types or only the point wise image type -(cf the on_the_same_line in the encode algorithm). How the properties are setted and to which values? -Depend on the PointWise property ?? + +*** Criteria needed by images to be encoded in runs + +Can we compress all the image types or only the piece wise image type +(cf the on_the_same_line in the encode algorithm). + +Currently, *encode function just work with image based on a **Point** set. +cf: +for (unsigned n = 0; same_line && n < dim - 1; ++n) + same_line = (p1[n] == p2[n]); + +Problem generalization: + A p_run take a Site as a parameter. + a p_run must be able to be casted into the Site (actually, const Site). + The only information that we have in the p_run is a len (integer). + So we must be able to refind all the site in the run from the site start, +and a addition with an integer. + + How do we name this property? + + + *** RLE Image @@ -470,7 +429,6 @@ *** Value morpher - *** Image_if ??? Index: sandbox/ballas/doc/draft.txt --- sandbox/ballas/doc/draft.txt (revision 0) +++ sandbox/ballas/doc/draft.txt (revision 0) @@ -0,0 +1,260 @@ +-*- outline -*- + +* Main image types to keep in mind + +- nD image + +- Runs image (RLE/ Sparse/ value encoded... + +- Image based on a Lut (to define) + +- graph image (see with Roland) + + +* New redefinition of the image property + + +** global Properties: + +*** Category (cf file milenadoc/tutorial/image_types.txt): + category: primary, + { domain_morpher, value_morpher, identity_morpher } < morpher + +*** Size: + + size: huge, + (large,) ? + regular + +**** Definition + The size property give us an indication about the size (in memory) taken by +the image. Indeed, an image can either huge or regular. +"huge" means that the image cannot be fully store in the ram. +"regular" means that the image can be fully store in the ram. + +**** Why does this property exist? +As we can imagine, access to a "huge" image is way slower than accessing +to "regular" image. Indeed, to access to "huge" image, we first need to load in +the ram the interesting part of the image from the hard drive. + +So, with this property, we can handle huge image in a different way in +the algorithms. (do inplace modification instead of create a new image for +instance, forbid the launch algorithm which create temporary image with huge +image...). + + + +FIXME define large? + + + +*** Border/Extended domain: + + border: none, + { stored, computed } <some + +FIXME: detail the problem, find example... + + The border property will probably change, because more information on +the border. +Actually a "border" is more extended domain of the image definition domain, +than a simple border. +The idea is an algorithm should have the same behavior for the same image with +or without an extended domain. + +problems: +If the image has stored data, but a computed border, we can do +(&ima(p)) if p is include in the image definition domain. +But we cannot do that if p is include int the image extended domain. + + + +*** Data + + data: stored, + linear < stored + raw < linear, + computed + +**** Definition + + The data property indicates how the image handle the value. +The value can be either computed (on the fly) or stored. +We say that a value is stored when we can takes a reference of it (&ima(p)) +Value stored can linear: + + | | | + v v v + -------- --------- ------------ +| | | | | | + -------- --------- ------------ + +Or they can be stored in a raw way (just in one block): + + | + v + ----------------------------------- + | | + ----------------------------------- + + +**** Note +computed data property implies that the image type has to be read_only. + + +*** io + io: read, + write, + read_only < read, + (write_only < write,) Is it useful? + read_write < both read'n write + +**** Definition + + Read/write Permissions for *ima(p)* + The name is probably bad, because image can also be accessed through +their value (ex: Run image encoded by value). + So, we should change the property name to access site_io/site_access or +something like that. + + + +*** Speed + + + speed: slow, + fast, + fastest + +**** Definition + + Give us information the time needed to access to a value from a site. + + slow: ima(p) is greater than O(1) + fast: ima(p) is in O(1) + fastest: ima(p) is in O(1) + Do we have an extended domain? + Can we access to the image value, directly with pointer in this +case? + + + + + + +** Properties related to I::value + +*** Kind + + kind: color, + gray, + label, + named < label, + logic, + mvlogic < logic && < label, + binary < mvlogic, + ternary < mvlogic + fuzzy < logic, + data, + map, + distance < map + data + +*** Quant + quant: low, + high + +*** Value + value: scalar, + vectorial, + structed, + pointer + + +** Properties related to I::pset + +*** access + + access: random, + browsing + +Does this property resume psite != site? + + + + +*** space + space: one_d, + two_d, + three_d + +Do we really need this property +If we look at the slice morpher, it must know that the input image is in 2d. +So, the answer tends to be yes. + +But in this case, how do we deal with no localized site? Specify the space +value in localization property. +Add a new value/property? + +*** support + + support: irregular, + regular + aligned < regular + +This property will certainly change. + +Support property represent the underlying support of an Image type. +It can be a grid, (regular, orthogonal, linear, isotrope, anysotrope). +It can be different from a grid (a mesh?). We also want to know the value of a +grid "stride (pas)", if this value is different from 1. Mask and convolution +directly depends on these properties. Most of the filter mask are developed to +work with grid with a stride of one. + +So, here we need to determine two things. +First, we need to know if the underling support of the image is a grid +or a mesh. Then, if the image is based on a grid we must provide an access +to the properties a grid. +We also have to choose which property of the grid are relevant, and which +property of the grid are useless. + + + + + +* Property that we need? +** has(p) + +Property that give the complexity of has(p). +-- O(1) when site = psite +-- depend on the implementation when site != psite + +** localization + + localization: none, + space, + grid < space, + regular_grid < grid + +This property should replace the support property. + +An image on a no-localized space has this property set to none. +An image on a localized space (each site can be convert into a site) has +this property set to space. +An image on a grid (which is not regular) has its property set to grid. +Does an image on a grid has always points as site? +An image based on the base grid (aligned, isotrope/anysotrop, othogonal) has +this property set to grid. + +When this property are set, you can access to new information from the image. +space -> box +regular_grid -> deltaX, deltaY of the grid + other grid property?? + + +Note: that a site can be localised without being point wise (site are localized +by vector). Do we need a property to say that the image is composed by point? +Furthermore, an image can be localized and have its type psite != site. +Ex: graph image<Point2d, T> We need a new property in order to define if the +image is localized. Is this property gives us the same information than the +support property? +

1 0

milena r1839: Add bounding box computation to p_image2d
by Matthieu Garrigues 04 Apr '08

04 Apr '08

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2008-04-04 Matthieu Garrigues <garrigues(a)lrde.epita.fr> Add bounding box computation to p_image2d. * mln/core/p_image2d.hh: (p_image2d::box_<mln_point(P)>& bbox()): New, give the exact boungind box. (p_image2d::insert(const p_image2d& set)): Rename as.. (p_image2d::insert(const Point_Set<S>& set)): ...this, and update it to be more generic. (p_image2d::remove(const p_image2d& set)): Rename as.. (p_image2d::remove(const Point_Set<S>& set)): ...this, likewise. --- p_image2d.hh | 59 ++++++++++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 46 insertions(+), 13 deletions(-) Index: trunk/milena/mln/core/p_image2d.hh =================================================================== --- trunk/milena/mln/core/p_image2d.hh (revision 1838) +++ trunk/milena/mln/core/p_image2d.hh (revision 1839) @@ -37,9 +37,7 @@ # include <mln/core/internal/point_set_base.hh> # include <mln/core/box2d.hh> # include <mln/core/image2d.hh> -# include <mln/accu/bbox.hh> # include <mln/level/fill.hh> -# include <mln/core/p_array.hh> namespace mln @@ -67,11 +65,13 @@ /// Insert a point \p p. p_image2d& insert(const P p); - p_image2d& insert(const p_image2d& set); + template <typename S> + p_image2d& insert(const Point_Set<S>& set); /// Remove a point \p p. p_image2d& remove(const P p); - p_image2d& remove(const p_image2d& set); + template <typename S> + p_image2d& remove(const Point_Set<S>& set); /// Give the number of points. unsigned npoints() const; @@ -82,6 +82,9 @@ /// Test if the set is empty. bool is_empty() const; + /// Give the exact bounding box. + const box_<mln_point(P)>& bbox();; + /// Clear this set. void clear(); @@ -90,6 +93,8 @@ private: image2d<bool> points_; unsigned npoints_; + accu::bbox bb_; + bool bb_need_update_; }; @@ -98,7 +103,8 @@ template <typename P> p_image2d::p_image2d(int nrows, int ncols) : points_(nrows, ncols), - npoints_(0) + npoints_(0), + bb_need_update_(false) { level::fill(points_, false); @@ -107,7 +113,8 @@ template <typename P> p_image2d::p_image2d(const box2d& b) : points_(b), - npoints_(0) + npoints_(0), + bb_need_update_(false) { level::fill(points_, false); } @@ -119,6 +126,7 @@ { if (points_(p) == false) { + bb_.take(p); points_(p) = true; npoints_++; } @@ -126,12 +134,12 @@ } template <typename P> + template <typename S> p_image2d& - p_image2d::insert(const p_image2d& set) + p_image2d::insert(const Point_Set<S>& set_) { - if (set->is_empty()) - return *this; - mln_fwd_piter(image2d<bool>) p(set.points_ | true); + const S& set = exact(set_); + mln_fwd_piter(S) p(set); for_all(p) if (this->points_.has(p)) this->insert(p); @@ -145,18 +153,21 @@ if (points_(p) == true) { points_(p) = false; + bb_need_update_ = true; npoints_--; } return *this; } template <typename P> + template <typename S> p_image2d& - p_image2d::remove(const p_image2d& set) + p_image2d::remove(const Point_Set<S>& set_) { - if (this->is_empty() || set->is_empty()) + const S& set = exact(set_); + if (this->is_empty()) return *this; - mln_fwd_piter(image2d<bool>) p(set.points_ | true); + mln_fwd_piter(S) p(set); for_all(p) if (this->points_.has(p)) this->remove(p); @@ -184,11 +195,33 @@ return npoints_ == 0; } + + template <typename P> + inline + const box_<mln_point(P)>& + p_image2d::bbox() + { + mln_precondition(npoints() != 0); + if (bb_need_update_) + { + bb_.init(); + + mln_fwd_piter(p_image2d) p(*this); + for_all(p) + bb_.take(p); + bb_need_update_ = false; + } + + return bb_.to_result(); + } + template <typename P> void p_image2d::clear() { level::fill(points_, false); + bb_.init(); + bb_need_update_ = false; } template <typename P>

1 0

1838: Sandboc : ICP: Update check script (gprof).
by Ugo Jardonnet 03 Apr '08

03 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandboc : ICP: Update check script (gprof). * sandbox/jardonnet/test/bench: Add info. * sandbox/jardonnet/test/icp_check.sh: Update testes. * sandbox/jardonnet/test/Makefile: Update rule for gprof. Makefile | 11 +++++++---- bench | 2 ++ icp_check.sh | 8 ++++---- 3 files changed, 13 insertions(+), 8 deletions(-) Index: sandbox/jardonnet/test/bench --- sandbox/jardonnet/test/bench (revision 1837) +++ sandbox/jardonnet/test/bench (working copy) @@ -30,3 +30,5 @@ 0m6.468s 0m1.212s +gprof icp_map 1 call chamfer 0.33s + icp_lazy 16 call memo 0.21 Index: sandbox/jardonnet/test/icp_check.sh --- sandbox/jardonnet/test/icp_check.sh (revision 1837) +++ sandbox/jardonnet/test/icp_check.sh (working copy) @@ -1,8 +1,8 @@ -#!/bin/sh +#!/bin/zsh for i in `\ls bin/` do - echo execute $i 01.pbm 02.pbm - time ./bin/$i 01.pbm 02.pbm > ./bin/log_$i - echo ./bin/log_$i + echo " " $i + \time -f "%U" ./bin/$i 01.pbm 02.pbm > ./bin/log_$i + echo "" done Index: sandbox/jardonnet/test/Makefile --- sandbox/jardonnet/test/Makefile (revision 1837) +++ sandbox/jardonnet/test/Makefile (working copy) @@ -17,12 +17,15 @@ time ./+sub.exe . . ; time ./+gsub.exe . . - icp: g++ icp.cc -I../../.. -g -o '+icp_map.exe' icp++: - g++ icp.cc -I../../.. -O3 -DNDEBUG -o '+icp_map.exe' + g++ icp.cc -I../../.. -O3 -DNDEBUG -pg -o '+icp_map.exe' + +icp_lazy++: + g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -pg -o '+icp_lazy.exe' + icp_sub: g++ icp_subsampled.cc -I../../.. -O3 -DNDEBUG -g -o '+icp_sub.exe' @@ -42,7 +45,7 @@ g++ icp.cc -I../../.. -O3 -o './bin/+icp_3' g++ icp.cc -I../../.. -O3 -ffloat-store -o './bin/+icp_3f' g++ icp.cc -I../../.. -O3 -DNDEBUG -o './bin/+icp_3D' - g++ icp.cc -I../../.. -O3 -DNDEBUG -ffloat-store -o './bin/+icp_3df' + g++ icp.cc -I../../.. -O3 -DNDEBUG -ffloat-store -o './bin/+icp_3Df' g++ icp_lazy.cc -I../../.. -O0 -o './bin/+icp_lazy_0' g++ icp_lazy.cc -I../../.. -O0 -ffloat-store -o './bin/+icp_lazy_0f' g++ icp_lazy.cc -I../../.. -O0 -DNDEBUG -o './bin/+icp_lazy_0D' @@ -50,7 +53,7 @@ g++ icp_lazy.cc -I../../.. -O3 -o './bin/+icp_lazy_3' g++ icp_lazy.cc -I../../.. -O3 -ffloat-store -o './bin/+icp_lazy_3f' g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -o './bin/+icp_lazy_3D' - g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -ffloat-store -o './bin/+icp_lazy_3df' + g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -ffloat-store -o './bin/+icp_lazy_3Df' ./icp_check.sh clean:

1 0

milena r1837: Fix documentation and guards in p_image2d
by Matthieu Garrigues 03 Apr '08

03 Apr '08

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2008-04-03 Matthieu Garrigues <garrigues(a)lrde.epita.fr> Fix documentation and guards in p_image2d. * mln/core/p_image2d.hh: Fix documentation and guards. * mln/core/p_image2d_piter.hh: Fix documentation. --- p_image2d.hh | 10 +++++----- p_image2d_piter.hh | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) Index: trunk/milena/mln/core/p_image2d.hh =================================================================== --- trunk/milena/mln/core/p_image2d.hh (revision 1836) +++ trunk/milena/mln/core/p_image2d.hh (revision 1837) @@ -25,13 +25,13 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_CORE_P_SET_HH -# define MLN_CORE_P_SET_HH +#ifndef MLN_CORE_P_IMAGE2D_HH +# define MLN_CORE_P_IMAGE2D_HH -/*! \file mln/core/p_set.hh +/*! \file mln/core/p_image2d.hh * - * \brief Definition of a point set class based on std::set. + * \brief Definition of a point set class based on a image of booleans. */ # include <mln/core/internal/point_set_base.hh> @@ -206,4 +206,4 @@ # include <mln/core/p_image2d_piter.hh> -#endif // ! MLN_CORE_P_SET_HH +#endif // ! MLN_CORE_P_IMAGE2D_HH Index: trunk/milena/mln/core/p_image2d_piter.hh =================================================================== --- trunk/milena/mln/core/p_image2d_piter.hh (revision 1836) +++ trunk/milena/mln/core/p_image2d_piter.hh (revision 1837) @@ -29,7 +29,7 @@ # define MLN_CORE_P_IMAGE2D_PITER_HH /// \file mln/core/p_image2d_piter.hh -/// \brief Definition of point iterators on mln::p_array. +/// \brief Definition of point iterators on mln::p_image2d. # include <mln/core/p_image2d.hh> # include <mln/core/box_piter.hh>

1 0

milena r1836: A fast point2d set : p_image2d
by Matthieu Garrigues 03 Apr '08

03 Apr '08

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2008-04-03 Matthieu Garrigues <garrigues(a)lrde.epita.fr> A fast point2d set : p_image2d. * mln/core/p_image2d.hh: The new set. * mln/core/p_image2d_piter.hh: To iter on it. * tests/core/p_image2d.cc: some small tests. --- mln/core/p_image2d.hh | 209 +++++++++++++++++++++++++++++ mln/core/p_image2d_piter.hh | 316 ++++++++++++++++++++++++++++++++++++++++++++ tests/core/p_image2d.cc | 56 +++++++ 3 files changed, 581 insertions(+) Index: trunk/milena/tests/core/p_image2d.cc =================================================================== --- trunk/milena/tests/core/p_image2d.cc (revision 0) +++ trunk/milena/tests/core/p_image2d.cc (revision 1836) @@ -0,0 +1,56 @@ +// 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. + +/*! \file tests/core/p_image2d.cc + * + * \brief Tests on mln::p_image2d. + */ + +#include <mln/core/p_image2d.hh> + +int main() +{ + using namespace mln; + + p_image2d<point2d> ps(20,20); + ps + .insert(make::point2d(6, 9)) + .insert(make::point2d(4, 2)) + .insert(make::point2d(4, 2)) + .insert(make::point2d(5, 1)); + mln_assertion(ps.npoints() == 3); + + ps.remove(make::point2d(5, 1)); + ps.remove(make::point2d(5, 1)); + + ps.remove(make::point2d(6, 9)); + ps.remove(make::point2d(4, 2)); + + mln_assertion(ps.npoints() == 0); + mln_assertion(ps.is_empty()); + +} Index: trunk/milena/mln/core/p_image2d.hh =================================================================== --- trunk/milena/mln/core/p_image2d.hh (revision 0) +++ trunk/milena/mln/core/p_image2d.hh (revision 1836) @@ -0,0 +1,209 @@ +// Copyright (C) 2007, 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_CORE_P_SET_HH +# define MLN_CORE_P_SET_HH + + +/*! \file mln/core/p_set.hh + * + * \brief Definition of a point set class based on std::set. + */ + +# include <mln/core/internal/point_set_base.hh> +# include <mln/core/box2d.hh> +# include <mln/core/image2d.hh> +# include <mln/accu/bbox.hh> +# include <mln/level/fill.hh> +# include <mln/core/p_array.hh> + + +namespace mln +{ + + + // Fwd decls. + template <typename P> struct p_image2d_fwd_piter_; + template <typename P> struct p_image2d_bkd_piter_; + + template <typename P> + class p_image2d : public internal::point_set_base_<P, p_image2d > + { + public: + + /// Forward Point_Iterator associated type. + typedef p_image2d_fwd_piter_ fwd_piter; + + /// Backward Point_Iterator associated type. + typedef p_image2d_bkd_piter_ bkd_piter; + + /// Constructor. + p_image2d(int nrows, int ncols); + p_image2d(const box2d& b); + + /// Insert a point \p p. + p_image2d& insert(const P p); + p_image2d& insert(const p_image2d& set); + + /// Remove a point \p p. + p_image2d& remove(const P p); + p_image2d& remove(const p_image2d& set); + + /// Give the number of points. + unsigned npoints() const; + + /// Test is \p p belongs to this point set. + bool has(const P& p) const; + + /// Test if the set is empty. + bool is_empty() const; + + /// Clear this set. + void clear(); + + /// Hook to the image2d + const image2d<bool>& image() const; + private: + image2d<bool> points_; + unsigned npoints_; + }; + + +# ifndef MLN_INCLUDE_ONLY + + template <typename P> + p_image2d::p_image2d(int nrows, int ncols) + : points_(nrows, ncols), + npoints_(0) + { + + level::fill(points_, false); + } + + template <typename P> + p_image2d::p_image2d(const box2d& b) + : points_(b), + npoints_(0) + { + level::fill(points_, false); + } + + + template <typename P> + p_image2d& + p_image2d::insert(const P p) + { + if (points_(p) == false) + { + points_(p) = true; + npoints_++; + } + return *this; + } + + template <typename P> + p_image2d& + p_image2d::insert(const p_image2d& set) + { + if (set->is_empty()) + return *this; + mln_fwd_piter(image2d<bool>) p(set.points_ | true); + for_all(p) + if (this->points_.has(p)) + this->insert(p); + return *this; + } + + template <typename P> + p_image2d& + p_image2d::remove(const P p) + { + if (points_(p) == true) + { + points_(p) = false; + npoints_--; + } + return *this; + } + + template <typename P> + p_image2d& + p_image2d::remove(const p_image2d& set) + { + if (this->is_empty() || set->is_empty()) + return *this; + mln_fwd_piter(image2d<bool>) p(set.points_ | true); + for_all(p) + if (this->points_.has(p)) + this->remove(p); + return *this; + } + + template <typename P> + unsigned + p_image2d::npoints() const + { + return npoints_; + } + + template <typename P> + bool + p_image2d::has(const P& p) const + { + return points_(p); + } + + template <typename P> + bool + p_image2d::is_empty() const + { + return npoints_ == 0; + } + + template <typename P> + void + p_image2d::clear() + { + level::fill(points_, false); + } + + template <typename P> + const image2d<bool>& + p_image2d::image() const + { + return points_; + } + +# endif // ! MLN_INCLUDE_ONLY + +} // end of namespace mln + + +# include <mln/core/p_image2d_piter.hh> + + +#endif // ! MLN_CORE_P_SET_HH Index: trunk/milena/mln/core/p_image2d_piter.hh =================================================================== --- trunk/milena/mln/core/p_image2d_piter.hh (revision 0) +++ trunk/milena/mln/core/p_image2d_piter.hh (revision 1836) @@ -0,0 +1,316 @@ +// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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_CORE_P_IMAGE2D_PITER_HH +# define MLN_CORE_P_IMAGE2D_PITER_HH + +/// \file mln/core/p_image2d_piter.hh +/// \brief Definition of point iterators on mln::p_array. + +# include <mln/core/p_image2d.hh> +# include <mln/core/box_piter.hh> + + +namespace mln +{ + + /// \brief Forward iterator on points of a p_array. + template <typename P> + struct p_image2d_fwd_piter_ + : public internal::point_iterator_base_< P, p_image2d_fwd_piter_ > + { + typedef p_image2d_fwd_piter_ self_; + typedef internal::point_iterator_base_< P, self_ > super_; + public: + /// The associated psite type. + typedef P psite; + + /// The associated point type. + typedef mln_point(P) point; + + enum { dim = super_::dim }; + + /// Coordinate associated type. + template <typename S> + p_image2d_fwd_piter_(const p_image2d<S>& s); + + /// Reference of the corresponding psite. + const psite& to_psite() const; + + /// Reference of the corresponding point. + const point& to_point() const; + + /// Read-only access to the \p i-th coordinate. + mln_coord(point) operator[](unsigned i) const; + + /// Test if the iterator is valid. + bool is_valid() const; + + /// Invalidate the iterator. + void invalidate(); + + /// Start an iteration. + void start(); + + /// Go to the next point. + void next_(); + + /// Convert the iterator into a psite. + operator psite() const; + + protected: + box_fwd_piter_ piter_; + const image2d<bool> ima_; + }; + + + + /// \brief Backward iterator on points of a p_array. + template <typename P> + struct p_image2d_bkd_piter_ + : public internal::point_iterator_base_< P, p_image2d_bkd_piter_ > + { + typedef p_image2d_bkd_piter_ self_; + typedef internal::point_iterator_base_< P, self_ > super_; + public: + /// The associated psite type. + typedef P psite; + + /// The associated point type. + typedef mln_point(P) point; + + enum { dim = super_::dim }; + + /// Coordinate associated type. + template <typename S> + p_image2d_bkd_piter_(const p_image2d<S>& s); + + /// Reference of the corresponding psite. + const psite& to_psite() const; + + /// Reference of the corresponding point. + const point& to_point() const; + + /// Read-only access to the \p i-th coordinate. + mln_coord(point) operator[](unsigned i) const; + + /// Test if the iterator is valid. + bool is_valid() const; + + /// Invalidate the iterator. + void invalidate(); + + /// Start an iteration. + void start(); + + /// Go to the next point. + void next_(); + + /// Convert the iterator into a psite. + operator psite() const; + + protected: + box_bkd_piter_ piter_; + image2d<bool> ima_; + + }; + + + +# ifndef MLN_INCLUDE_ONLY + + /*------------------------. + | p_image2d_fwd_piter_. | + `------------------------*/ + + template <typename P> + template <typename S> + inline + p_image2d_fwd_piter_::p_image2d_fwd_piter_(const p_image2d<S>& s) + : piter_(exact(s).image().domain()), + ima_(exact(s).image()) + { + invalidate(); + } + + template <typename P> + inline + const P& + p_image2d_fwd_piter_::to_psite() const + { + return piter_.to_psite(); + } + + template <typename P> + inline + const mln_point(P)& + p_image2d_fwd_piter_::to_point() const + { + return piter_.to_point(); + } + + template <typename P> + inline + mln_coord(mln_point_(P)) + p_image2d_fwd_piter_::operator[](unsigned i) const + { + return piter_[i]; + } + + template <typename P> + inline + bool + p_image2d_fwd_piter_::is_valid() const + { + return piter_.is_valid(); + } + + template <typename P> + inline + void + p_image2d_fwd_piter_::invalidate() + { + return piter_.invalidate(); + } + + template <typename P> + inline + void + p_image2d_fwd_piter_::start() + { + piter_.start(); + while(is_valid() && !ima_(piter_)) + piter_.next(); + } + + template <typename P> + inline + void + p_image2d_fwd_piter_::next_() + { + piter_.next(); + while(is_valid() && !ima_(piter_)) + piter_.next(); + } + + template <typename P> + inline + p_image2d_fwd_piter_::operator P() const + { + return P(piter_); + } + + + /*------------------------. + | p_image2d_bkd_piter_. | + `------------------------*/ + + template <typename P> + template <typename S> + inline + p_image2d_bkd_piter_::p_image2d_bkd_piter_(const p_image2d<S>& s) + : piter_(exact(s).image().domain()), + ima_(exact(s).image()) + { + invalidate(); + } + + template <typename P> + inline + const P& + p_image2d_bkd_piter_::to_psite() const + { + return piter_.to_psite(); + } + + template <typename P> + inline + const mln_point(P)& + p_image2d_bkd_piter_::to_point() const + { + return piter_.to_point(); + } + + template <typename P> + inline + mln_coord(mln_point_(P)) + p_image2d_bkd_piter_::operator[](unsigned i) const + { + mln_precondition(is_valid()); + return piter_[i]; + } + + template <typename P> + inline + bool + p_image2d_bkd_piter_::is_valid() const + { + return piter_.is_valid(); + } + + template <typename P> + inline + void + p_image2d_bkd_piter_::invalidate() + { + return piter_.invalidate(); + } + + template <typename P> + inline + void + p_image2d_bkd_piter_::start() + { + piter_.start(); + while(is_valid() && !ima_(piter_)) + piter_.next(); + } + + template <typename P> + inline + void + p_image2d_bkd_piter_::next_() + { + piter_.next(); + while(is_valid() && !ima_(piter_)) + piter_.next(); + } + + template <typename P> + inline + p_image2d_bkd_piter_::operator P() const + { + mln_precondition(is_valid()); + return P(piter_); + } + +# endif // ! MLN_INCLUDE_ONLY + +} // end of namespace mln + + +#endif // ! MLN_CORE_P_IMAGE2D_PITER_HH

1 0

1835: Fix a bug in basic_najman.
by Alexandre Abraham 03 Apr '08

03 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Alexandre Abraham <abraham(a)lrde.epita.fr> Fix a bug in basic_najman. * sandbox/abraham/morpho/basic_najman.hh: (mln::morpho::basic_najman<I,N>::BuildComponentTree): Fix a bug, the algorithm is still bugged but now gives a nice tree. * sandbox/abraham/morpho/test.cc: Add picture inversion for test purpose. basic_najman.hh | 103 +++++++++++++++++++++++++++++++++++--------------------- test.cc | 10 ++++- 2 files changed, 75 insertions(+), 38 deletions(-) Index: sandbox/abraham/morpho/basic_najman.hh --- sandbox/abraham/morpho/basic_najman.hh (revision 1834) +++ sandbox/abraham/morpho/basic_najman.hh (working copy) @@ -37,13 +37,13 @@ const Image& ima; const Neighborhood<N>& nbh; - mln_ch_value(I, psite) par; + mln_ch_value(I, psite) Par_node; mln_ch_value(I, psite) Par_tree; // Par_node is handled by par (from ap_maxtree) mln_ch_value(I, int) Rnk_tree; mln_ch_value(I, int) Rnk_node; - mln_ch_value(I, psite) lowestNode; + mln_ch_value(I, psite) subtreeRoot; mln_ch_value(I, node) nodes; mln_ch_value(I, bool) isproc; psite Root; @@ -54,11 +54,11 @@ const Neighborhood<N>& n) : ima(i), nbh(n), - par(exact(i).domain(), exact(i).border()), + Par_node(exact(i).domain(), exact(i).border()), Par_tree(exact(i).domain(), exact(i).border()), Rnk_tree(exact(i).domain(), exact(i).border()), Rnk_node(exact(i).domain(), exact(i).border()), - lowestNode(exact(i).domain(), exact(i).border()), + subtreeRoot(exact(i).domain(), exact(i).border()), nodes(exact(i).domain(), exact(i).border()), isproc(exact(i).domain(), exact(i).border()) { @@ -73,7 +73,7 @@ void MakeSet_node(psite x) { - par(x) = x; // was Par_node(x) = x; + Par_node(x) = x; // was Par_node(x) = x; Rnk_node(x) = 0; } @@ -91,29 +91,29 @@ y = memo; } - bool is_root(const psite& x) const - { - return par(x) == x; - } +// bool is_root(const psite& x) const +// { +// return Par_node(x) == x; +// } - bool is_level_root(const psite& x) const - { - return is_root(x) || exact(ima)(par(x)) < exact(ima)(x); - } +// bool is_level_root(const psite& x) const +// { +// return is_root(x) || exact(ima)(Par_node(x)) < exact(ima)(x); +// } - psite find_level_root(const psite& x) - { - if (is_level_root(x)) - return x; - else - return par(x) = find_level_root(par(x)); - } +// psite find_level_root(const psite& x) +// { +// if (is_level_root(x)) +// return x; +// else +// return Par_node(x) = find_level_root(Par_node(x)); +// } psite Find_node(psite x) { - if (par(x) != x) - par(x) = Find_node(par(x)); - return par(x); + if (Par_node(x) != x) + Par_node(x) = Find_node(Par_node(x)); + return Par_node(x); } void go() @@ -135,8 +135,8 @@ psite p = S[ip]; MakeSet_node(p); MakeSet_tree(p); - // if (lowestNode.hold(p)) - lowestNode(p) = p; + // if (subtreeRoot.hold(p)) + subtreeRoot(p) = p; // if (nodes.hold(p)) nodes(p) = MakeNode(exact(ima)(p)); } @@ -148,22 +148,32 @@ for (int ip = 0; ip < int(S.npoints()); ++ip) { psite p = S[ip]; - isproc(p) = true; + // isproc(p) = true; psite curTree = Find_tree(p); - psite curNode = Find_node(lowestNode(curTree)); + psite curNode = Find_node(subtreeRoot(curTree)); mln_niter(N) q(nbh, p); for_all(q) - if (exact(ima).has(q) and isproc(q) and exact(ima)(q) >= exact(ima)(p)) // f[q] >= f[p] parce qu'on prend l'ordre de l'histogramme + if (exact(ima).has(q) and isproc(q) and exact(ima)(q) <= exact(ima)(p)) { psite adjTree = Find_tree(q); - psite adjNode = Find_node(lowestNode(adjTree)); + psite adjNode = Find_node(subtreeRoot(adjTree)); if (curNode != adjNode) { + std::cout << "curNode != adjNode; " << nodes(curNode).level << " vs. " << nodes(adjNode).level << std::endl; if (nodes(curNode).level == nodes(adjNode).level) - curNode = MergeNode(adjNode, curNode); + { + // curNode = MergeNode(adjNode, curNode); + psite tmpNode = MergeNode(adjNode, curNode); + /*if (tmpNode == curNode) + nodes(curNode).addChildren(nodes(adjNode)); + else + nodes(adjNode).addChildren(nodes(curNode)); + nodes(adjNode) = nodes(curNode);*/ + curNode = tmpNode; + } else { // we have nodes[curNode].level < nodes[adjNode].level @@ -172,12 +182,16 @@ nodes(curNode).area += nodes(adjNode).area; nodes(curNode).highest += nodes(adjNode).highest; } + // curTree = Link_tree(adjTree, curTree); + // subtreeRoot(curTree) = curNode; + } curTree = Link_tree(adjTree, curTree); - lowestNode(curTree) = curNode; + subtreeRoot(curTree) = curNode; } + isproc(p) = true; } - Root = lowestNode(Find_tree(Find_node(psite(0, 0)))); + Root = subtreeRoot(Find_tree(Find_node(psite(0, 0)))); // Pour garder une map de correspondance point <-> local_root // for (int ip = 0; ip < int(S.size()); ++ip) // { @@ -185,9 +199,10 @@ // M(p) = Find_node(p); // } - mln_piter(I) r(par.domain()); + mln_piter(I) r(Par_node.domain()); for_all(r) - par(r) = Find_node(r); + Par_node(r) = Find_node(r); + } @@ -228,7 +243,7 @@ else if (Rnk_node(x) == Rnk_node(y)) Rnk_node(y) += 1; - par(x) = y; // was Par_node(x) = y; + Par_node(x) = y; // was Par_node(x) = y; return y; } @@ -260,7 +275,7 @@ void print_tree(psite p) { node& n = nodes(p); - std::cout << "psite " << p << "=" << (p.row() * exact(lowestNode).domain().len(1) + p.col()) << " : "; + std::cout << "psite " <::fwd_piter it(n.children); for (it.start(); @@ -268,7 +283,7 @@ it.next()) { psite q = it.to_psite(); - std::cout << q << "=" << (q.row() * lowestNode.domain().len(2) + q.col()) << " "; + std::cout << q << "=" << (q.row() * subtreeRoot.domain().len(1) + q.col()) << " "; } std::cout << std::endl; @@ -280,6 +295,20 @@ } } + void print_sup_tree() + { + psite max = psite(0, 0); + + for(unsigned int i = 0; i < nodes.domain().len(0); ++i) + for(unsigned int j = 0; j < nodes.domain().len(1); ++j) + { + if (nodes(psite(i, j)).area > nodes(max).area) + max = psite(i, j); + } + + print_tree(max); + } + }; // struct basic_najman }; // namespace morpho Index: sandbox/abraham/morpho/test.cc --- sandbox/abraham/morpho/test.cc (revision 1834) +++ sandbox/abraham/morpho/test.cc (working copy) @@ -25,16 +25,24 @@ image2d<int_u8> input; io::pgm::load(input, "./images/test_component_tree.pgm"); + + for(unsigned int i = 0; i < input.domain().len(0); ++i) + for(unsigned int j = 0; j < input.domain().len(1); ++j) + input.at(i, j) = 255 - input.at(i, j); + + morpho::basic_najman< image2d<int_u8>, neighb2d> n(input, c8()); n.image_output(input); n.go(); n.image_output(n.Par_tree); n.print_tree(n.Root); - n.image_output(n.par); + n.image_output(n.Par_node); std::cout << " --------------------------- " << std::endl; n.image_output(n.Rnk_tree); std::cout << " --------------------------- " << std::endl; n.image_output(n.Rnk_node); + std::cout << " --------------------------- " << std::endl; + n.print_sup_tree(); // image2d<int_u8> output = morpho::topo_wst(input, c4()); // io::pgm::save(output, "out.pgm");

1 0