Olena-patches February 2009

olena-patches@lrde.epita.fr

12 participants
266 discussions

by Guillaume Lazzara

* doc/tutorial/Makefile.am: fix wrong include path. * mln/canvas/morpho/algebraic_union_find.hh: add missing MLN_INCLUDE_ONLY guards. * mln/core/clock_neighb.hh: add missing is_valid(). * mln/labeling/mean_values.hh, * mln/fun/v2v/rgb_to_hsl.hh: add missing includes. * mln/fun/x2x/rotation.hh: add more comments. * mln/io/dump/load.hh: check if the loaded image is valid. * mln/labeling/level.hh: Fix warnings. * mln/registration/all.hh: update comments. * mln/registration/icp2.hh: rename as... * mln/registration/icp.hh: ... this. * mln/registration/registration.hh: fix wrong call to registration_tests. * mln/util/graph.hh: add NDEBUG guards. * mln/value/rgb.hh: add missing forward declaration. --- milena/ChangeLog | 32 + milena/doc/tutorial/Makefile.am | 22 +- milena/mln/canvas/morpho/algebraic_union_find.hh | 12 + milena/mln/core/clock_neighb.hh | 18 +- milena/mln/fun/v2v/rgb_to_hsl.hh | 3 +- milena/mln/fun/x2x/rotation.hh | 3 + milena/mln/io/dump/load.hh | 2 + milena/mln/labeling/level.hh | 4 +- milena/mln/labeling/mean_values.hh | 9 +- milena/mln/registration/all.hh | 10 +- milena/mln/registration/icp.hh | 723 ++++++++++++++++++---- milena/mln/registration/icp2.hh | 693 --------------------- milena/mln/registration/registration.hh | 9 +- milena/mln/util/graph.hh | 2 + milena/mln/value/rgb.hh | 5 + 15 files changed, 692 insertions(+), 855 deletions(-) delete mode 100644 milena/mln/registration/icp2.hh diff --git a/milena/ChangeLog b/milena/ChangeLog index 1c446cd..6531f76 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,37 @@ 2009-02-18 Guillaume Lazzara <z(a)lrde.epita.fr> + Various small fixes. + + * doc/tutorial/Makefile.am: fix wrong include path. + + * mln/canvas/morpho/algebraic_union_find.hh: add missing + MLN_INCLUDE_ONLY guards. + + * mln/core/clock_neighb.hh: add missing is_valid(). + + * mln/labeling/mean_values.hh, + * mln/fun/v2v/rgb_to_hsl.hh: add missing includes. + + * mln/fun/x2x/rotation.hh: add more comments. + + * mln/io/dump/load.hh: check if the loaded image is valid. + + * mln/labeling/level.hh: Fix warnings. + + * mln/registration/all.hh: update comments. + + * mln/registration/icp2.hh: rename as... + * mln/registration/icp.hh: ... this. + + * mln/registration/registration.hh: fix wrong call to + registration_tests. + + * mln/util/graph.hh: add NDEBUG guards. + + * mln/value/rgb.hh: add missing forward declaration. + +2009-02-18 Guillaume Lazzara <z(a)lrde.epita.fr> + Add support for 3D images in fun::x2v::bilinear. * mln/fun/x2v/bilinear.hh: Add a new specialization for 3d Images. diff --git a/milena/doc/tutorial/Makefile.am b/milena/doc/tutorial/Makefile.am index ceb4534..447dcbe 100644 --- a/milena/doc/tutorial/Makefile.am +++ b/milena/doc/tutorial/Makefile.am @@ -41,24 +41,24 @@ tutorial.tex \ tutorial.pdf .PHONY: regen-dist -regen-dist: $(srcdir)/headers.stamp +regen-dist: $(top_srcdir)/headers.stamp -$(srcdir)/examples/examples.mk: $(srcdir)/headers.stamp -$(srcdir)/figures/figures.mk: $(srcdir)/headers.stamp -$(srcdir)/outputs/outputs.mk: $(srcdir)/headers.stamp -$(srcdir)/samples/samples.mk: $(srcdir)/headers.stamp +$(top_srcdir)/milena/doc/tutorial/examples/examples.mk: $(top_srcdir)/headers.stamp +$(top_srcdir)/milena/doc/tutorial/figures/figures.mk: $(top_srcdir)/headers.stamp +$(top_srcdir)/milena/doc/tutorial/outputs/outputs.mk: $(top_srcdir)/headers.stamp +$(top_srcdir)/milena/doc/tutorial/samples/samples.mk: $(top_srcdir)/headers.stamp -EXTRA_DIST += $(srcdir)/headers.stamp -$(srcdir)/headers.stamp: $(srcdir)/generate_dist_files.sh +EXTRA_DIST += $(top_srcdir)/headers.stamp +$(top_srcdir)/headers.stamp: $(top_srcdir)/milena/doc/tutorial/generate_dist_files.sh @rm -f $@.tmp @touch $@.tmp cd $(srcdir) && ./generate_dist_files.sh @mv -f $@.tmp $@ -include $(srcdir)/examples/examples.mk -include $(srcdir)/figures/figures.mk -include $(srcdir)/outputs/outputs.mk -include $(srcdir)/samples/samples.mk +include $(top_srcdir)/milena/doc/tutorial/examples/examples.mk +include $(top_srcdir)/milena/doc/tutorial/figures/figures.mk +include $(top_srcdir)/milena/doc/tutorial/outputs/outputs.mk +include $(top_srcdir)/milena/doc/tutorial/samples/samples.mk EXTRA_DIST += \ tools/sample_utils.hh \ diff --git a/milena/mln/canvas/morpho/algebraic_union_find.hh b/milena/mln/canvas/morpho/algebraic_union_find.hh index 34a922c..a19ace8 100644 --- a/milena/mln/canvas/morpho/algebraic_union_find.hh +++ b/milena/mln/canvas/morpho/algebraic_union_find.hh @@ -51,6 +51,16 @@ namespace mln namespace morpho { + template <typename I, typename N, typename F> + inline + mln_concrete(I) + algebraic_union_find(const Image& input, + const Neighborhood<N>& nbh, + F& f); + + +# ifndef MLN_INCLUDE_ONLY + namespace impl { @@ -393,6 +403,8 @@ namespace mln } +# endif // ! MLN_INCLUDE_ONLY + } // end of namespace mln::canvas::morpho } // end of namespace mln::canvas diff --git a/milena/mln/core/clock_neighb.hh b/milena/mln/core/clock_neighb.hh index 2c72506..1474978 100644 --- a/milena/mln/core/clock_neighb.hh +++ b/milena/mln/core/clock_neighb.hh @@ -1,5 +1,5 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory -// (LRDE) +// Copyright (C) 2007, 2008, 2009 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 @@ -105,6 +105,9 @@ namespace mln /// FIXME: not in constant time! mln::window<D> win() const; + /// Return whether this neighborhood is valid. + bool is_valid() const; + private: std::vector<D> vec_; }; @@ -145,6 +148,17 @@ namespace mln result.insert(vec_[i]); return result; } + + template <typename D> + inline + bool + clock_neighb<D>::is_valid() const + { + //FIXME: correct? + return true; + } + + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln diff --git a/milena/mln/fun/v2v/rgb_to_hsl.hh b/milena/mln/fun/v2v/rgb_to_hsl.hh index 73a68f2..39433f2 100644 --- a/milena/mln/fun/v2v/rgb_to_hsl.hh +++ b/milena/mln/fun/v2v/rgb_to_hsl.hh @@ -36,6 +36,8 @@ # include <mln/trait/value_.hh> +# include <mln/value/rgb8.hh> +# include <mln/value/rgb16.hh> namespace mln { @@ -43,7 +45,6 @@ namespace mln // Forward declaration namespace value { - template <unsigned n> struct rgb; template <typename H, typename S, typename L> class hsl_; typedef hsl_<float, float, float> hsl_f; } diff --git a/milena/mln/fun/x2x/rotation.hh b/milena/mln/fun/x2x/rotation.hh index 9628c26..9addb50 100644 --- a/milena/mln/fun/x2x/rotation.hh +++ b/milena/mln/fun/x2x/rotation.hh @@ -35,6 +35,9 @@ /// /// \todo store the quaternion instead of (axis, alpha) /// => better precision while composing two rotation matrices. +/// +/// Conversion from Quaternion to (angle,axis), Source: +/// http://jeux.developpez.com/faq/matquat/?page=quaternions#Q56 # include <cmath> # include <mln/core/concept/function.hh> diff --git a/milena/mln/io/dump/load.hh b/milena/mln/io/dump/load.hh index 7bbce9e..21d6534 100644 --- a/milena/mln/io/dump/load.hh +++ b/milena/mln/io/dump/load.hh @@ -129,6 +129,8 @@ namespace mln internal::load_header(ima, file); internal::load_data(ima, file); + mln_postcondition(exact(ima).is_valid()); + trace::exiting("mln::io::dump::load"); } diff --git a/milena/mln/labeling/level.hh b/milena/mln/labeling/level.hh index 050b2d0..11b3a06 100644 --- a/milena/mln/labeling/level.hh +++ b/milena/mln/labeling/level.hh @@ -121,9 +121,9 @@ namespace mln bool handles_(unsigned p) const { return input.element(p) == val; } bool equiv_(unsigned n, unsigned) const { return input.element(n) == val; } bool labels_(unsigned) const { return true; } - void do_no_union_(unsigned n, unsigned p) {} + void do_no_union_(unsigned, unsigned) {} void init_attr_(unsigned) {} - void merge_attr_(unsigned r, unsigned p) {} + void merge_attr_(unsigned, unsigned) {} // end of Requirements. diff --git a/milena/mln/labeling/mean_values.hh b/milena/mln/labeling/mean_values.hh index d7c131a..2ea06e4 100644 --- a/milena/mln/labeling/mean_values.hh +++ b/milena/mln/labeling/mean_values.hh @@ -29,6 +29,13 @@ # define MLN_LABELING_MEAN_VALUES_HH # include <mln/core/concept/image.hh> +# include <mln/core/alias/vec3d.hh> + +# include <mln/accu/mean.hh> + +# include <mln/labeling/compute.hh> + +# include <mln/literal/colors.hh> /// \file mln/labeling/mean_values.hh /// @@ -116,7 +123,7 @@ namespace mln util::array<mln_value(I)> m; convert::from_to(m_3f, m); - m[0] = literal::white; //FIXME: handle label 0 correctly. + m[0] = literal::yellow; //FIXME: handle label 0 correctly. mln_concrete(I) output = level::transform(lbl, convert::to< fun::i2v::array<mln_value(I)> >(m)); diff --git a/milena/mln/registration/all.hh b/milena/mln/registration/all.hh index dd29f26..bb4e18e 100644 --- a/milena/mln/registration/all.hh +++ b/milena/mln/registration/all.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory +// Copyright (C) 2008, 2009 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 @@ -28,10 +29,9 @@ #ifndef MLN_REGISTRATION_ALL_HH # define MLN_REGISTRATION_ALL_HH -/*! \file mln/registration/all.hh - * - * \brief File that includes all "point-wise" expression tools. - */ +/// \file mln/registration/all.hh +/// +/// File that includes all "point-wise" expression tools. namespace mln diff --git a/milena/mln/registration/icp.hh b/milena/mln/registration/icp.hh index b116b0e..de763ea 100644 --- a/milena/mln/registration/icp.hh +++ b/milena/mln/registration/icp.hh @@ -1,4 +1,4 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory +// Copyright (C) 2009 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 @@ -28,217 +28,666 @@ #ifndef MLN_REGISTRATION_ICP_HH # define MLN_REGISTRATION_ICP_HH -/*! \file mln/registration/icp.hh - * - * \brief Register an image over an another using the ICP algorithm. - */ +/// \file mln/registration/icp.hh +/// +/// Register an image over an another using the ICP algorithm. +/// +/// \todo encode distances on 12 bits. +# include <cmath> +# include <algorithm> + +# include <mln/core/alias/vec3d.hh> +# include <mln/math/jacobi.hh> # include <mln/fun/x2x/all.hh> # include <mln/fun/x2v/all.hh> -# include <mln/registration/get_rtransf.hh> -# include <mln/registration/internal/rms.hh> +# include <mln/convert/to.hh> +# include <mln/accu/compute.hh> +# include <mln/accu/center.hh> +# include <mln/accu/rms.hh> +# include <mln/trait/image_from_grid.hh> +# include <mln/set/compute.hh> + +//Should be removed when closest_point functors are moved. +# include <mln/core/alias/neighb3d.hh> +# include <mln/transform/internal/closest_point_functor.hh> +# include <mln/canvas/distance_geodesic.hh> +# include <mln/pw/all.hh> + +# include <mln/io/ppm/save.hh> +# include <mln/io/pbm/save.hh> +# include <mln/debug/colorize.hh> + +# include <mln/literal/black.hh> +# include <mln/literal/white.hh> +# include <mln/literal/colors.hh> + +# include <mln/core/image/slice_image.hh> + +#include <mln/core/image/tr_image.hh> +#include <mln/core/image/extension_fun.hh> + +#include <mln/accu/histo.hh> +#include <mln/accu/sum.hh> +#include <mln/debug/histo.hh> + +# include <mln/util/timer.hh> namespace mln { - namespace util + + namespace registration { - // FIXME: Remove use of this class - template <unsigned int n, typename T> - struct buffer + //FIXME: used for debug purpose. Should be removed later. + + using namespace fun::x2x; + + /*! Register point in \p c using a function of closest points + * \p closest_point. + * + * \param[in] P_ The cloud of points. + * \param[in] X the reference surface. + * \param[in] closest_point The function of closest points. + * \param[out] qk The rigid transformation obtained. + * + * \return the rigid transformation which may be use later to create + * a registered image. + * + * WARNING: the function \p closest_point *MUST* take float/double + * vector as arguments. Otherwise the resulting transformation may be + * wrong due to the truncation of the vector coordinate values. + * + * \pre \p P_ and \p X must not be empty. + * + * Reference article: "A Method for Registration of 3-D Shapes", Paul J. + * Besl and Neil D. McKay, IEEE, 2, February 1992. + * + */ + template <typename P, typename F> + std::pair<algebra::quat,mln_vec(P)> + icp(const p_array& P_, + const p_array& X, + const F& closest_point, + const algebra::quat& initial_rot, + const mln_vec(P)& initial_translation); + + + template <typename P, typename F> + composed< translation<P::dim,float>,rotation<P::dim,float> > + icp(const p_array& P_, + const p_array& X, + const F& closest_point); + + +# ifndef MLN_INCLUDE_ONLY + + + /// Closest point functor based on map distance. + template <typename P> + class closest_point_with_map { - buffer() - : setted(0) + typedef mln_image_from_grid(mln_grid(P), P) I; + typedef mln_ch_value(I, P) cp_ima_t; + typedef mln_ch_value(I,value::int_u16) dmap_t; + + public: + + closest_point_with_map(const p_array& X) { + box3d box = geom::bbox(X); + box.enlarge(1, box.nrows()); + box.enlarge(2, box.ncols()); + + mln_postcondition(box.is_valid()); + + std::cout << "Map image defined on " << box << std::endl; + + init(X, box); } - void store(T e) + closest_point_with_map(const p_array& X, const box& box) { - for (unsigned i = 0; i < n-1; i++) - buf[i+1] = buf[i]; - buf[0] = e; - - setted++; + init(X, box); } - T& operator[](unsigned int i) + void init(const p_array& X, const box& box) { - assert(i < n && i < setted); - return buf[i]; + typedef mln_ch_value(I, bool) model_t; + model_t model(box); + data::fill(model, false); + data::fill((model | X).rw(), true); + + transform::internal::closest_point_functor<model_t> f; + util::timer t; + t.start(); + dmap_X_ = canvas::distance_geodesic(model, c6(), + mln_max(value::int_u16), + f); + std::cout << "canvas::distance_geodesic - " << t << "s" << std::endl; + + initialize(this->cp_ima_, f.cp_ima); + this->cp_ima_ = f.cp_ima; + + mln_postcondition(cp_ima_.is_valid()); + mln_postcondition(cp_ima_.domain().is_valid()); + std::cout << "pmin = " << cp_ima_.domain().pmin() << std::endl;; + std::cout << "pmax = " << cp_ima_.domain().pmax() << std::endl;; + +#ifndef NDEBUG +// mln_ch_value(I, bool) debug2(box); +// data::fill(debug2, false); +// mln_ch_value(I, value::rgb8) debug(box); +// mln_piter(p_array) p(X); +// for_all(p) +// { +// debug(p) = debug::internal::random_color(value::rgb8()); +// debug2(p) = true; +// } +// io::pbm::save(slice(debug2,0), "debug2-a.ppm"); +// +// mln_piter(I) pi(cp_ima_.domain()); +// for_all(pi) +// { +// debug(pi) = debug(cp_ima_(pi)); +// debug2(pi) = debug2(cp_ima_(pi)); +// } +// +// io::pbm::save(slice(debug2,0), "debug2-b.ppm"); +// io::ppm::save(slice(debug,0), "debug.ppm"); +// std::cout << "map saved" << std::endl; +#endif } - const T * get_array() + mln_site(I) operator()(const mln_site(I)& p) const { - return buf; + return cp_ima_(p); } + // Distance map + dmap_t dmap_X_; + // Closest point image. + cp_ima_t cp_ima_; private: - T buf[n]; - unsigned int setted; }; - } // end of namespace mln::util - namespace registration - { + /// Closest point functor based on map distance. + template <typename P> + class closest_point_basic + { + typedef mln_image_from_grid(mln_grid(P), P) I; + typedef p_array X_t; - using namespace fun::x2x; + public: - /*! Register point in \p c using a map of closest points \p map - * - * \param[in] c The cloud of points. - * \param[in] map The map of closest points. - * \param[in] epsilon ICP stops if sqr_norm(qk - qk-1) / - * sqr_norm(qk) > epsilon - * \param[out] qk The rigid transformation obtained. - * - * \pre \p ima has to be initialized. - */ - template <typename P, typename M> - inline - composed<rotation<P::dim, float>, translation<P::dim, float> > - icp(const p_array& c, const M& map, - const float epsilon = 1e-3); + closest_point_basic(const p_array& X) + : X_(X) + { + } + mln_site(I) operator()(const vec3d_f& v) const + { + vec3d_f best_x = convert::to<vec3d_f>(X_[0].to_vec()); + + float best_d = norm::l2_distance(v, best_x); + mln_piter(X_t) X_i(X_); + for_all(X_i) + { + float d = norm::l2_distance(v, convert::to<vec3d_f>(X_i)); + if (d < best_d) + { + best_d = d; + best_x = X_i.to_vec(); + } + } + return best_x; + } - /*! Register point in \p c using a map of closest points \p map - * - * \param[in] c The cloud of points. - * \param[in] c_length points from cloud to use. - * \param[in] map The map of closest points. - * \param[in] qk The initial rigid transformation applied. - * \param[in] epsilon ICP stops if sqr_norm(qk - qk-1) / - * sqr_norm(qk) > epsilon - * \param[out] qk The rigid transformation obtained. - * - * \pre \p ima has to be initialized. - */ - template <typename P, typename M, typename T> - inline + private: + const p_array& X_; + }; + + + template <typename P> void - icp_subset(const p_array& c, - const unsigned int c_length, - const M& map, - T& qk, - float epsilon = 1e-3); + draw_last_run(const box3d& box, const p_array& kept, + const p_array& removed, const p_array& X, + const algebra::quat& qR, const vec3d_f qT) + { + typedef image3d<value::rgb8> result_t; + result_t result(box); + typedef extension_fun<result_t,pw::cst_<mln_value(result_t)> > ext_result_t; + ext_result_t ext_result(result, pw::cst(value::rgb8(0,0,0))); + data::fill(ext_result, literal::black); + data::fill((ext_result | X).rw(), literal::white); -# ifndef MLN_INCLUDE_ONLY + mln_piter(p_array) p(kept); + for_all(p) + ext_result(qR.rotate(p.to_vec()) + qT) = literal::green; + + mln_piter(p_array) p2(removed); + for_all(p2) + ext_result(qR.rotate(p2.to_vec()) + qT) = literal::red; - namespace impl + io::ppm::save(slice(ext_result,0), "registered-2.ppm"); + } + + + + template <typename P, typename F> + float compute_standard_deviation(const p_array& P_, + const std::pair<algebra::quat,mln_vec(P)>& pair, + const F& closest_point) { + accu::rms<vec3d_f,float> e_k_accu; - // FIXME: Move elsewhere - template <typename P> - algebra::vec<P::dim,float> - center(const p_array& a) + // Standard Deviation + float sd; + mln_piter(p_array) p(P_); + for_all(p) { - algebra::vec<P::dim,float> c(literal::zero); - for (unsigned i = 0; i < a.nsites(); ++i) - c += convert::to< algebra::vec<P::dim,float> > (a[i]); + vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; + vec3d_f Yk_i = closest_point(Pk_i).to_vec(); + // yk_i - pk_i + e_k_accu.take(Yk_i - Pk_i); + } + + float d = e_k_accu.to_result(); + sd = std::sqrt(e_k_accu.hook_value_() / P_.nsites() - d * d); + return sd; + } + - return c / a.nsites(); + template <typename P, typename F> + void + remove_too_far_sites_debug(image3d<value::rgb8>& out, const p_array& P_, + const F& closest_point, + const std::pair<algebra::quat,mln_vec(P)>& pair, + const p_array& X, + unsigned r, int d_min, int d_max, unsigned prefix) + { + unsigned removed = 0; + accu::histo<value::int_u8> h; + mln_piter(p_array) p(P_); + data::fill(out, literal::black); + data::fill((out | X).rw(), literal::white); + + for_all(p) + { + vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; + vec3d_f Yk_i = closest_point(Pk_i).to_vec(); + + int d_i = closest_point.dmap_X_(Pk_i); + if (d_i >= d_min && d_i <= d_max) + out(Pk_i) = literal::green; + else + { + ++removed; + out(Pk_i) = literal::red; + } } - // FIXME: Make use of level::apply - template <typename P, typename F> - void apply(F& f, - const p_array& a, - p_array& b) +#ifndef NDEBUG + std::ostringstream ss1; + ss1 << "histo_" << prefix << r << ".dat"; + std::cout << h << std::endl; + + std::ostringstream ss2; + ss2 << "out_" << prefix << r << ".ppm"; + io::ppm::save(mln::slice(out,0), ss2.str()); +#endif + std::cout << "Points removed with the whole set and current d_min/d_max: " << removed << std::endl; + + } + + + template <typename P, typename F> + void + compute_distance_criteria(const p_array& P_, + const F& closest_point, + const std::pair<algebra::quat,mln_vec(P)>& pair, + unsigned r, int& d_min, int& d_max) + { + mln_piter(p_array) p(P_); + accu::histo<value::int_u8> h; + + float sd; { - for (unsigned i = 0; i < a.nsites(); i++) - b[i] = f(convert::to< algebra::vec<P::dim,float> >(a[i])); + accu::sum<float> s, s2; + for_all(p) + { + vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; + unsigned d_i = closest_point.dmap_X_(Pk_i); + h.take(d_i); + s.take(d_i); + s2.take(d_i * d_i); + } + float mean = s / P_.nsites(); + sd = std::sqrt(s2 / P_.nsites() - mean * mean); + d_min = int(mean - sd); + d_max = int(mean + sd); } - template <typename P, typename M, typename T> - inline - void - icp_(const p_array& c_, - const unsigned c_length, - const M& map, - T& qk, - const float epsilon = 1e-3) + std::cout << "Standard deviation = " << sd << std::endl; + std::ostringstream ss1; + ss1 << "histo_" << r << ".dat"; + std::cout << h << std::endl; + std::cout << "d thresholds = " << d_min << ' ' << d_max << std::endl; + } + + template <typename P, typename F> + p_array + remove_too_far_sites(image3d<value::rgb8>& out, const p_array& P_, + const F& closest_point, + const std::pair<algebra::quat,mln_vec(P)>& pair, + const p_array& X, p_array& removed_set, + unsigned r, int d_min, int d_max, + const std::string& method) + { + p_array tmp; + unsigned removed = 0; + +# ifndef NDEBUG + data::fill(out, literal::black); + data::fill((out | X).rw(), literal::white); +# endif // ! NDEBUG + + mln_piter(p_array) p(P_); + for_all(p) { - trace::entering("registration::impl::icp_"); + vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; + vec3d_f Yk_i = closest_point(Pk_i).to_vec(); + + int d_i = closest_point.dmap_X_(Pk_i); + if (d_i >= d_min && d_i <= d_max) + { + tmp.append(p); + out(Pk_i) = literal::green; + } + else + { + ++removed; + removed_set.append(p); + out(Pk_i) = literal::red; + } + } + - util::buffer<4,T> buf_qk; - util::buffer<3,float> buf_dk; +# ifndef NDEBUG + std::ostringstream ss2; + ss2 << method << "_" << r << "_removed_sites" << ".ppm"; + io::ppm::save(mln::slice(out,0), ss2.str()); + + std::cout << "Points removed: " << removed << std::endl; +# endif // ! NDEBUG + // They are used for debug purpose only. + // When NDEBUG is set, they are unused. + (void) X; + (void) r; + (void) method; + + return tmp; + } - float d_k = 10000; + template <typename P> + void + display_sites_used_in_icp(image3d<value::rgb8>& out, const p_array& P_sub, + const p_array& P_, const p_array& X, + unsigned r, const std::string& prefix, + const std::pair<algebra::quat,mln_vec(P)>& pair, + const std::string& period, const value::rgb8& c) + { +# ifndef NDEBUG + data::fill(out, literal::black); + data::fill((out | X).rw(), literal::white); +# endif // !NDEBUG - //work on a reduced version of c_ - p_array c = c_; - //FIXME: Use c_ and specify c_length to function using c_ - c.hook_std_vector_().resize(c_length); + mln_piter(p_array) p1(P_); + for_all(p1) + { + vec3d_f Pk_i = pair.first.rotate(p1.to_vec()) + pair.second; + out(Pk_i) = literal::red; + } - //FIXME: avoid use of ck - p_array ck(c); - algebra::vec<P::dim,float> mu_c = center(c); + mln_piter(p_array) p2(P_sub); + for_all(p2) + { + vec3d_f Pk_i = pair.first.rotate(p2.to_vec()) + pair.second; + out(Pk_i) = c; + } - buf_qk.store(qk); +# ifndef NDEBUG + std::ostringstream ss; + ss << prefix << "_" << r << "_" << period << ".ppm"; - apply(qk, c, ck); + io::ppm::save(slice(out,0), ss.str()); +# endif // ! NDEBUG + } - do { - std::cout << "*------" << std::endl; - buf_dk.store(d_k); - //compute qk - qk = get_rtransf(c, mu_c, ck, map); - buf_qk.store(qk); + template <typename P, typename F> + inline + float + compute_d_k(const p_array& P_, + const F& closest_point, + const algebra::quat& qR, + const algebra::quat& qR_old, + const vec3d_f& qT, + const vec3d_f& qT_old) + { + accu::rms<vec3d_f, float> accu; + mln_piter(p_array) p(P_); + for_all(p) + { + // yk_i - pk+1_i + vec3d_f Pk_i = qR_old.rotate(p.to_vec()) + qT_old; + vec3d_f Pk_1_i = qR.rotate(p.to_vec()) + qT; + accu.take(closest_point(Pk_i).to_vec() - Pk_1_i); + } + return accu.to_result(); + } - //Ck+1 = qk(c) - apply(qk, c, ck); - // d_k = d(Yk, Pk+1) - // = d(closest(qk-1(P)), qk(P)) - d_k = internal::rms(c, map, buf_qk[1], qk); + /// FIXME: work only for 3d images. + template <typename P, typename F> + algebra::quat + get_rot(const p_array& P_, + const vec3d_f& mu_P, + const vec3d_f& mu_Yk, + const F& closest_point, + const algebra::quat& qR, + const vec3d_f& qT) + { + /// Spx: cross-covariance matrix. + algebra::mat<3u,3u,float> Spx; + mln_piter(p_array) p(P_); - //FIXME: Add matrix norm - //} while ((qk - buf_qk[1]).sqr_norm() / qk.sqr_norm() > epsilon); + // FIXME: could we use an accu? + for_all(p) + { + vec3d_f P_i = p.to_vec(); + vec3d_f Pk_i = qR.rotate(p.to_vec()) + qT; + vec3d_f Yk_i = closest_point(Pk_i).to_vec(); + Spx += make::mat(P_i - mu_P) * trans(make::mat(Yk_i - mu_Yk)); + } + Spx /= P_.nsites(); - std::cout << d_k << std::endl; + vec3d_f A; + A[0] = Spx(1,2) - Spx(2,1); + A[1] = Spx(2,0) - Spx(0,2); + A[2] = Spx(0,1) - Spx(1,0); - } while (d_k > 9.2 + epsilon); + algebra::mat<4u,4u,float> Qk; + float t = tr(Spx); - trace::exiting("registration::impl::icp_"); + Qk(0,0) = t; + for (int i = 1; i < 4; ++i) + { + Qk(i,0) = A[i - 1]; + Qk(0,i) = A[i - 1]; + for (int j = 1; j < 4; ++j) + if (i == j) + Qk(i,j) = 2 * Spx(i - 1,i - 1) - t; } - } // end of namespace mln::registration::impl + Qk(1,2) = Spx(0,1) + Spx(1,0); + Qk(2,1) = Spx(0,1) + Spx(1,0); + + Qk(1,3) = Spx(0,2) + Spx(2,0); + Qk(3,1) = Spx(0,2) + Spx(2,0); + + Qk(2,3) = Spx(1,2) + Spx(2,1); + Qk(3,2) = Spx(1,2) + Spx(2,1); + + return math::jacobi(Qk); + } + - template <typename P, typename M> + // Compute mu_Yk, mass center of Yk. + template <typename P, typename F> inline - composed<rotation<P::dim, float>, translation<P::dim, float> > - icp(const p_array& c, const M& map, - const float epsilon = 1e-3) + vec3d_f + get_mu_yk(const p_array& P_, + const F& closest_point, + const algebra::quat& qR, + const vec3d_f& qT, + float& e_k) { - composed<rotation<P::dim, float>, translation<P::dim, float> > qk; + accu::rms<vec3d_f,float> e_k_accu; + accu::center<P,vec3d_f> mu_yk; - std::cout << " icp " << std::endl; + mln_piter(p_array) p(P_); + for_all(p) + { + // yk_i - pk_i + vec3d_f Pk_i = qR.rotate(p.to_vec()) + qT; + vec3d_f Yk_i = closest_point(Pk_i).to_vec(); + mu_yk.take(Yk_i); + e_k_accu.take(Yk_i - Pk_i); + } - impl::icp_(c, c.nsites(), map, qk, epsilon); - return qk; + e_k = e_k_accu.to_result(); + return mu_yk.to_result(); } - template <typename P, typename M, typename T> + + + /// Base version of the ICP algorithm. It is called in other variants. + template <typename P, typename F> inline - void - icp_subset(const p_array& c, - const unsigned int c_length, - const M& map, - T& qk, - float epsilon = 1e-3) + std::pair<algebra::quat,mln_vec(P)> + icp(const p_array& P_, + const p_array& X, + const F& closest_point, + const algebra::quat& initial_rot, + const mln_vec(P)& initial_translation) { - impl::icp_(c, c_length, map, qk, epsilon); + trace::entering("registration::icp"); + + (void) X; + mln_precondition(P::dim == 3); + mln_precondition(!P_.is_empty()); + mln_precondition(!X.is_empty()); + mln_precondition(!initial_rot.is_null()); + + typedef p_array cloud_t; + + vec3d_f mu_P = set::compute(accu::center<P,vec3d_f>(), P_); + + vec3d_f qT_old, qT = initial_translation; + algebra::quat qR_old, qR = initial_rot; + float e_k, e_k_old = mln_max(float); + float d_k, d_k_old = mln_max(float); + unsigned k = 0; + +# ifndef NDEBUG + box3d box = geom::bbox(X); + //FIXME: too large? + box.enlarge(1, box.nrows() / 2); + box.enlarge(2, box.ncols() / 2); + image3d<value::rgb8> debug(box); + data::fill(debug, literal::black); + data::fill((debug | X).rw(), literal::white); +# endif + + do + { + qT_old = qT; + qR_old = qR; + + /// Compute transformation + /// + // mu_Yk - Pk's mass center. + // + Compute error ek = d(Pk,Yk) + vec3d_f mu_Yk = get_mu_yk(P_, closest_point, qR_old, qT_old, e_k); + + // quaternion qR - rotation + qR = get_rot(P_, mu_P, mu_Yk, closest_point, qR_old, qT_old); + vec3d_f tmp = qR.v(); + + // vector qT - translation + qT = mu_Yk - qR.rotate(mu_P); + /// + /// End of "compute transformation" + + // Distance dk = d(Pk+1, Yk) + d_k = compute_d_k(P_, closest_point, qR, qR_old, qT, qT_old); + + +#ifndef NDEBUG + image3d<value::rgb8> tmp_ = duplicate(debug); + mln_piter(p_array) p_dbg(P_); + for_all(p_dbg) + tmp_(qR_old.rotate(p_dbg.to_vec()) + qT_old) = literal::green; + std::ostringstream ss; + ss << "tmp_0"; + if (k < 10) + ss << "0"; + ss << k << ".ppm"; + io::ppm::save(mln::slice(tmp_,0), ss.str()); +#endif + + std::cout << "e_" << k << "=" << e_k << std::endl; + std::cout << "d_" << k << "=" << d_k << std::endl; + + // Check distance and error according to the related paper. + // Disabled because of the following 'if' +// mln_assertion(0 <= d_k); +// mln_assertion(d_k <= e_k); + +// mln_assertion(e_k <= d_k_old); +// mln_assertion(d_k_old <= e_k_old); + + // During the first runs, d_k may be higher than e_k. + // Hence, there we test k > 3. + if (k > 3 && (d_k > e_k || d_k > d_k_old || e_k > e_k_old)) + { + qR = qR_old; + qT = qT_old; + break; + } + + // Backing up results. + d_k_old = d_k; + e_k_old = e_k; + + ++k; + + } while ((k < 3) + || norm::l2((qT - qT_old)) + norm::l2((qR - qR_old).to_vec()) > 1e-3); + + trace::exiting("registration::icp"); + return std::make_pair(qR, qT); } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::registration } // end of namespace mln - - #endif // ! MLN_REGISTRATION_ICP_HH diff --git a/milena/mln/registration/icp2.hh b/milena/mln/registration/icp2.hh deleted file mode 100644 index de763ea..0000000 --- a/milena/mln/registration/icp2.hh +++ /dev/null @@ -1,693 +0,0 @@ -// Copyright (C) 2009 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_REGISTRATION_ICP_HH -# define MLN_REGISTRATION_ICP_HH - -/// \file mln/registration/icp.hh -/// -/// Register an image over an another using the ICP algorithm. -/// -/// \todo encode distances on 12 bits. - -# include <cmath> -# include <algorithm> - -# include <mln/core/alias/vec3d.hh> -# include <mln/math/jacobi.hh> -# include <mln/fun/x2x/all.hh> -# include <mln/fun/x2v/all.hh> -# include <mln/convert/to.hh> -# include <mln/accu/compute.hh> -# include <mln/accu/center.hh> -# include <mln/accu/rms.hh> -# include <mln/trait/image_from_grid.hh> -# include <mln/set/compute.hh> - -//Should be removed when closest_point functors are moved. -# include <mln/core/alias/neighb3d.hh> -# include <mln/transform/internal/closest_point_functor.hh> -# include <mln/canvas/distance_geodesic.hh> -# include <mln/pw/all.hh> - -# include <mln/io/ppm/save.hh> -# include <mln/io/pbm/save.hh> -# include <mln/debug/colorize.hh> - -# include <mln/literal/black.hh> -# include <mln/literal/white.hh> -# include <mln/literal/colors.hh> - -# include <mln/core/image/slice_image.hh> - -#include <mln/core/image/tr_image.hh> -#include <mln/core/image/extension_fun.hh> - -#include <mln/accu/histo.hh> -#include <mln/accu/sum.hh> -#include <mln/debug/histo.hh> - -# include <mln/util/timer.hh> - -namespace mln -{ - - - namespace registration - { - - //FIXME: used for debug purpose. Should be removed later. - - using namespace fun::x2x; - - /*! Register point in \p c using a function of closest points - * \p closest_point. - * - * \param[in] P_ The cloud of points. - * \param[in] X the reference surface. - * \param[in] closest_point The function of closest points. - * \param[out] qk The rigid transformation obtained. - * - * \return the rigid transformation which may be use later to create - * a registered image. - * - * WARNING: the function \p closest_point *MUST* take float/double - * vector as arguments. Otherwise the resulting transformation may be - * wrong due to the truncation of the vector coordinate values. - * - * \pre \p P_ and \p X must not be empty. - * - * Reference article: "A Method for Registration of 3-D Shapes", Paul J. - * Besl and Neil D. McKay, IEEE, 2, February 1992. - * - */ - template <typename P, typename F> - std::pair<algebra::quat,mln_vec(P)> - icp(const p_array& P_, - const p_array& X, - const F& closest_point, - const algebra::quat& initial_rot, - const mln_vec(P)& initial_translation); - - - template <typename P, typename F> - composed< translation<P::dim,float>,rotation<P::dim,float> > - icp(const p_array& P_, - const p_array& X, - const F& closest_point); - - -# ifndef MLN_INCLUDE_ONLY - - - /// Closest point functor based on map distance. - template <typename P> - class closest_point_with_map - { - typedef mln_image_from_grid(mln_grid(P), P) I; - typedef mln_ch_value(I, P) cp_ima_t; - typedef mln_ch_value(I,value::int_u16) dmap_t; - - public: - - closest_point_with_map(const p_array& X) - { - box3d box = geom::bbox(X); - box.enlarge(1, box.nrows()); - box.enlarge(2, box.ncols()); - - mln_postcondition(box.is_valid()); - - std::cout << "Map image defined on " << box << std::endl; - - init(X, box); - } - - closest_point_with_map(const p_array& X, const box& box) - { - init(X, box); - } - - void init(const p_array& X, const box& box) - { - typedef mln_ch_value(I, bool) model_t; - model_t model(box); - data::fill(model, false); - data::fill((model | X).rw(), true); - - transform::internal::closest_point_functor<model_t> f; - util::timer t; - t.start(); - dmap_X_ = canvas::distance_geodesic(model, c6(), - mln_max(value::int_u16), - f); - std::cout << "canvas::distance_geodesic - " << t << "s" << std::endl; - - initialize(this->cp_ima_, f.cp_ima); - this->cp_ima_ = f.cp_ima; - - mln_postcondition(cp_ima_.is_valid()); - mln_postcondition(cp_ima_.domain().is_valid()); - std::cout << "pmin = " << cp_ima_.domain().pmin() << std::endl;; - std::cout << "pmax = " << cp_ima_.domain().pmax() << std::endl;; - -#ifndef NDEBUG -// mln_ch_value(I, bool) debug2(box); -// data::fill(debug2, false); -// mln_ch_value(I, value::rgb8) debug(box); -// mln_piter(p_array) p(X); -// for_all(p) -// { -// debug(p) = debug::internal::random_color(value::rgb8()); -// debug2(p) = true; -// } -// io::pbm::save(slice(debug2,0), "debug2-a.ppm"); -// -// mln_piter(I) pi(cp_ima_.domain()); -// for_all(pi) -// { -// debug(pi) = debug(cp_ima_(pi)); -// debug2(pi) = debug2(cp_ima_(pi)); -// } -// -// io::pbm::save(slice(debug2,0), "debug2-b.ppm"); -// io::ppm::save(slice(debug,0), "debug.ppm"); -// std::cout << "map saved" << std::endl; -#endif - } - - mln_site(I) operator()(const mln_site(I)& p) const - { - return cp_ima_(p); - } - - // Distance map - dmap_t dmap_X_; - // Closest point image. - cp_ima_t cp_ima_; - private: - }; - - - /// Closest point functor based on map distance. - template <typename P> - class closest_point_basic - { - typedef mln_image_from_grid(mln_grid(P), P) I; - typedef p_array X_t; - - public: - - closest_point_basic(const p_array& X) - : X_(X) - { - } - - mln_site(I) operator()(const vec3d_f& v) const - { - vec3d_f best_x = convert::to<vec3d_f>(X_[0].to_vec()); - - float best_d = norm::l2_distance(v, best_x); - mln_piter(X_t) X_i(X_); - for_all(X_i) - { - float d = norm::l2_distance(v, convert::to<vec3d_f>(X_i)); - if (d < best_d) - { - best_d = d; - best_x = X_i.to_vec(); - } - } - return best_x; - } - - private: - const p_array& X_; - }; - - - template <typename P> - void - draw_last_run(const box3d& box, const p_array& kept, - const p_array& removed, const p_array& X, - const algebra::quat& qR, const vec3d_f qT) - { - typedef image3d<value::rgb8> result_t; - result_t result(box); - typedef extension_fun<result_t,pw::cst_<mln_value(result_t)> > ext_result_t; - ext_result_t ext_result(result, pw::cst(value::rgb8(0,0,0))); - - data::fill(ext_result, literal::black); - data::fill((ext_result | X).rw(), literal::white); - - mln_piter(p_array) p(kept); - for_all(p) - ext_result(qR.rotate(p.to_vec()) + qT) = literal::green; - - mln_piter(p_array) p2(removed); - for_all(p2) - ext_result(qR.rotate(p2.to_vec()) + qT) = literal::red; - - io::ppm::save(slice(ext_result,0), "registered-2.ppm"); - } - - - - template <typename P, typename F> - float compute_standard_deviation(const p_array& P_, - const std::pair<algebra::quat,mln_vec(P)>& pair, - const F& closest_point) - { - accu::rms<vec3d_f,float> e_k_accu; - - // Standard Deviation - float sd; - mln_piter(p_array) p(P_); - for_all(p) - { - vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; - vec3d_f Yk_i = closest_point(Pk_i).to_vec(); - // yk_i - pk_i - e_k_accu.take(Yk_i - Pk_i); - } - - float d = e_k_accu.to_result(); - sd = std::sqrt(e_k_accu.hook_value_() / P_.nsites() - d * d); - return sd; - } - - - template <typename P, typename F> - void - remove_too_far_sites_debug(image3d<value::rgb8>& out, const p_array& P_, - const F& closest_point, - const std::pair<algebra::quat,mln_vec(P)>& pair, - const p_array& X, - unsigned r, int d_min, int d_max, unsigned prefix) - { - unsigned removed = 0; - accu::histo<value::int_u8> h; - mln_piter(p_array) p(P_); - data::fill(out, literal::black); - data::fill((out | X).rw(), literal::white); - - for_all(p) - { - vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; - vec3d_f Yk_i = closest_point(Pk_i).to_vec(); - - int d_i = closest_point.dmap_X_(Pk_i); - if (d_i >= d_min && d_i <= d_max) - out(Pk_i) = literal::green; - else - { - ++removed; - out(Pk_i) = literal::red; - } - } - -#ifndef NDEBUG - std::ostringstream ss1; - ss1 << "histo_" << prefix << r << ".dat"; - std::cout << h << std::endl; - - std::ostringstream ss2; - ss2 << "out_" << prefix << r << ".ppm"; - io::ppm::save(mln::slice(out,0), ss2.str()); -#endif - std::cout << "Points removed with the whole set and current d_min/d_max: " << removed << std::endl; - - } - - - template <typename P, typename F> - void - compute_distance_criteria(const p_array& P_, - const F& closest_point, - const std::pair<algebra::quat,mln_vec(P)>& pair, - unsigned r, int& d_min, int& d_max) - { - mln_piter(p_array) p(P_); - accu::histo<value::int_u8> h; - - float sd; - { - accu::sum<float> s, s2; - for_all(p) - { - vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; - unsigned d_i = closest_point.dmap_X_(Pk_i); - h.take(d_i); - s.take(d_i); - s2.take(d_i * d_i); - } - float mean = s / P_.nsites(); - sd = std::sqrt(s2 / P_.nsites() - mean * mean); - d_min = int(mean - sd); - d_max = int(mean + sd); - } - - std::cout << "Standard deviation = " << sd << std::endl; - std::ostringstream ss1; - ss1 << "histo_" << r << ".dat"; - std::cout << h << std::endl; - std::cout << "d thresholds = " << d_min << ' ' << d_max << std::endl; - } - - template <typename P, typename F> - p_array - remove_too_far_sites(image3d<value::rgb8>& out, const p_array& P_, - const F& closest_point, - const std::pair<algebra::quat,mln_vec(P)>& pair, - const p_array& X, p_array& removed_set, - unsigned r, int d_min, int d_max, - const std::string& method) - { - p_array tmp; - unsigned removed = 0; - -# ifndef NDEBUG - data::fill(out, literal::black); - data::fill((out | X).rw(), literal::white); -# endif // ! NDEBUG - - mln_piter(p_array) p(P_); - for_all(p) - { - vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; - vec3d_f Yk_i = closest_point(Pk_i).to_vec(); - - int d_i = closest_point.dmap_X_(Pk_i); - if (d_i >= d_min && d_i <= d_max) - { - tmp.append(p); - out(Pk_i) = literal::green; - } - else - { - ++removed; - removed_set.append(p); - out(Pk_i) = literal::red; - } - } - - -# ifndef NDEBUG - std::ostringstream ss2; - ss2 << method << "_" << r << "_removed_sites" << ".ppm"; - io::ppm::save(mln::slice(out,0), ss2.str()); - - std::cout << "Points removed: " << removed << std::endl; -# endif // ! NDEBUG - // They are used for debug purpose only. - // When NDEBUG is set, they are unused. - (void) X; - (void) r; - (void) method; - - return tmp; - } - - template <typename P> - void - display_sites_used_in_icp(image3d<value::rgb8>& out, const p_array& P_sub, - const p_array& P_, const p_array& X, - unsigned r, const std::string& prefix, - const std::pair<algebra::quat,mln_vec(P)>& pair, - const std::string& period, const value::rgb8& c) - { -# ifndef NDEBUG - data::fill(out, literal::black); - data::fill((out | X).rw(), literal::white); -# endif // !NDEBUG - - mln_piter(p_array) p1(P_); - for_all(p1) - { - vec3d_f Pk_i = pair.first.rotate(p1.to_vec()) + pair.second; - out(Pk_i) = literal::red; - } - - mln_piter(p_array) p2(P_sub); - for_all(p2) - { - vec3d_f Pk_i = pair.first.rotate(p2.to_vec()) + pair.second; - out(Pk_i) = c; - } - -# ifndef NDEBUG - std::ostringstream ss; - ss << prefix << "_" << r << "_" << period << ".ppm"; - - io::ppm::save(slice(out,0), ss.str()); -# endif // ! NDEBUG - } - - - template <typename P, typename F> - inline - float - compute_d_k(const p_array& P_, - const F& closest_point, - const algebra::quat& qR, - const algebra::quat& qR_old, - const vec3d_f& qT, - const vec3d_f& qT_old) - { - accu::rms<vec3d_f, float> accu; - mln_piter(p_array) p(P_); - for_all(p) - { - // yk_i - pk+1_i - vec3d_f Pk_i = qR_old.rotate(p.to_vec()) + qT_old; - vec3d_f Pk_1_i = qR.rotate(p.to_vec()) + qT; - accu.take(closest_point(Pk_i).to_vec() - Pk_1_i); - } - return accu.to_result(); - } - - - /// FIXME: work only for 3d images. - template <typename P, typename F> - algebra::quat - get_rot(const p_array& P_, - const vec3d_f& mu_P, - const vec3d_f& mu_Yk, - const F& closest_point, - const algebra::quat& qR, - const vec3d_f& qT) - { - /// Spx: cross-covariance matrix. - algebra::mat<3u,3u,float> Spx; - mln_piter(p_array) p(P_); - - // FIXME: could we use an accu? - for_all(p) - { - vec3d_f P_i = p.to_vec(); - vec3d_f Pk_i = qR.rotate(p.to_vec()) + qT; - vec3d_f Yk_i = closest_point(Pk_i).to_vec(); - Spx += make::mat(P_i - mu_P) * trans(make::mat(Yk_i - mu_Yk)); - } - Spx /= P_.nsites(); - - vec3d_f A; - A[0] = Spx(1,2) - Spx(2,1); - A[1] = Spx(2,0) - Spx(0,2); - A[2] = Spx(0,1) - Spx(1,0); - - algebra::mat<4u,4u,float> Qk; - float t = tr(Spx); - - Qk(0,0) = t; - for (int i = 1; i < 4; ++i) - { - Qk(i,0) = A[i - 1]; - Qk(0,i) = A[i - 1]; - for (int j = 1; j < 4; ++j) - if (i == j) - Qk(i,j) = 2 * Spx(i - 1,i - 1) - t; - } - - Qk(1,2) = Spx(0,1) + Spx(1,0); - Qk(2,1) = Spx(0,1) + Spx(1,0); - - Qk(1,3) = Spx(0,2) + Spx(2,0); - Qk(3,1) = Spx(0,2) + Spx(2,0); - - Qk(2,3) = Spx(1,2) + Spx(2,1); - Qk(3,2) = Spx(1,2) + Spx(2,1); - - return math::jacobi(Qk); - } - - - // Compute mu_Yk, mass center of Yk. - template <typename P, typename F> - inline - vec3d_f - get_mu_yk(const p_array& P_, - const F& closest_point, - const algebra::quat& qR, - const vec3d_f& qT, - float& e_k) - { - accu::rms<vec3d_f,float> e_k_accu; - accu::center<P,vec3d_f> mu_yk; - - mln_piter(p_array) p(P_); - for_all(p) - { - // yk_i - pk_i - vec3d_f Pk_i = qR.rotate(p.to_vec()) + qT; - vec3d_f Yk_i = closest_point(Pk_i).to_vec(); - mu_yk.take(Yk_i); - e_k_accu.take(Yk_i - Pk_i); - } - - e_k = e_k_accu.to_result(); - return mu_yk.to_result(); - } - - - - /// Base version of the ICP algorithm. It is called in other variants. - template <typename P, typename F> - inline - std::pair<algebra::quat,mln_vec(P)> - icp(const p_array& P_, - const p_array& X, - const F& closest_point, - const algebra::quat& initial_rot, - const mln_vec(P)& initial_translation) - { - trace::entering("registration::icp"); - - (void) X; - mln_precondition(P::dim == 3); - mln_precondition(!P_.is_empty()); - mln_precondition(!X.is_empty()); - mln_precondition(!initial_rot.is_null()); - - typedef p_array cloud_t; - - vec3d_f mu_P = set::compute(accu::center<P,vec3d_f>(), P_); - - vec3d_f qT_old, qT = initial_translation; - algebra::quat qR_old, qR = initial_rot; - float e_k, e_k_old = mln_max(float); - float d_k, d_k_old = mln_max(float); - unsigned k = 0; - -# ifndef NDEBUG - box3d box = geom::bbox(X); - //FIXME: too large? - box.enlarge(1, box.nrows() / 2); - box.enlarge(2, box.ncols() / 2); - image3d<value::rgb8> debug(box); - data::fill(debug, literal::black); - data::fill((debug | X).rw(), literal::white); -# endif - - do - { - qT_old = qT; - qR_old = qR; - - /// Compute transformation - /// - // mu_Yk - Pk's mass center. - // + Compute error ek = d(Pk,Yk) - vec3d_f mu_Yk = get_mu_yk(P_, closest_point, qR_old, qT_old, e_k); - - // quaternion qR - rotation - qR = get_rot(P_, mu_P, mu_Yk, closest_point, qR_old, qT_old); - vec3d_f tmp = qR.v(); - - // vector qT - translation - qT = mu_Yk - qR.rotate(mu_P); - /// - /// End of "compute transformation" - - // Distance dk = d(Pk+1, Yk) - d_k = compute_d_k(P_, closest_point, qR, qR_old, qT, qT_old); - - -#ifndef NDEBUG - image3d<value::rgb8> tmp_ = duplicate(debug); - mln_piter(p_array) p_dbg(P_); - for_all(p_dbg) - tmp_(qR_old.rotate(p_dbg.to_vec()) + qT_old) = literal::green; - std::ostringstream ss; - ss << "tmp_0"; - if (k < 10) - ss << "0"; - ss << k << ".ppm"; - io::ppm::save(mln::slice(tmp_,0), ss.str()); -#endif - - std::cout << "e_" << k << "=" << e_k << std::endl; - std::cout << "d_" << k << "=" << d_k << std::endl; - - // Check distance and error according to the related paper. - // Disabled because of the following 'if' -// mln_assertion(0 <= d_k); -// mln_assertion(d_k <= e_k); - -// mln_assertion(e_k <= d_k_old); -// mln_assertion(d_k_old <= e_k_old); - - // During the first runs, d_k may be higher than e_k. - // Hence, there we test k > 3. - if (k > 3 && (d_k > e_k || d_k > d_k_old || e_k > e_k_old)) - { - qR = qR_old; - qT = qT_old; - break; - } - - // Backing up results. - d_k_old = d_k; - e_k_old = e_k; - - ++k; - - } while ((k < 3) - || norm::l2((qT - qT_old)) + norm::l2((qR - qR_old).to_vec()) > 1e-3); - - trace::exiting("registration::icp"); - return std::make_pair(qR, qT); - } - - -# endif // ! MLN_INCLUDE_ONLY - - } // end of namespace mln::registration - -} // end of namespace mln - -#endif // ! MLN_REGISTRATION_ICP_HH diff --git a/milena/mln/registration/registration.hh b/milena/mln/registration/registration.hh index f96787a..21d6ea5 100644 --- a/milena/mln/registration/registration.hh +++ b/milena/mln/registration/registration.hh @@ -34,11 +34,14 @@ /// \sa registration::icp # include <mln/core/image/image3d.hh> -# include <mln/registration/icp2.hh> +# include <mln/registration/icp.hh> # include <mln/fun/x2x/all.hh> # include <mln/fun/x2p/closest_point.hh> # include <mln/convert/to_p_array.hh> +//FIXME: to be removed. +# include <mln/util/timer.hh> + namespace mln { @@ -331,7 +334,7 @@ namespace mln { trace::entering("registration::registration1"); - registration_tests(cloud, surface); + internal::registration_tests(cloud, surface); composed< translation<P::dim,float>, rotation<P::dim,float> > qk = impl::registration1(cloud, surface); @@ -350,7 +353,7 @@ namespace mln { trace::entering("registration::registration2"); - registration_tests(cloud, surface); + internal::registration_tests(cloud, surface); composed< translation<P::dim,float>, rotation<P::dim,float> > qk = impl::registration2(cloud, surface); diff --git a/milena/mln/util/graph.hh b/milena/mln/util/graph.hh index f95e1e8..10e900f 100644 --- a/milena/mln/util/graph.hh +++ b/milena/mln/util/graph.hh @@ -405,7 +405,9 @@ namespace mln unsigned id = data_->edges_.size() - 1; // Update the set of edges. +# ifndef NDEBUG data_->edges_set_.insert(edge); +# endif // ! NDEBUG data_->vertices_[edge.first()].push_back(id); data_->vertices_[edge.second()].push_back(id); diff --git a/milena/mln/value/rgb.hh b/milena/mln/value/rgb.hh index 9c64936..839d89c 100644 --- a/milena/mln/value/rgb.hh +++ b/milena/mln/value/rgb.hh @@ -51,6 +51,11 @@ namespace mln { + // Forward declaration. + namespace value { template <unsigned n> struct rgb; } + + + namespace fun { -- 1.5.6.5

15 years, 10 months

3389: Add support for 3D images in fun::x2v::bilinear.

by Guillaume Lazzara

* mln/fun/x2v/bilinear.hh: Add a new specialization for 3d Images. --- milena/ChangeLog | 6 ++++ milena/mln/fun/x2v/bilinear.hh | 65 ++++++++++++++++++++++++++++++++++++---- 2 files changed, 65 insertions(+), 6 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 93b5737..1c446cd 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,11 @@ 2009-02-18 Guillaume Lazzara <z(a)lrde.epita.fr> + Add support for 3D images in fun::x2v::bilinear. + + * mln/fun/x2v/bilinear.hh: Add a new specialization for 3d Images. + +2009-02-18 Guillaume Lazzara <z(a)lrde.epita.fr> + Fix fun::x2v::rotation. * mln/fun/x2x/rotation.hh: diff --git a/milena/mln/fun/x2v/bilinear.hh b/milena/mln/fun/x2v/bilinear.hh index b0add18..e264462 100644 --- a/milena/mln/fun/x2v/bilinear.hh +++ b/milena/mln/fun/x2v/bilinear.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// Copyright (C) 2008, 2009 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 @@ -59,9 +60,15 @@ namespace mln bilinear(const I& ima); - template <unsigned n, typename T> + /// Bilinear filtering on 2d images. + template <typename T> mln_value(I) - operator()(const algebra::vec<n,T>& v) const; + operator()(const algebra::vec<2,T>& v) const; + + /// Bilinear filtering on 3d images. Work on slices. + template <typename T> + mln_value(I) + operator()(const algebra::vec<3,T>& v) const; const I& ima; }; @@ -72,13 +79,13 @@ namespace mln template <typename I> bilinear::bilinear(const I& ima) : ima(ima) { - mlc_bool(I::psite::dim == 2)::check(); + //mlc_bool(I::psite::dim == 2)::check(); } template <typename I> - template <unsigned n, typename T> + template <typename T> mln_value(I) - bilinear::operator()(const algebra::vec<n,T>& v) const + bilinear::operator()(const algebra::vec<2,T>& v) const { typedef mln_sum(mln_value(I)) vsum; @@ -121,6 +128,52 @@ namespace mln } + template <typename I> + template <typename T> + mln_value(I) + bilinear::operator()(const algebra::vec<3,T>& v) const + { + typedef mln_sum(mln_value(I)) vsum; + + // q12----r2----q22 + // | | | + // | x | + // | | | + // q11----r1----q21 + + double x = v[1]; + double y = v[2]; + + double x1 = std::floor(x); + double x2 = std::floor(x) + 1; + double y1 = std::floor(y); + double y2 = std::floor(y) + 1; + def::coord z = math::round<float>()(v[0]); + + //Following access are supposed valid. + vsum q11 = ima(point3d(z, static_cast<unsigned>(x1), static_cast<unsigned>(y1))); + vsum q12 = ima(point3d(z, static_cast<unsigned>(x1), static_cast<unsigned>(y2))); + vsum q21 = ima(point3d(z, static_cast<unsigned>(x2), static_cast<unsigned>(y1))); + vsum q22 = ima(point3d(z, static_cast<unsigned>(x2), static_cast<unsigned>(y2))); + + double x2_x1 = x2 - x1; + double y2_y1 = y2 - y1; + + // linear interpolation #1 + vsum img_r1 = q11 * (x2 - x) / (x2_x1) + + q21 * (x - x1) / (x2_x1); + + // linear interpolation #2 + vsum img_r2 = q12 * (x2 - x) / (x2_x1) + q22 * (x - x1) / (x2_x1); + + // interpolating in y direction + vsum res = (img_r1 * (y2 - y) / (y2_y1) + + img_r2 * (y - y1) / (y2_y1)); + + return convert::to<mln_value(I)>(res); + } + + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::fun::x2v -- 1.5.6.5

15 years, 10 months

3388: Fix fun::x2v::rotation.

by Guillaume Lazzara

* mln/fun/x2x/rotation.hh: (get_rot_h_mat 3d): fix matrix construction. (rotation(algebra::quat)): update attributes in order to produce a valid inverted rotation. * tests/fun/x2x/rotation.cc: improve test. --- milena/ChangeLog | 19 +++++++++++++ milena/mln/fun/x2x/rotation.hh | 36 ++++++++++++++++++------- milena/tests/fun/x2x/rotation.cc | 55 ++++++++++++++----------------------- 3 files changed, 66 insertions(+), 44 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 62c0696..93b5737 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,24 @@ 2009-02-18 Guillaume Lazzara <z(a)lrde.epita.fr> + Fix fun::x2v::rotation. + + * mln/fun/x2x/rotation.hh: + (get_rot_h_mat 3d): fix matrix construction. + (rotation(algebra::quat)): update attributes in order to produce a + valid inverted rotation. + + * tests/fun/x2x/rotation.cc: improve test. + +2009-02-18 Guillaume Lazzara <z(a)lrde.epita.fr> + + Add fun::x2v::trilinear. + + * mln/fun/x2v/trilinear.hh: new file. New interpolation algorithm. + + * mln/fun/x2v/all.hh: include new file. + +2009-02-18 Guillaume Lazzara <z(a)lrde.epita.fr> + Add aliases for 2D algebra::vec. * mln/core/alias/vec2d.hh: new. new aliases for float and double 2d diff --git a/milena/mln/fun/x2x/rotation.hh b/milena/mln/fun/x2x/rotation.hh index dce37ed..9628c26 100644 --- a/milena/mln/fun/x2x/rotation.hh +++ b/milena/mln/fun/x2x/rotation.hh @@ -32,6 +32,9 @@ /// \file mln/fun/x2x/rotation.hh /// /// Define a rotation function. +/// +/// \todo store the quaternion instead of (axis, alpha) +/// => better precision while composing two rotation matrices. # include <cmath> # include <mln/core/concept/function.hh> @@ -70,6 +73,10 @@ namespace mln //FIXME: cannot use '!=' operator. mln_precondition(!(axis_ == vec_t(literal::zero))); + /// Be sure that axis is normalized. + algebra::vec<3,C> axis = axis_; + axis.normalize(); + const float cos_a = cos(alpha_); const float sin_a = sin(alpha_); const float u = axis_[0]; @@ -78,23 +85,22 @@ namespace mln const float u2 = u * u; const float v2 = v * v; const float w2 = w * w; - const float uvw2 = u2 + v2 + w2; algebra::h_mat<3, C> m_; - m_(0,0) = (u2 + (v2 + w2) * cos_a) / uvw2; - m_(0,1) = (u*v * (1 - cos_a) - u * std::sqrt(uvw2) * sin_a) / uvw2; - m_(0,2) = (u*w * (1 - cos_a) + v * std::sqrt(uvw2) * sin_a) / uvw2; + m_(0,0) = u2 + (1 - u2) * cos_a; + m_(0,1) = u*v * (1 - cos_a) - w * sin_a; + m_(0,2) = u*w * (1 - cos_a) + v * sin_a; m_(0,3) = 0; - m_(1,0) = (u*v * (1 - cos_a) + w * std::sqrt(uvw2) * sin_a) / uvw2; - m_(1,1) = (v2 + (u2 + w2) * cos_a) / uvw2; - m_(1,2) = (v*w * (1 - cos_a) - u * std::sqrt(uvw2) * sin_a) / uvw2; + m_(1,0) = u*v * (1 - cos_a) + w * sin_a; + m_(1,1) = v2 + (1 - v2) * cos_a; + m_(1,2) = v * w * (1 - cos_a) - u * sin_a; m_(1,3) = 0; - m_(2,0) = (u*w * (1 - cos_a) - v * std::sqrt(uvw2) * sin_a) / uvw2; - m_(2,1) = (v*w * (1 - cos_a) + u * std::sqrt(uvw2) * sin_a) / uvw2; - m_(2,2) = (u2 + (u2 + v2) * cos_a) / uvw2; + m_(2,0) = u * w * (1 - cos_a) - v * sin_a; + m_(2,1) = u * w * (1 - cos_a) + u * sin_a; + m_(2,2) = w2 + (1 - w2) * cos_a; m_(2,3) = 0; m_(3,0) = 0; @@ -189,19 +195,29 @@ namespace mln rotation<n,C>::rotation(const algebra::quat& q) { mln_precondition(q.is_unit()); + // FIXME: Should also work for 2d. mln_precondition(n == 3); + float w = q.to_vec()[0], x = q.to_vec()[1], x2 = 2*x*x, xw = 2*x*w, y = q.to_vec()[2], y2 = 2*y*y, xy = 2*x*y, yw = 2*y*w, z = q.to_vec()[3], z2 = 2*z*z, xz = 2*x*z, yz = 2*y*z, zw = 2*z*w; + float t[9] = {1.f - y2 - z2, xy - zw, xz + yw, xy + zw, 1.f - x2 - z2, yz - xw, xz - yw, yz + xw, 1.f - x2 - y2}; this->m_ = mln::make::h_mat(t); mln_assertion(check_rotation(q)); + + /// Update attributes + + alpha_ = acos(w) * 2; + axis_[0] = x; + axis_[1] = y; + axis_[2] = z; } diff --git a/milena/tests/fun/x2x/rotation.cc b/milena/tests/fun/x2x/rotation.cc index b96c551..fc36a17 100644 --- a/milena/tests/fun/x2x/rotation.cc +++ b/milena/tests/fun/x2x/rotation.cc @@ -33,14 +33,9 @@ #include <iostream> #include <mln/fun/x2x/rotation.hh> -#include <mln/core/image/image2d.hh> -#include <mln/value/int_u8.hh> -#include <mln/io/pgm/load.hh> -#include <mln/io/pgm/save.hh> -#include <mln/core/image/interpolated.hh> +#include <mln/core/alias/vec2d.hh> +#include <mln/core/alias/vec3d.hh> #include <mln/make/vec.hh> -#include <mln/fun/x2v/bilinear.hh> -#include <mln/extension/adjust.hh> #include "tests/data.hh" @@ -48,36 +43,28 @@ int main() { using namespace mln; - using value::int_u8; - algebra::vec<2,float> axis; - axis[0] = 0; - axis[1] = 1; - typedef image2d<int_u8> ima_t; - ima_t lena; - io::pgm::load(lena, MLN_IMG_DIR "/lena.pgm"); + { + algebra::vec<2,float> axis; + axis[0] = 0; + axis[1] = 1; + fun::x2x::rotation<2,float> rot(0.1f, axis); - ima_t out; - initialize(out, lena); + vec2d_f v = make::vec(4, 4); + vec2d_f res = rot(v); + mln_assertion(rot.inv()(res) == v); + } - interpolated<ima_t, fun::x2v::bilinear> inter(lena); + { + algebra::vec<3,float> axis; + axis[0] = 1; + axis[1] = 0; + axis[2] = 0; + fun::x2x::rotation<3,float> rot(0.1f, axis); - fun::x2x::rotation<2,float> rot1(0.1, axis); + vec3d_f v = make::vec(4, 1, 2); + vec3d_f res = rot(v); + mln_assertion(rot.inv()(res) == v); + } - mln_piter_(ima_t) p(out.domain()); - for_all(p) - { - algebra::vec<2,float> v = rot1.inv()(p.to_site().to_vec()); - if (inter.domain().has(v)) - out(p) = inter(v); - else - out(p) = 255; - } - io::pgm::save(out, "out.pgm"); - - fun::x2x::rotation<2,float> rot2(3.14116, axis); - mln_assertion(fabs(rot2(make::vec(0.0, 1.0))[0] - - make::vec(0.0, -1.0)[0]) <= 0.125); - mln_assertion(fabs(rot2(make::vec(0.0, 1.0))[1] - - make::vec(0.0, -1.0)[1]) <= 0.125); } -- 1.5.6.5

15 years, 10 months

3387: Add fun::x2v::trilinear.

by Guillaume Lazzara

* mln/fun/x2v/trilinear.hh: new file. New interpolation algorithm. * mln/fun/x2v/all.hh: include new file. --- milena/mln/fun/x2v/all.hh | 12 ++-- milena/mln/fun/x2v/trilinear.hh | 139 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 146 insertions(+), 5 deletions(-) create mode 100644 milena/mln/fun/x2v/trilinear.hh diff --git a/milena/mln/fun/x2v/all.hh b/milena/mln/fun/x2v/all.hh index 91c4562..d8ba3ae 100644 --- a/milena/mln/fun/x2v/all.hh +++ b/milena/mln/fun/x2v/all.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory +// Copyright (C) 2008, 2009 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 @@ -28,10 +29,10 @@ #ifndef MLN_FUN_X2V_ALL_HH # define MLN_FUN_X2V_ALL_HH -/*! \file mln/fun/x2v/all.hh - * - * \brief File that includes all functions from vector to value. - */ +/// \file mln/fun/x2v/all.hh +/// +/// File that includes all functions from vector to value. + namespace mln @@ -57,5 +58,6 @@ namespace mln # include <mln/fun/x2v/linear.hh> # include <mln/fun/x2v/bilinear.hh> # include <mln/fun/x2v/nneighbor.hh> +# include <mln/fun/x2v/trilinear.hh> #endif // ! MLN_FUN_X2V_ALL_HH diff --git a/milena/mln/fun/x2v/trilinear.hh b/milena/mln/fun/x2v/trilinear.hh new file mode 100644 index 0000000..3015502 --- /dev/null +++ b/milena/mln/fun/x2v/trilinear.hh @@ -0,0 +1,139 @@ +// Copyright (C) 2008, 2009 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_FUN_X2V_TRILINEAR_HH +# define MLN_FUN_X2V_TRILINEAR_HH + +# include <mln/core/image/image2d.hh> +# include <mln/core/concept/function.hh> +# include <mln/fun/internal/selector.hh> +# include <mln/convert/to.hh> +# include <mln/algebra/vec.hh> + +/// \file mln/fun/x2v/trilinear.hh +/// +/// Define a trilinear interpolation of values from an underlying image +/// +/// Source: http://en.wikipedia.org/wiki/Rotation_matrix + +namespace mln +{ + + namespace fun + { + + namespace x2v + { + + /// Represent a trilinear interolation of values from an underlying image + /// + template < typename I > + struct trilinear + : public fun::internal::selector_<const algebra::vec<3,float>, + // 3,float is a dummy parameter (real is n,T) + mln_value(I), trilinear >::ret + { + typedef mln_value(I) result; + + trilinear(const I& ima); + + template <typename T> + mln_value(I) + operator()(const algebra::vec<3,T>& v) const; + + const I& ima; + }; + + +# ifndef MLN_INCLUDE_ONLY + + template <typename I> + trilinear::trilinear(const I& ima) : ima(ima) + { + mlc_bool(I::psite::dim == 3)::check(); + } + + + template <typename I> + template <typename T> + mln_value(I) + trilinear::operator()(const algebra::vec<3,T>& v) const + { + typedef mln_sum(mln_value(I)) vsum; + + double x = v[1]; // row + double y = v[2]; // col + double z = v[0]; // sli + + math::round<double> f; + unsigned x1 = f(std::floor(x)); + unsigned x2 = f(std::floor(x) + 1); + unsigned y1 = f(std::floor(y)); + unsigned y2 = f(std::floor(y) + 1); + unsigned z1 = f(std::floor(z)); + unsigned z2 = f(std::floor(z) + 1); + + double xd = x - x1; + double yd = y - y1; + double zd = z - z1; + + // interpolating in z direction + // Following access are supposed valid. + vsum i1 = ima(point3d(z1,x1,y1)) * (1 - zd) + + ima(point3d(z2,x1,y1)) * zd; + + vsum i2 = ima(point3d(z1,x1,y2)) * (1 - zd) + + ima(point3d(z2,x1,y2)) * zd; + + vsum j1 = ima(point3d(z1,x2,y1)) * (1 - zd) + + ima(point3d(z2,x2,y1)) * zd; + + vsum j2 = ima(point3d(z1,x2,y2)) * (1 - zd) + + ima(point3d(z2,x2,y2)) * zd; + + // interpolating in y direction + vsum w1 = i1 * (1 - yd) + i2 * yd; + vsum w2 = j1 * (1 - yd) + j2 * yd; + + // interpolating in x direction + vsum res = w1 * (1 - xd) + w2 * xd; + + return convert::to<mln_value(I)>(res); + } + + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::fun::x2v + + } // end of namespace mln::fun + +} // end of namespace mln + + +#endif // ! MLN_FUN_X2V_TRILINEAR_HH -- 1.5.6.5

15 years, 10 months

3386: Add aliases for 2D algebra::vec.

by Guillaume Lazzara

* mln/core/alias/vec2d.hh: new. new aliases for float and double 2d vectors. * mln/core/alias/all.hh, * mln/essential/2d.hh: include the new header. --- milena/ChangeLog | 10 ++++++++++ 1 files changed, 10 insertions(+), 0 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index f9e905f..62c0696 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -8,6 +8,16 @@ * mln/core/alias/all.hh, * mln/essential/2d.hh: include the new header. +2009-02-18 Guillaume Lazzara <z(a)lrde.epita.fr> + + Add aliases for 2D algebra::vec. + + * mln/core/alias/vec2d.hh: new. new aliases for float and double 2d + vectors. + + * mln/core/alias/all.hh, + * mln/essential/2d.hh: include the new header. + 2009-02-18 Edwin Carlinet <carlinet(a)lrde.epita.fr> Fix bugs in connected filters tests and add Makefile rules. -- 1.5.6.5

15 years, 10 months

milena r3383: Fix bugs in connected filters tests and add Makefile rules

by Edwin Carlinet

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2009-02-18 Edwin Carlinet <carlinet(a)lrde.epita.fr> Fix bugs in connected filters tests and add Makefile rules. * mln/canvas/morpho/connected_filter.hh: Fix bugs to avoid warning at compile. * tests/canvas/morpho/Makefile.am: Add check rule. * tests/canvas/morpho/connected_filter.cc: Fix include mistakes. --- mln/canvas/morpho/connected_filter.hh | 6 +++++- tests/canvas/morpho/Makefile.am | 4 ++++ tests/canvas/morpho/connected_filter.cc | 5 ++++- 3 files changed, 13 insertions(+), 2 deletions(-) Index: trunk/milena/tests/canvas/morpho/Makefile.am =================================================================== --- trunk/milena/tests/canvas/morpho/Makefile.am (revision 3382) +++ trunk/milena/tests/canvas/morpho/Makefile.am (revision 3383) @@ -1,3 +1,7 @@ ## Process this file through Automake to create Makefile.in -*- Makefile -*- include $(top_srcdir)/milena/tests/tests.mk + +check_PROGRAMS=connected_filter + +connected_filter_SOURCES=connected_filter.cc Index: trunk/milena/tests/canvas/morpho/connected_filter.cc =================================================================== --- trunk/milena/tests/canvas/morpho/connected_filter.cc (revision 3382) +++ trunk/milena/tests/canvas/morpho/connected_filter.cc (revision 3383) @@ -6,13 +6,16 @@ #include <mln/core/alias/neighb2d.hh> #include <mln/morpho/attribute/volume.hh> #include <mln/morpho/attribute/card.hh> -#include "connected_filter.hh" +#include <mln/canvas/morpho/connected_filter.hh> int main(int argc, char** argv) { using namespace mln; using value::int_u8; + if (argc < 2) + return 1; + std::cout << "Leveling filter test" << std::endl; typedef mln::image2d<int_u8> I; Index: trunk/milena/mln/canvas/morpho/connected_filter.hh =================================================================== --- trunk/milena/mln/canvas/morpho/connected_filter.hh (revision 3382) +++ trunk/milena/mln/canvas/morpho/connected_filter.hh (revision 3383) @@ -74,6 +74,8 @@ void take_as_init_fastest (trait::accumulator::when_pix::use_none, A& accu, const I& input, const unsigned p) { + (void)input; + (void)p; accu.take_as_init (); } @@ -81,6 +83,7 @@ void take_as_init (trait::accumulator::when_pix::use_p, A& accu, const I& input, const P& p) { + (void)input; accu.take_as_init (p); } @@ -88,6 +91,7 @@ void take_as_init (trait::accumulator::when_pix::use_none, A& accu, const I& input, const P& p) { + (void)input; accu.take_as_init (p); } @@ -279,11 +283,11 @@ const typename A::result& lambda) { trace::entering("canvas::morpho::impl::connected_filter_fastest"); - // FIXME: Tests? const I& input = exact(input_); const N& nbh = exact(nbh_); + (void)a_; mln_concrete(I) output; initialize(output, input);

15 years, 10 months

milena r3382: Merge leveling and algebraic filters into connected filter

by Edwin Carlinet

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2009-02-18 Edwin Carlinet <carlinet(a)lrde.epita.fr> Merge leveling and algebraic filters into connected filter. * mln/canvas/morpho/connected_filter.hh: File merging the two filters. * sandbox/edwin/filtres/connectes/connected_filter.hh: Add static tests about accumaltors properties when calling directly algebraic and leveling filters. Add leveling and algebraic facade. * sandbox/edwin/filtres/connectes/filter.cc: Fix bugs. * tests/canvas/morpho/connected_filter.cc: Test file for connected filters. --- mln/canvas/morpho/connected_filter.hh | 601 ++++++++++++++++++++ sandbox/edwin/filtres/connectes/connected_filter.hh | 45 + sandbox/edwin/filtres/connectes/filter.cc | 40 + tests/canvas/morpho/connected_filter.cc | 75 ++ 4 files changed, 755 insertions(+), 6 deletions(-) Index: trunk/milena/tests/canvas/morpho/connected_filter.cc =================================================================== --- trunk/milena/tests/canvas/morpho/connected_filter.cc (revision 0) +++ trunk/milena/tests/canvas/morpho/connected_filter.cc (revision 3382) @@ -0,0 +1,75 @@ + + +#include <mln/core/image/image2d.hh> +#include <mln/io/pgm/all.hh> +#include <mln/util/timer.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/morpho/attribute/volume.hh> +#include <mln/morpho/attribute/card.hh> +#include "connected_filter.hh" + +int main(int argc, char** argv) +{ + using namespace mln; + using value::int_u8; + + std::cout << "Leveling filter test" << std::endl; + + typedef mln::image2d<int_u8> I; + I lena; + + float elapsed; + mln::util::timer chrono; + mln::morpho::attribute::volume c; + mln::morpho::attribute::card c2; + int lambda = atoi(argv[1]); + + mln::io::pgm::load(lena, "../../../../img/lena.pgm"); + I out; + + std::cout << "Test algebraic" << std::endl; + chrono.start(); + out = mln::canvas::morpho::connected_filter(lena, c4(), c2, lambda, true); + elapsed = chrono.stop(); + std::cout << "(auto) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "alg_auto.pgm"); + + chrono.start(); + out = mln::canvas::morpho::algebraic_filter(lena, c4(), c2, lambda, true); + elapsed = chrono.stop(); + std::cout << "(fast) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "alg_fast.pgm"); + + chrono.start(); + out = mln::canvas::morpho::algebraic_filter(lena, c4(), c2, lambda, true); + elapsed = chrono.stop(); + std::cout << "(slow) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "alg_slow.pgm"); + +// Try force algebraic dispatch with wrong accu (volume). +// out = mln::canvas::morpho::algebraic_filter(lena, c4(), c, lambda, true); + + std::cout << "Test leveling" << std::endl; + chrono.start(); + out = mln::canvas::morpho::connected_filter(lena, c4(), c, lambda, true); + elapsed = chrono.stop(); + std::cout << "(auto) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "lev_auto.pgm"); + + chrono.start(); + out = mln::canvas::morpho::leveling_filter(lena, c4(), c, lambda, true); + elapsed = chrono.stop(); + std::cout << "(fast) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "lev_fast.pgm"); + + chrono.start(); + out = mln::canvas::morpho::leveling_filter(lena, c4(), c, lambda, true); + elapsed = chrono.stop(); + std::cout << "(slow) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "lev_slow.pgm"); + +// Try force leveling dispatch with wrong accu (card). +// mln::canvas::morpho::leveling_filter(lena, c4(), c2, lambda, true); + + +} Index: trunk/milena/mln/canvas/morpho/connected_filter.hh =================================================================== --- trunk/milena/mln/canvas/morpho/connected_filter.hh (revision 0) +++ trunk/milena/mln/canvas/morpho/connected_filter.hh (revision 3382) @@ -0,0 +1,601 @@ +// Copyright (C) 2007, 2008, 2009 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_CANVAS_MORPHO_CONNECTED_FILTER_HH +# define MLN_CANVAS_MORPHO_CONNECTED_FILTER_HH + +/// \file mln/canvas/morpho/connected_filter.hh +/// +/// Connected filters dispatch (algebraic & leveling filters). +/// +# include <mln/core/concept/image.hh> +# include <mln/core/concept/neighborhood.hh> +# include <mln/core/concept/accumulator.hh> + +# include <mln/trait/accumulators.hh> + +# include <mln/level/sort_psites.hh> +# include <mln/level/sort_offsets.hh> + + +namespace mln { + namespace canvas { + namespace morpho { + + // Facade Fwd Declaration + template <typename I, typename N, typename A> + mln_concrete(I) + connected_filter(const Image& input, const Neighborhood<N>& nbh, + const Accumulator<A>& a, const typename A::result& lambda, + bool increasing); + + template <typename I, typename N, typename A> + mln_concrete(I) + leveling_filter(const Image& input, const Neighborhood<N>& nbh, + const Accumulator<A>& a, const typename A::result& lambda, + bool increasing); + + template <typename I, typename N, typename A> + mln_concrete(I) + algebraic_filter(const Image& input, const Neighborhood<N>& nbh, + const Accumulator<A>& a, const typename A::result& lambda, + bool increasing); + +# ifndef MLN_INCLUDE_ONLY + + namespace impl { + + template <typename A, typename I> + void take_as_init_fastest (trait::accumulator::when_pix::use_none, A& accu, + const I& input, const unsigned p) + { + accu.take_as_init (); + } + + template <typename A, typename I, typename P> + void take_as_init (trait::accumulator::when_pix::use_p, A& accu, + const I& input, const P& p) + { + accu.take_as_init (p); + } + + template <typename A, typename I, typename P> + void take_as_init (trait::accumulator::when_pix::use_none, A& accu, + const I& input, const P& p) + { + accu.take_as_init (p); + } + + template <typename A, typename I, typename P> + void take_as_init (trait::accumulator::when_pix::use_pix, A& accu, + const I& input, const P& p) + { + accu.take_as_init (make::pix(input, p)); + } + + template <typename A, typename I, typename P> + void take_as_init (trait::accumulator::when_pix::use_v, A& accu, + const I& input, const P& p) + { + accu.take_as_init (make::pix(input, p)); + } + + template <typename A, typename I> + void take_as_init_fastest (trait::accumulator::when_pix::use_v, A& accu, + const I& input, const unsigned p) + { + accu.take_as_init (input.element (p)); + } + + + template <typename A, typename I, typename P> + void take_as_init (A& accu, const I& input, const P& p) + { + take_as_init (mln_trait_accumulator_when_pix(A)(), accu, input, p); + } + + template <typename A, typename I, typename P> + void take_as_init_fastest (A& accu, const I& input, const P& p) + { + take_as_init_fastest (mln_trait_accumulator_when_pix(A)(), accu, input, p); + } + + + namespace generic { + + + //////////////////////// + /// Generic version. /// + //////////////////////// + + + template <typename I> + static inline + mln_psite(I) + find_root(I& parent, const mln_psite(I) & x) + { + if (parent(x) == x) + return x; + else + return parent(x) = find_root(parent, parent(x)); + } + + template <typename I, typename N, typename S, typename A> + mln_concrete(I) + connected_filter (const Image& input_, + const Neighborhood<N>& nbh_, + const Site_Set<S>& s_, + const Accumulator<A>& a_, + const typename A::result& lambda) + { + trace::entering("canvas::morpho::impl::generic::connected_filter"); + // FIXME: Test?! + + const I& input = exact(input_); + const N& nbh = exact(nbh_); + const S& s = exact(s_); + (void)a_; // To avoid warning at compilation + + mln_concrete(I) output; + initialize(output, input); + + // Local type. + typedef mln_psite(I) P; + + + // Auxiliary data. + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, bool) activity; + mln_ch_value(I, P) parent; + mln_ch_value(I, A) data; + + // Initialization. + { + initialize(deja_vu, input); + mln::data::fill(deja_vu, false); + initialize(activity, input); + mln::data::fill(activity, true); + initialize(parent, input); + initialize(data, input); + //a.init(); // init required. + } + + // First pass. + { + mln_fwd_piter(S) p(s); // s required. + mln_niter(N) n(nbh, p); + + for_all(p) + { + // Make set. + { + parent(p) = p; + + // Check accumulator trait to handle argument type (Pix or Site). + take_as_init (data(p), input, p); + } + + for_all(n) + if (input.domain().has(n) && deja_vu(n)) + { + //do_union(n, p); + P r = find_root(parent, n); + if (r != p) + { + if (input(r) == input(p) || (activity(r) && (data(r) < lambda))) // Equiv(r, p) + // Either a flat zone or the component of r is still growing. + { + /* FIXME: Same remark as above concerning the + initialization of data(p); instead of + + data(p).take(data(r)); + + we should (or could) have + + unite_data(p, r); + + so as to keep the generic aspect of this canvas + (as long as the set of acceptable types for the + template parameter A is not bound). */ + + data(p).take(data(r)); + parent(r) = p; + if (activity(r) == false) + activity(p) = false; + } + else + { + activity(p) = false; + } + } + } + deja_vu(p) = true; + } + } + + // Second pass. + { + mln_bkd_piter(S) p(s); + for_all(p) + if (parent(p) == p) // p is root. + output(p) = input(p); + else + output(p) = output(parent(p)); + } + + trace::exiting("canvas::morpho::impl::generic::connected_filter"); + return output; + } + + } // end of namespace mln::canvas::morpho::impl::generic + + + //////////////////////// + /// Fastest version. /// + //////////////////////// + + template <typename I> + inline + unsigned + find_root_fastest(I& parent, unsigned x) + { + if (parent.element(x) == x) + return x; + else + return parent.element(x) = find_root_fastest(parent, parent.element(x)); + } + + template <typename I, typename N, typename A> + mln_concrete(I) + connected_filter_fastest(const Image& input_, + const Neighborhood<N>& nbh_, + const util::array<unsigned>& s, + const Accumulator<A>& a_, + const typename A::result& lambda) + { + trace::entering("canvas::morpho::impl::connected_filter_fastest"); + + // FIXME: Tests? + + const I& input = exact(input_); + const N& nbh = exact(nbh_); + + mln_concrete(I) output; + initialize(output, input); + + // Local type. + typedef mln_psite(I) P; + + // Auxiliary data. + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, bool) activity; + mln_ch_value(I, unsigned) parent; + mln_ch_value(I, A) data; + + // Initialization. + { + initialize(deja_vu, input); + mln::data::fill(deja_vu, false); + initialize(activity, input); + mln::data::fill(activity, true); + initialize(parent, input); + mln::data::fill(parent, 0); + initialize(data, input); + } + + util::array<int> dp = offsets_wrt(input, nbh); + const unsigned n_nbhs = dp.nelements(); + const unsigned n_points = s.nelements(); + + // First pass. + { + for (unsigned i = 0; i < n_points; ++i) + { + unsigned p = s[i]; // An offset. + + // Make set. + parent.element(p) = p; + + // Check accumulator trait to handle argument type (Value or None). + take_as_init_fastest (data.element(p), input, p); + + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p + dp[j]; + if (!deja_vu.element(n)) + continue; + + unsigned r = find_root_fastest(parent, n); + if (r != p) + { + if (input.element(r) == input.element(p) + || (activity.element(r) + && (data.element(r) < lambda))) + { + data.element(p).take(data.element(r)); + parent.element(r) = p; + if (activity.element(r) == false) + activity.element(p) = false; + } + else + activity.element(p) = false; + } + } + + deja_vu.element(p) = true; + } + } + + + // Second pass. + { + for (int i = n_points - 1; i >= 0 ; --i) + { + unsigned p = s[i]; + if (parent.element(p) == p) // p is root. + output.element(p) = input.element(p); + else + output.element(p) = output.element(parent.element(p)); + } + } + + trace::exiting("canvas::morpho::impl::connected_filter_fastest"); + return output; + } + + + } // end of namespace mln::canvas::morpho::impl + + + + + // Dispatch. + + + namespace internal + { + // Leveling + template <typename I, typename N, typename A> + mln_concrete(I) + leveling_filter_dispatch(metal::false_, + const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + p_array < mln_psite(I) > s = + increasing ? + level::sort_psites_increasing(input) : + level::sort_psites_decreasing(input); + return impl::generic::connected_filter(input, nbh, s, a, lambda); + } + + template <typename I, typename N, typename A> + mln_concrete(I) + leveling_filter_dispatch(metal::true_, + const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + util::array<unsigned> s = + increasing ? + level::sort_offsets_increasing(input) : + level::sort_offsets_decreasing(input); + return impl::connected_filter_fastest(input, nbh, s, a, lambda); + } + + template <typename I, typename N, typename A> + inline + mln_concrete(I) + leveling_filter_dispatch(const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + mlc_or(mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_pix), + mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_v))::check(); + + enum + { + test = mlc_equal(mln_trait_image_speed(I), + trait::image::speed::fastest)::value + && mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_v)::value + && mln_is_simple_neighborhood(N)::value + }; + return leveling_filter_dispatch(metal::bool_<test>(), + input, nbh, a, lambda, increasing); + } + + // Alegebraic + template <typename I, typename N, typename A> + inline + mln_concrete(I) + algebraic_filter_dispatch(metal::false_, + const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& accu, + const mln_result(A)& lambda, + bool increasing) + { + p_array<mln_psite(I)> s = increasing ? + level::sort_psites_increasing(input) : + level::sort_psites_decreasing(input); + return impl::generic::connected_filter(input, nbh, s, accu, lambda); + } + + template <typename I, typename N, typename A> + inline + mln_concrete(I) + algebraic_filter_dispatch(metal::true_, + const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& accu, + const mln_result(A)& lambda, + bool increasing) + { + util::array<unsigned> s = increasing ? + level::sort_offsets_increasing(input) : + level::sort_offsets_decreasing(input); + return impl::connected_filter_fastest(input, nbh, s, accu, lambda); + } + + template <typename I, typename N, typename A> + inline + mln_concrete(I) + algebraic_filter_dispatch(const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& accu, + const mln_result(A)& lambda, + bool increasing) + { + mlc_or(mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_none), + mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_p))::check(); + + enum { + test = (mlc_equal(mln_trait_image_speed(I), + trait::image::speed::fastest)::value && + mln_is_simple_neighborhood(N)::value && + mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_none)::value) + }; + return algebraic_filter_dispatch(metal::bool_<test>(), input, nbh, + accu, lambda, increasing); + } + + + + //connected + template <typename I, typename N, typename A> + inline + mln_concrete(I) + connected_filter_dispatch(mln::trait::accumulator::when_pix::use_none, + const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + return algebraic_filter_dispatch(input, nbh, a, lambda, increasing); + } + + template <typename I, typename N, typename A> + inline + mln_concrete(I) + connected_filter_dispatch(mln::trait::accumulator::when_pix::use_p, + const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + return algebraic_filter_dispatch(input, nbh, a, lambda, increasing); + } + + template <typename I, typename N, typename A> + inline + mln_concrete(I) + connected_filter_dispatch(mln::trait::accumulator::when_pix::use_v, + const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + return leveling_filter_dispatch(input, nbh, a, lambda, increasing); + } + + template <typename I, typename N, typename A> + inline + mln_concrete(I) + connected_filter_dispatch(mln::trait::accumulator::when_pix::use_pix, + const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + return leveling_filter_dispatch(input, nbh, a, lambda, increasing); + } + } // end of namespace mln::canvas::morpho::internal + + +# endif // ! MLN_INCLUDE_ONLY + + // Facade. + + template <typename I, typename N, typename A> + mln_concrete(I) + connected_filter(const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + return internal::connected_filter_dispatch(mln_trait_accumulator_when_pix(A)(), + input, nbh, a, lambda, increasing); + } + + template <typename I, typename N, typename A> + mln_concrete(I) + algebraic_filter(const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + return internal::algebraic_filter_dispatch(input, nbh, a, lambda, increasing); + } + + template <typename I, typename N, typename A> + mln_concrete(I) + leveling_filter(const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + return internal::leveling_filter_dispatch(input, nbh, a, lambda, increasing); + } + + + + } // end of namespace mln::canvas::morpho + } // end of namespace mln::canvas +} // end of namespace mln + + +#endif // ! MLN_CANVAS_MORPHO_CONNECTED_FILTER_HH Index: trunk/milena/sandbox/edwin/filtres/connectes/connected_filter.hh =================================================================== --- trunk/milena/sandbox/edwin/filtres/connectes/connected_filter.hh (revision 3381) +++ trunk/milena/sandbox/edwin/filtres/connectes/connected_filter.hh (revision 3382) @@ -54,6 +54,18 @@ const Accumulator<A>& a, const typename A::result& lambda, bool increasing); + template <typename I, typename N, typename A> + mln_concrete(I) + leveling_filter(const Image& input, const Neighborhood<N>& nbh, + const Accumulator<A>& a, const typename A::result& lambda, + bool increasing); + + template <typename I, typename N, typename A> + mln_concrete(I) + algebraic_filter(const Image& input, const Neighborhood<N>& nbh, + const Accumulator<A>& a, const typename A::result& lambda, + bool increasing); + # ifndef MLN_INCLUDE_ONLY namespace impl { @@ -409,6 +421,11 @@ const typename A::result& lambda, bool increasing) { + mlc_or(mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_pix), + mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_v))::check(); + enum { test = mlc_equal(mln_trait_image_speed(I), @@ -463,6 +480,11 @@ const mln_result(A)& lambda, bool increasing) { + mlc_or(mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_none), + mlc_equal(mln_trait_accumulator_when_pix(A), + trait::accumulator::when_pix::use_p))::check(); + enum { test = (mlc_equal(mln_trait_image_speed(I), trait::image::speed::fastest)::value && @@ -547,6 +569,29 @@ input, nbh, a, lambda, increasing); } + template <typename I, typename N, typename A> + mln_concrete(I) + algebraic_filter(const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + return internal::algebraic_filter_dispatch(input, nbh, a, lambda, increasing); + } + + template <typename I, typename N, typename A> + mln_concrete(I) + leveling_filter(const Image& input, + const Neighborhood<N>& nbh, + const Accumulator<A>& a, + const typename A::result& lambda, + bool increasing) + { + return internal::leveling_filter_dispatch(input, nbh, a, lambda, increasing); + } + + } // end of namespace mln::canvas::morpho } // end of namespace mln::canvas Index: trunk/milena/sandbox/edwin/filtres/connectes/filter.cc =================================================================== --- trunk/milena/sandbox/edwin/filtres/connectes/filter.cc (revision 3381) +++ trunk/milena/sandbox/edwin/filtres/connectes/filter.cc (revision 3382) @@ -5,6 +5,7 @@ #include <mln/util/timer.hh> #include <mln/core/alias/neighb2d.hh> #include <mln/morpho/attribute/volume.hh> +#include <mln/morpho/attribute/card.hh> #include "connected_filter.hh" int main(int argc, char** argv) @@ -20,28 +21,55 @@ float elapsed; mln::util::timer chrono; mln::morpho::attribute::volume c; + mln::morpho::attribute::card c2; int lambda = atoi(argv[1]); mln::io::pgm::load(lena, "../../../../img/lena.pgm"); I out; + std::cout << "Test algebraic" << std::endl; chrono.start(); - out = mln::canvas::morpho::connected_filter(lena, c4(), c, lambda, true); + out = mln::canvas::morpho::connected_filter(lena, c4(), c2, lambda, true); elapsed = chrono.stop(); std::cout << "(auto) " << elapsed << "s" << std::endl; - mln::io::pgm::save(out, "auto.pgm"); + mln::io::pgm::save(out, "alg_auto.pgm"); chrono.start(); - out = mln::canvas::morpho::internal::leveling_filter_dispatch(mln::metal::true_(), lena, c4(), c, lambda, true); + out = mln::canvas::morpho::algebraic_filter(lena, c4(), c2, lambda, true); elapsed = chrono.stop(); std::cout << "(fast) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "alg_fast.pgm"); + + chrono.start(); + out = mln::canvas::morpho::algebraic_filter(lena, c4(), c2, lambda, true); + elapsed = chrono.stop(); + std::cout << "(slow) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "alg_slow.pgm"); + +// Try force algebraic dispatch with wrong accu (volume). +// out = mln::canvas::morpho::algebraic_filter(lena, c4(), c, lambda, true); + + std::cout << "Test leveling" << std::endl; + chrono.start(); + out = mln::canvas::morpho::connected_filter(lena, c4(), c, lambda, true); + elapsed = chrono.stop(); + std::cout << "(auto) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "lev_auto.pgm"); - mln::io::pgm::save(out, "fast.pgm"); + chrono.start(); + out = mln::canvas::morpho::leveling_filter(lena, c4(), c, lambda, true); + elapsed = chrono.stop(); + std::cout << "(fast) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "lev_fast.pgm"); chrono.start(); - out = mln::canvas::morpho::internal::leveling_filter_dispatch(mln::metal::false_(), lena, c4(), c, lambda, true); + out = mln::canvas::morpho::leveling_filter(lena, c4(), c, lambda, true); elapsed = chrono.stop(); std::cout << "(slow) " << elapsed << "s" << std::endl; + mln::io::pgm::save(out, "lev_slow.pgm"); + +// Try force leveling dispatch with wrong accu (card). +// mln::canvas::morpho::leveling_filter(lena, c4(), c2, lambda, true); + - mln::io::pgm::save(out, "slow.pgm"); }

15 years, 10 months

r3381: Add sandbox for ImageMagick support

by Fabien Freling

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-02-17 Fabien Freling <fabien.freling(a)lrde.epita.fr> Add sandbox for ImageMagick support. * fabien/dicom/load.hh: Minor update of dcmtk support. * fabien/dicom/save.hh: Minor update of dcmtk support. * fabien/magick/Makefile: New file to help testing. * fabien/magick/load.hh: ImageMagick load support. * fabien/magick/magick.cc: Test file. --- dicom/load.hh | 591 ++++++++++++++++++++----------------------------------- dicom/save.hh | 220 -------------------- magick/Makefile | 2 magick/load.hh | 142 +++++++++++++ magick/magick.cc | 14 + 5 files changed, 374 insertions(+), 595 deletions(-) Index: trunk/milena/sandbox/fabien/dicom/save.hh =================================================================== --- trunk/milena/sandbox/fabien/dicom/save.hh (revision 3380) +++ trunk/milena/sandbox/fabien/dicom/save.hh (revision 3381) @@ -949,223 +949,3 @@ return 0; } - -//******************************* - -/* -** CVS/RCS Log: -** $Log: dump2dcm.cc,v $ -** Revision 1.51 2005/12/16 09:07:03 onken -** - Added variable initialization to avoid compiler warning -** -** Revision 1.50 2005/12/08 15:40:50 meichel -** Changed include path schema for all DCMTK header files -** -** Revision 1.49 2004/07/13 09:43:10 meichel -** Fixed memory leak occuring when raw data is read from file. -** -** Revision 1.48 2004/03/05 09:59:00 joergr -** Avoid wrong warning for LUTData (0028,3006) having a VR of US or SS. -** Added initial "hooks" for (compressed) pixel items. -** Added "ignore errors" option (similar to dcmdump). -** -** Revision 1.47 2004/01/16 10:53:16 joergr -** Adapted type casts to new-style typecast operators defined in ofcast.h. -** Removed acknowledgements with e-mail addresses from CVS log. -** -** Revision 1.46 2003/11/05 16:15:27 meichel -** Removed useless "--write-xfer-same" command line option -** -** Revision 1.45 2002/12/05 13:59:29 joergr -** Fixed typo. -** -** Revision 1.44 2002/11/27 12:07:18 meichel -** Adapted module dcmdata to use of new header file ofstdinc.h -** -** Revision 1.43 2002/11/26 08:43:02 meichel -** Replaced all includes for "zlib.h" with <zlib.h> -** to avoid inclusion of zlib.h in the makefile dependencies. -** -** Revision 1.42 2002/09/23 17:52:04 joergr -** Prepared code for future support of 'config.guess' host identifiers. -** -** Revision 1.41 2002/09/23 13:50:42 joergr -** Added new command line option "--version" which prints the name and version -** number of external libraries used. -** -** Revision 1.40 2002/08/21 10:14:16 meichel -** Adapted code to new loadFile and saveFile methods, thus removing direct -** use of the DICOM stream classes. -** -** Revision 1.39 2002/04/16 13:38:55 joergr -** Added configurable support for C++ ANSI standard includes (e.g. streams). -** -** Revision 1.38 2001/12/11 14:00:39 joergr -** Fixed bug in 'dump2dcm' parser causing AT attribute values to be ignored. -** -** Revision 1.37 2001/11/09 15:50:53 joergr -** Renamed some of the getValue/getParam methods to avoid ambiguities reported -** by certain compilers. -** -** Revision 1.36 2001/09/25 17:21:01 meichel -** Adapted dcmdata to class OFCondition -** -** Revision 1.35 2001/06/01 15:48:30 meichel -** Updated copyright header -** -** Revision 1.34 2000/04/14 15:42:54 meichel -** Global VR generation flags are now derived from OFGlobal and, thus, -** safe for use in multi-thread applications. -** -** Revision 1.33 2000/03/08 16:26:06 meichel -** Updated copyright header. -** -** Revision 1.32 2000/03/06 18:09:38 joergr -** Avoid empty statement in the body of if-statements (MSVC6 reports warnings). -** -** Revision 1.31 2000/03/03 14:05:16 meichel -** Implemented library support for redirecting error messages into memory -** instead of printing them to stdout/stderr for GUI applications. -** -** Revision 1.30 2000/02/29 11:48:51 meichel -** Removed support for VS value representation. This was proposed in CP 101 -** but never became part of the standard. -** -** Revision 1.29 2000/02/23 15:11:36 meichel -** Corrected macro for Borland C++ Builder 4 workaround. -** -** Revision 1.28 2000/02/10 16:02:51 joergr -** Enhanced handling of PixelData/Item element. Externally stored raw data is -** now always imported as little endian and swapped if necessary. This change -** reflects the new 'export' feature of dcmdump. -** -** Revision 1.27 2000/02/01 10:11:58 meichel -** Avoiding to include <stdlib.h> as extern "C" on Borland C++ Builder 4, -** workaround for bug in compiler header files. -** -** Revision 1.26 1999/05/03 14:13:40 joergr -** Minor code purifications to keep Sun CC 2.0.1 quiet. -** -** Revision 1.25 1999/04/27 17:50:53 joergr -** Adapted console applications to new OFCommandLine and OFConsoleApplication -** functionality. -** -** Revision 1.24 1999/04/27 12:23:27 meichel -** Prevented dcmdata applications from opening a file with empty filename, -** leads to application crash on Win32. -** -** Revision 1.23 1999/03/31 09:24:23 meichel -** Updated copyright header in module dcmdata -** -** Revision 1.22 1999/03/29 10:14:15 meichel -** Adapted command line options of dcmdata applications to new scheme. -** -** Revision 1.21 1999/03/22 16:16:01 meichel -** dump2dcm now allows to include the contents of binary files -** as OB/OW values while converting a dump to a DICOM file. -** -** Revision 1.20 1999/01/07 14:13:12 meichel -** Corrected bug in dump2dcm that prevented the correct processing of -** dumps created with dcmdump if they contained the "internal" VR markers -** "xs" (US or SS) or "ox" (OB or OW). -** -** Revision 1.19 1998/01/27 10:51:27 meichel -** Removed some unused variables, meaningless const modifiers -** and unreached statements. -** -** Revision 1.18 1998/01/14 14:41:15 hewett -** Modified existing -u command line option to also disable generation -** of UT and VS (previously just disabled generation of UN). -** -** Revision 1.17 1997/08/05 07:34:54 andreas -** Corrected Error handling of SQ in dump2dcm -** -** Revision 1.16 1997/07/21 07:59:02 andreas -** - Deleted support for DcmPixelItems and DcmPixelSequences in dump2dcm -** ToDo: New support should be added in the future compatible to -** the new DcmPixel class. -** - Replace all boolean types (BOOLEAN, CTNBOOLEAN, DICOM_BOOL, BOOL) -** with one unique boolean type OFBool. -** -** Revision 1.15 1997/07/03 15:09:40 andreas -** - removed debugging functions Bdebug() and Edebug() since -** they write a static array and are not very useful at all. -** Cdebug and Vdebug are merged since they have the same semantics. -** The debugging functions in dcmdata changed their interfaces -** (see dcmdata/include/dcdebug.h) -** -** Revision 1.14 1997/05/30 06:44:57 andreas -** - fixed scanf format problem leading to warnings on 64 bit machines. -** -** Revision 1.13 1997/05/29 15:52:52 meichel -** Added constant for dcmtk release date in dcuid.h. -** All dcmtk applications now contain a version string -** which is displayed with the command line options ("usage" message) -** and which can be queried in the binary with the "ident" command. -** -** Revision 1.12 1997/05/22 13:26:25 hewett -** Modified the test for presence of a data dictionary to use the -** method DcmDataDictionary::isDictionaryLoaded(). -** -** Revision 1.11 1997/05/20 07:57:12 andreas -** - Removed obsolete applications file2ds and ds2file. The functionality of these -** applications is now peformed by dcmconv. Unified calling parameters -** are implemented in dump2dcm, dcmdump and dcmconv. -** -** Revision 1.10 1997/05/16 08:31:06 andreas -** - Revised handling of GroupLength elements and support of -** DataSetTrailingPadding elements. The enumeratio E_GrpLenEncoding -** got additional enumeration values (for a description see dctypes.h). -** addGroupLength and removeGroupLength methods are replaced by -** computeGroupLengthAndPadding. To support Padding, the parameters of -** element and sequence write functions changed. -** -** Revision 1.9 1997/04/18 08:06:56 andreas -** - Minor corrections: correct some warnings of the SUN-C++ Compiler -** concerning the assignments of wrong types and inline compiler -** errors -** - The put/get-methods for all VRs did not conform to the C++-Standard -** draft. Some Compilers (e.g. SUN-C++ Compiler, Metroworks -** CodeWarrier, etc.) create many warnings concerning the hiding of -** overloaded get methods in all derived classes of DcmElement. -** So the interface of all value representation classes in the -** library are changed rapidly, e.g. -** OFCondition get(Uint16 & value, const unsigned long pos); -** becomes -** OFCondition getUint16(Uint16 & value, const unsigned long pos); -** All (retired) "returntype get(...)" methods are deleted. -** For more information see dcmdata/include/dcelem.h -** -** Revision 1.8 1997/03/27 15:47:25 hewett -** Added command line switche to allow generation of UN to be -** disabled (it is enabled by default). -** -** Revision 1.7 1996/09/24 16:13:51 hewett -** Added preliminary support for the Macintosh environment (GUSI library). -** -** Revision 1.6 1996/05/02 17:00:23 hewett -** Corrected program name in usage description. -** -** Revision 1.5 1996/05/02 15:55:11 hewett -** Stopped whitespace being stripped from inside value strings when -** no [] delimiter present. Now only leading and trailing whitespace -** is stripped. -** -** Revision 1.4 1996/04/27 12:13:01 hewett -** Corrected bug in last bug-fix. A tag value [some text] was being -** parsed as an empty string. Now both [] and [some text] appear to -** work as intended. -** -** Revision 1.3 1996/03/22 12:38:44 andreas -** Correct some mistakes: handling [] as empty string (no value field) -** handling =Name correct if Name is not correct -** -** Revision 1.2 1996/03/12 15:11:39 hewett -** Added call to prepareCmdLineArgs to enable command line arguments -** in environments which do not provide them. -** -** Revision 1.1 1996/01/29 13:36:38 andreas -** dump2dcm added convert ASCII descriptions into DICOM files -** -** -*/ Index: trunk/milena/sandbox/fabien/dicom/load.hh =================================================================== --- trunk/milena/sandbox/fabien/dicom/load.hh (revision 3380) +++ trunk/milena/sandbox/fabien/dicom/load.hh (revision 3381) @@ -1,3 +1,30 @@ +// Copyright (C) 2009 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. + /* * * Copyright (C) 1994-2005, OFFIS @@ -30,45 +57,109 @@ * */ -#include "dcmtk/config/osconfig.h" /* make sure OS specific configuration is included first */ -#include "dcmtk/ofstd/ofstream.h" -#include "dcmtk/dcmdata/dctk.h" -#include "dcmtk/dcmdata/dcdebug.h" -#include "dcmtk/dcmdata/cmdlnarg.h" -#include "dcmtk/ofstd/ofconapp.h" -#include "dcmtk/dcmdata/dcuid.h" /* for dcmtk version name */ -#include "dcmtk/dcmdata/dcistrmz.h" /* for dcmZlibExpectRFC1950Encoding */ +#ifndef MLN_IO_DICOM_LOAD_HH +# define MLN_IO_DICOM_LOAD_HH + +/*! + * \file mln/io/dicom/load.hh + * + * \brief Define a function which loads an image of kind dicom with + * given path. + * + */ + +# include <mln/io/pnm/load.hh> + +# include <dcmtk/config/osconfig.h> /* make sure OS specific configuration is included first */ +# include <dcmtk/ofstd/ofstream.h> +# include <dcmtk/dcmdata/dctk.h> +# include <dcmtk/dcmdata/dcdebug.h> +# include <dcmtk/dcmdata/cmdlnarg.h> +# include <dcmtk/ofstd/ofconapp.h> +# include <dcmtk/dcmdata/dcuid.h> /* for dcmtk version name */ +# include <dcmtk/dcmdata/dcistrmz.h> /* for dcmZlibExpectRFC1950Encoding */ #define INCLUDE_CSTDLIB #define INCLUDE_CSTRING -#include "dcmtk/ofstd/ofstdinc.h" +# include <dcmtk/ofstd/ofstdinc.h> #ifdef WITH_ZLIB #include <zlib.h> /* for zlibVersion() */ #endif -#define OFFIS_CONSOLE_APPLICATION "dcmdump" +# define SHORTCOL 3 +# define LONGCOL 20 +# define PATH_SEPARATOR '/' + + + +namespace mln +{ + + namespace io + { + + namespace dicom + { + + /*! Load a dicom image in a Milena image. + * + * \param[out] ima A reference to the image which will receive + * data. + * \param[in] filename The source. + */ + template <typename I> + void load(Image& ima, + const std::string& filename); + + /*! Load a dicom image in a Milena image. To use this routine, you + * should specialize the template whith the value type of the + * image loaded. (ex : load<value::int_u8>("...") ) + * + * \param[in] filename The image source. + * + * \return An image2d which contains loaded data. + */ + template <typename V> + image2d<V> load(const std::string& filename); + + /*! Load a dicom image in a Milena image. To use this routine, you + * should specialize the template whith the value type of the + * image loaded. (ex : load<value::int_u8>("...") ) + * + * \param[in] filename The image source. + * + * \return An image2d which contains loaded data. + */ + template <typename V> + image3d<V> load(const std::string& filename); + +# ifndef MLN_INCLUDE_ONLY + + template <typename V> + inline + image2d<V> load(const std::string& filename) + { + trace::entering("mln::io::dicom::load"); + image2d<V> ima;// = io::pnm::load<V>(DICOM, filename); + trace::exiting("mln::io::dicom::load"); + return ima; + } + + template <typename V> + inline + image3d<V> load(const std::string& filename) + { + trace::entering("mln::io::dicom::load"); + image2d<V> ima;// = io::pnm::load<V>(DICOM, filename); + trace::exiting("mln::io::dicom::load"); + return ima; + } + -static char rcsid[] = "$dcmtk: " OFFIS_CONSOLE_APPLICATION " v" - OFFIS_DCMTK_VERSION " " OFFIS_DCMTK_RELEASEDATE " $"; -#ifdef HAVE_GUSI_H - /* needed for Macintosh */ -#include <GUSI.h> -#include <SIOUX.h> -#endif -static int dumpFile(ostream & out, - const char *ifname, - const E_FileReadMode readMode, - const E_TransferSyntax xfer, - const size_t printFlags, - const OFBool loadIntoMemory, - const OFBool stopOnErrors, - const OFBool writePixelData, - const char *pixelDirectory); -// ******************************************** static OFBool printAllInstances = OFTrue; static OFBool prependSequenceHierarchy = OFFalse; @@ -81,7 +172,7 @@ static OFBool addPrintTagName(const char* tagName) { if (printTagCount >= MAX_PRINT_TAG_NAMES) { - CERR << "error: too many print Tag options (max: " << MAX_PRINT_TAG_NAMES << ")" << endl; + std::cerr << "error: too many print Tag options (max: " << MAX_PRINT_TAG_NAMES << ")" << endl; return OFFalse; } @@ -92,16 +183,21 @@ /* it is a name */ const DcmDataDictionary& globalDataDict = dcmDataDict.rdlock(); const DcmDictEntry *dicent = globalDataDict.findEntry(tagName); - if( dicent == NULL ) { - CERR << "error: unrecognised tag name: '" << tagName << "'" << endl; + if (dicent == NULL) + { + std::cerr << "error: unrecognised tag name: '" << tagName << "'" << endl; dcmDataDict.unlock(); return OFFalse; - } else { + } + else + { /* note for later */ printTagKeys[printTagCount] = new DcmTagKey(dicent->getKey()); } dcmDataDict.unlock(); - } else { + } + else + { /* tag name has format xxxx,xxxx */ /* do not lookup in dictionary, tag could be private */ printTagKeys[printTagCount] = NULL; @@ -112,327 +208,13 @@ return OFTrue; } -#define SHORTCOL 3 -#define LONGCOL 20 - -int main(int argc, char *argv[]) -{ - int opt_debugMode = 0; - OFBool loadIntoMemory = OFTrue; - size_t printFlags = DCMTypes::PF_shortenLongTagValues /*| DCMTypes::PF_showTreeStructure*/; - OFBool printFilename = OFFalse; - OFBool writePixelData = OFFalse; - E_FileReadMode readMode = ERM_autoDetect; - E_TransferSyntax xfer = EXS_Unknown; - OFBool stopOnErrors = OFTrue; - const char *current = NULL; - const char *pixelDirectory = NULL; - - -#ifdef HAVE_GUSI_H - /* needed for Macintosh */ - /* set options for the Metrowerks CodeWarrior SIOUX console */ - SIOUXSettings.autocloseonquit = OFFalse; - SIOUXSettings.asktosaveonclose = OFFalse; - SIOUXSettings.showstatusline = OFTrue; - SIOUXSettings.setupmenus = OFTrue; - /* set options for the GUSI sockets library */ - GUSISetup(GUSIwithSIOUXSockets); - GUSISetup(GUSIwithInternetSockets); -#endif - - SetDebugLevel(( 0 )); - - OFConsoleApplication app(OFFIS_CONSOLE_APPLICATION, "Dump DICOM file and data set", rcsid); - OFCommandLine cmd; - cmd.setOptionColumns(LONGCOL, SHORTCOL); - cmd.setParamColumn(LONGCOL + SHORTCOL + 4); - - cmd.addParam("dcmfile-in", "DICOM input filename to be dumped", OFCmdParam::PM_MultiMandatory); - - cmd.addGroup("general options:", LONGCOL, SHORTCOL + 2); - cmd.addOption("--help", "-h", "print this help text and exit"); - cmd.addOption("--version", "print version information and exit", OFTrue /* exclusive */); -#ifdef USE_EXPERIMENTAL_QUIET_MODE - cmd.addOption("--quiet", "-q", "quiet mode, print no warnings and errors"); -#endif - cmd.addOption("--debug", "-d", "debug mode, print debug information"); - - cmd.addGroup("input options:"); - cmd.addSubGroup("input file format:"); - cmd.addOption("--read-file", "+f", "read file format or data set (default)"); - cmd.addOption("--read-file-only", "+fo", "read file format only"); - cmd.addOption("--read-dataset", "-f", "read data set without file meta information"); - cmd.addSubGroup("input transfer syntax:"); - cmd.addOption("--read-xfer-auto", "-t=", "use TS recognition (default)"); - cmd.addOption("--read-xfer-detect", "-td", "ignore TS specified in the file meta header"); - cmd.addOption("--read-xfer-little", "-te", "read with explicit VR little endian TS"); - cmd.addOption("--read-xfer-big", "-tb", "read with explicit VR big endian TS"); - cmd.addOption("--read-xfer-implicit", "-ti", "read with implicit VR little endian TS"); - cmd.addSubGroup("parsing of odd-length attributes:"); - cmd.addOption("--accept-odd-length", "+ao", "accept odd length attributes (default)"); - cmd.addOption("--assume-even-length", "+ae", "assume real length is one byte larger"); - cmd.addSubGroup("handling of undefined length UN elements:"); - cmd.addOption("--enable-cp246", "+ui", "read undefined len UN as implicit VR (default)"); - cmd.addOption("--disable-cp246", "-ui", "read undefined len UN as explicit VR"); - cmd.addSubGroup("handling of defined length UN elements:"); - cmd.addOption("--retain-un", "-uc", "retain elements as UN (default)"); - cmd.addOption("--convert-un", "+uc", "convert to real VR if known"); - cmd.addSubGroup("automatic data correction:"); - cmd.addOption("--enable-correction", "+dc", "enable automatic data correction (default)"); - cmd.addOption("--disable-correction", "-dc", "disable automatic data correction"); -#ifdef WITH_ZLIB - cmd.addSubGroup("bitstream format of deflated input:"); - cmd.addOption("--bitstream-deflated", "+bd", "expect deflated bitstream (default)"); - cmd.addOption("--bitstream-zlib", "+bz", "expect deflated zlib bitstream"); -#endif - - cmd.addGroup("output options:"); - cmd.addSubGroup("printing:"); - cmd.addOption("--load-all", "+M", "load very long tag values (default)"); - cmd.addOption("--load-short", "-M", "do not load very long values (e.g. pixel data)"); - cmd.addOption("--max-read-length", "+R", 1, "[k]bytes: integer [4..4194302] (default: 4)", - "set threshold for long values to k kbytes"); - cmd.addOption("--print-all", "+L", "print long tag values completely"); - cmd.addOption("--print-short", "-L", "print long tag values shortened (default)"); - cmd.addOption("--print-filename", "+F", "print header with filename for each input file"); - - cmd.addSubGroup("error handling:"); - cmd.addOption("--stop-on-error", "-E", "do not print if file is damaged (default)"); - cmd.addOption("--ignore-errors", "+E", "attempt to print even if file is damaged"); - - cmd.addSubGroup("searching:"); - cmd.addOption("--search", "+P", 1, "[t]ag: \"xxxx,xxxx\" or a data dictionary name", - "print the value of tag t\nthis option can be specified multiple times\n(default: the complete file is printed)"); - - cmd.addOption("--search-all", "+s", "print all instances of searched tags (default)"); - cmd.addOption("--search-first", "-s", "only print first instance of searched tags"); - - cmd.addOption("--prepend", "+p", "prepend sequence hierarchy to printed tag,\ndenoted by: (xxxx,xxxx).(xxxx,xxxx).*\n(only with --search-all or --search-first)"); - cmd.addOption("--no-prepend", "-p", "do not prepend hierarchy to tag (default)"); - - cmd.addSubGroup("writing:"); - cmd.addOption("--write-pixel", "+W", 1, "[d]irectory : string", - "write pixel data to a .raw file stored in d\n(little endian, filename created automatically)"); - - /* evaluate command line */ - prepareCmdLineArgs(argc, argv, OFFIS_CONSOLE_APPLICATION); - if (app.parseCommandLine(cmd, argc, argv, OFCommandLine::ExpandWildcards)) - { - /* check exclusive options first */ - - if (cmd.getParamCount() == 0) - { - if (cmd.findOption("--version")) - { - app.printHeader(OFTrue /*print host identifier*/); // uses ofConsole.lockCerr() - CERR << endl << "External libraries used:"; -#ifdef WITH_ZLIB - CERR << endl << "- ZLIB, Version " << zlibVersion() << endl; -#else - CERR << " none" << endl; -#endif - return 0; - } - } - - /* options */ - -#ifdef USE_EXPERIMENTAL_QUIET_MODE - if (cmd.findOption("--quiet")) - { - // tbd: disable ofConsole output! - app.setQuietMode(); - } -#endif - if (cmd.findOption("--debug")) opt_debugMode = 5; - - cmd.beginOptionBlock(); - if (cmd.findOption("--read-file")) readMode = ERM_autoDetect; - if (cmd.findOption("--read-file-only")) readMode = ERM_fileOnly; - if (cmd.findOption("--read-dataset")) readMode = ERM_dataset; - cmd.endOptionBlock(); - - cmd.beginOptionBlock(); - if (cmd.findOption("--read-xfer-auto")) - xfer = EXS_Unknown; - if (cmd.findOption("--read-xfer-detect")) - dcmAutoDetectDatasetXfer.set(OFTrue); - if (cmd.findOption("--read-xfer-little")) - { - app.checkDependence("--read-xfer-little", "--read-dataset", readMode == ERM_dataset); - xfer = EXS_LittleEndianExplicit; - } - if (cmd.findOption("--read-xfer-big")) - { - app.checkDependence("--read-xfer-big", "--read-dataset", readMode == ERM_dataset); - xfer = EXS_BigEndianExplicit; - } - if (cmd.findOption("--read-xfer-implicit")) - { - app.checkDependence("--read-xfer-implicit", "--read-dataset", readMode == ERM_dataset); - xfer = EXS_LittleEndianImplicit; - } - cmd.endOptionBlock(); - - cmd.beginOptionBlock(); - if (cmd.findOption("--accept-odd-length")) - { - dcmAcceptOddAttributeLength.set(OFTrue); - } - if (cmd.findOption("--assume-even-length")) - { - dcmAcceptOddAttributeLength.set(OFFalse); - } - cmd.endOptionBlock(); - - cmd.beginOptionBlock(); - if (cmd.findOption("--enable-cp246")) - { - dcmEnableCP246Support.set(OFTrue); - } - if (cmd.findOption("--disable-cp246")) - { - dcmEnableCP246Support.set(OFFalse); - } - cmd.endOptionBlock(); - - cmd.beginOptionBlock(); - if (cmd.findOption("--retain-un")) - { - dcmEnableUnknownVRConversion.set(OFFalse); - } - if (cmd.findOption("--convert-un")) - { - dcmEnableUnknownVRConversion.set(OFTrue); - } - cmd.endOptionBlock(); - - cmd.beginOptionBlock(); - if (cmd.findOption("--enable-correction")) - { - dcmEnableAutomaticInputDataCorrection.set(OFTrue); - } - if (cmd.findOption("--disable-correction")) - { - dcmEnableAutomaticInputDataCorrection.set(OFFalse); - } - cmd.endOptionBlock(); - -#ifdef WITH_ZLIB - cmd.beginOptionBlock(); - if (cmd.findOption("--bitstream-deflated")) - { - dcmZlibExpectRFC1950Encoding.set(OFFalse); - } - if (cmd.findOption("--bitstream-zlib")) - { - dcmZlibExpectRFC1950Encoding.set(OFTrue); - } - cmd.endOptionBlock(); -#endif - - if (cmd.findOption("--max-read-length")) - { - app.checkValue(cmd.getValueAndCheckMinMax(maxReadLength, 4, 4194302)); - maxReadLength *= 1024; // convert kbytes to bytes - } - cmd.beginOptionBlock(); - if (cmd.findOption("--load-all")) loadIntoMemory = OFTrue; - if (cmd.findOption("--load-short")) loadIntoMemory = OFFalse; - cmd.endOptionBlock(); - - cmd.beginOptionBlock(); - if (cmd.findOption("--print-all")) printFlags &= ~DCMTypes::PF_shortenLongTagValues; - if (cmd.findOption("--print-short")) printFlags |= DCMTypes::PF_shortenLongTagValues; - cmd.endOptionBlock(); - - if (cmd.findOption("--print-filename")) - printFilename = OFTrue; - - cmd.beginOptionBlock(); - if (cmd.findOption("--stop-on-error")) stopOnErrors = OFTrue; - if (cmd.findOption("--ignore-errors")) stopOnErrors = OFFalse; - cmd.endOptionBlock(); - - if (cmd.findOption("--search", 0, OFCommandLine::FOM_First)) - { - do - { - app.checkValue(cmd.getValue(current)); - if (!addPrintTagName(current)) return 1; - } while (cmd.findOption("--search", 0, OFCommandLine::FOM_Next)); - } - - cmd.beginOptionBlock(); - if (cmd.findOption("--search-all")) - { - app.checkDependence("--search-all", "--search", printTagCount>0); - printAllInstances = OFTrue; - } - if (cmd.findOption("--search-first")) - { - app.checkDependence("--search-first", "--search", printTagCount>0); - printAllInstances = OFFalse; - } - cmd.endOptionBlock(); - - cmd.beginOptionBlock(); - if (cmd.findOption("--prepend")) - { - app.checkDependence("--prepend", "--search", printTagCount>0); - prependSequenceHierarchy = OFTrue; - } - if (cmd.findOption("--no-prepend")) - { - app.checkDependence("--no-prepend", "--search", printTagCount>0); - prependSequenceHierarchy = OFFalse; - } - cmd.endOptionBlock(); - - if (cmd.findOption("--write-pixel")) - { - app.checkValue(cmd.getValue(pixelDirectory)); - writePixelData = OFTrue; - } - } - - SetDebugLevel((opt_debugMode)); - - /* make sure data dictionary is loaded */ - if (!dcmDataDict.isDictionaryLoaded()) - { - CERR << "Warning: no data dictionary loaded, " - << "check environment variable: " - << DCM_DICT_ENVIRONMENT_VARIABLE; - } - - int errorCount = 0; - int count = cmd.getParamCount(); - for (int i=1; i<=count; i++) - { - cmd.getParam(i, current); - if (printFilename) - { - /* a newline separates two consecutive "dumps" */ - if (i > 1) - COUT << endl; - /* print header with filename */ - COUT << "# " << OFFIS_CONSOLE_APPLICATION << " (" << i << "/" << count << "): " << current << endl; - } - errorCount += dumpFile(COUT, current, readMode, xfer, printFlags, loadIntoMemory, stopOnErrors, - writePixelData, pixelDirectory); - } - - return errorCount; -} -static void printResult(ostream& out, DcmStack& stack, size_t printFlags) + template <typename I> + static void printResult(Image& ima, DcmStack& stack, size_t printFlags) { unsigned long n = stack.card(); - if (n == 0) { + if (n == 0) return; - } if (prependSequenceHierarchy) { /* print the path leading up to the top stack elem */ @@ -446,17 +228,18 @@ sprintf(buf, "(%x,%x).", OFstatic_cast(unsigned, dobj->getGTag()), OFstatic_cast(unsigned, dobj->getETag())); - out << buf; + std::cout << buf; } } } /* print the tag and its value */ DcmObject *dobj = stack.top(); - dobj->print(out, printFlags); + dobj->print(std::cout, printFlags); } -static int dumpFile(ostream & out, + template <typename I> + static int dumpFile(Image& ima, const char *ifname, const E_FileReadMode readMode, const E_TransferSyntax xfer, @@ -468,16 +251,11 @@ { int result = 0; - if ((ifname == NULL) || (strlen(ifname) == 0)) - { - CERR << OFFIS_CONSOLE_APPLICATION << ": invalid filename: <empty string>" << endl; - return 1; - } - DcmFileFormat dfile; DcmObject *dset = &dfile; - if (readMode == ERM_dataset) dset = dfile.getDataset(); + if (readMode == ERM_dataset) + dset = dfile.getDataset(); // Load file @@ -487,14 +265,15 @@ if (! cond.good()) { - CERR << OFFIS_CONSOLE_APPLICATION << ": error: " << dfile.error().text() + std::cerr << "error: " << dfile.error().text() << ": reading file: "<< ifname << endl; - result = 1; - if (stopOnErrors) return result; + if (stopOnErrors) + return result; } - if (loadIntoMemory) dfile.loadAllDataIntoMemory(); + if (loadIntoMemory) + dfile.loadAllDataIntoMemory(); if (printTagCount == 0) { @@ -510,10 +289,13 @@ else rname += str.substr(pos + 1); size_t counter = 0; - dset->print(out, printFlags, 0 /*level*/, rname.c_str(), &counter); - } else - dset->print(out, printFlags); - } else { + dset->print(std::cout, printFlags, 0 /*level*/, rname.c_str(), &counter); + } + else + dset->print(std::cout, printFlags); + } + else + { /* only print specified tags */ for (int i=0; i<printTagCount; i++) { @@ -521,23 +303,29 @@ unsigned int elem = 0xffff; DcmTagKey searchKey; const char* tagName = printTagNames[i]; - if (printTagKeys[i]) searchKey = *printTagKeys[i]; - else if (sscanf( tagName, "%x,%x", &group, &elem ) == 2 ) searchKey.set(group, elem); - else { - CERR << "Internal ERROR in File " << __FILE__ << ", Line " - << __LINE__ << endl - << "-- Named tag inconsistency" << endl; + if (printTagKeys[i]) + searchKey = *printTagKeys[i]; + else + { + if (sscanf(tagName, "%x,%x", &group, &elem) == 2) + searchKey.set(group, elem); + else + { + std::cerr << "Internal ERROR in File " << __FILE__ << ", Line " + << __LINE__ << std::endl + << "-- Named tag inconsistency" << std::endl; abort(); } + } DcmStack stack; if (dset->search(searchKey, stack, ESM_fromHere, OFTrue) == EC_Normal) { - printResult(out, stack, printFlags); + printResult(ima, stack, printFlags); if (printAllInstances) { while (dset->search(searchKey, stack, ESM_afterStackTop, OFTrue) == EC_Normal) - printResult(out, stack, printFlags); + printResult(ima, stack, printFlags); } } } @@ -545,3 +333,56 @@ return result; } + + + template <typename I> + inline + void load(Image& ima, + const std::string& filename) + { + trace::entering("mln::io::dicom::load"); + + std::ifstream file(filename.c_str()); + if (! file) + { + std::cerr << "error: cannot open file '" << filename << "'!"; + abort(); + } + + int opt_debugMode = 0; + OFBool loadIntoMemory = OFTrue; + size_t printFlags = DCMTypes::PF_shortenLongTagValues /*| DCMTypes::PF_showTreeStructure*/; + OFBool printFilename = OFFalse; + OFBool writePixelData = OFFalse; + E_FileReadMode readMode = ERM_autoDetect; + E_TransferSyntax xfer = EXS_Unknown; + OFBool stopOnErrors = OFTrue; + const char *current = NULL; + const char *pixelDirectory = NULL; + + + /* make sure data dictionary is loaded */ + if (!dcmDataDict.isDictionaryLoaded()) + { + std::cerr << "Warning: no data dictionary loaded, " + << "check environment variable: " + << DCM_DICT_ENVIRONMENT_VARIABLE; + } + + int errorCount = 0; + dumpFile(ima, current, readMode, xfer, printFlags, loadIntoMemory, stopOnErrors, + writePixelData, pixelDirectory); + + trace::exiting("mln::io::dicom::load"); + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::io::dicom + + } // end of namespace mln::io + +} // end of namespace mln + + +#endif // ! MLN_IO_DICOM_LOAD_HH Index: trunk/milena/sandbox/fabien/magick/magick.cc =================================================================== --- trunk/milena/sandbox/fabien/magick/magick.cc (revision 0) +++ trunk/milena/sandbox/fabien/magick/magick.cc (revision 3381) @@ -0,0 +1,14 @@ +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> + +#include "load.hh" + +int main() +{ + using namespace mln; + using value::int_u8; + + image2d<int_u8> lena; + + io::magick::load(lena, "/Users/HiSoKa/Work/LRDE/Olena/resources/CardiacCT/178562160.dcm"); +} Index: trunk/milena/sandbox/fabien/magick/Makefile =================================================================== --- trunk/milena/sandbox/fabien/magick/Makefile (revision 0) +++ trunk/milena/sandbox/fabien/magick/Makefile (revision 3381) @@ -0,0 +1,2 @@ +all: magick.cc + g++ -I../../../ `Magick++-config --cppflags --cxxflags --ldflags --libs` magick.cc -o mln_magick Index: trunk/milena/sandbox/fabien/magick/load.hh =================================================================== --- trunk/milena/sandbox/fabien/magick/load.hh (revision 0) +++ trunk/milena/sandbox/fabien/magick/load.hh (revision 3381) @@ -0,0 +1,142 @@ +// Copyright (C) 2009 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_IO_MAGICK_LOAD_HH +# define MLN_IO_MAGICK_LOAD_HH + +/*! + * \file mln/io/magick/load.hh + * + * \brief Define a function which loads an image of kind magick with + * given path. + * + */ + +# include <mln/core/image/image2d.hh> +# include <Magick++.h> + + +namespace mln +{ + + namespace io + { + + namespace magick + { + + /*! Load a magick image in a Milena image. + * + * \param[out] ima A reference to the image which will receive + * data. + * \param[in] filename The source. + */ + template <typename I> + void load(Image& ima, + const std::string& filename); + + /*! Load a magick image in a Milena image. To use this routine, you + * should specialize the template whith the value type of the + * image loaded. (ex : load<value::int_u8>("...") ) + * + * \param[in] filename The image source. + * + * \return An image2d which contains loaded data. + */ + template <typename V> + image2d<V> load(const std::string& filename); + + /*! Load a magick image in a Milena image. To use this routine, you + * should specialize the template whith the value type of the + * image loaded. (ex : load<value::int_u8>("...") ) + * + * \param[in] filename The image source. + * + * \return An image2d which contains loaded data. + */ + template <typename V> + image3d<V> load(const std::string& filename); + +# ifndef MLN_INCLUDE_ONLY + + template <typename V> + inline + image2d<V> load(const std::string& filename) + { + trace::entering("mln::io::magick::load"); + image2d<V> ima;// = io::pnm::load<V>(MAGICK, filename); + trace::exiting("mln::io::magick::load"); + return ima; + } + + template <typename V> + inline + image3d<V> load(const std::string& filename) + { + trace::entering("mln::io::magick::load"); + image2d<V> ima;// = io::pnm::load<V>(MAGICK, filename); + trace::exiting("mln::io::magick::load"); + return ima; + } + + + template <typename I> + inline + void load(Image& ima, + const std::string& filename) + { + trace::entering("mln::io::magick::load"); + + //std::ifstream file(filename.c_str()); + //if (! file) + //{ + // std::cerr << "error: cannot open file '" << filename << "'!"; + // abort(); + //} + + Magick::Image file(filename); + //std::cout << "file attribute: " << file.attribute() << std::endl; + std::cout << "width: " << file.columns() << std::endl; + std::cout << "height: " << file.rows() << std::endl; + std::cout << "x resolution: " << file.xResolution() << std::endl; + std::cout << "y resolution: " << file.yResolution() << std::endl; + std::cout << "comment: " << file.comment() << std::endl; + std::cout << "format: " << file.format() << std::endl; + + trace::exiting("mln::io::magick::load"); + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::io::magick + + } // end of namespace mln::io + +} // end of namespace mln + + +#endif // ! MLN_IO_MAGICK_LOAD_HH

