April 2008 - Olena-patches

1855: Speed up iterations on line graph neighborhoods.
by Roland Levillain 09 Apr '08

09 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Roland Levillain <roland(a)lrde.epita.fr> Speed up iterations on line graph neighborhoods. * mln/core/line_graph_neighborhood_piter.hh (mln::line_graph_neighborhood_fwd_piter::neighbors_t) (mln::line_graph_neighborhood_bkd_piter::neighbors_t): New typedefs. (mln::line_graph_neighborhood_fwd_piter::first_) (mln::line_graph_neighborhood_fwd_piter::step_) (mln::line_graph_neighborhood_fwd_piter::adjacent_or_equal_to_p_ref_) (mln::line_graph_neighborhood_bkd_piter::first_) (mln::line_graph_neighborhood_bkd_piter::step_) (mln::line_graph_neighborhood_bkd_piter::adjacent_or_equal_to_p_ref_): Remove methods. (mln::line_graph_neighborhood_fwd_piter::id_) (mln::line_graph_neighborhood_bkd_piter::id_): Remove attributes. (mln::line_graph_neighborhood_fwd_piter::saved_p_ref_) (mln::line_graph_neighborhood_fwd_piter::neighbors_) (mln::line_graph_neighborhood_fwd_piter::i_) (mln::line_graph_neighborhood_bkd_piter::saved_p_ref_) (mln::line_graph_neighborhood_bkd_piter::neighbors_) (mln::line_graph_neighborhood_bkd_piter::i_): New attributes. (line_graph_neighborhood_fwd_piter) (line_graph_neighborhood_bkd_piter): Adjust ctors. (mln::line_graph_neighborhood_fwd_piter::p_ref) (mln::line_graph_neighborhood_fwd_piter::plg) (mln::line_graph_neighborhood_fwd_piter::neighbors) (mln::line_graph_neighborhood_bkd_piter::p_ref) (mln::line_graph_neighborhood_bkd_piter::plg) (mln::line_graph_neighborhood_bkd_piter::neighbors): New accessors. (mln::line_graph_neighborhood_fwd_piter::is_valid) (mln::line_graph_neighborhood_fwd_piter::invalidate) (mln::line_graph_neighborhood_fwd_piter::start) (mln::line_graph_neighborhood_fwd_piter::next_) (mln::line_graph_neighborhood_bkd_piter::is_valid) (mln::line_graph_neighborhood_bkd_piter::invalidate) (mln::line_graph_neighborhood_bkd_piter::start) (mln::line_graph_neighborhood_bkd_piter::next_): Change the algorithms used in these routines, so that they use a list of computed neighbors (new member neighbors_) instead of iterating over the whole set of edges. (mln::line_graph_neighborhood_fwd_piter::update_) (mln::line_graph_neighborhood_bkd_piter::update_): Adjust. * mln/core/line_graph_elt_neighborhood.hh: Catch up with the new interface of piters on line graph neighborhoods. (mln::line_graph_elt_neighborhood::neighbors_t): New typedef. (mln::line_graph_elt_neighborhood::start) (mln::line_graph_elt_neighborhood::next_): Remove methods. (mln::line_graph_elt_neighborhood::compute_neighbors_): New method. line_graph_elt_neighborhood.hh | 58 ++++---- line_graph_neighborhood_piter.hh | 262 ++++++++++++++++++--------------------- 2 files changed, 158 insertions(+), 162 deletions(-) Index: mln/core/line_graph_neighborhood_piter.hh --- mln/core/line_graph_neighborhood_piter.hh (revision 1854) +++ mln/core/line_graph_neighborhood_piter.hh (working copy) @@ -42,6 +42,8 @@ - mln::line_graph_window_bkd_piter - mln::line_graph_neighborhood_bkd_piter. */ +# include <set> + # include <mln/core/concept/point_iterator.hh> # include <mln/core/p_line_graph.hh> # include <mln/core/line_graph_psite.hh> @@ -80,6 +82,9 @@ // FIXME: Dummy typedef. typedef void mesh; + // The type of the set of neighbors (adjacent edge ids). + typedef std::set<util::edge_id> neighbors_t; + public: /// Construction. /// \{ @@ -114,28 +119,31 @@ /// Convert the iterator into a line graph psite. operator psite() const; + /// Return the reference psite. + const psite& p_ref() const; + /// Return the mln::p_line_graph corresponding to this piter. + const p_line_graph& plg() const; + /// Return the set of neighbors (adjacent edge ids). + neighbors_t& neighbors(); + /// Read-only access to the \a i-th coordinate. // FIXME: Dummy. coord operator[](unsigned i) const; /// \} - /// Internals, used by the neighborhood. - /// \{ - public: - /// Set the iterator to the first site of the graph. - void first_(); - /// Advance the position of the iterator by one step. - void step_(); - - /// An internal iterator on the set of edges of the underlying graph. - util::edge_id id_; - /// \} - private: /// The neighborhood. const N& nbh_; /// The ``central'' psite of the neighborhood. const psite& p_ref_; + + /// The last reference psite whose ajacent psites have been computed. + psite saved_p_ref_; + /// The set of edge ids adjacent to the reference psite. + neighbors_t neighbors_; + /// An iterator on the set of adjacent edges. + neighbors_t::const_iterator i_; + /// The psite corresponding to this iterator. psite psite_; /// The point corresponding to this iterator. @@ -178,6 +186,9 @@ // FIXME: Dummy typedef. typedef void mesh; + // The type of the set of neighbors (adjacent edge ids). + typedef std::set<util::edge_id> neighbors_t; + public: /// Construction. /// \{ @@ -212,28 +223,31 @@ /// Convert the iterator into a line graph psite. operator psite() const; + /// Return the reference psite. + const psite& p_ref() const; + /// Return the mln::p_line_graph corresponding to this piter. + const p_line_graph& plg() const; + /// Return the set of neighbors (adjacent edge ids). + neighbors_t& neighbors(); + /// Read-only access to the \a i-th coordinate. // FIXME: Dummy. coord operator[](unsigned i) const; /// \} - /// Internals, used by the neighborhood. - /// \{ - public: - /// Set the iterator to the first site of the graph. - void first_(); - /// Advance the position of the iterator by one step. - void step_(); - - /// An internal iterator on the set of edges of the underlying graph. - util::edge_id id_; - /// \} - private: /// The neighborhood. const N& nbh_; /// The ``central'' psite of the neighborhood. const psite& p_ref_; + + /// The last reference psite whose ajacent psites have been computed. + psite saved_p_ref_; + /// The set of edge ids adjacent to the reference psite. + neighbors_t neighbors_; + /// An iterator on the set of adjacent edges. + neighbors_t::const_reverse_iterator i_; + /// The psite corresponding to this iterator. psite psite_; /// The point corresponding to this iterator. @@ -260,7 +274,6 @@ | line_graph_neighborhood_fwd_piter<P, N>. | `------------------------------------------*/ - // FIXME: Currently, argument nbh is ignored. template <typename P, typename N> template <typename Pref> inline @@ -271,7 +284,7 @@ // Initialize psite_ to a dummy value. psite_() { - // Invalidate id_. + // Invalidate i_. invalidate(); } @@ -280,10 +293,14 @@ bool line_graph_neighborhood_fwd_piter<P, N>::is_valid() const { - // FIXME: We depend too much on the implementation of util::graph - // here. The util::graph should provide the service to abstract - // these manipulations. - return p_ref_.is_valid() && id_ < p_ref_.plg().gr_->nedges(); + return + // The reference point must be valid... + p_ref_.is_valid() + // ...and must not have changed since the neighborhood has been + // computed... + && p_ref_ == saved_p_ref_ + // ...and the iterator i_ must point a valid value. + && i_ != neighbors_.end(); } template <typename P, typename N> @@ -291,7 +308,7 @@ void line_graph_neighborhood_fwd_piter<P, N>::invalidate() { - id_ = -1; + i_ = neighbors_.end(); } template <typename P, typename N> @@ -299,7 +316,15 @@ void line_graph_neighborhood_fwd_piter<P, N>::start() { - nbh_.start(*this); + mln_precondition(p_ref_.is_valid()); + // Update the neighbors, if needed. + if (!saved_p_ref_.is_valid() || p_ref_ != saved_p_ref_) + { + saved_p_ref_ = p_ref_; + nbh_.compute_neighbors_(*this); + } + i_ = neighbors_.begin(); + // FIXME: We might move the is_valid condition within update_. if (is_valid()) update_(); } @@ -309,7 +334,10 @@ void line_graph_neighborhood_fwd_piter<P, N>::next_() { - nbh_.next_(*this); + // Ensure the p_ref_ has not changed. + mln_precondition(p_ref_ == saved_p_ref_); + ++i_; + // FIXME: We might move the is_valid condition within update_. if (is_valid()) update_(); } @@ -317,84 +345,58 @@ template <typename P, typename N> inline void - line_graph_neighborhood_fwd_piter<P, N>::first_() + line_graph_neighborhood_fwd_piter<P, N>::update_() { - id_ = 0; + // Update psite_. + psite_ = line_graph_psite(plg(), *i_); } template <typename P, typename N> inline - void - line_graph_neighborhood_fwd_piter<P, N>::step_() + const P& + line_graph_neighborhood_fwd_piter<P, N>::to_point() const { - ++id_; + return p_; } - template <typename P, typename N> inline - bool - line_graph_neighborhood_fwd_piter<P, N>::adjacent_or_equal_to_p_ref_() const + const line_graph_psite& + line_graph_neighborhood_fwd_piter<P, N>::to_psite() const { - // Check wether the iterator points to P_REF_. - if (id_ == p_ref_.id()) - return true; - - // Check whether the iterator is among the neighbors of P_REF_. - { - // Paranoid assertion. - assert (p_ref_.id() < p_ref_.plg().gr_->nedges()); - /* FIXME: We should have a convenient shortcut for these - repetitive accesses to p_ref_.plg().gr_ (e.g., a method gr() - or g() in this class. */ - const util::tracked_ptr<typename p_line_graph::graph>& gr = - p_ref_.plg().gr_; - // Check whether the edge this iterator points to is adjacent to - // the one p_ref points to, by comparing their ajacent nodes. - /* FIXME: This is too low-level. Yet another service the graph - // should provide. */ - if (gr->edge(id_).n1() == gr->edge(p_ref_.id()).n1() || - gr->edge(id_).n1() == gr->edge(p_ref_.id()).n2() || - gr->edge(id_).n2() == gr->edge(p_ref_.id()).n1() || - gr->edge(id_).n2() == gr->edge(p_ref_.id()).n2()) - return true; - } - - // Otherwise, the iterator is not adjacent to P_REF_. - return false; + return psite_; } template <typename P, typename N> inline - void - line_graph_neighborhood_fwd_piter<P, N>::update_() + line_graph_neighborhood_fwd_piter<P, N>::operator line_graph_psite () const { - // Update psite_. - psite_ = line_graph_psite(p_ref_.plg(), id_); + mln_precondition(is_valid()); + return psite_; } template <typename P, typename N> inline - const P& - line_graph_neighborhood_fwd_piter<P, N>::to_point() const + const line_graph_psite& + line_graph_neighborhood_fwd_piter<P, N>::p_ref() const { - return p_; + return p_ref_; } template <typename P, typename N> inline - const line_graph_psite& - line_graph_neighborhood_fwd_piter<P, N>::to_psite() const + const p_line_graph& + line_graph_neighborhood_fwd_piter<P, N>::plg() const { - return psite_; + return p_ref_.plg(); } template <typename P, typename N> inline - line_graph_neighborhood_fwd_piter<P, N>::operator line_graph_psite () const + std::set<util::edge_id>& + line_graph_neighborhood_fwd_piter<P, N>::neighbors() { - mln_precondition(is_valid()); - return psite_; + return neighbors_; } template <typename P, typename N> @@ -420,7 +422,6 @@ | line_graph_neighborhood_bkd_piter<P, N>. | `------------------------------------------*/ - // FIXME: Currently, argument nbh is ignored. template <typename P, typename N> template <typename Pref> inline @@ -431,7 +432,7 @@ // Initialize psite_ to a dummy value. psite_() { - // Invalidate id_. + // Invalidate i_. invalidate(); } @@ -440,10 +441,14 @@ bool line_graph_neighborhood_bkd_piter<P, N>::is_valid() const { - // FIXME: We depend too much on the implementation of util::graph - // here. The util::graph should provide the service to abstract - // these manipulations. - return p_ref_.is_valid() && id_ < p_ref_.plg().gr_->nedges(); + return + // The reference point must be valid... + p_ref_.is_valid() + // ...and must not have changed since the neighborhood has been + // computed... + && p_ref_ == saved_p_ref_ + // ...and the iterator i_ must point a valid value. + && i_ != neighbors_.rend(); } template <typename P, typename N> @@ -451,7 +456,7 @@ void line_graph_neighborhood_bkd_piter<P, N>::invalidate() { - id_ = -1; + i_ = neighbors_.rend(); } template <typename P, typename N> @@ -459,7 +464,15 @@ void line_graph_neighborhood_bkd_piter<P, N>::start() { - nbh_.start(*this); + mln_precondition(p_ref_.is_valid()); + // Update the neighbors, if needed. + if (!saved_p_ref_.is_valid() || p_ref_ != saved_p_ref_) + { + saved_p_ref_ = p_ref_; + nbh_.compute_neighbors_(*this); + } + i_ = neighbors_.rbegin(); + // FIXME: We might move the is_valid condition within update_. if (is_valid()) update_(); } @@ -469,7 +482,10 @@ void line_graph_neighborhood_bkd_piter<P, N>::next_() { - nbh_.next_(*this); + // Ensure the p_ref_ has not changed. + mln_precondition(p_ref_ == saved_p_ref_); + ++i_; + // FIXME: We might move the is_valid condition within update_. if (is_valid()) update_(); } @@ -477,84 +493,58 @@ template <typename P, typename N> inline void - line_graph_neighborhood_bkd_piter<P, N>::first_() + line_graph_neighborhood_bkd_piter<P, N>::update_() { - id_ = p_ref_.plg().gr_->nedges() - 1; + // Update psite_. + psite_ = line_graph_psite(plg(), *i_); } template <typename P, typename N> inline - void - line_graph_neighborhood_bkd_piter<P, N>::step_() + const P& + line_graph_neighborhood_bkd_piter<P, N>::to_point() const { - --id_; + return p_; } - template <typename P, typename N> inline - bool - line_graph_neighborhood_bkd_piter<P, N>::adjacent_or_equal_to_p_ref_() const + const line_graph_psite& + line_graph_neighborhood_bkd_piter<P, N>::to_psite() const { - // Check wether the iterator points to P_REF_. - if (id_ == p_ref_.id()) - return true; - - // Check whether the iterator is among the neighbors of P_REF_. - { - // Paranoid assertion. - assert (p_ref_.id() < p_ref_.plg().gr_->nedges()); - /* FIXME: We should have a convenient shortcut for these - repetitive accesses to p_ref_.plg().gr_ (e.g., a method gr() - or g() in this class. */ - const util::tracked_ptr<typename p_line_graph::graph>& gr = - p_ref_.plg().gr_; - // Check whether the edge this iterator points to is adjacent to - // the one p_ref points to, by comparing their ajacent nodes. - /* FIXME: This is too low-level. Yet another service the graph - // should provide. */ - if (gr->edge(id_).n1() == gr->edge(p_ref_.id()).n1() || - gr->edge(id_).n1() == gr->edge(p_ref_.id()).n2() || - gr->edge(id_).n2() == gr->edge(p_ref_.id()).n1() || - gr->edge(id_).n2() == gr->edge(p_ref_.id()).n2()) - return true; - } - - // Otherwise, the iterator is not adjacent to P_REF_. - return false; + return psite_; } template <typename P, typename N> inline - void - line_graph_neighborhood_bkd_piter<P, N>::update_() + line_graph_neighborhood_bkd_piter<P, N>::operator line_graph_psite () const { - // Update psite_. - psite_ = line_graph_psite(p_ref_.plg(), id_); + mln_precondition(is_valid()); + return psite_; } template <typename P, typename N> inline - const P& - line_graph_neighborhood_bkd_piter<P, N>::to_point() const + const line_graph_psite& + line_graph_neighborhood_bkd_piter<P, N>::p_ref() const { - return p_; + return p_ref_; } template <typename P, typename N> inline - const line_graph_psite& - line_graph_neighborhood_bkd_piter<P, N>::to_psite() const + const p_line_graph& + line_graph_neighborhood_bkd_piter<P, N>::plg() const { - return psite_; + return p_ref_.plg(); } template <typename P, typename N> inline - line_graph_neighborhood_bkd_piter<P, N>::operator line_graph_psite () const + std::set<util::edge_id>& + line_graph_neighborhood_bkd_piter<P, N>::neighbors() { - mln_precondition(is_valid()); - return psite_; + return neighbors_; } template <typename P, typename N> Index: mln/core/line_graph_elt_neighborhood.hh --- mln/core/line_graph_elt_neighborhood.hh (revision 1854) +++ mln/core/line_graph_elt_neighborhood.hh (working copy) @@ -29,7 +29,7 @@ # define MLN_CORE_LINE_GRAPH_ELT_NEIGHBORHOOD_HH /// \file mln/core/line_graph_elt_neighborhood.hh -/// \brief Definition of the elementary ``window'' on a line graph. +/// \brief Definition of the elementary ``neighborhood'' on a line graph. /* FIXME: Have a consistent naming: we have window (without '_') but point_, neighb_, etc. */ @@ -40,6 +40,8 @@ - mln::line_graph_elt_window - mln::line_graph_elt_neighborhood. */ +# include <set> + # include <mln/core/concept/neighborhood.hh> # include <mln/core/line_graph_psite.hh> # include <mln/core/line_graph_neighborhood_piter.hh> @@ -66,6 +68,9 @@ typedef P point; /// The type of psite corresponding to the neighborhood. typedef line_graph_psite psite; + // The type of the set of neighbors (edge ids adjacent to the + // reference psite). + typedef std::set<util::edge_id> neighbors_t; // FIXME: This is a dummy value. typedef void dpoint; @@ -84,43 +89,44 @@ /// Services for iterators. /// \{ - /// Move \a piter to the beginning of the iteration on this neighborhood. - template <typename Piter> - void start(Point_Iterator<Piter>& piter) const; - /// Move \a piter to the next site on this neighborhood. template <typename Piter> - void next_(Point_Iterator<Piter>& piter) const; + void compute_neighbors_(Point_Iterator<Piter>& piter) const; /// \} }; -# ifndef MLN_INCLUDE_ONLY - template <typename P> - template <typename Piter> - inline - void - line_graph_elt_neighborhood::start(Point_Iterator<Piter>& piter_) const - { - Piter& piter = exact(piter_); - piter.first_(); - if (!piter.adjacent_or_equal_to_p_ref_()) - next_(piter); - } + +# ifndef MLN_INCLUDE_ONLY template <typename P> template <typename Piter> inline void - line_graph_elt_neighborhood::next_(Point_Iterator<Piter>& piter_) const + line_graph_elt_neighborhood::compute_neighbors_(Point_Iterator<Piter>& piter_) const { Piter& piter = exact(piter_); - /* FIXME: This is inefficient. The graph structure should be able - to produce the set of adjacent nodes fast. Boost Graphs - probably provides adequates structures to fetch these - neighbors in constant time. */ - do - piter.step_(); - while (piter.is_valid() && !piter.adjacent_or_equal_to_p_ref_()); + neighbors_t& neighbors = piter.neighbors(); + neighbors.clear(); + // Add the reference piter itself. We might reconsider this, or + // create an elementary neighbors without the reference point. + neighbors.insert(piter.p_ref().id()); + /* FIXME: Move this computation out of the window. In fact, + this should be a service of the graph, also proposed by the + p_line_graph. */ + // Ajacent edges connected through node 1. + // FIXME: Far too low-level. + util::node_id id1 = piter.p_ref().first_id(); + const util::node& node1 = piter.plg().gr_->node(id1); + for (std::vector<util::edge_id>::const_iterator e = + node1.edges.begin(); e != node1.edges.end(); ++e) + neighbors.insert(*e); + // Ajacent edges connected through node 2. + // FIXME: Likewise. + util::node_id id2 = piter.p_ref().second_id(); + const util::node& node2 = piter.plg().gr_->node(id2); + for (std::vector<util::edge_id>::const_iterator e = + node2.edges.begin(); e != node2.edges.end(); ++e) + neighbors.insert(*e); } # endif // ! MLN_INCLUDE_ONLY

1 0

1854: Sandbox: ICP: Add multi-scale registration.
by Ugo Jardonnet 09 Apr '08

09 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: Add multi-scale registration. * sandbox/jardonnet/test/icp.cc: Add nb_it as argument for a registration. * sandbox/jardonnet/test/Makefile: Remove icp_lazy test. Add parameter for multi-scale icp: * sandbox/jardonnet/registration/tools.hh: Update. * sandbox/jardonnet/registration/quat7.hh: Use virtual length (c_length) for p_array. * sandbox/jardonnet/registration/cloud.hh: Use virtual length. * sandbox/jardonnet/registration/icp.hh: Use virtual length. Following files now exist eleswhere (registration/) * sandbox/jardonnet/registration/quat: Remove. * sandbox/jardonnet/registration/quat/all.hh: Remove. * sandbox/jardonnet/registration/quat/misc.hh: Remove. * sandbox/jardonnet/registration/quat/rotation.hh: Remove. * sandbox/jardonnet/registration/projection.hh: Use virtual length. registration/cloud.hh | 15 ++++++---- registration/icp.hh | 62 +++++++++++++++++++++++++++++++-------------- registration/projection.hh | 7 ++--- registration/quat7.hh | 24 ++++++++++------- registration/tools.hh | 26 +----------------- test/Makefile | 8 ----- test/icp.cc | 6 ++-- 7 files changed, 77 insertions(+), 71 deletions(-) Index: sandbox/jardonnet/test/icp.cc --- sandbox/jardonnet/test/icp.cc (revision 1852) +++ sandbox/jardonnet/test/icp.cc (working copy) @@ -41,18 +41,20 @@ closest_point<mln_point_(I3d)> fun(x, box_<point3d>(1000,1000,1)); lazy_image< closest_point<mln_point_(I3d)> > map(fun); - quat7<3> q = registration::icp(c, map); + quat7<3> q = registration::icp(c, map, 2); #ifndef NDEBUG std::cout << "closest_point(Ck[i]) = " << fun.i << std::endl; std::cout << "Pts processed = " << registration::pts << std::endl; #endif + //FIXME: register img1 + //init output image //FIXME: Should be //mln_concrete(I) output(res.bbox()) = convert::to_image(res)? - q.apply_on(c, c); + q.apply_on(c, c, c.npoints()); const box_<point2d> box2d(1000,1000); image2d<bool> output(box2d, 1); Index: sandbox/jardonnet/test/Makefile --- sandbox/jardonnet/test/Makefile (revision 1852) +++ sandbox/jardonnet/test/Makefile (working copy) @@ -48,14 +48,6 @@ g++ icp.cc -I../../.. -O3 -ffloat-store -o './bin/+icp_3f' g++ icp.cc -I../../.. -O3 -DNDEBUG -o './bin/+icp_3D' g++ icp.cc -I../../.. -O3 -DNDEBUG -ffloat-store -o './bin/+icp_3Df' - g++ icp_lazy.cc -I../../.. -O0 -o './bin/+icp_lazy_0' - g++ icp_lazy.cc -I../../.. -O0 -ffloat-store -o './bin/+icp_lazy_0f' - g++ icp_lazy.cc -I../../.. -O0 -DNDEBUG -o './bin/+icp_lazy_0D' - g++ icp_lazy.cc -I../../.. -O0 -DNDEBUG -ffloat-store -o './bin/+icp_lazy_0Df' - g++ icp_lazy.cc -I../../.. -O3 -o './bin/+icp_lazy_3' - g++ icp_lazy.cc -I../../.. -O3 -ffloat-store -o './bin/+icp_lazy_3f' - g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -o './bin/+icp_lazy_3D' - g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -ffloat-store -o './bin/+icp_lazy_3Df' ./icp_check.sh clean: Index: sandbox/jardonnet/registration/tools.hh --- sandbox/jardonnet/registration/tools.hh (revision 1852) +++ sandbox/jardonnet/registration/tools.hh (working copy) @@ -10,6 +10,7 @@ namespace mln { + //FIXME: groe length template <typename P> struct closest_point { @@ -18,11 +19,9 @@ closest_point(const p_array& X, const box_& box) : X(X), box(box) - #ifndef NDEBUG , i(0) #endif - { } result @@ -80,6 +79,7 @@ : value(nrows, ncols,1), is_known(nrows,ncols,1), fun(fun) { } + // FIXME: gore length const mln_result(F) //inline operator () (const typename F::input& p) const @@ -254,28 +254,6 @@ } // end of namespace convert - - namespace fun - { - //FIXME: temporary - template <typename C, typename T= float> - struct cham : public Function_p2v< cham<C,T> > - { - typedef T result; - //bad - T operator()(dpoints_fwd_piter<mln::dpoint_<mln::grid::cube, int> >& v) const - { - C o = C::origin; - if (v < o) - return 1.; - else - return 0.; - } - }; - } // end of namespace fun - - - namespace algebra { Index: sandbox/jardonnet/registration/quat7.hh --- sandbox/jardonnet/registration/quat7.hh (revision 1852) +++ sandbox/jardonnet/registration/quat7.hh (working copy) @@ -59,12 +59,15 @@ } template <typename P> - void apply_on(const p_array& input, p_array& output) const + void apply_on(const p_array& input, + p_array& output, + const size_t c_length) const { - assert(input.npoints() == output.npoints()); + assert(c_length <= input.npoints()); + assert(c_length <= output.npoints()); assert(_qR.is_unit()); - for (size_t i = 0; i < input.npoints(); i++) + for (size_t i = 0; i < c_length; i++) output.hook_()[i] = algebra::to_point((*this)(input[i])); } @@ -95,27 +98,30 @@ quat7<P::dim> match(const p_array& C, const algebra::vec<P::dim,float>& mu_C, const p_array& Ck, - M& map) + M& map, + size_t c_length) { + //FIXME: mix mu_Xk and qk loop + + //mu_Xk = center map(Ck) algebra::vec<P::dim,float> mu_Xk(literal::zero); - for (size_t i = 0; i < Ck.npoints(); ++i) + for (size_t i = 0; i < c_length; ++i) { algebra::vec<P::dim,float> xki = map(Ck[i]); mu_Xk += xki; } - mu_Xk /= Ck.npoints(); - + mu_Xk /= c_length; // qR algebra::mat<P::dim,P::dim,float> Mk(literal::zero); - for (size_t i = 0; i < C.npoints(); ++i) + for (size_t i = 0; i < c_length; ++i) { algebra::vec<P::dim,float> Ci = C[i]; algebra::vec<P::dim,float> Xki = map(Ck[i]); Mk += make::mat(Ci - mu_C) * trans(make::mat(Xki - mu_Xk)); } - Mk /= C.npoints(); + Mk /= c_length; /* Index: sandbox/jardonnet/registration/cloud.hh --- sandbox/jardonnet/registration/cloud.hh (revision 1852) +++ sandbox/jardonnet/registration/cloud.hh (working copy) @@ -18,16 +18,16 @@ { template <typename P> - P center(const p_array& a) + P center(const p_array& a, size_t length) { algebra::vec<P::dim,float> c(literal::zero); - for (size_t i = 0; i < a.npoints(); ++i) + for (size_t i = 0; i < length; ++i) { algebra::vec<P::dim,float> ai = a[i]; c += ai; } - return algebra::to_point(c / a.npoints()); + return algebra::to_point(c / length); } @@ -43,16 +43,19 @@ template <typename P> float rms(const p_array& a1, - const p_array& a2) + const p_array& a2, + const size_t length) { - assert(a1.npoints() == a2.npoints()); + assert(length <= a1.npoints()); + assert(length <= a2.npoints()); /* float f = 0.f; for (size_t i = 0; i < a1.npoints(); ++i) f += sqr_norm(a1[i] - a2[i]); return f / a1.npoints();*/ + float f = 0.f; - for (size_t i = 0; i < a1.npoints(); ++i) + for (size_t i = 0; i < length; ++i) { algebra::vec<3,float> a1f = a1[i]; algebra::vec<3,float> a2f = a2[i]; Index: sandbox/jardonnet/registration/icp.hh --- sandbox/jardonnet/registration/icp.hh (revision 1852) +++ sandbox/jardonnet/registration/icp.hh (working copy) @@ -35,6 +35,7 @@ # include <iostream> # include <string> +# include <cmath> # include <mln/algebra/quat.hh> # include <mln/algebra/vec.hh> @@ -78,48 +79,56 @@ p_array icp_(const p_array& C, M& map, - quat7<P::dim>& qk) + quat7<P::dim>& qk, + const size_t c_length, + const float epsilon = 1e-3) { trace::entering("registration::impl::icp_"); #ifndef NDEBUG //display registred points std::cout << "Register " - << C.npoints() << " points" << std::endl; + << c_length << " points" << std::endl; std::cout << "k\terror\tdqk" << std::endl; #endif - - unsigned int k; quat7<P::dim> old_qk; float err; //float err_bis; p_array Ck(C), Xk(C); //FIXME: Xk copy C - algebra::vec<P::dim,float> mu_C = center(C), mu_Xk; + algebra::vec<P::dim,float> mu_C = center(C, c_length), mu_Xk; - const float epsilon = 1;//1e-3; + qk.apply_on(C, Ck, c_length); - //// step 1 - k = 0; + unsigned int k = 0; do { //compute qk old_qk = qk; - qk = match(C, mu_C, Ck, map); + qk = match(C, mu_C, Ck, map, c_length); //Ck+1 = qk(C) - qk.apply_on(C, Ck); + qk.apply_on(C, Ck, c_length); //err = d(Ck+1,Xk) - err = rms(Ck, Xk); + err = rms(Ck, Xk, c_length); #ifndef NDEBUG - { using namespace std; - image2d<bool> img = convert::to_image2d(Ck); + + image2d<bool> img(box2d(500,800), 0); + for (size_t i = 0; i < c_length; i++) + { + point2d p = convert::to_point2d(Ck[i]); + if (img.has(p)) + img(p) = true; + } + + //image2d<bool> img = convert::to_image2d(Ck); stringstream oss; - oss << "reg" << k << ".pbm"; + static int pimp = 0; + oss << "i_" << pimp++ << ".pbm"; io::pbm::save(img, oss.str()); } @@ -128,10 +137,10 @@ << (qk - old_qk).sqr_norm() << '\t' << std::endl; //count the number of points processed - pts += Ck.npoints(); + pts += c_length; #endif - ++k; + k++; } while (k < 3 || (qk - old_qk).sqr_norm() > epsilon); trace::exiting("registration::impl::icp_"); @@ -145,19 +154,34 @@ template <typename P, typename M> inline quat7<P::dim> - icp(const p_array& cloud, - M& map) + icp(p_array cloud, + M& map, + const unsigned nb_it) { trace::entering("registration::icp"); mln_precondition(P::dim == 3); mln_precondition(cloud.npoints() != 0); + // Shuffle cloud + for (size_t i = 0; i < cloud.npoints(); i++) + { + size_t r = rand() % cloud.npoints(); + P tmp; + tmp = cloud[i]; + cloud.hook_()[i] = cloud[r]; + cloud.hook_()[r] = tmp; + } + //init rigid transform qk quat7<P::dim> qk; //run icp - p_array res = impl::icp_(cloud, map, qk); + for (int i = nb_it; i >= 0; i--) + { + size_t l = cloud.npoints() / std::pow(2., i); // 3 is arbitrary fixed + impl::icp_(cloud, map, qk, l, i + 1e-3); + } trace::exiting("registration::icp"); Index: sandbox/jardonnet/registration/projection.hh --- sandbox/jardonnet/registration/projection.hh (revision 1852) +++ sandbox/jardonnet/registration/projection.hh (working copy) @@ -70,11 +70,12 @@ template <typename P, typename T> void memo(const p_array& Ck, p_array& Xk, - lazy_image<T>& map) // first: closest points, second: is_computed + lazy_image<T>& map, + const size_t c_length) // first: closest points, second: is_computed { - assert(Ck.npoints() == Xk.npoints()); + //assert(Ck.npoints() == Xk.npoints()); //FIXME: - for (size_t i = 0; i < Ck.npoints(); ++i) + for (size_t i = 0; i < c_length; ++i) Xk.hook_()[i] = map(Ck[i]); }