15 years, 10 months

3380: Update Igr's code.

by Guillaume Lazzara

* lazzara/irm/hsl_grad_and_wst.cc, * lazzara/irm/irm_seg_with_mm_and_rag.cc, * lazzara/irm/wst_rag.cc, * lazzara/irm/wst_rag_hsl.cc: revamp and make use of make::graph. --- milena/sandbox/ChangeLog | 10 ++ milena/sandbox/lazzara/irm/hsl_grad_and_wst.cc | 79 +++++++-- .../sandbox/lazzara/irm/irm_seg_with_mm_and_rag.cc | 93 +++++++++++ milena/sandbox/lazzara/irm/wst_rag.cc | 100 +++--------- milena/sandbox/lazzara/irm/wst_rag_hsl.cc | 170 ++++++++------------ 5 files changed, 262 insertions(+), 190 deletions(-) create mode 100644 milena/sandbox/lazzara/irm/irm_seg_with_mm_and_rag.cc diff --git a/milena/sandbox/ChangeLog b/milena/sandbox/ChangeLog index ec838b2..79e7d3c 100644 --- a/milena/sandbox/ChangeLog +++ b/milena/sandbox/ChangeLog @@ -1,3 +1,13 @@ +2009-02-16 Guillaume Lazzara <z(a)lrde.epita.fr> + + Update Igr's code. + + * lazzara/irm/hsl_grad_and_wst.cc, + * lazzara/irm/irm_seg_with_mm_and_rag.cc, + * lazzara/irm/wst_rag.cc, + * lazzara/irm/wst_rag_hsl.cc: revamp and make use of make::graph. + + 2009-02-13 Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Finalize ICIP 2009 code. diff --git a/milena/sandbox/lazzara/irm/hsl_grad_and_wst.cc b/milena/sandbox/lazzara/irm/hsl_grad_and_wst.cc index f7b0335..ac29dd6 100644 --- a/milena/sandbox/lazzara/irm/hsl_grad_and_wst.cc +++ b/milena/sandbox/lazzara/irm/hsl_grad_and_wst.cc @@ -1,18 +1,22 @@ #include <mln/essential/2d.hh> -#include <mln/value/hsl.hh> + #include <mln/core/image/violent_cast_image.hh> + #include <mln/extract/all.hh> + #include <mln/value/label_16.hh> #include <mln/value/int_u8.hh> +#include <mln/value/hsl.hh> + #include <mln/morpho/closing_area.hh> +#include <mln/morpho/closing_volume.hh> + #include <mln/io/dump/save.hh> +#include <mln/io/ppm/all.hh> -mln::value::int_u8 foo(unsigned u) -{ - return u == 0 ? - 0 : // wshed line - 1 + (u - 1) % 255; // basin -} +#include <mln/labeling/mean_values.hh> + +#include <mln/fun/l2l/wrap.hh> namespace mln @@ -33,9 +37,7 @@ namespace mln } - - -int main(int, char *argv[]) +int main(int argc, char *argv[]) { using namespace mln; @@ -48,16 +50,48 @@ int main(int, char *argv[]) typedef image2d<rgb8> I; typedef image2d<hsl_f> J; + if (argc < 3) + { + std::cout << "Usage: " << argv[0] << " <ima.ppm> <closure_lambda>" + << std::endl; + } + I input_; io::ppm::load(input_, argv[1]); J input = level::convert(hsl_f(), input_); + io::ppm::save(level::transform(input, hsl2rgb()), "input-hsl.ppm"); - std::cout << "input_rgb(p) = " << input_(point2d(310,0)) << " - input_hsl(p) = " << input(point2d(310,0)) - << " hsl2rgb(p) = " << hsl2rgb()(input(point2d(310,0))) - << std::endl; - std::cout << "input_rgb(p) = " << input_(point2d(311,0)) << " - input_hsl(p) = " << input(point2d(311,0)) - << " hsl2rgb(p) = " << hsl2rgb()(input(point2d(311,0))) - << std::endl; + { + I tmp = level::transform(input, hsl2rgb()); + + image2d<int_u8> grad_hue; + initialize(grad_hue, input); + data::fill(grad_hue, extract::blue(tmp)); + image2d<int_u8> grad_sat; + initialize(grad_sat, input); + data::fill(grad_sat, extract::green(tmp)); + image2d<int_u8> grad_lum; + initialize(grad_lum, input); + data::fill(grad_lum, extract::blue(tmp)); + + + grad_hue = morpho::gradient(grad_hue, win_c4p()); + io::pgm::save(grad_hue, "hsl2rgb_grad_hue.pgm"); + + grad_sat = morpho::gradient(grad_sat, win_c4p()); + io::pgm::save(grad_sat, "hsl2rgb_grad_sat.pgm"); + + grad_lum = morpho::gradient(grad_lum, win_c4p()); + io::pgm::save(grad_lum, "hsl2rgb_grad_lum.pgm"); + + image2d<value::int_u8> grad; + initialize(grad, grad_hue); + mln_piter_(image2d<value::int_u8>) p(grad.domain()); + for_all(p) + grad(p) = math::max(grad_hue(p), math::max(grad_sat(p), grad_lum(p))); + io::pgm::save(grad, "hsl2rgb_grad_rgb.pgm"); + + } /// Compute the gradient on each component of the HSL image. /// Then keep the max values in the three images. @@ -70,9 +104,15 @@ int main(int, char *argv[]) image2d<float> grad_lum; initialize(grad_lum, input); data::fill(grad_lum, extract::lum(input)); + grad_hue = morpho::gradient(grad_hue, win_c4p()); +// io::pgm::save(grad_hue, "hsl2rgb_grad_hue.pgm"); + grad_sat = morpho::gradient(grad_sat, win_c4p()); +// io::pgm::save(grad_sat, "hsl2rgb_grad_sat.pgm"); + grad_lum = morpho::gradient(grad_lum, win_c4p()); +// io::pgm::save(grad_lum, "hsl2rgb_grad_lum.pgm"); image2d<value::int_u8> grad; initialize(grad, grad_hue); @@ -82,14 +122,19 @@ int main(int, char *argv[]) io::pgm::save(grad, "hsl2rgb_grad.pgm"); image2d<int_u8> clo = morpho::closing_area(grad, c4(), 10); +// image2d<int_u8> clo = morpho::closing_volume(grad, c4(), atoi(argv[2])); io::pgm::save(clo, "hsl2rgb_clo_a100.pgm"); label_16 nbasins; image2d<label_16> wshed = morpho::meyer_wst(clo, c4(), nbasins); io::dump::save(wshed, "hsl2rgb_wshed.dump"); std::cout << "nbasins = " << nbasins << std::endl; - io::pgm::save(level::transform(wshed, convert::to_fun(foo)), + + io::pgm::save(level::transform(wshed, fun::l2l::wrap<int_u8>()), "hsl2rgb_wshed.pgm"); + io::ppm::save(labeling::mean_values(level::transform(input, hsl2rgb()), wshed, nbasins), + "hsl2rgb_wshed_mean_colors.ppm"); + I out = level::transform(input, hsl2rgb()); io::ppm::save(out, "hsl2rgb_out.ppm"); diff --git a/milena/sandbox/lazzara/irm/irm_seg_with_mm_and_rag.cc b/milena/sandbox/lazzara/irm/irm_seg_with_mm_and_rag.cc new file mode 100644 index 0000000..9647986 --- /dev/null +++ b/milena/sandbox/lazzara/irm/irm_seg_with_mm_and_rag.cc @@ -0,0 +1,93 @@ +#include <iostream> +#include <mln/core/image/image2d.hh> + +#include <mln/core/alias/neighb2d.hh> +#include <mln/core/alias/window2d.hh> +#include <mln/core/image/image_if.hh> + +#include <mln/io/ppm/save.hh> +#include <mln/io/ppm/load.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/io/dump/save.hh> + +#include <mln/value/rgb8.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/label_16.hh> + +#include <mln/level/transform.hh> + +#include <mln/labeling/mean_values.hh> + +#include <mln/convert/to_fun.hh> + +#include <mln/morpho/gradient.hh> +#include <mln/morpho/closing_area.hh> +#include <mln/morpho/meyer_wst.hh> + +#include <mln/fun/l2l/wrap.hh> + +#include <mln/core/var.hh> +#include <mln/morpho/elementary/dilation.hh> + +#include <mln/core/routine/extend.hh> + +mln::value::int_u8 foo(unsigned u) +{ + return u == 0 ? + 0 : // wshed line + 1 + (u - 1) % 255; // basin +} + + +int main(int argc, char *argv[]) +{ + using namespace mln; + using value::int_u8; + using value::rgb8; + using value::label_16; + + if (argc < 4) + { + std::cout << "Usage: " << argv[0] << " <ima.pgm> <ima.ppm> <closure_lambda>" + << std::endl; + return 1; + } + + image2d<int_u8> irm; + io::pgm::load(irm, argv[1]); + + image2d<rgb8> irm_color; + io::ppm::load(irm_color, argv[2]); + + + image2d<int_u8> grad = morpho::gradient(irm, win_c4p()); + io::pgm::save(grad, "tmp_grad_c4p.pgm"); + +// image2d<int_u8> clo = morpho::closing_area(grad, c4(), atoi(argv[3])); +// io::pgm::save(clo, "tmp_clo_a100.pgm"); + + label_16 nbasins; + image2d<label_16> wshed = morpho::meyer_wst(grad, c4(), nbasins); + std::cout << "nbasins = " << nbasins << std::endl; + + io::pgm::save(level::transform(wshed, fun::l2l::wrap<int_u8>()), + "tmp_wshed.pgm"); + + + mln_VAR(vol2_, + morpho::elementary::dilation(extend(wshed | (pw::value(wshed) == 0u), + wshed), + c8())); + + data::fill((wshed | (pw::value(wshed) == 0u)).rw(), vol2_); + + io::ppm::save(labeling::mean_values(irm_color, wshed, nbasins), + "tmp_wshed_mean_colors.ppm"); + + data::fill((irm_color | (pw::value(wshed) == 0u)).rw(), literal::yellow); + + io::ppm::save(irm_color, "tmp_wshed_color.ppm"); + io::dump::save(wshed, "tmp_wshed.dump"); + +} diff --git a/milena/sandbox/lazzara/irm/wst_rag.cc b/milena/sandbox/lazzara/irm/wst_rag.cc index e397ba7..b88b8ec 100644 --- a/milena/sandbox/lazzara/irm/wst_rag.cc +++ b/milena/sandbox/lazzara/irm/wst_rag.cc @@ -17,9 +17,9 @@ #include <mln/core/image/line_graph_elt_neighborhood.hh> #include <mln/morpho/elementary/dilation.hh> #include <mln/labeling/mean_values.hh> - +#include <mln/extension/adjust_fill.hh> #include <mln/extract/all.hh> - +#include <mln/make/region_adjacency_graph.hh> // Given a color image and a wshed image, computes the component graph. // Vertex values are computed thanks to a RGB image. @@ -27,71 +27,20 @@ namespace mln { - - template <typename I, typename N> - util::graph - make_graph(const I& w, const N& nbh, const mln_value(I)& nbasins) + template <typename V> + value::int_u8 dist(const V& c1, const V& c2) { + unsigned d = math::diff_abs(c1.red(), c2.red()); + unsigned d_; + d_ = math::diff_abs(c1.green(), c2.green()); - image2d<bool> adj(box2d(nbasins.next(), nbasins.next())); - data::fill(adj, false); - - typedef mln_value(I) L; - L l1, l2; - mln_piter(I) p(w.domain()); - mln_niter(N) n(nbh, p); - for_all(p) - { - if (w(p) != 0u) - continue; - // p is in the watershed line. - l1 = l2 = 0; - for_all(n) - if (w.has(n) && w(n) != 0u) - { - if (l1 == 0u) // First label to be stored. - l1 = w(n); - else - if (w(n) != l1) // Useless: && l2 == 0) - { // Second label to be stored. - mln_invariant(l2 == 0u); - l2 = w(n); - break; - } - } - if (l2 == 0u || l1 == 0u) - continue; - if (l2 < l1) - std::swap(l1, l2); - - // adjacency l1 l2 - adj(point2d(l1,l2)) = true; - adj(point2d(l2,l1)) = true; - } + if (d_ > d) + d = d_; - // Construct graph. - util::graph g; - g.add_vertices(nbasins.next()); - for (unsigned i = 1; i < geom::nrows(adj); ++i) - for (unsigned j = 1; j < i; ++j) - if (adj(point2d(i,j))) - g.add_edge(i, j); + d_ = math::diff_abs(c1.blue(), c2.blue()); - return g; - } - - - - - template <typename V> - value::int_u8 dist(const V& c1, const V& c2) - { - unsigned d = 0; - d += (math::diff_abs(c1.red(), c2.red()) + 2) / 3; - d += (math::diff_abs(c1.green(), c2.green()) + 2) / 3; - d += (math::diff_abs(c1.blue(), c2.blue()) + 2) / 3; - if (d > 255) - d = 255; + if (d_ > d) + d = d_; return d; } @@ -218,12 +167,12 @@ namespace mln image2d<mln_value(I)> make_debug_graph_image(const I& input, const V& ima_v, const E& ima_e, - unsigned box_size) + unsigned box_size, const value::rgb8& bg) { image2d<mln_value(I)> ima; initialize(ima, input); - data::fill(ima, literal::white); + data::fill(ima, bg); debug::draw_graph(ima, ima_v.domain(), pw::cst(mln_value(I)(literal::green)), edge_to_color<E, mln_value(I)>(ima_e)); @@ -258,7 +207,7 @@ int main(int argc, char *argv[]) typedef image2d<rgb8> I; typedef image2d<label_16> J; - if (argc < 5) + if (argc < 6) { std::cout << argv[0] << " <input.ppm> <wsh-label16.dump> <nbasins> <box_size> <dist_max>" << std::endl << std::endl; @@ -282,15 +231,17 @@ int main(int argc, char *argv[]) label_16 nbasins = atoi(argv[3]); std::cout << "nbasins = " << nbasins << std::endl; - util::graph g = make_graph(vol, c4(), nbasins); + util::graph g = make::graph(vol, c4(), nbasins); // Compute value distances with a RGB image. mln_VAR(ima_v, make_vertex_graph_image(g, input, vol, nbasins)); mln_VAR(ima_e, make_edge_graph_image(ima_v, g)); //DEBUG -// io::ppm::save(make_debug_graph_image(input, ima_v, ima_e, box_size), -// "wst_rag_graph_image.ppm"); + io::ppm::save(make_debug_graph_image(input, ima_v, ima_e, box_size, literal::white), + "wst_rag_graph_image_white.ppm"); + io::ppm::save(make_debug_graph_image(input, ima_v, ima_e, box_size, literal::black), + "wst_rag_graph_image_black.ppm"); mln_piter_(ima_e_t) e(ima_e.domain()); @@ -313,16 +264,17 @@ int main(int argc, char *argv[]) for (unsigned i = 0; i < parent.nelements(); ++i) { unsigned p = find_root(parent, i); + mln_assertion(parent[p] == find_root(parent, i)); if (new_label[p] == 0) new_label[p] = nbasins2++; f(i) = new_label[p]; } - mln_invariant(f(0) == 0); + mln_invariant(f(0) == 0u); --nbasins2; //nbasins2 does not count the basin with label 0. std::cout << "nbasins2 = " << nbasins2 << std::endl; J vol2 = level::transform(vol, f); - util::graph g2 = make_graph(vol2, c4(), nbasins2); + util::graph g2 = make::graph(vol2, c4(), nbasins2); // Compute value distances with a RGB image. mln_VAR(ima_v2, make_vertex_graph_image(g2, input, vol2, nbasins2)); mln_VAR(ima_e2, make_edge_graph_image(ima_v2, g2)); @@ -336,6 +288,8 @@ int main(int argc, char *argv[]) io::ppm::save(labeling::mean_values(input, vol2, nbasins2), "wst_rag_mean_colors.ppm"); - io::ppm::save(make_debug_graph_image(input, ima_v2, ima_e2, box_size), - "wst_rag_graph_image2.ppm"); + io::ppm::save(make_debug_graph_image(input, ima_v2, ima_e2, box_size, literal::white), + "wst_rag_graph_image2_white.ppm"); + io::ppm::save(make_debug_graph_image(input, ima_v2, ima_e2, box_size, literal::black), + "wst_rag_graph_image2_black.ppm"); } diff --git a/milena/sandbox/lazzara/irm/wst_rag_hsl.cc b/milena/sandbox/lazzara/irm/wst_rag_hsl.cc index 156f1b3..48cca5a 100644 --- a/milena/sandbox/lazzara/irm/wst_rag_hsl.cc +++ b/milena/sandbox/lazzara/irm/wst_rag_hsl.cc @@ -1,6 +1,5 @@ #include <mln/essential/2d.hh> - #include <mln/core/alias/dpoint2d.hh> #include <mln/core/alias/vec3d.hh> #include <mln/core/image/line_graph_elt_neighborhood.hh> @@ -32,9 +31,10 @@ #include <mln/value/label_8.hh> #include <mln/value/rgb16.hh> #include <mln/value/rgb8.hh> - #include <mln/value/hsl.hh> +#include <mln/make/region_adjacency_graph.hh> + // Given a color image and a wshed image, computes the component graph. // Vertex values are computed thanks to a HSL image. @@ -44,71 +44,33 @@ namespace mln { - template <typename I, typename N> - util::graph - make_graph(const I& w, const N& nbh, const mln_value(I)& nbasins) + struct hsl2rgb : Function_v2v<hsl2rgb> { + typedef value::rgb8 result; - image2d<bool> adj(box2d(nbasins.next(), nbasins.next())); - data::fill(adj, false); - - typedef mln_value(I) L; - L l1, l2; - mln_piter(I) p(w.domain()); - mln_niter(N) n(nbh, p); - for_all(p) + value::rgb8 operator()(const value::hsl_f& v) const { - if (w(p) != 0u) - continue; - // p is in the watershed line. - l1 = l2 = 0; - for_all(n) - if (w.has(n) && w(n) != 0u) - { - if (l1 == 0u) // First label to be stored. - l1 = w(n); - else - if (w(n) != l1) // Useless: && l2 == 0) - { // Second label to be stored. - mln_invariant(l2 == 0u); - l2 = w(n); - break; - } - } - if (l2 == 0u || l1 == 0u) - continue; - if (l2 < l1) - std::swap(l1, l2); - - // adjacency l1 l2 - adj(point2d(l1,l2)) = true; - adj(point2d(l2,l1)) = true; + value::rgb8 res((v.hue() / 360.0f) * 255.0f, + v.sat() * 255.0f, + v.lum()); + return res; } - - // Construct graph. - util::graph g; - g.add_vertices(nbasins.next()); - for (unsigned i = 1; i < geom::nrows(adj); ++i) - for (unsigned j = 1; j < i; ++j) - if (adj(point2d(i,j))) - g.add_edge(i, j); - - return g; - } + }; - value::int_u8 dist(const value::hsl_f& c1, const value::hsl_f& c2) + value::int_u8 dist(const value::rgb8& c1, const value::rgb8& c2) { unsigned d = 0; - d += (math::diff_abs(c1.hue(), c2.hue()) + 2) / 3; - d += (math::diff_abs(c1.sat(), c2.sat()) + 2) / 3; - d += (math::diff_abs(c1.lum(), c2.lum()) + 2) / 3; + d += (math::diff_abs(c1.red(), c2.red()) + 2) / 3; + d += (math::diff_abs(c1.green(), c2.green()) + 2) / 3; + d += (math::diff_abs(c1.blue(), c2.blue()) + 2) / 3; if (d > 255) d = 255; return d; } + // ima_v, image on graph vertices; value = mean color per vertex (watershed basin) template <typename I> inline @@ -134,7 +96,7 @@ namespace mln mln_piter(ima_e_t) e(ima_e.domain()); for_all(e) // in ima_e ima_e(e) = dist(ima_v.function()(e.element().v1()), - ima_v.function()(e.element().v2())); + ima_v.function()(e.element().v2())); return ima_e; } @@ -147,15 +109,15 @@ namespace mln { // Cf. sandbox/theo/color/segment_rgb_pixels.cc - util::array<value::hsl_f> m_3f = labeling::compute(accu::mean<mln_value(I)>(), + util::array<vec3d_f> m_3f = labeling::compute(accu::mean<mln_value(I)>(), input, // input color image w, // watershed labeling nbasins); - m_3f[0] = value::hsl_f(0,0,0);; + m_3f[0] = literal::zero; util::array<mln_value(I)> m; convert::from_to(m_3f, m); - m[0] = value::hsl_f(360,0.5,255); + m[0] = literal::yellow; return m; } @@ -269,7 +231,7 @@ int main(int argc, char *argv[]) typedef image2d<hsl_f> I; typedef image2d<label_16> J; - if (argc < 5) + if (argc < 6) { std::cout << argv[0] << " <input.ppm> <wsh-label16.dump> <nbasins> <box_size> <dist_max>" << std::endl << std::endl; @@ -283,9 +245,10 @@ int main(int argc, char *argv[]) unsigned box_size = atoi(argv[4]); - K input_; - io::ppm::load(input_, argv[1]); - I input = level::convert(hsl_f(), input_); + K input; + io::ppm::load(input, argv[1]); + I input_ = level::convert(hsl_f(), input); + input = level::transform(input_, hsl2rgb()); J vol; io::dump::load(vol, argv[2]); @@ -293,46 +256,53 @@ int main(int argc, char *argv[]) label_16 nbasins = atoi(argv[3]); std::cout << "nbasins = " << nbasins << std::endl; - util::graph g = make_graph(vol, c4(), nbasins); + J vol2; + unsigned nbasins2; + fun::i2v::array<label_16> f; + { + util::graph g = make::graph(vol, c4(), nbasins); - mln_VAR(ima_v, make_vertex_graph_image(g, input, vol, nbasins)); - mln_VAR(ima_e, make_edge_graph_image(ima_v, g)); + mln_VAR(ima_v, make_vertex_graph_image(g, input, vol, nbasins)); + mln_VAR(ima_e, make_edge_graph_image(ima_v, g)); - //DEBUG - io::ppm::save(make_debug_graph_image(input, ima_v, ima_e, box_size), - "wst_rag_graph_image.ppm"); + //DEBUG + io::ppm::save(make_debug_graph_image(input, ima_v, ima_e, box_size), + "wst_rag_graph_image.ppm"); - mln_piter_(ima_e_t) e(ima_e.domain()); - util::array<label_16> parent(g.v_nmax()); - for (unsigned i = 0; i < parent.nelements(); ++i) - parent[i] = i; + mln_piter_(ima_e_t) e(ima_e.domain()); + util::array<label_16> parent(g.v_nmax()); + for (unsigned i = 0; i < parent.nelements(); ++i) + parent[i] = i; - for_all(e) - { - unsigned v1 = e.element().v1(); - unsigned v2 = e.element().v2(); - if (ima_e(e) <= (unsigned)(atoi(argv[5])) - && find_root(parent, v1) != find_root(parent, v2)) - parent[find_root(parent, v1)] = find_root(parent, v2); - } + for_all(e) + { + unsigned v1 = e.element().v1(); + unsigned v2 = e.element().v2(); + if (ima_e(e) <= (unsigned)(atoi(argv[5])) + && find_root(parent, v1) != find_root(parent, v2)) + parent[find_root(parent, v1)] = find_root(parent, v2); + } - fun::i2v::array<label_16> f(parent.nelements()); - std::vector<unsigned> new_label(parent.nelements(), 0); - unsigned nbasins2 = 0; - for (unsigned i = 0; i < parent.nelements(); ++i) - { - unsigned p = find_root(parent, i); - if (new_label[p] == 0) - new_label[p] = nbasins2++; - f(i) = new_label[p]; - } - mln_invariant(f(0) == 0); - --nbasins2; //nbasins2 does not count the basin with label 0. - std::cout << "nbasins2 = " << nbasins2 << std::endl; + f = fun::i2v::array<label_16>(parent.nelements()); + std::vector<unsigned> new_label(parent.nelements(), 0); + nbasins2 = 0; + for (unsigned i = 0; i < parent.nelements(); ++i) + { + unsigned p = find_root(parent, i); + if (new_label[p] == 0) + new_label[p] = nbasins2++; + f(i) = new_label[p]; + } + mln_invariant(f(0) == 0); + --nbasins2; //nbasins2 does not count the basin with label 0. - J vol2 = level::transform(vol, f); - util::graph g2 = make_graph(vol2, c4(), nbasins2); + std::cout << "nbasins2 = " << nbasins2 << std::endl; + } + // debug::println(vol); + vol2 = level::transform(vol, f); + // debug::println(vol2); + util::graph g2 = make::graph(vol2, c4(), nbasins2); // Compute values distance on the HSL Image. mln_VAR(ima_v2, make_vertex_graph_image(g2, input, vol2, nbasins2)); @@ -340,14 +310,14 @@ int main(int argc, char *argv[]) mln_VAR(vol2_, morpho::elementary::dilation(extend(vol2 | (pw::value(vol2) == 0u), - vol2), - c8())); + vol2), + c8())); data::fill((vol2 | (pw::value(vol2) == 0u)).rw(), vol2_); - io::ppm::save(labeling::mean_values(input_, vol2, nbasins2), - "wst_rag_mean_colors.ppm"); - io::ppm::save(make_debug_graph_image(input_, ima_v2, ima_e2, box_size), - "wst_rag_graph_image2.ppm"); + io::ppm::save(labeling::mean_values(input, vol2, nbasins2), + "wst_rag_mean_colors.ppm"); + io::ppm::save(make_debug_graph_image(input, ima_v2, ima_e2, box_size), + "wst_rag_graph_image2.ppm"); } -- 1.5.6.5

15 years, 10 months

3379: Revamp icp2.hh.

by Guillaume Lazzara

* registration/icp2.hh: move registration*... * registration/registration.hh: ... here. --- milena/ChangeLog | 7 + milena/mln/registration/icp2.hh | 301 ++++++--------------------- milena/mln/registration/registration.hh | 346 +++++++++++++++++++++++++++---- 3 files changed, 384 insertions(+), 270 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 2282bcd..682d896 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,12 @@ 2009-02-16 Guillaume Lazzara <z(a)lrde.epita.fr> + Revamp icp2.hh. + + * registration/icp2.hh: move registration*... + * registration/registration.hh: ... here. + +2009-02-16 Guillaume Lazzara <z(a)lrde.epita.fr> + Fix make::graph and add make::region_adjacency_graph. * mln/make/graph.hh: use an adjacency matrix and a neighborhood. diff --git a/milena/mln/registration/icp2.hh b/milena/mln/registration/icp2.hh index 916f8c0..de763ea 100644 --- a/milena/mln/registration/icp2.hh +++ b/milena/mln/registration/icp2.hh @@ -113,6 +113,7 @@ namespace mln const algebra::quat& initial_rot, const mln_vec(P)& initial_translation); + template <typename P, typename F> composed< translation<P::dim,float>,rotation<P::dim,float> > icp(const p_array& P_, @@ -120,40 +121,6 @@ namespace mln const F& closest_point); - //FIXME: move to registration.hh - /// Call ICP once and return the resulting transformation. - template <typename P, typename F> - inline - composed< translation<P::dim,float>,rotation<P::dim,float> > - registration1(const p_array& P_, - const p_array& X); - - //FIXME: move to registration.hh - /// Call ICP 10 times. - /// Do the first call to ICP with all sites then work on a subset of - /// which size is decreasing. - /// For each call, sites part of the subset which are too far or too - /// close are removed. - /// Removed sites are *NOT* reused later in the subset. - template <typename P> - inline - composed< translation<P::dim,float>,rotation<P::dim,float> > - registration2(const p_array& P_, - const p_array& X); - - //FIXME: move to registration.hh - /// Call ICP 10 times. - /// Do the first call to ICP with all sites then work on a subset. - /// For each call, a distance criterion is computed on a subset. - /// A new subset is computed from the whole set of points according - /// to this distance. It will be used in the next call. - /// Sites which have been removed can be reintegrated. - template <typename P> - inline - composed< translation<P::dim,float>,rotation<P::dim,float> > - registration3(const p_array& P_, - const p_array& X); - # ifndef MLN_INCLUDE_ONLY @@ -172,8 +139,21 @@ namespace mln box3d box = geom::bbox(X); box.enlarge(1, box.nrows()); box.enlarge(2, box.ncols()); + + mln_postcondition(box.is_valid()); + std::cout << "Map image defined on " << box << std::endl; + init(X, box); + } + + closest_point_with_map(const p_array& X, const box& box) + { + init(X, box); + } + + void init(const p_array& X, const box& box) + { typedef mln_ch_value(I, bool) model_t; model_t model(box); data::fill(model, false); @@ -187,30 +167,36 @@ namespace mln f); std::cout << "canvas::distance_geodesic - " << t << "s" << std::endl; - cp_ima_ = f.cp_ima; - -#ifndef NDEBUG - mln_ch_value(I, bool) debug2(box); - data::fill(debug2, false); - mln_ch_value(I, value::rgb8) debug(box); - mln_piter(p_array) p(X); - for_all(p) - { - debug(p) = debug::internal::random_color(value::rgb8()); - debug2(p) = true; - } - io::pbm::save(slice(debug2,0), "debug2-a.ppm"); + initialize(this->cp_ima_, f.cp_ima); + this->cp_ima_ = f.cp_ima; - mln_piter(I) pi(cp_ima_.domain()); - for_all(pi) - { - debug(pi) = debug(cp_ima_(pi)); - debug2(pi) = debug2(cp_ima_(pi)); - } + mln_postcondition(cp_ima_.is_valid()); + mln_postcondition(cp_ima_.domain().is_valid()); + std::cout << "pmin = " << cp_ima_.domain().pmin() << std::endl;; + std::cout << "pmax = " << cp_ima_.domain().pmax() << std::endl;; - io::pbm::save(slice(debug2,0), "debug2-b.ppm"); - io::ppm::save(slice(debug,0), "debug.ppm"); - std::cout << "map saved" << std::endl; +#ifndef NDEBUG +// mln_ch_value(I, bool) debug2(box); +// data::fill(debug2, false); +// mln_ch_value(I, value::rgb8) debug(box); +// mln_piter(p_array) p(X); +// for_all(p) +// { +// debug(p) = debug::internal::random_color(value::rgb8()); +// debug2(p) = true; +// } +// io::pbm::save(slice(debug2,0), "debug2-a.ppm"); +// +// mln_piter(I) pi(cp_ima_.domain()); +// for_all(pi) +// { +// debug(pi) = debug(cp_ima_(pi)); +// debug2(pi) = debug2(cp_ima_(pi)); +// } +// +// io::pbm::save(slice(debug2,0), "debug2-b.ppm"); +// io::ppm::save(slice(debug,0), "debug.ppm"); +// std::cout << "map saved" << std::endl; #endif } @@ -219,9 +205,11 @@ namespace mln return cp_ima_(p); } + // Distance map dmap_t dmap_X_; - private: + // Closest point image. cp_ima_t cp_ima_; + private: }; @@ -355,6 +343,7 @@ namespace mln } + template <typename P, typename F> void compute_distance_criteria(const p_array& P_, @@ -401,8 +390,10 @@ namespace mln p_array tmp; unsigned removed = 0; +# ifndef NDEBUG data::fill(out, literal::black); data::fill((out | X).rw(), literal::white); +# endif // ! NDEBUG mln_piter(p_array) p(P_); for_all(p) @@ -424,11 +415,19 @@ namespace mln } } + +# ifndef NDEBUG std::ostringstream ss2; ss2 << method << "_" << r << "_removed_sites" << ".ppm"; io::ppm::save(mln::slice(out,0), ss2.str()); std::cout << "Points removed: " << removed << std::endl; +# endif // ! NDEBUG + // They are used for debug purpose only. + // When NDEBUG is set, they are unused. + (void) X; + (void) r; + (void) method; return tmp; } @@ -441,8 +440,10 @@ namespace mln const std::pair<algebra::quat,mln_vec(P)>& pair, const std::string& period, const value::rgb8& c) { +# ifndef NDEBUG data::fill(out, literal::black); data::fill((out | X).rw(), literal::white); +# endif // !NDEBUG mln_piter(p_array) p1(P_); for_all(p1) @@ -458,10 +459,12 @@ namespace mln out(Pk_i) = c; } +# ifndef NDEBUG std::ostringstream ss; ss << prefix << "_" << r << "_" << period << ".ppm"; io::ppm::save(slice(out,0), ss.str()); +# endif // ! NDEBUG } @@ -487,6 +490,7 @@ namespace mln return accu.to_result(); } + /// FIXME: work only for 3d images. template <typename P, typename F> algebra::quat @@ -639,7 +643,7 @@ namespace mln for_all(p_dbg) tmp_(qR_old.rotate(p_dbg.to_vec()) + qT_old) = literal::green; std::ostringstream ss; - ss << "reg/tmp_0"; + ss << "tmp_0"; if (k < 10) ss << "0"; ss << k << ".ppm"; @@ -650,14 +654,16 @@ namespace mln std::cout << "d_" << k << "=" << d_k << std::endl; // Check distance and error according to the related paper. - mln_assertion(0 <= d_k); - mln_assertion(d_k <= e_k); - // Disabled because of the following 'if' +// mln_assertion(0 <= d_k); +// mln_assertion(d_k <= e_k); + // mln_assertion(e_k <= d_k_old); // mln_assertion(d_k_old <= e_k_old); - if (d_k > d_k_old) + // During the first runs, d_k may be higher than e_k. + // Hence, there we test k > 3. + if (k > 3 && (d_k > e_k || d_k > d_k_old || e_k > e_k_old)) { qR = qR_old; qT = qT_old; @@ -678,177 +684,6 @@ namespace mln } - - template <typename P, typename F> - inline - composed< translation<P::dim,float>,rotation<P::dim,float> > - registration1(const p_array& P_, - const p_array& X) - { - trace::entering("mln::registration::registration1"); - - util::timer t; - t.start(); - - registration::closest_point_with_map closest_point(X); - - std::pair<algebra::quat,mln_vec(P)> pair = icp(P_, X, closest_point, - algebra::quat(1,0,0,0), - literal::zero); - std::cout << "icp = " << t << std::endl; - - typedef rotation<3u,float> rot_t; - rot_t tqR(pair.first); - typedef translation<3u,float> trans_t; - trans_t tqT(pair.second); - composed<trans_t, rot_t> result(tqT, tqR); - - trace::exiting("mln::registration::registration1"); - - return result; - } - - - template <typename P> - inline - composed< translation<P::dim,float>,rotation<P::dim,float> > - registration2(const p_array& P_, - const p_array& X) - { - trace::entering("mln::registration::registration2"); - - std::string method = "compute_on_subsets"; - - registration::closest_point_with_map closest_point(X); - - util::timer t; - t.start(); - - // P_bak is shuffled. - p_array P_bak = P_; - - unsigned r = 0; - std::pair<algebra::quat,mln_vec(P)> pair; - pair.first = algebra::quat(1,0,0,0); - pair.second = literal::zero; - box3d box = geom::bbox(X); - box.enlarge(1, 60); - box.enlarge(2, 60); - image3d<value::rgb8> out(box); - p_array removed_set; - - do - { - std::cout << std::endl << std::endl << "==== New run - " << r << std::endl; - pair = icp(P_bak, X, closest_point, - pair.first, - pair.second); - - display_sites_used_in_icp(out, P_bak, P_, X, r, method, pair, "final", literal::blue); - - int d_min, d_max; - compute_distance_criteria(P_bak, closest_point, pair, r, d_min, d_max); - - P_bak = remove_too_far_sites(out, P_bak, - closest_point, pair, X, removed_set, - r, d_min, d_max, method); - -#ifndef NDEBUG - display_sites_used_in_icp(out, P_bak, P_, X, r, method, pair, "schanges", literal::green); -#endif - - ++r; - - std::cout << "==== End of run" << std::endl; - } while (r < 10); - std::cout << "icp = " << t << std::endl; - -#ifndef NDEBUG - draw_last_run(box, P_bak, removed_set, X, pair.first, pair.second); -#endif - - typedef rotation<3u,float> rot_t; - rot_t tqR(pair.first); - typedef translation<3u,float> trans_t; - trans_t tqT(pair.second); - composed<trans_t,rot_t> result(tqT, tqR); - - trace::exiting("mln::registration::registration2"); - - return result; - } - - - template <typename P> - inline - composed< translation<P::dim,float>,rotation<P::dim,float> > - registration3(const p_array& P_, - const p_array& X) - { - trace::entering("mln::registration::registration3"); - - registration::closest_point_with_map closest_point(X); - - std::string method = "compute_on_subset_and_p"; - - util::timer t; - t.start(); - - // P_bak is shuffled. - p_array P_bak = P_; - - unsigned r = 0; - std::pair<algebra::quat,mln_vec(P)> pair; - pair.first = algebra::quat(1,0,0,0); - pair.second = literal::zero; - box3d box = geom::bbox(X); - box.enlarge(1, 60); - box.enlarge(2, 60); - image3d<value::rgb8> out(box); - p_array removed_set; - - do - { - std::cout << std::endl << std::endl << "==== New run - " << r << std::endl; - pair = icp(P_bak, X, closest_point, - pair.first, - pair.second); - - display_sites_used_in_icp(out, P_bak, P_, X, r, method, pair, "final", literal::blue); - - int d_min, d_max; - compute_distance_criteria(P_bak, closest_point, pair, r, d_min, d_max); - - P_bak = remove_too_far_sites(out, P_, - closest_point, pair, X, removed_set, - r, d_min, d_max, method); - -#ifndef NDEBUG - display_sites_used_in_icp(out, P_bak, P_, X, r, method, pair, "schanges", literal::green); -#endif - - ++r; - - std::cout << "==== End of run" << std::endl; - } while (r < 10); - std::cout << "icp = " << t << std::endl; - -#ifndef NDEBUG - draw_last_run(box, P_bak, removed_set, X, pair.first, pair.second); -#endif - - typedef rotation<3u,float> rot_t; - rot_t tqR(pair.first); - typedef translation<3u,float> trans_t; - trans_t tqT(pair.second); - composed<trans_t,rot_t> result(tqT, tqR); - - trace::exiting("mln::registration::registration3"); - - return result; - } - - # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::registration diff --git a/milena/mln/registration/registration.hh b/milena/mln/registration/registration.hh index 7ad7b00..f96787a 100644 --- a/milena/mln/registration/registration.hh +++ b/milena/mln/registration/registration.hh @@ -1,4 +1,4 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory +// Copyright (C) 2008, 2009 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 @@ -30,12 +30,13 @@ /// \file mln/registration/registration.hh /// -/// image registration +/// Image registration +/// \sa registration::icp -# include <mln/registration/icp.hh> +# include <mln/core/image/image3d.hh> +# include <mln/registration/icp2.hh> # include <mln/fun/x2x/all.hh> # include <mln/fun/x2p/closest_point.hh> -# include <mln/core/image/lazy_image.hh> # include <mln/convert/to_p_array.hh> namespace mln @@ -46,63 +47,334 @@ namespace mln using namespace mln::fun::x2x; - /// Register an image \p cloud over the image \p surface. - template <typename I, typename J> + + //FIXME: move to registration.hh + /// Call ICP once and return the resulting transformation. + template <typename P> + inline + composed< translation<P::dim,float>,rotation<P::dim,float> > + registration1(const p_array& P_, + const p_array& X); + + //FIXME: move to registration.hh + /// Call ICP 10 times. + /// Do the first call to ICP with all sites then work on a subset of + /// which size is decreasing. + /// For each call, a distance criterion is computed on a subset. + /// Sites part of the subset which are too far or too + /// close are removed. + /// Removed sites are *NOT* reused later in the subset. + template <typename P> inline - composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > - registration(const Image& cloud, - const Image<J>& surface); + composed< translation<P::dim,float>,rotation<P::dim,float> > + registration2(const p_array& P_, + const p_array& X); + + //FIXME: move to registration.hh + /// Call ICP 10 times. + /// Do the first call to ICP with all sites then work on a subset. + /// For each call, a distance criterion is computed on a subset. + /// A new subset is computed from the whole set of points according + /// to this distance. It will be used in the next call. + /// Removed Sites *MAY* be reintegrated. + template <typename P> + inline + composed< translation<P::dim,float>,rotation<P::dim,float> > + registration3(const p_array& P_, + const p_array& X); + # ifndef MLN_INCLUDE_ONLY + + namespace internal + { + + template <typename P> + inline + void + registration_tests(const p_array& P_, const p_array& X) + { + mln_assertion(P_.is_valid()); + mln_assertion(X.is_valid()); + mln_assertion(!X.is_empty()); + mln_assertion(!P_.is_empty()); + + // FIXME: Work only in 3D for now... + mln_precondition(P::dim == 3); + (void) P_; + (void) X; + } + + } // end of namespace mln::registration::internal + + namespace impl { - template <typename I, typename J> + template <typename P> + inline + composed< translation<P::dim,float>,rotation<P::dim,float> > + registration1(const p_array& P_, + const p_array& X) + { + trace::entering("mln::registration::registration1"); + +# ifndef NDEBUG + util::timer t; + t.start(); +# endif // ! NDEBUG + + registration::closest_point_with_map closest_point(X); + + std::pair<algebra::quat,mln_vec(P)> pair = icp(P_, X, closest_point, + algebra::quat(1,0,0,0), + literal::zero); +# ifndef NDEBUG + std::cout << "icp = " << t << std::endl; +# endif // ! NDEBUG + + typedef rotation<3u,float> rot_t; + rot_t tqR(pair.first); + typedef translation<3u,float> trans_t; + trans_t tqT(pair.second); + composed<trans_t, rot_t> result(tqT, tqR); + + trace::exiting("mln::registration::registration1"); + + return result; + } + + + template <typename P> + inline + composed< translation<P::dim,float>,rotation<P::dim,float> > + registration2(const p_array& P_, + const p_array& X) + { + trace::entering("mln::registration::registration2"); + + // Used for debug. + std::string method = "registration2"; + + registration::closest_point_with_map closest_point(X); + +# ifndef NDEBUG + util::timer t; + t.start(); +# endif // ! NDEBUG + + // P_bak is shuffled. + p_array P_bak = P_; + + unsigned r = 0; + std::pair<algebra::quat,mln_vec(P)> pair; + pair.first = algebra::quat(1,0,0,0); + pair.second = literal::zero; + box3d box = geom::bbox(X); + box.enlarge(1, 60); + box.enlarge(2, 60); + + // Used for debug. + image3d<value::rgb8> out(box); + + p_array removed_set; + + do + { + +# ifndef NDEBUG + std::cout << std::endl << std::endl << "==== New run - " << r << std::endl; +# endif // ! NDEBUG + + pair = icp(P_bak, X, closest_point, + pair.first, + pair.second); + +# ifndef NDEBUG + display_sites_used_in_icp(out, P_bak, P_, X, r, method, pair, + "final", literal::blue); +# endif // ! NDEBUG + + int d_min, d_max; + compute_distance_criteria(P_bak, closest_point, pair, r, d_min, d_max); + + P_bak = remove_too_far_sites(out, P_bak, + closest_point, pair, X, removed_set, + r, d_min, d_max, method); + +# ifndef NDEBUG + display_sites_used_in_icp(out, P_bak, P_, X, r, method, pair, + "schanges", literal::green); + std::cout << "==== End of run" << std::endl; +# endif + + ++r; + + } while (r < 10); + +# ifndef NDEBUG + std::cout << "icp = " << t << std::endl; + draw_last_run(box, P_bak, removed_set, X, pair.first, pair.second); +# endif + + typedef rotation<3u,float> rot_t; + rot_t tqR(pair.first); + typedef translation<3u,float> trans_t; + trans_t tqT(pair.second); + composed<trans_t,rot_t> result(tqT, tqR); + + trace::exiting("mln::registration::registration2"); + + return result; + } + + + template <typename P> inline - composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > - registration_(const I& cloud, - const J& surface) + composed< translation<P::dim,float>,rotation<P::dim,float> > + registration3(const p_array& P_, + const p_array& X) { - //FIXME: Use convert::to< p_array<mln_psite(I)> >() - 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); + trace::entering("mln::registration::registration3"); + + registration::closest_point_with_map closest_point(X); + std::cout << " pmin and pmax: " << std::endl; + std::cout << closest_point.cp_ima_.domain().pmin() << std::endl; + std::cout << closest_point.cp_ima_.domain().pmax() << std::endl; + + // Used for debug. + std::string method = "registration3"; + +# ifndef NDEBUG + util::timer t; + t.start(); +# endif // ! NDEBUG + + // P_bak is shuffled. + p_array P_bak = P_; + + unsigned r = 0; + std::pair<algebra::quat,mln_vec(P)> pair; + pair.first = algebra::quat(1,0,0,0); + pair.second = literal::zero; + box3d box = geom::bbox(X); + box.enlarge(1, 60); + box.enlarge(2, 60); + + // Used for debug. + image3d<value::rgb8> out(box); + + p_array removed_set; + + do + { +# ifndef NDEBUG + std::cout << std::endl << std::endl << "==== New run - " + << r << std::endl; +# endif // ! NDEBUG + + pair = icp(P_bak, X, closest_point, + pair.first, + pair.second); + +# ifndef NDEBUG + display_sites_used_in_icp(out, P_bak, P_, X, r, method, pair, + "final", literal::blue); +# endif // ! NDEBUG + + int d_min, d_max; + compute_distance_criteria(P_bak, closest_point, pair, r, d_min, d_max); + + P_bak = remove_too_far_sites(out, P_, + closest_point, pair, X, removed_set, + r, d_min, d_max, method); + +# ifndef NDEBUG + display_sites_used_in_icp(out, P_bak, P_, X, r, method, pair, + "schanges", literal::green); + std::cout << "==== End of run" << std::endl; +# endif // ! NDEBUG + + ++r; - //init rigid transform qk - composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > qk; + } while (r < 10); - //working box - const box<mln_psite(I)> working_box = - larger_than(geom::bbox(c), geom::bbox(x)).to_larger(100); +# ifndef NDEBUG + std::cout << "icp = " << t << std::endl; + draw_last_run(box, P_bak, removed_set, X, pair.first, pair.second); +# endif // ! NDEBUG - //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)>, box<mln_psite(I)> > - map(fun, fun.domain()); + typedef rotation<3u,float> rot_t; + rot_t tqR(pair.first); + typedef translation<3u,float> trans_t; + trans_t tqT(pair.second); + composed<trans_t,rot_t> result(tqT, tqR); - //run registration - return registration::icp(c, map, 1e-3); + trace::exiting("mln::registration::registration3"); + return result; } + } // end of namespace mln::registration::impl + + + + // Facade + + template <typename P> + inline + composed< translation<P::dim,float>,rotation<P::dim,float> > + registration1(const p_array& cloud, + const p_array& surface) + { + trace::entering("registration::registration1"); + + registration_tests(cloud, surface); + + composed< translation<P::dim,float>, rotation<P::dim,float> > + qk = impl::registration1(cloud, surface); + + trace::exiting("registration::registration1"); + + return qk; } - template <typename I, typename J> + + template <typename P> + inline + composed< translation<P::dim,float>,rotation<P::dim,float> > + registration2(const p_array& cloud, + const p_array& surface) + { + trace::entering("registration::registration2"); + + registration_tests(cloud, surface); + + composed< translation<P::dim,float>, rotation<P::dim,float> > + qk = impl::registration2(cloud, surface); + + trace::exiting("registration::registration2"); + + return qk; + } + + + template <typename P> inline - composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > - registration(const Image& cloud, - const Image<J>& surface) + composed< translation<P::dim,float>,rotation<P::dim,float> > + registration3(const p_array& cloud, + const p_array& surface) { - trace::entering("registration::registration"); + trace::entering("registration::registration3"); - mln_precondition(I::psite::dim == J::psite::dim); - mln_precondition(I::psite::dim == 3 || J::psite::dim == 2); + internal::registration_tests(cloud, surface); - composed< rotation<I::psite::dim, float>, translation<I::psite::dim, float> > - qk = impl::registration_(exact(cloud), exact(surface)); + composed< translation<P::dim,float>, rotation<P::dim,float> > + qk = impl::registration3(cloud, surface); - trace::exiting("registration::registration"); + trace::exiting("registration::registration3"); return qk; } -- 1.5.6.5

15 years, 10 months

← Newer
1
...
8
9
10
11
12
13
14
...
27
Older →

Jump to page:

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

2006

2005

2004

Olena-patches February 2009