1 0

1853: Add a log version of front-based distance computation.
by Thierry Geraud 09 Apr '08

09 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk Index: ChangeLog from Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Add a log version of front-based distance computation. * milena/sandbox/folio/psn_log.cc: New. Derived from psn.cc to trace the behavior of the algorithm in the non-convex case. psn_log.cc | 290 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 290 insertions(+) Index: milena/sandbox/folio/psn_log.cc --- milena/sandbox/folio/psn_log.cc (revision 0) +++ milena/sandbox/folio/psn_log.cc (revision 0) @@ -0,0 +1,290 @@ + +// Copyright (C) 2007 EPITA Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +#ifndef MLN_DT_CHAMFER_HH +# define MLN_DT_CHAMFER_HH + +# include <vector> +# include <queue> +# include <map> +# include <cmath> + +# include <mln/core/concept/image.hh> +# include <mln/core/concept/neighborhood.hh> +# include <mln/literal/zero.hh> + +#include <mln/core/image2d.hh> +# include <mln/debug/println.hh> + +namespace mln +{ + + namespace dt + { + + /*! Propagation using a single neighborhood (PSN). + * + * \param[in] input_ The input image. + * \param[in] chamfer The chamfer window to use for distance calcul. + * \return A pair <distance map, closest point map>. + * + * \pre \p img has to be initialized. + */ + template<typename I, typename N> + mln_ch_value(I, unsigned) + psn(const Image& input_, const N& nbh); + + +# ifndef MLN_INCLUDE_ONLY + + + namespace impl + { + + template <typename BP> + class compare + { + public: + bool + operator()(const BP& lhs, const BP& rhs) const + { + return lhs.second > rhs.second; + } + }; + + template <typename D> + unsigned + sq(const D& dp) + { + unsigned res = 0; + + for (unsigned i = 0; i < D::dim; ++i) + res += std::abs(dp[i]); // FIXME: dp[i] * dp[i]; + + return res; + } + + template <typename BP> + unsigned + size(const std::map<unsigned, std::queue<BP> >& bucket, + int d) + { + unsigned s = 0; + typename std::map<unsigned, std::queue<BP> >::const_iterator i; + for (i = bucket.begin(); i != bucket.end(); ++i) + if (i->first >= d) + s += i->second.size(); + return s; + } + + template <typename BP> + unsigned + size(const std::map<unsigned, std::queue<BP> >& bucket) + { + unsigned s = 0; + typename std::map<unsigned, std::queue<BP> >::const_iterator i; + for (i = bucket.begin(); i != bucket.end(); ++i) + s += i->second.size(); + return s; + } + + template <typename BP> + void + print(const std::map<unsigned, std::queue<BP> >& bucket) + { + typename std::map<unsigned, std::queue<BP> >::const_iterator i; + int d = -1; + for (i = bucket.begin(); i != bucket.end(); ++i) + { + if (i->first != d) + { + d = i->first; + std::cout << std::endl << d << ": "; + } + std::queue<BP> qu = i->second; // copy + std::vector<BP> v; + v.push_back(qu.front()); qu.pop(); + typename std::vector<BP>::const_iterator j; + for (j = v.begin(); j != v.end(); ++j) + std::cout << j->first << " "; // point + } + std::cout << std::endl; + } + + } // end of namespace mln::dt::impl + + + + // Facade. + + template<typename I, typename N> + inline + mln_ch_value(I, unsigned) + psn(const Image& input_, const N& nbh) + { + const I& input = exact(input_); + mln_precondition(input.has_data()); + + mln_ch_value(I, unsigned) D; + initialize(D, input); + + static const unsigned M = 666; // FIXME + + // Initialization. + typedef mln_point(I) point; + typedef mln_dpoint(I) dpoint; + typedef std::pair<point, dpoint> BP; + + std::map<unsigned, std::queue<BP> > bucket; + unsigned bucket_size = 0; + + mln_fwd_piter(I) p(input.domain()); + for_all(p) + if (input(p)) + { + D(p) = literal::zero; + bucket[0].push(BP(p, literal::zero)); + ++bucket_size; + } + else + D(p) = M; + + debug::println(input); + debug::println(D); + impl::print(bucket); + + unsigned d = 0; + + while (impl::size(bucket, d) != 0) + { + + mln_invariant(impl::size(bucket) == bucket_size); + + std::cout << "PROCESSING d = " << d << std::endl; + + std::queue<BP>& bucket_d = bucket[d]; + + while (! bucket_d.empty()) + { + point p = bucket_d.front().first; + dpoint dp = bucket_d.front().second; + + bucket_d.pop(); + --bucket_size; + + std::cout << "pop " << p << " at D=" << D(p) << std::endl; + + if (D(p) == d) + { + mln_niter(N) n_(nbh, p); + + for_all(n_) + if (D.has(n_)) + { + dpoint n = n_ - p; + unsigned newD = impl::sq(dp + n); + + std::cout << " n=" << n_ << " at D=" << D(n_) + << " and newD=" << newD + << (newD < D(p + n) ? " push" : "") + << std::endl; + + if (newD < D(p + n)) // p + n = p + n_ - p = n_ + { + D(n_) = newD; + bucket[newD].push(BP(p + n, dp + n)); + ++bucket_size; + } + } + } + } + ++d; + } + + return D; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::dt + +} // end of namespace mln + +#endif // ! MLN_DT_CHAMFER_HH + +#include <iostream> +#include <mln/debug/println.hh> +#include <mln/make/win_chamfer.hh> +#include <mln/level/fill.hh> +#include <mln/core/neighb2d.hh> + +#include <mln/io/pbm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/level/stretch.hh> +#include <mln/value/int_u8.hh> + +#include <mln/core/sub_image.hh> +#include <mln/core/image_if.hh> +#include <mln/pw/value.hh> + +int main() +{ + using namespace mln; + +// image2d<bool> ima(5,5); +// bool vals[] = { 1, 1, 1, 0, 0, +// 0, 0, 1, 0, 0, +// 0, 0, 1, 0, 0, +// 0, 0, 0, 0, 0, +// 0, 0, 0, 0, 0 }; + + image2d<bool> ima(3,3); + bool vals[] = { 1, 0, 0, + 0, 0, 0, + 0, 0, 0}; + level::fill(ima, vals); + + image2d<bool> msk(3,3); + bool rest[] = { 1, 0, 1, + 1, 0, 1, + 1, 1, 1}; + level::fill(msk, rest); + + image2d<unsigned> out; + out = dt::psn(ima | pw::value(msk), c4()); + + debug::println(ima | pw::value(msk)); + debug::println(out); + +// image2d<bool> ima = io::pbm::load("../../img/c01.pbm"); + +// image2d<value::int_u8> out2(out.domain()); +// level::stretch(out, out2); + +// io::pgm::save(out2, "out.pgm"); +}

1 0

1852: Add img/quite-small.pgm.
by Roland Levillain 09 Apr '08

09 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Roland Levillain <roland(a)lrde.epita.fr> * img/quite-small.pgm: New image. quite-small.pgm | 4 ++++ 1 file changed, 4 insertions(+) Index: img/quite-small.pgm --- img/quite-small.pgm (revision 0) +++ img/quite-small.pgm (revision 0) @@ -0,0 +1,4 @@ +P5 +128 128 +255 +��v_`fijkkkkioux{}��zxrjp��ҝiktxyyzyyy{||z|}}zuy��t^]dhhhiiiintx{}~��~zxqjn��֮pjswxyzyy{||{{{}}|y||f��s][afghhgginsx{|~��~{wrliu��ˏintwxyyxy{z{{|}}�~d<.��r][afijigfhnsvz|~��|yspkl��ٸyjpuvxyxyyz{|}~��b;/0��s]Zbfiihgfgnswz|}��~��}zuqliv��Ҟklsuxzzz{|}��]5.//��w^\agiggfghnrvx{|��~~�|{xqlip��Ájpuwyzyy~��]6/1/1��rZX`fhhhffilqwy{~~��~}}~~~~��~~~~~|{ypljr��׫omtyzzxy|~��[8.,-14��q[Yafiihggikqwy{}~~��|zz|}}}~��~~~}~~~ywsnks��Ώkqvyzyz|��Z1,-.125��q[[aehihffhlqvyz|}�~~��|{yz|}|{{~��~}~}~~~zwqnks��ܽzjruyz}�[2+223325��pWW`egffffglquxx{|}}}|~��~��~}{z{{{xxy|��~}}~~~}zvupkq��֞hlrwz�Y1*-575310��oWW`fffgffgkqwwyz{{}}|}��~}{zyyxwvx|~��~~~�~~~~~{yusrmq}��~gmrz}\2*-39640--��pXW`eghggfimrwxz{||~}}��~��}xx{|}~��~~~~~|yurpko{��اmmtxY3,-16973/-.��rYW^dhhfgfinqvyy|}}}~��zvy{��~}~}{xvqokn{��͌pvX0*,.5686.,-.��qUU]cgfgghjotwyy|~}|��~~��wvz}}}|{ywsqojn|��߼[2)-/254201/--��z��oSR\`cdedgimsvwz{{{|}~~~~��}��sz||z|zzvtqolio{��ێ4),/4553012/-/��on��pSPX^bacdfhlpuuxxz{{|~}}��}z��ŷ�uqvxyxxuspmjgm{��t&*-1554-/4/.03��zZi��pQMV]_`cdegilrtuvzz{|}{}�syzyzy��Ŀ�ykrstutroligk{��ޠ7'+0566/,10,03*��dNi��nOMW]cddcefimqqtuwx{|||�zovvy{y��Ƿ�fgkpqpmkjfj{��B)+-1351132.2;.!��uSOl��mNMX_cddcceilprrvwwy{z��pqoqv|}��ž�rbhmljjjehz��>*+,342/-5403=6'Q��aRSl��lONW`cbcddfhmqsrvxxyy{�}nknuuy}}��¿��[gifggdi{��qp~��H++,/44/-4501=@*L��pTSUm��kONY_bccddfjnqrtwyzyy~|olouvvyyz|~��b[deddahy��fWi��X,*,0441,15//7G>H��}\TWWm��mROY^acdfdeimqrtvxy{||ojnsvxxwvx~��ݡWY`bb_fx��d@Ov��o0')0554,04/0:CFP}��gTXXWl��lQOZ_`defdfjlpstxyyz{pkmpswuurv|��؊PW]_^dx��i:8a~��@2-/555/071.5@IUz��pXX[XWj��kQP[abbcddfhlqtuvxz|ukopruuqpr|��oMRWZcy��k8.Pr��T)..454/,63137CWw��[UZ[YWj��mQNZ`cccceghlrutvx|{qqoortstux~��ƊQMT`x��i:,<_|��o0&+5861+35.46;It��WY[[ZXj��mRMX_bcdcefimqtvvx�sjoqoqrsvx��הHK\v��i>++Cj��>'*/6760070255Aj��W[\\ZZl��mSMX^`ceccehlptur��ggkqppqqw{��nCWs��l@,)/Ou��W++.3664/500427_��Z\\[YXj��nRNV]abeedfjmqsrn��ejlnortvz}��~v��޷VMp��lB-,)4Ux��p2+.3556144162+P��\[ZZWUh��pRMV]^acddijloqly�vfhkopoux|��{v��߉Fm��lB-,($S��ȧB*/1667754062-B{��ZYXYYTf��mNJTY^```bfhjlmg��oeimqqquw��~px��߸Qf��j@,''Y��M)-467635338/4m��YZ[]\Xg��kLHQ[_bbbcehinnj��odikrqtu|�~{z��}sz��؄_��g:$-t��{(-2641110542Y��\\]_^Zg��kNJR\acddcehloot��gdhlmmt{{w|��|q|��q��_.8�Ʒ��*-2410.,261D��^_aba_i��mQMV^bcdedhjmpk~Ùhdhfjsy{yxz��r~��Գ��yTc�ƻ��)063---13,0i��bbcca_i��nROW_cefgeflmpj�ŉhafjpptuw|~{r~��ƻ��Ծ��¿��*0351-052#J��cddcbal��pQOW^cefecdhkme��ndhhglrs||xvs��ƺ��ͩ��ž��q*2030+/3/7t��dccccbo��oQOZ_ccedcdhjic��qghggkmuwvxu��¶��¿��Z.431,-15<_��ddcbcco��pRNY]bdcbbdfheg��vhihkorvv{u��B,142,/3:Q��cbbabcm��pQMW^acbbbcehcpŵ��vnhjnpxzv~��ܨ.+.30-43;o��abbabal��rSOZ^adbcbcfgbtǲ��tjhkqto��Ʃ��x),030.33K��a````al��qSO[_cccbbcefauı��}l`jro}��|��˰��F,43/-.2Ht��a_^^_am��qSOZ`bcabbdec\y¯��kdik��|��Ͽ��**43+1/;b��a_^_^ao��sTPZ`aaabbcccY�ĵ��vjh|��|z��¾��̨��k(/11232H~��`a``aer��tRMX_a`_aabbaU�˼��fu��~{��|��ü��һ��x��=(221966^��aaabcgr��rQKW_baabbbb`T��r{~�|w��~��xvklr��z��kx��ɢ�a^�Ɍ%)0./64C{��bbbbdhv��sPKU^`bccccc`T��ø��svzwt��|��zn_PY��rg_heXJB>C��ʷ�jV\|��C(//-031]��cecbekv��uPKV^`bbcbccaXr��ç��ouum��}��x�|qlo^Yo|fKLZKD3))A|��ŧ�f]Xb��k+/21/4-6x��efeaekv��tPLW_`bbdcdb`Za��Ʃ~��{rm{��z��x�~rsbSCAWm_\bQ@K-$C��˰w^VXft��T,-2//36.O��hhebdiu��tRLW^`abdbcba_W��Ƭ��kv��z~��u{�gUTK:4Kl_AV^S:I%O��Ƭ��eai{��W*/44--69:p��ihgedgu��vRKU]_abcbbbcbT��ù��tr��{t~�~{��ytWFMB717fb;D\Q</O]��Ű�ei��W(/563//51L��jhgeefs��wRKW]_`bdabbeeZm�Ⱦ��o�rz}|z��ltmMFK=66+KV;AbcC>;y��·��YRrn��_(,3652-234e��jjheddo��xPJV\_`aaaabghc\�ҽ��wz�}s|~yy��}\MLKM>7P.3R94FMMMQ[��þ��MN_ozw��Q)*05861-33Cy��kjhfcal��wPJS[_aaabadgih`j�Ƿ��wy��|z�zm]PKT9*EW*<H5;8DD_po��fVde|x��?),25472./48Y��kjgfc_g��xQHRY_``abbeillg[��ç�yrt�{v|~zpPHMYP6.HM9;37BAHL]v��Ż�YY_c��u3*-06777/-23An��kigec]e��zQIRZ_abaabehmnj^��̭zpwxx��sqbG;8MJ@7@VM..01*<ZWw��ßfCHTs��V**-57995-,33R~��ihheb[d��yQIU]abccbbeioqma��ώowxus��~dCBE0.@JN;;Jh;1,&1`h|��ŪqD1Jc��K--/59;;2-.07j��ghgec]e��zSJS]`cddbdgiorpgx��tkpsv}wroZ1F4-/5GG@B4fQ8#*k|��͸L+:V��>-/268<9/-1/H��fhgeb^d��|VJS\beffdddt|rrmk��qtns�{eTHA:L,5617=>QCSC#&b��ĒX/,D��E2368:>6-13.]��hhedb]e��{WMV_dfgfddgxustsm��wqu|�~gJ:.2G2953468WVRW"Y��ˢb6%9s�O357:=>3.218t��ggdc_[c��{VOZ_cefgffhssuwwt{�}s��fE8I.-F:<627<3EB6YP��ѯqB)1c�U069:;;204.N��ffdba\b��zUMX_ceghffiotssqszxw}��rF=/I703GG//;<99F!L{��һO-,V�\36:<>61324j��fecaa\b��zULV_degffgiotvwx�{y�{�zaAF,V21EP8+2=6G?:1hx��K+*G�j;>:>=325/C}��ddccb]a��zUMV]dgfeefipuww��J6m�xdDQ:<[*2;8./4=;H=.^y��ljpr��¾��ŧn:)(@�y;><A<1240W��feedc]`��{TOW^cgeeefhnsx��|4,t�^HEP4O[)6=>9,*CMB.H~��cPVv��~i^I0*8��=>@>51366i��dfgeb[_~��}VOX_efeccdgjo��l[ofF>WGLCM7;1<N)(.C70s��qwum�lg`^z��f^ceS<,6{�D:?=3244Cz��dedc`X\~��WLU]cfdcbabx��v?9FV^J\6==A/+RC,,(([��ö�jLHA:8?4Lt|}��ó�YD99>A9)6v�N:@8014-O��a`_^\SVy��}UJT[`cba_Yp��Q=*OPUMWg578L8'5ZC/$:|��®f<331-FtK+O~��t;7OT2571*5o�Q??4/240b��b`_]WORy��|ULX_bbbc_Up��o;0IED9Ae`1.,BI++364/c��ȺxM57LA:X��[Woy��ԖB6Ao�F/2/-9j�WG?213/=z��cba]XNQx��|XO[cfffh_DM��L6,IG;CASie1..4I/,/''Fo��ż�rnT@fXO�Ǌmou��zeVT~�R/--->d�YO;124.L��dcb]YMPx��~\T^fjljHA^{��T<5SY5<MPahrB1,,33//#1m��Ĺ�w��eky��|y��}||��G8514@^�^P8/041`��ffc\VMNw��\U`hnoqBi��y>0V`ZK*CPL]tu]6+,./1)!T{��r��w{��}~��͡��rQIK8/7@W�gM5/25Aw��ed`XSJJt��[R^glkx��ZQL/0>nva6:KXXhx_6...//$6w��c{��Ү��vccW:-6>T�oF4/50J��b]YSOCCq��\S^fllq��fUT:3GV`XQ+6MOoac]4.-,1/+a{��ÅYp��չ��xr`?18<O�tA2.2-Y��]ZUQK@Co��[Q\dmw��vnH:EZje^9'4TYelr_L6/-/,Dx��ÓIbv��׼��yf?2;;G�|91/14n��WWPOI?Bo��ZP[jz��bKG;AP^opH0,,6N\P5gs@--,-h~��CGky��ؾ��}h>4><A��9-0.>~��QNKKG>An��[Q\intp[[XMBXIfqO;.0019=9*h�N((%A{��\.Ppz��}b93==;��=+6.N��MJIHD;?m��YP`diqmpd^[\K;VnL8340/4:42dx{6#&b{�ȅ.6Ooz��Ĥ��{Z41=<4��L+72d��LJKHB8=m��XMogw�wJ\s]CJM`y]:;4))7J:C}qra!9v�Ű<,<Mjy��}}��ʣ��xN/6B</��Z+38v��MKKG@69k��VJj|�}aGoaQC\^kyg>50&+F]Fa��FFX|��Z!1CLhx��{z��ԣ��qA-7F<+}�c*-F��PMJG>24i��XI[rt|_L\^OO_U`|]<+2.-GPlo[X.(;o��&$6MMfw~��|�{��ڤ��e4-9F;(p�r.*\��SOLG>23h��ZVqqmrQN``KI^B`w^L371*G]XlH,#;��S#':ULgw~��w~��֟��{Q+/9F9(`�5+p��SOMIA7<o��c[acmoNPm^8]Y>^sn]96=(39G\/& c��l+&*:YRhw}��{�qcsw��n8+2<C8*Q��;6��OMLIC@Ly��]KXcmgRSaM>jQI_n~e:0KD$1eH25:z��*+-+<`Ykv|��~qknm��V/.3@E9-B��DM��NLKJDGZ��\MYbpn@KI;\LNOcgy|L(X|,4O@Xv�o��8);5+?f[nw|��v=.06GJ=18��M`��MKMJCEa��^MYau~@FOaL4GQcjf~m=Qe:4Ew��vo�Y/4890?g_oz}��¿��^2/0<KJ?61��Zu��IHHG?Cb��]KVg�qEOiJ<7CHNaXu|W_C2;m��^h,016=08ebny~��Ǵ��}@//2?LG@3+w�g��FDDD@F^��\JYwskG^]4CCDKETLc{yV,Qiz��kns08516B23`_hw~��ɼſ��d.114@LF?4*o�p��@=BC@E\��[Iknki<RR9??L[RMRJd~MF{��xf�=*H308F52ZZ`t}��D,326BLJC7&d�y��9:>A>F^��Z^q^j9BA?B9MreYM[`\\|��ˈb�[*-<2/8C91TUYp}��rchqutnikz}iqu�k0/11:BIL=<)X��79=?>Md��nuXj�`A;2:CBS\d`]`N`��i�~42-4419>>3MQUkx||��yqpu|��C-2/1<?KK<E.M��78;=<Od��zTS|rXKS-3KG^�JYf[Jb��~��D15./508=?4HLLcqvy|��}}��_,.2,35?NM>G2J��5588:Pa��[Ek{jHYU-@LLs}`YSimo��R35943459=A7GJFYintz��~~|��z=,/2255BMOCG6I��3353:W_}��VGyj_;`J+F?@�vkwV\~po|��c23693348;<>8?IBN_gmu|��}}��R2*.4:64FPQBG:G��~~��221/Acg~��YGrhSDT:,G<=�|o�mbin`r��z22238334469;86A>DR`iq{��j:1,18:41IRQFJ<C��||��/011Iei��[FbhOQQ5.I5?pgk}xdYvs`��>1426:144206985<;?@M^hu��I6//5;971OUSHK@C��{��.1.1Hah��ZG]]AEU73E4@oVCi�zk{�ag\<3528<1635027879<>>4DWhy��~78337<462TWSJIDC��z��01-/>Yl��WB_N3FT:<A14tj@Umz�{�q^G53429@2946214899<;E61=Rcmy��m37777?646WXRNNHE��x��k-/,,6Td}��Rc|c5GE<>9/.m�WKbnnw��s6/34@E/;434115;<>9DG<BKS_jt~��e45888>738V\UQRHM��y��v��zK1-34,)2KVv��~~�]d3?G>>90/f�pGFcqo��^�z9,2BE.<2141128::5:LUgt|��yS:559329T[VPRGS��soqwz|��{w��ɹ�V.*3>@=3-0BJn��TNT):M<:FB5Uxz]Fs�]gz]a�s73DC0@4323246:957Mbw~��Ƹ�[:311<XYSOSKZ��z{tsqonptz}�~xz��ȼ�`0.;LOiXD88<;[��XPoC&.FA0DWC9Ztch�gXcpd�k^2@?/C4106214;;85Kfw|}}��˿�U-&7TVPLNF]��}xrlhffihknhn��W45?KSU�|fOG=4Q��b|j)(*@X'=cT>2ildlQhisx}~NhJ=;,B41/67159A=4Fcrx|��пz22MTKIJCa��zrlea\draZ��X4CILTVV��p]F1O��yF)('MhO_[L?ILXkNI]r}�|_`Ok>6/E8516<64;BF6Ebrw}}��Δ@<MDFI=g��~yrlbnfc��ʚN5KPNX]WS��sQ2M��d<.--=>TP?BOS6dF@G[d}��RhhkM/0D8646>66@FV9?bu{x}��ԜFB;DH9m��}yqv�m��ͺo2JZUWdbXYo��e6L��i.-36C7:7E=175VJ?GV`o��uf�`d.1D;636@86CLfD9dwuz��֒<1?C9v��}z|�y��ǞUDZ[\embZ`Ep��u:J��I)*-8<6.5YA0:B9UCHL\lz��rXp~0*@:404>50=QrO6^r|��s-896~��~��Ҽ|LW]^dng^af*R��>C��R)+0@320=`=6*G?INU]S`l}�x3e�9&==4/3A628YtT2Zu��ݽC116��|��Ϫ_R[_fmg[_fg#>��:=��{P+,67.0,Dd6=(3KKBTgaRe��b4h�:'?D6/4<13=alV5]x��})'9��|{��ƌZ\]dlk^Zcfe#9{��7<��L.061-1-Yf-=3675=OcaO_��t5a�1-EF824;03Fdn_9^z��ݸ8!@��z��аp[\ckqe[`fdg/t��99��}8260,-*4r]+2:7014DWjDZ��vby�W.C>853659Tq{b=b{��`J��~z��Ɠc[\gpn^[egec+m��?7}��m88.,,-1S]M+2?0)53=EjYk��{��]B7::97636Eh|{\Be{��ݏS��{|��ҳ}d_^jpf\aheaY)j��I7{��s@0,-/6B7JQ15=/,74F:@k��n`|�iZOC@A761:Pq}|VEj~��۷0Z��|��Țwibckna\dgb[S&a��L7y��w=./-1372CI5654-49P41_|{uP9DGK[ceQ843C\w~}RMp��Qb��~}��Է�tjbcjeZ]gc^Y[!%W��Q:x��t82.,/5<2;>7068)39O>L\nrmyJA-0:<@K</6Ofy�yKWx��zi��ϧzqj`abV[ee`^\_%$M��[?v��D3)/?E5.;:574;06CBOCjuuQdu^@347733.?Xnz~nEa}��ةy��ʚskc]bXPdg_Z^_\&#E��eBt��N6+8>2.2>;9=792;E0UHAOmMEw�G1062/.1Gcs{{[Gm��ȏzkt|��oe\_[OZgbZ[_]O#"C��hCr��yP8+21.19<7C9/2*99:ak>.HKJcg12-2-)/<Rku|tFQ{��֚YQU^hr{��԰}eVY_NSdb^[c_S@!!9��sCo��rV:-2-.98:5@<0,/1;Id{_@@FWPH0-/-+*6IauuxZ<e��ۤPQMMOTW]el|��ǖiQN[TNab^\`e]H8%"6��Gm��iY8.,+6;3:3?>8/242/Lp�g^cJ-23--))->Vmvwh7Lx��ڼ\\]XRJC?=>Ol��ֳwM@NQQ_f[[[ccW?58=E��Lm��kP1//355176>;84-.531Q��kaY-0,-+()3JcqvqC8i��t^ge_VL?3*)I�ũ��͖T>DNJZi`VU]gaM;5VU[��Vm��qQ-284.2239<>:/,/2?52W~fd5.+*'(.CVnvrI3Y}��ؒYgjg`XM?0'G�Ǫ�ά^:?HHRfbROTdl]N@:W^l��`m��~nM/////5401<;=:-6<:;>EXjpH/--*&,;NbtkE1Su��ذ[bkic^VK@7N��h99DDNbdTKN]lgYL?@Mgr��dm��}kD-124365-7K<EJ.5L:F@HaiKF,-,''4JZ]O8=[q|��h\fgea_URR^��_@;EEM`gYIJWimbVDBA=[e{��hm��\;+14745=.8O@H_>2H3ML?^�P:+,,+,5>;9?Shu{}��ӀVaeged_bce��ZKEIIMaj\LBMcplaK@C<3EXq��qn��I2.11:89C08L8GlU/OA7WD^tg?,.-+,3<JWfrx|��מR[bhkhhoolvy_STUOJZg^MCC\ntlZE=<57<Ie��|p��zB7010;::F25O/A`f2M_,IRJi[@+/+-8L\hmtz|}��ַVVahmnrvsqke_df[PU`YMJFPkvqcN>:7357<]��s��l=43/399<?<8[0EZlS7]97[AP]A,*08K`jruu{}~��gR_gksxwtrmikrk]Z`\LGHOdwwiUA886=257\��šv��wB0138:>G<?GY3MWXj3VJ7LE?[bE'5IZgnswx{~~~��҂O\bgrvzytsu|yi^_YSSNOax}o[D53:EU023R��v��e>./48?BHBAcO6NQE_R[M.BD3HP_CGT^ipuy{|~~��ԡ][\dlw�~z{��vbYQLUXRXo|eG6/3BYd \ No newline at end of file

1 0

1851: Speed up edge insertion into graphs.
by Roland Levillain 09 Apr '08

09 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Roland Levillain <roland(a)lrde.epita.fr> Speed up edge insertion into graphs. * mln/util/internal/graph_base.hh: Stick to the coding style w.r.t. spaces. (mln::util::internal::graph_base<N, E>::edges_set_t): New typedef. (mln::util::internal::graph_base<N, E>::edges_set_): New member. (mln::util::internal::graph_base<N, E>::graph_base): Adjust ctor. (mln::util::internal::graph_base<N, E>::add_edge_): Use edges_set_ instead of edges_ to test the presence of an edge in the graph. * mln/util/graph.hh (mln::util::graph<void, void>::add_edge): Add missing preconditions. graph.hh | 2 ++ internal/graph_base.hh | 29 +++++++++++++++++++++-------- 2 files changed, 23 insertions(+), 8 deletions(-) Index: mln/util/graph.hh --- mln/util/graph.hh (revision 1850) +++ mln/util/graph.hh (working copy) @@ -162,6 +162,8 @@ void graph<void, void>::add_edge(node_id n1, node_id n2) { + mln_assertion(n1 < this->nnodes()); + mln_assertion(n2 < this->nnodes()); super::add_edge_(new util::edge<void>(n1, n2)); } Index: mln/util/internal/graph_base.hh --- mln/util/internal/graph_base.hh (revision 1850) +++ mln/util/internal/graph_base.hh (working copy) @@ -31,12 +31,15 @@ /// \file mln/util/internal/graph.hh /// \brief Factored implementation of undirected graphs. -# include <mln/core/concept/object.hh> # include <cstddef> -# include <ostream> -# include <vector> -# include <list> + # include <algorithm> + +# include <vector> +# include <set> +# include <ostream> + +# include <mln/core/concept/object.hh> # include <mln/util/ordpair.hh> // FIXME: Rename node(s) as vertex (vertices). @@ -156,6 +159,9 @@ typedef std::vector< node<N>* > nodes_t; /// The type of the set of edges. typedef std::vector< edge<E>* > edges_t; + /// A set to test the presence of a given edge. + typedef std::set< edge<E>* > edges_set_t; + /// Constructor. graph_base(); @@ -210,6 +216,8 @@ nodes_t nodes_; /// The edges. edges_t edges_; + /// An index of the set of edges, for fast-access purpose. + edges_set_t edges_set_; }; } // end of namespace mln::util::internal @@ -242,7 +250,7 @@ template<typename N, typename E> inline graph_base<N, E>::graph_base() - : nodes_(), edges_() + : nodes_(), edges_(), edges_set_() { } @@ -354,22 +362,27 @@ nodes_.push_back (node); } + // FIXME: Return the (new) edge id. template<typename N, typename E> inline void graph_base<N, E>::add_edge_(util::edge<E>* edge) { // Does this edge already exist in the graph? - if (std::find(edges_.begin(), edges_.end(), edge) != edges_.end ()) + if (edges_set_.find(edge) != edges_set_.end ()) { + // Remove the previously allocated data for EDGE. delete edge; edge = 0; } else { + // Otherwise insert it into the graph. edges_.push_back (edge); - nodes_[edge->n1()]->edges.push_back(edge->n2()); - nodes_[edge->n2()]->edges.push_back(edge->n1()); + edge_id id = edges_.size() - 1; + edges_set_.insert(edge); + nodes_[edge->n1()]->edges.push_back(id); + nodes_[edge->n2()]->edges.push_back(id); } }

1 0

1850: Delay construction of pw::value.
by Thierry Geraud 09 Apr '08

09 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk Index: ChangeLog from Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Delay construction of pw::value. * milena/mln/pw/value.hh (value): New overload. Such function shall have a delayed construction mechanism so that 'initialize' can work. value.hh | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) Index: milena/mln/pw/value.hh --- milena/mln/pw/value.hh (revision 1849) +++ milena/mln/pw/value.hh (working copy) @@ -54,8 +54,10 @@ typedef mln_value(I) result; value_(const I& ima); mln_rvalue(I) operator()(const mln_psite(I)& p) const; + + value_(); // A function should have a ctor without arg so that "initialize" can work. protected: - const I& ima_; + const I* ima_; }; @@ -74,7 +76,14 @@ template <typename I> inline value_::value_(const I& ima) - : ima_(ima) + : ima_(&ima) + { + } + + template <typename I> + inline + value_::value_() + : ima_(0) { } @@ -83,8 +92,9 @@ mln_rvalue(I) value_::operator()(const mln_psite(I)& p) const { - mln_precondition(ima_.owns_(p)); - return ima_(p); + mln_precondition(ima_ != 0); + mln_precondition(ima_->owns_(p)); + return (*ima_)(p); } // pw::value(ima)

1 0

milena r1849: Some stuff for tests
by Etienne FOLIO 09 Apr '08

09 Apr '08

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2008-04-09 Etienne FOLIO <folio(a)lrde.epita.fr> Some stuff for tests. * mln/trait/ch_value.hh: New specialization for functions. * sandbox/folio/psn.cc: Some test stuff. --- mln/trait/ch_value.hh | 7 +++++++ sandbox/folio/psn.cc | 33 +++++++++++++++++++++++---------- 2 files changed, 30 insertions(+), 10 deletions(-) Index: trunk/milena/mln/trait/ch_value.hh =================================================================== --- trunk/milena/mln/trait/ch_value.hh (revision 1848) +++ trunk/milena/mln/trait/ch_value.hh (revision 1849) @@ -133,6 +133,13 @@ typedef typename image_from_mesh< mesh, V >::ret ret; }; + template < template <class, class> class M, typename I, typename F, + typename V > + struct ch_value_< M< tag::image_, tag::function_<F> >, V > + { + typedef M< mln_ch_value(I, V), F > ret; + }; + } // end of namespace mln::trait::impl Index: trunk/milena/sandbox/folio/psn.cc =================================================================== --- trunk/milena/sandbox/folio/psn.cc (revision 1848) +++ trunk/milena/sandbox/folio/psn.cc (revision 1849) @@ -182,6 +182,10 @@ #include <mln/level/stretch.hh> #include <mln/value/int_u8.hh> +#include <mln/core/sub_image.hh> +#include <mln/core/image_if.hh> +#include <mln/pw/value.hh> + int main() { using namespace mln; @@ -193,19 +197,28 @@ // 0, 0, 0, 0, 0, // 0, 0, 0, 0, 0 }; -// level::fill(ima, vals); -// debug::println(ima); - - image2d<bool> ima = io::pbm::load("../../img/c01.pbm"); + image2d<bool> ima(3,3); + bool vals[] = { 1, 0, 0, + 0, 0, 0, + 0, 0, 0}; + level::fill(ima, vals); + + image2d<bool> msk(3,3); + bool rest[] = { 1, 0, 1, + 1, 0, 1, + 1, 1, 1}; + level::fill(msk, rest); image2d<unsigned> out; - out = dt::psn(ima, c4()); + out = dt::psn(ima | pw::value(msk), c4()); + + debug::println(ima | pw::value(msk)); + debug::println(out); - image2d<value::int_u8> out2(out.domain()); - level::stretch(out, out2); +// image2d<bool> ima = io::pbm::load("../../img/c01.pbm"); - io::pgm::save(out2, "out.pgm"); +// image2d<value::int_u8> out2(out.domain()); +// level::stretch(out, out2); -// std::cerr << "Distance:" << std::endl; -// debug::println(out); +// io::pgm::save(out2, "out.pgm"); }

1 0

[Olena] #123: Revamp general graph-based images (aka badly named « mesh images »)
by Olena 09 Apr '08

09 Apr '08

#123: Revamp general graph-based images (aka badly named « mesh images ») ----------------------+----------------------------------------------------- Reporter: levill_r | Owner: levill_r Type: task | Status: new Priority: major | Milestone: Olena 1.0ß Component: Milena | Version: 1.0 Keywords: | ----------------------+----------------------------------------------------- For starters, * Use `graph_based` prefix instead of `mesh`. * Improve the graph structure; possibly provide a BGL-based version of the graph-based image. Then later: * Create a « add_neighborhood » morpher, and adjust window-based morpho algorithms so that they can work on images with neighborhoods. * Rename « mesh elementary window » as « graph elementary neighborhood » (or adjacency? Check formal definitions to choose the right term.). * Revamp this elementary neighborhood too (currently, the graph-based piters ignore them, and assume they always work on elementary neighborhoods). These latest tasks might deserve their own ticket, as their scope is beyond graph-based images. -- Ticket URL: <https://trac.lrde.org/olena/ticket/123> Olena <http://olena.lrde.epita.fr> Olena, a generic and efficient C++ image library.

1 7

1848: Sandbox: ICP: various tools and improvments.
by Ugo Jardonnet 08 Apr '08

08 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: various tools and improvments. * sandbox/jardonnet/test/icp.cc: Fix regression (output registred.pbm). * sandbox/jardonnet/registration/tools.hh: Add multiple 2d to 3d and 3d to 2d conversion . * sandbox/jardonnet/registration/quat7.hh: Rewrite Covariance matrix computation. * sandbox/jardonnet/registration/frankel_young.hh: New, resolve linear equation Ax = b. Useless for icp. * sandbox/jardonnet/registration/power_it.hh: Find greatest eigen value, should be more efficient than jacobi ?. * sandbox/jardonnet/registration/icp.hh: Add intermediate image outputs. mln/algebra/quat.hh | 2 sandbox/jardonnet/registration/cloud.hh | 2 sandbox/jardonnet/registration/frankel_young.hh | 46 ++++++++++ sandbox/jardonnet/registration/icp.hh | 14 ++- sandbox/jardonnet/registration/jacobi.hh | 6 - sandbox/jardonnet/registration/power_it.hh | 41 +++++++++ sandbox/jardonnet/registration/quat7.hh | 38 +++++++- sandbox/jardonnet/registration/tools.hh | 109 ++++++++++++++++-------- sandbox/jardonnet/test/icp.cc | 33 +++---- 9 files changed, 233 insertions(+), 58 deletions(-) Index: mln/algebra/quat.hh --- mln/algebra/quat.hh (revision 1847) +++ mln/algebra/quat.hh (working copy) @@ -186,7 +186,7 @@ quat conj() const; /// Give the invert. - quat inv() const; // FIXME: Rename inv as invert. + quat inv() const; //FIXME: rename invert. /// Transform into unit quaternion. quat& set_unit(); Index: sandbox/jardonnet/test/icp.cc --- sandbox/jardonnet/test/icp.cc (revision 1847) +++ sandbox/jardonnet/test/icp.cc (working copy) @@ -41,27 +41,30 @@ closest_point<mln_point_(I3d)> fun(x, box_<point3d>(1000,1000,1)); lazy_image< closest_point<mln_point_(I3d)> > map(fun); - registration::icp(c, map); + quat7<3> q = registration::icp(c, map); #ifndef NDEBUG std::cout << "closest_point(Ck[i]) = " << fun.i << std::endl; std::cout << "Pts processed = " << registration::pts << std::endl; #endif - // //init output image - // //FIXME: Should be - // //mln_concrete(I) output(res.bbox()) = convert::to_image(res)? - // const box_<point2d> box2d(1000,1000,0); - // image2d<bool> output(box2d, 1); - - // //to 2d : projection (FIXME:if 3d) - // for (size_t i = 0; i < res.npoints(); i++) - // { - // point2d p(res[i][0], res[i][1]); - // if (output.has(p)) - // output(p) = true; - // } + //init output image + //FIXME: Should be + //mln_concrete(I) output(res.bbox()) = convert::to_image(res)? + + q.apply_on(c, c); + + const box_<point2d> box2d(1000,1000); + image2d<bool> output(box2d, 1); + + //to 2d : projection (FIXME:if 3d) + for (size_t i = 0; i < c.npoints(); i++) + { + point2d p(c[i][0], c[i][1]); + if (output.has(p)) + output(p) = true; + } - //io::pbm::save("./+registred.pbm"); + io::pbm::save(output, "registred.pbm"); } Index: sandbox/jardonnet/registration/jacobi.hh --- sandbox/jardonnet/registration/jacobi.hh (revision 1847) +++ sandbox/jardonnet/registration/jacobi.hh (working copy) @@ -17,7 +17,7 @@ #define rotateJacobi(a,i,j,k,l) g=a(i,j);h=a(k,l);a(i,j)=g-s*(h+g*tau); \ a(k,l)=h+s*(g-h*tau); - void jacobi(algebra::mat<4,4,float> a, algebra::quat& q) + algebra::quat jacobi(algebra::mat<4,4,float> a) { float dd, d[4]; algebra::mat<4,4,float> v(literal::zero); @@ -46,12 +46,12 @@ iq = ip; dd = d[ip]; } - q = algebra::quat(v(0,iq), + algebra::quat q(v(0,iq), v(1,iq), v(2,iq), v(3,iq)); q.set_unit(); - return; + return q; } if (i < 4) { i++; Index: sandbox/jardonnet/registration/power_it.hh --- sandbox/jardonnet/registration/power_it.hh (revision 0) +++ sandbox/jardonnet/registration/power_it.hh (revision 0) @@ -0,0 +1,41 @@ +#ifndef POWER_IT_HH +# define POWER_IT_HH + + +#include <mln/algebra/mat.hh> +#include <cmath> +#include "misc.hh" + +// from num. rec. in C + + +namespace mln +{ + + algebra::quat power_it(algebra::mat<4,4,float> A) + { + float normex = 1.; + + algebra::vec<4,float> x0(literal::zero); + algebra::vec<4,float> b(literal::zero); + + algebra::vec<4,float> x(literal::zero); + for (int i = 0; i < 4; i++) + x[i] = 1.; + + while (fabs(norm::l2(x) - norm::l2(x0)) > 1e-6) + { + x0 = x; + normex = norm::l2(x); + b = x / normex; + x = A * b; + } + normex = norm::l2(x); + + algebra::quat q(x / normex); + q.set_unit(); + return q; + } +} + +#endif /* POWER_IT_HH */ Index: sandbox/jardonnet/registration/tools.hh --- sandbox/jardonnet/registration/tools.hh (revision 1847) +++ sandbox/jardonnet/registration/tools.hh (working copy) @@ -17,12 +17,17 @@ typedef P result; closest_point(const p_array& X, const box_& box) - : X(X), box(box), i(0) + : X(X), box(box) + +#ifndef NDEBUG + , i(0) +#endif + { } result //inline - operator () (const P& Ck) + operator () (const P& Ck) const { #ifndef NDEBUG @@ -45,14 +50,17 @@ return algebra::to_point(best_x); } - const box_& domain() + const box_& domain() const { return box; } const p_array& X; const box_ box; + +#ifndef NDEBUG mutable unsigned i; +#endif }; @@ -74,7 +82,7 @@ const mln_result(F) //inline - operator() (const typename F::input& p) + operator () (const typename F::input& p) const { mln_precondition(fun.domain().has(p)); //FIXME: What about domain? @@ -86,31 +94,12 @@ } //FIXME: 3d -> //mln_dim(ml_input(input)) - image3d<mln_result(F)> value; - image3d<bool> is_known; + mutable image3d<mln_result(F)> value; + mutable image3d<bool> is_known; F& fun; }; -// // Point - -// template <typename P> -// P min(const P& a, const P& b) -// { -// if (a > b) -// return b; -// return a; -// } - -// template <typename P> -// P max(const P& a, const P& b) -// { -// if (a < b) -// return b; -// return a; -// } - - // Box template <typename P> @@ -142,6 +131,7 @@ namespace convert { + // to_p_array template <typename I> inline p_array<mln_point(I)> @@ -158,23 +148,74 @@ return a; } + template < typename P > + inline + box_<point2d> + to_box2d(const box_& b) + { + point2d pmin(b.pmin()[0], b.pmin()[1]); + point2d pmax(b.pmax()[0], b.pmax()[1]); + + return box_<point2d>(pmin, pmax); + } + + //FIXME: to write - /* - template <typename A> + template <typename P> inline - image2d<mln_value(A)> - to_image2d(const A& a) + image2d<bool> + to_image2d(const p_array& a) { - image2d<mln_value(A)> img(a.bbox()); + image2d<bool> output (to_box2d(a.bbox()), 0); for (size_t i = 0; i < a.npoints(); i++) { - point2d p(res[i][0], res[i][1]); - //FIXME: BOF - //if (output.has(p)) + point2d p(a[i][0], a[i][1]); + if (output.has(p)) output(p) = true; } + return output; + } + + // to_pointNd + + //FIXME: Should we really provide this + //point3d -> point2d + template <typename T> + inline + point_<grid::square, T> + to_point2d(const point_<grid::cube, T>& p) + { + return point_<grid::square, T>(p[0], p[1]); + } + //point2d -> point2d + template <typename T> + inline + point_<grid::square, T>& + to_point2d(const point_<grid::square, T>& p) + { + return p; + } + + //point2d -> point3d + template <typename T> + inline + point_<grid::cube, T> + to_point3d(const point_<grid::square, T>& p) + { + return point_<grid::cube, T>(p[0], p[1], 0); + } + //point3d -> point3d + template <typename T> + inline + point_<grid::cube, T>& + to_point3d(const point_<grid::cube, T>& p) + { + return p; } - */ + + + // to_imageNd + //const return a const template <typename T> Index: sandbox/jardonnet/registration/quat7.hh --- sandbox/jardonnet/registration/quat7.hh (revision 1847) +++ sandbox/jardonnet/registration/quat7.hh (working copy) @@ -10,6 +10,8 @@ # include "rotation.hh" # include "jacobi.hh" +# include "power_it.hh" +# include "frankel_young.hh" # include <mln/norm/l2.hh> @@ -115,6 +117,8 @@ } Mk /= C.npoints(); + + /* const algebra::mat<P::dim,P::dim,float> Ak = Mk - trans(Mk); const float v[3] = {Ak(1,2), Ak(2,0), Ak(0,1)}; @@ -127,14 +131,42 @@ put(D, 1,0, Qk); put(Mk + trans(Mk) - algebra::mat<P::dim,P::dim,float>::identity() * tr(Mk), 1,1, Qk); + */ - algebra::quat qR(literal::zero); + algebra::vec<3,float> a; + a[0] = Mk(1,2) - Mk(2,1); + a[1] = Mk(2,0) - Mk(0,2); + a[2] = Mk(0,1) - Mk(1,0); + + algebra::mat<4,4,float> Qk(literal::zero); + float t = tr(Mk); + + 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 * Mk(i - 1,i - 1) - t; + } + Qk(1,2) = Mk(0,1) + Mk(1,0); + Qk(2,1) = Mk(0,1) + Mk(1,0); + + Qk(1,3) = Mk(0,2) + Mk(2,0); + Qk(3,1) = Mk(0,2) + Mk(2,0); + + Qk(2,3) = Mk(1,2) + Mk(2,1); + Qk(3,2) = Mk(1,2) + Mk(2,1); + + algebra::quat qR(literal::zero); + qR = jacobi(Qk); + //std::cout << qR << std::endl; + //qR = power_it(Qk); //std::cout << qR << std::endl; - jacobi(Qk, qR); // qT - const algebra::vec<P::dim,float> qT = mu_Xk - rotate(qR, mu_C); return quat7<P::dim>(qR, qT); Index: sandbox/jardonnet/registration/frankel_young.hh --- sandbox/jardonnet/registration/frankel_young.hh (revision 0) +++ sandbox/jardonnet/registration/frankel_young.hh (revision 0) @@ -0,0 +1,46 @@ +#ifndef _FRANKEL_YOUNG_HH +# define _FRANKEL_YOUNG_HH + +//Successive over-relaxation (SOR) is a numerical method used to speed up +//convergence of the Gauss-Seidel method for solving a linear system of +//equations. + +namespace mln +{ + + algebra::quat frankel_young(algebra::mat<4,4,float> a, float w) + { + //Inputs: A , b, w + //Output: phi + + algebra::vec<4,float> phi(literal::zero); + algebra::vec<4,float> old_phi(literal::zero); + + algebra::vec<4,float> b; + for (int i = 0; i < 4; i++) + b[i] = phi[i] = a(i,i); + + //Choose an initial guess phi to the solution + while (fabs(norm::l2(old_phi) - norm::l2(phi)) > 1e-6) + { + old_phi = phi; + for (int i = 1; i < 4; i++) + { + float sigma = 0; + + for (int j = 1; j < i-1; j++) + sigma += a(i,j) * phi[i]; + + for (int j = i + 1; j < 4; j++) + sigma += a(i,j) * phi[i]; + + phi[i] = (1 - w) * phi[i] + w / a(i,i) * (b[i] - sigma); + } + } + algebra::quat q(phi); + q.set_unit(); + return q; + } +} + +#endif /* _FRANKEL_YOUNG_HH */ Index: sandbox/jardonnet/registration/cloud.hh --- sandbox/jardonnet/registration/cloud.hh (revision 1847) +++ sandbox/jardonnet/registration/cloud.hh (working copy) @@ -56,7 +56,7 @@ { algebra::vec<3,float> a1f = a1[i]; algebra::vec<3,float> a2f = a2[i]; - f += norm::l2_distance(a1f,a2f); + f += norm::l2(a1f - a2f); } return f / a1.npoints(); } Index: sandbox/jardonnet/registration/icp.hh --- sandbox/jardonnet/registration/icp.hh (revision 1847) +++ sandbox/jardonnet/registration/icp.hh (working copy) @@ -33,6 +33,9 @@ * \brief image registration */ +# include <iostream> +# include <string> + # include <mln/algebra/quat.hh> # include <mln/algebra/vec.hh> # include <mln/make/w_window.hh> @@ -111,6 +114,15 @@ err = rms(Ck, Xk); #ifndef NDEBUG + + { + using namespace std; + image2d<bool> img = convert::to_image2d(Ck); + stringstream oss; + oss << "reg" << k << ".pbm"; + io::pbm::save(img, oss.str()); + } + //plot file std::cout << k << '\t' << err << '\t' << (qk - old_qk).sqr_norm() << '\t' @@ -138,8 +150,8 @@ { trace::entering("registration::icp"); - mln_precondition(cloud.npoints() != 0); mln_precondition(P::dim == 3); + mln_precondition(cloud.npoints() != 0); //init rigid transform qk quat7<P::dim> qk;

1 0

milena r1846: Add pixter for p_image2d
by Matthieu Garrigues 07 Apr '08

07 Apr '08

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2008-04-07 Matthieu Garrigues <garrigues(a)lrde.epita.fr> Add pixter for p_image2d. * mln/core/p_image2d_pixter.hh: New, a pixter to optimize iteration on p_image2d. * mln/core/p_image2d.hh: Update, optimize clear and the number of update of bb_. * tests/core/p_image2d.cc: More tests on the clean method. --- mln/core/p_image2d.hh | 53 +++++++++- mln/core/p_image2d_pixter.hh | 208 +++++++++++++++++++++++++++++++++++++++++++ tests/core/p_image2d.cc | 14 ++ 3 files changed, 270 insertions(+), 5 deletions(-) Index: trunk/milena/tests/core/p_image2d.cc =================================================================== --- trunk/milena/tests/core/p_image2d.cc (revision 1845) +++ trunk/milena/tests/core/p_image2d.cc (revision 1846) @@ -36,6 +36,7 @@ { using namespace mln; + trace::quiet = false; p_image2d<point2d> ps(20,20); ps .insert(make::point2d(6, 9)) @@ -53,4 +54,17 @@ mln_assertion(ps.npoints() == 0); mln_assertion(ps.is_empty()); + std::cout << ps << std::endl; + + mln_fwd_piter_(box2d) p(inplace(make::box2d(13,13,19,15))); + for_all(p) + { + ps.insert(p); + } + ps.clear(); + for_all(p) + { + mln_assertion(!ps.has(p)); + } + } Index: trunk/milena/mln/core/p_image2d_pixter.hh =================================================================== --- trunk/milena/mln/core/p_image2d_pixter.hh (revision 0) +++ trunk/milena/mln/core/p_image2d_pixter.hh (revision 1846) @@ -0,0 +1,208 @@ +// Copyright (C) 2007 EPITA Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +#ifndef MLN_CORE_P_IMAGE2D_PIXTER_HH +# define MLN_CORE_P_IMAGE2D_PIXTER_HH + +/*! \file mln/core/p_image2d_pixter.hh + * + * \brief Pixel iterator class on a image 2d with border. + */ + +# include <mln/core/internal/pixel_iterator_base.hh> +# include <mln/core/point2d.hh> +# include <mln/geom/size2d.hh> + + + +namespace mln +{ + + template <typename P> + class p_image2d_fwd_pixter : public internal::pixel_iterator_base_< typename p_image2d::image_type, p_image2d_fwd_pixter > + { + typedef internal::pixel_iterator_base_< typename p_image2d::image_type, p_image2d_fwd_pixter > super_; + + public: + + /// Image type. + typedef typename p_image2d::image_type image; + + /// Point type. + typedef point2d point; + + /*! \brief Constructor. + * + * \param[in] image Image to iterate over its pixels. + */ + p_image2d_fwd_pixter(p_image2d& s); + + /// Reference of the corresponding point. + const point to_point() const; + + /// Go to the next pixel. + void next_(); + + void start(); + + private: + + /// Row offset. + unsigned row_offset_; + /// End of the current row. + mln_qlf_value(image)* eor_; + }; + + + + template <typename P> + class p_image2d_bkd_pixter : public internal::pixel_iterator_base_< typename image2d::image_type, p_image2d_bkd_pixter > + { + typedef internal::pixel_iterator_base_< typename image2d::image_type, p_image2d_bkd_pixter > super_; + + public: + + /// Image type. + typedef typename p_image2d::image_type image; + + /// Point type. + typedef point2d point; + + /*! \brief Constructor. + * + * \param[in] image Image to iterate over its pixels. + */ + p_image2d_bkd_pixter(p_image2d& s); + + /// Go to the next pixel. + void next_(); + + void start(); + + /// Reference of the corresponding point. + const point to_point() const; + + private: + + /// Row offset. + unsigned row_offset_; + + /// Beginning of the current row. + mln_qlf_value(image)* bor_; + }; + + + + +#ifndef MLN_INCLUDE_ONLY + + // Fwd. + + template <typename P> + inline + p_image2d_fwd_pixter::p_image2d_fwd_pixter(p_image2d& s) : + super_(s.image_non_const()) + { + mln_precondition(this->image_.has_data()); + row_offset_ = geom::max_col(this->image_) - geom::min_col(this->image_) + 1; + eor_ = & this->image_.at(geom::min_row(this->image_), geom::max_col(this->image_)) + 1; + } + + template <typename P> + inline + const typename p_image2d_fwd_pixter::point + p_image2d_fwd_pixter::to_point() const + { + return this->image_.point_at_offset(*this); + } + + template <typename P> + inline + void + p_image2d_fwd_pixter::next_() + { + ++this->value_ptr_; + while(this->is_valid() && !(*this->value_ptr_)) + ++this->value_ptr_; + } + + + template <typename P> + inline + void + p_image2d_fwd_pixter::start() + { + this->value_ptr_ = this->boi_ + 1; + while(this->is_valid() && !(*this->value_ptr_)) + ++this->value_ptr_; + } + + // Bkd. + + template <typename P> + inline + p_image2d_bkd_pixter::p_image2d_bkd_pixter(p_image2d& s) : + super_(s.image_non_const()) + { + mln_precondition(this->image_.has_data()); + row_offset_ = geom::max_col(this->image_) - geom::min_col(this->image_) + 1; + bor_ = & this->image_.at(geom::max_row(this->image_), geom::min_col(this->image_)) - 1; + } + + template <typename P> + inline + const typename p_image2d_bkd_pixter::point + p_image2d_bkd_pixter::to_point() const + { + return this->image_.point_at_offset(*this); + } + + template <typename P> + inline + void + p_image2d_bkd_pixter::next_() + { + --this->value_ptr_; + while(this->is_valid() && !(*this->value_ptr_)) + --this->value_ptr_; + } + + template <typename P> + inline + void + p_image2d_bkd_pixter::start() + { + this->value_ptr_ = this->eoi_ - 1; + while(this->is_valid() && !(*this->value_ptr_)) + --this->value_ptr_; + } + +#endif // ! MLN_INCLUDE_ONLY + +} // end of namespace mln + +#endif // ! MLN_CORE_P_IMAGE2D_PIXTER_HH Index: trunk/milena/mln/core/p_image2d.hh =================================================================== --- trunk/milena/mln/core/p_image2d.hh (revision 1845) +++ trunk/milena/mln/core/p_image2d.hh (revision 1846) @@ -37,8 +37,14 @@ # include <mln/core/internal/point_set_base.hh> # include <mln/core/box2d.hh> # include <mln/core/image2d.hh> +# include <mln/core/sub_image.hh> + +# include <mln/accu/bbox.hh> +# include <mln/geom/ncols.hh> +# include <mln/accu/bbox.hh> # include <mln/level/fill.hh> +# include <mln/level/memset_.hh> namespace mln { @@ -53,6 +59,8 @@ { public: + typedef image2d<bool> image_type; + /// Forward Point_Iterator associated type. typedef p_image2d_fwd_piter_ fwd_piter; @@ -90,6 +98,7 @@ /// Hook to the image2d const image2d<bool>& image() const; + image2d<bool>& image_non_const(); private: image2d<bool> points_; unsigned npoints_; @@ -102,11 +111,10 @@ template <typename P> p_image2d::p_image2d(int nrows, int ncols) - : points_(nrows, ncols), + : points_(nrows, ncols, 0), npoints_(0), bb_need_update_(false) { - level::fill(points_, false); } @@ -153,8 +161,16 @@ if (points_(p) == true) { points_(p) = false; - bb_need_update_ = true; npoints_--; + + if (npoints_ == 0) + { + bb_.init(); + bb_need_update_ = false; + } + else + bb_need_update_ = true; + } return *this; } @@ -171,6 +187,7 @@ for_all(p) if (this->points_.has(p)) this->remove(p); + return *this; } @@ -205,7 +222,6 @@ if (bb_need_update_) { bb_.init(); - mln_fwd_piter(p_image2d) p(*this); for_all(p) bb_.take(p); @@ -215,11 +231,30 @@ return bb_.to_result(); } + template <typename P> void p_image2d::clear() { - level::fill(points_, false); + if (npoints_ == 0) + return; + + unsigned bb_nrows = geom::nrows(bb_.to_result()); + unsigned ima_nrows = geom::nrows(points_); + + if (bb_nrows * 3 < ima_nrows * 2) + { + unsigned bb_ncols = geom::ncols(bb_.to_result()); + mln_line_piter_(image2d<bool>) p(bb_.to_result()); + for_all(p) + { + level::memset_(points_, p, false, bb_ncols); + } + } + else + level::fill(inplace(points_), false); + + npoints_ = 0; bb_.init(); bb_need_update_ = false; } @@ -231,12 +266,20 @@ return points_; } + template <typename P> + image2d<bool>& + p_image2d::image_non_const() + { + return points_; + } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln # include <mln/core/p_image2d_piter.hh> +# include <mln/core/p_image2d_pixter.hh> #endif // ! MLN_CORE_P_IMAGE2D_HH

1 0