- Olena-patches

1946: Sandbox: ICP: icp_ref + cleanup.
by Ugo Jardonnet 14 May '08

14 May '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Ugo Jardonnet <jardonnet(a)lrde.epita.fr> Sandbox: ICP: icp_ref + cleanup. icp_ref uses no rounding. * jardonnet/test/icp_ref.cc: New. * jardonnet/registration/icp_ref.hh: New. * jardonnet/test/icp.cc: clean up. * jardonnet/test/Makefile: Update. * jardonnet/registration/final_qk.hh: Put final_qk function in a dedicated file. * jardonnet/registration/tools.hh: Add buffer and shuffle. * jardonnet/registration/icp.hh: Minor fix. registration/cloud.hh | 2 registration/final_qk.hh | 107 +++++++++++++++++++ registration/icp.hh | 9 - registration/icp_ref.hh | 254 +++++++++++++++++++++++++++++++++++++++++++++++ registration/quat7.hh | 2 registration/tools.hh | 116 +++++---------------- test/Makefile | 4 test/icp.cc | 12 +- test/icp_ref.cc | 87 ++++++++++++++++ 9 files changed, 495 insertions(+), 98 deletions(-) Index: jardonnet/test/icp_ref.cc --- jardonnet/test/icp_ref.cc (revision 0) +++ jardonnet/test/icp_ref.cc (revision 0) @@ -0,0 +1,87 @@ +#include <mln/core/image3d.hh> + +#include <mln/io/pbm/load.hh> +#include <mln/io/pbm/save.hh> +#include <mln/io/ppm/save.hh> +#include <mln/norm/l2.hh> + +#include <sandbox/jardonnet/registration/icp_ref.hh> +#include <sandbox/jardonnet/registration/tools.hh> + +void usage(char *argv[]) +{ + std::cout << "usage : " << argv[0] + << " cloud.pbm surface.pbm q e" << std::endl + << " q >= 1 and e >= 1" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 5) + usage(argv); + + float q = std::atof(argv[3]); + int e = std::atoi(argv[4]); + + if (q < 1 or e < 1) + usage(argv); + + using namespace mln; + typedef image3d< bool > I3d; + + image2d< bool > img1; + image2d< bool > img2; + + //load images + io::pbm::load(img1, argv[1]); + io::pbm::load(img2, argv[2]); + + //convert to image3d + I3d cloud = convert::to_image3d(img1); + I3d surface = convert::to_image3d(img2); + + //build p_arrays. + p_array<mln_point_(I3d)> c = convert::to_p_array(cloud); + p_array<mln_point_(I3d)> x = convert::to_p_array(surface); + + // FIXME : TODO : map : vec<3,float> -> point + + quat7<3> qk = registration::icp(c, map, q, e, x); + +#ifndef NDEBUG + std::cout << "closest_point(Ck[i]) = " << fun.i << std::endl; + std::cout << "Pts processed = " << registration::pts << std::endl; +#endif + + qk.apply_on(c, c, c.npoints()); + + //init output image + //FIXME: Should be like + //mln_concrete(I) output(res.bbox()) = convert::to_image(res) ? + + const box_<point2d> box2d(400,700); + image2d<value::rgb8> output(box2d, 1); + + //to 2d : projection (FIXME:if 3d) + for (size_t i = 0; i < c.npoints(); i++) + { + //Ck points + point2d p(c[i][0], c[i][1]); + if (output.has(p)) + output(p) = literal::white; + //Xk points + } + + for (size_t i = 0; i < x.npoints(); i++) + { + point2d x(map(c[i])[0], map(c[i])[1]); + if (output.has(x)) + output(x) = literal::green; + } + + io::ppm::save(output, "registred.ppm"); + +} + Index: jardonnet/test/icp.cc --- jardonnet/test/icp.cc (revision 1945) +++ jardonnet/test/icp.cc (working copy) @@ -7,6 +7,7 @@ #include <sandbox/jardonnet/registration/icp.hh> #include <sandbox/jardonnet/registration/tools.hh> +#include <sandbox/jardonnet/registration/final_qk.hh> void usage(char *argv[]) { @@ -57,18 +58,16 @@ std::cout << "Pts processed = " << registration::pts << std::endl; #endif - //FIXME: register img1 + qk.apply_on(c, c, c.npoints()); //init output image //FIXME: Should be like //mln_concrete(I) output(res.bbox()) = convert::to_image(res) ? - qk.apply_on(c, c, c.npoints()); const box_<point2d> box2d(400,700); image2d<value::rgb8> output(box2d, 1); - float stddev, mean; registration::mean_stddev(c, map, mean, stddev); @@ -103,10 +102,15 @@ else output(p) = literal::white; } + } + + // print surface + for (size_t i = 0; i < c.npoints(); i++) + { //Xk points point2d x(map(c[i])[0], map(c[i])[1]); if (output.has(x)) - output(x) = value::rgb8(0,255,0); + output(x) = literal::green; } io::ppm::save(output, "registred.ppm"); Index: jardonnet/test/Makefile --- jardonnet/test/Makefile (revision 1945) +++ jardonnet/test/Makefile (working copy) @@ -19,6 +19,9 @@ icp: icp.cc +depend icp.o g++ icp.o -I../../.. -O3 -DNDEBUG -o 'icp' +icp_ref: + g++ icp_ref.cc -I../../.. -g -o 'icp_ref' + icp.o: g++ -c icp.cc -I../../.. -O3 -DNDEBUG @@ -30,6 +33,7 @@ rm -f sub gsub gau icpD icp rm -f log.dat registred.pbm rm -f i_* + rm -f tmp.ppm registred.ppm semantic.cache .PHONY: check bench icpD Index: jardonnet/registration/final_qk.hh --- jardonnet/registration/final_qk.hh (revision 0) +++ jardonnet/registration/final_qk.hh (revision 0) @@ -0,0 +1,107 @@ +#ifndef MLN_REGISTRATION_FINAL_QK_HH +# define MLN_REGISTRATION_FINAL_QK_HH + +namespace mln +{ + + namespace registration + { + + template <typename P, typename M> + void mean_stddev(const p_array& c, + const M& map, + float& mean, + float& stddev) + { + //mean + length + std::vector<float> length(c.npoints()); + mean = 0; + + for (size_t i = 0; i < c.npoints(); i++) + { + float f = norm::l2(algebra::vec<3,int> (c[i] - map(c[i]))); + length[i] = f; + mean += f; + } + mean /= c.npoints(); + + //standar variation + stddev = 0; + for (size_t i = 0; i < c.npoints(); i++) + stddev += (length[i] - mean) * (length[i] - mean); + stddev /= c.npoints(); + stddev = math::sqrt(stddev); + } + + + //final qk : only use point less than nstddev (2*stddev); + template <typename P, typename M> + quat7<P::dim> + final_qk(const p_array& c, + const M& map, + float nstddev) + { + p_array newc; + algebra::vec<3,float> mu_newc(literal::zero); + + for (size_t i = 0; i < c.npoints(); ++i) + { + algebra::vec<3,float> xki = map(c[i]); + algebra::vec<3,float> ci = c[i]; + + if (norm::l2(ci - xki) > nstddev) + { + newc.append(c[i]); + mu_newc += ci; + } + } + mu_newc /= newc.npoints(); + + quat7<P::dim> qk = match(newc, mu_newc, newc, map, newc.npoints()); + + return qk; + } + + //final > nstddev for translation; < nstddev for rotation + template <typename P, typename M> + quat7<P::dim> + final_qk2(const p_array& c, + const M& map, + float nstddev) + { + //mu_Xk = center map(Ck) + algebra::vec<3,float> mu_Xk(literal::zero); + algebra::vec<3,float> mu_C(literal::zero); + + float nb_point = 0; + for (size_t i = 0; i < c.npoints(); ++i) + { + algebra::vec<3,float> xki = map(c[i]); + algebra::vec<3,float> ci = c[i]; + + if (norm::l2(ci - xki) > nstddev) + { + mu_C += ci; + mu_Xk += xki; + nb_point++; + } + } + mu_C /= nb_point; + mu_Xk /= nb_point; + + // qT + const algebra::vec<3,float> qT = mu_C - mu_Xk; + + // qR + quat7<P::dim> qk = final_qk(c, map, nstddev); + qk._qT = qT; + + return qk; + } + + + } // end of namespace mln::registration + +} + +#endif // MLN_REGISTRATION_FINAL_QK_HH Index: jardonnet/registration/tools.hh --- jardonnet/registration/tools.hh (revision 1945) +++ jardonnet/registration/tools.hh (working copy) @@ -12,104 +12,46 @@ namespace mln { - namespace registration - { - - template <typename P, typename M> - void mean_stddev(const p_array& c, - const M& map, - float& mean, - float& stddev) + template <typename P> + void shuffle(p_array& a) { - //mean + length - std::vector<float> length(c.npoints()); - mean = 0; - - for (size_t i = 0; i < c.npoints(); i++) + for (unsigned int i = 0; i < a.npoints(); i++) { - float f = norm::l2(algebra::vec<3,int> (c[i] - map(c[i]))); - length[i] = f; - mean += f; + unsigned int r = rand() % a.npoints(); + P tmp; + tmp = a[i]; + a.hook_()[i] = a[r]; + a.hook_()[r] = tmp; } - mean /= c.npoints(); - - //standar variation - stddev = 0; - for (size_t i = 0; i < c.npoints(); i++) - stddev += (length[i] - mean) * (length[i] - mean); - stddev /= c.npoints(); - stddev = math::sqrt(stddev); } - - //final qk : only use point less than nstddev (2*stddev); - template <typename P, typename M> - quat7<P::dim> - final_qk(const p_array& c, - const M& map, - float nstddev) + template <unsigned int n, typename T> + struct buffer { - p_array newc; - algebra::vec<3,float> mu_newc(literal::zero); - - for (size_t i = 0; i < c.npoints(); ++i) - { - algebra::vec<3,float> xki = map(c[i]); - algebra::vec<3,float> ci = c[i]; - - if (norm::l2(ci - xki) > nstddev) + buffer() + : setted(0) { - newc.append(c[i]); - mu_newc += ci; } - } - mu_newc /= newc.npoints(); - - quat7<P::dim> qk = match(newc, mu_newc, newc, map, newc.npoints()); - return qk; - } - - //final > nstddev for translation; < nstddev for rotation - template <typename P, typename M> - quat7<P::dim> - final_qk2(const p_array& c, - const M& map, - float nstddev) - { - //mu_Xk = center map(Ck) - algebra::vec<3,float> mu_Xk(literal::zero); - algebra::vec<3,float> mu_C(literal::zero); - - float nb_point = 0; - for (size_t i = 0; i < c.npoints(); ++i) + void add(T e) { - algebra::vec<3,float> xki = map(c[i]); - algebra::vec<3,float> ci = c[i]; + for (int i = 0; i < n-1; i++) + buf[i+1] = buf[i]; + buf[0] = e; - if (norm::l2(ci - xki) > nstddev) - { - mu_C += ci; - mu_Xk += xki; - nb_point++; - } + setted++; } - mu_C /= nb_point; - mu_Xk /= nb_point; - - // qT - const algebra::vec<3,float> qT = mu_C - mu_Xk; - // qR - quat7<P::dim> qk = final_qk(c, map, nstddev); - qk._qT = qT; - - return qk; + T& operator[](unsigned int i) + { + assert(i < n && i < setted); + return buf[i]; } - - } // end of namespace mln::registration - + private: + T buf[n]; + unsigned int setted; + }; //FIXME: groe length template <typename P> @@ -127,7 +69,7 @@ result //inline - operator () (const P& Ck) const + operator () (const input& Ck) const { #ifndef NDEBUG @@ -170,13 +112,13 @@ struct lazy_image { // Fun is potentially an image. - lazy_image(F& fun) + lazy_image(const F& fun) : value(fun.domain()), is_known(fun.domain()), fun(fun) { } // FIXME: hack, remove this constructor - lazy_image(F& fun, int nrows, int ncols, int nslis) + lazy_image(const F& fun, int nrows, int ncols, int nslis) : value(nrows, ncols,1), is_known(nrows,ncols,1), fun(fun) { } @@ -197,7 +139,7 @@ //FIXME: 3d -> //mln_dim(ml_input(input)) mutable image3d<mln_result(F)> value; mutable image3d<bool> is_known; - F& fun; + const F& fun; }; Index: jardonnet/registration/quat7.hh --- jardonnet/registration/quat7.hh (revision 1945) +++ jardonnet/registration/quat7.hh (working copy) @@ -126,7 +126,7 @@ quat7<P::dim> match(const p_array& C, const algebra::vec<P::dim,float>& mu_C, const p_array& Ck, - M& map, + const M& map, size_t c_length) { //mu_Xk = center map(Ck) Index: jardonnet/registration/cloud.hh --- jardonnet/registration/cloud.hh (revision 1945) +++ jardonnet/registration/cloud.hh (working copy) @@ -46,7 +46,7 @@ template <typename P, typename M> float rms(const p_array& a1, quat7<P::dim>& qk, - M& map, + const M& map, const size_t length) { assert(length <= a1.npoints()); Index: jardonnet/registration/icp_ref.hh --- jardonnet/registration/icp_ref.hh (revision 0) +++ jardonnet/registration/icp_ref.hh (revision 0) @@ -0,0 +1,254 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +#ifndef MLN_REGISTRATION_ICP_HH +# define MLN_REGISTRATION_ICP_HH + +/*! \file mln/registration/icp.hh + * + * \brief image registration + */ + +# include <iostream> +# include <string> +# include <cmath> + +# include <mln/algebra/quat.hh> +# include <mln/algebra/vec.hh> +# include <mln/make/w_window.hh> +# include <mln/make/w_window3d.hh> + +# include <mln/value/rgb8.hh> +# include <mln/literal/colors.hh> +# include <mln/literal/black.hh> +# include <mln/level/fill.hh> +# include <mln/io/ppm/save.hh> + + +# include "tools.hh" + +# include "cloud.hh" +# include "quat7.hh" +# include "update_qk.hh" +# include "chamfer.hh" + +# include "save.hh" + +namespace mln +{ + + namespace registration + { + +#ifndef NDEBUG + static unsigned pts = 0; +#endif + + /*! Registration FIXME : doxy + * + * + */ + template <typename I, typename J> + inline + mln_concrete(I) + icp(const Image& cloud, + const Image<J>& surface); + +# ifndef MLN_INCLUDE_ONLY + + namespace impl + { + + template <typename P, typename M> + inline + void + icp_(const p_array& C, + const M& map, + 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_length << " points" << std::endl; + std::cout << "k\t\te_k >=\td_k\tdqk" << std::endl; +#endif + + quat7<P::dim> buf_qk[4]; + float buf_dk[3]; + + float err = 100; + //float err_bis; + p_array Ck(C); + + algebra::vec<P::dim,float> mu_C = center(C, c_length), mu_Xk; + + buf_qk[0] = qk; + + qk.apply_on(C, Ck, c_length); + + unsigned int k = 0; + + do { + //update buff dk //FIXME: rewrite it + buf_dk[2] = buf_dk[1]; + buf_dk[1] = buf_dk[0]; + buf_dk[0] = err; + + //compute qk + qk = match(C, mu_C, Ck, map, c_length); + + //update buff qk //FIXME: rewrite it + buf_qk[3] = buf_qk[2]; + buf_qk[2] = buf_qk[1]; + buf_qk[1] = buf_qk[0]; + buf_qk[0] = qk; + + //update qk + /* + if (k > 3) + qk = update_qk(buf_qk, buf_dk); + qk._qR.set_unit(); + buf_qk[0] = qk; + */ + + //Ck+1 = qk(C) + qk.apply_on(C, Ck, c_length); + + // e_k = d(Yk, Pk) + // = d(closest(Pk), Pk) + // = d(closest(qk-1(P)), qk-1(P)) + float e_k = rms(C, map, c_length, buf_qk[1], buf_qk[1]); + + // d_k = d(Yk, Pk+1) + // = d(closest(qk-1(P)), qk(P)) + float d_k = rms(C, map, c_length, buf_qk[1], qk); + + + //err = d(Ck+1,Xk) = d_k + err = rms(C, qk, map, c_length); + + //err = d(Ck,Xk) = e_k + float err_bis = rms(C, buf_qk[1], map, c_length); + +#ifndef NDEBUG + //plot file + std::cout << k << '\t' << (e_k >= d_k ? ' ' : '-') << '\t' << e_k << '\t' << d_k << '\t' + << ((qk - buf_qk[1]).sqr_norm() / qk.sqr_norm()) << '\t' + << std::endl; + //count the number of points processed + pts += c_length; +#endif + k++; + } while ((qk - buf_qk[1]).sqr_norm() / qk.sqr_norm() > epsilon); + + trace::exiting("registration::impl::icp_"); + } + + } // end of namespace mln::registration::impl + + + // Only for 3d images + template <typename P, typename M> + inline + quat7<P::dim> + icp(p_array cloud, // FIXME : is almost const (shuffled) + const M& map, + const float q, + const unsigned nb_it, + const p_array& x) + { + trace::entering("registration::icp"); + + mln_precondition(P::dim == 3); + mln_precondition(cloud.npoints() != 0); + + shuffle(cloud); + + //init rigid transform qk + quat7<P::dim> qk; + +#ifndef NDEBUG // FIXME: theo + image2d<value::rgb8> tmp(500,800); + level::fill(tmp, literal::black); + //write X + mln_piter(p_array) p(x); + for_all(p) + { + point2d qp = make::point2d(p[0], p[1]); + if (tmp.has(qp)) + tmp(qp) = literal::white; + } +#endif + + + //run icp + for (int e = nb_it-1; e >= 0; e--) + { + + size_t l = cloud.npoints() / std::pow(q, e); + l = (l<1) ? 1 : l; + impl::icp_(cloud, map, qk, l, 1e-3); + +#ifndef NDEBUG + { + value::rgb8 c; + switch (e) { + case 2: c = literal::green; break; + case 1: c = literal::blue; break; + case 0: c = literal::red; break; + } + mln_piter(p_array) p(cloud); + for_all(p) + { + algebra::vec<3,float> p_ = point3d(p), qp_ = qk(p_); + point2d qp = make::point2d(int(qp_[0]), int(qp_[1])); + if (tmp.has(qp)) + tmp(qp) = c; + } + if (e == 0) + io::ppm::save(tmp, "tmp.ppm"); + } +#endif + } + + trace::exiting("registration::icp"); + + return qk; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::registration + +} // end of namespace mln + + +#endif // ! MLN_REGISTRATION_ICP_HH Index: jardonnet/registration/icp.hh --- jardonnet/registration/icp.hh (revision 1945) +++ jardonnet/registration/icp.hh (working copy) @@ -85,9 +85,9 @@ template <typename P, typename M> inline - p_array + void icp_(const p_array& C, - M& map, + const M& map, quat7<P::dim>& qk, const size_t c_length, const float epsilon = 1e-3) @@ -171,7 +171,6 @@ } while (/*k < 3 ||*/ (qk - buf_qk[1]).sqr_norm() / qk.sqr_norm() > epsilon); trace::exiting("registration::impl::icp_"); - return Ck; } } // end of namespace mln::registration::impl @@ -181,8 +180,8 @@ template <typename P, typename M> inline quat7<P::dim> - icp(p_array cloud, - M& map, + icp(p_array cloud, //here reference implies low efficiency (FIXME:check again) + const M& map, const float q, const unsigned nb_it, const p_array& x)

1 0

Improve FLLT : Make the two level set trees and start to merge it.
by Matthieu Garrigues 14 May '08

14 May '08

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2008-05-14 Matthieu Garrigues <garrigues(a)lrde.epita.fr> * geraud/fllt.cc: Make the two level set trees and start to merge it. --- fllt.cc | 313 +++++++++++++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 272 insertions(+), 41 deletions(-) Index: trunk/milena/sandbox/geraud/fllt.cc =================================================================== --- trunk/milena/sandbox/geraud/fllt.cc (revision 1944) +++ trunk/milena/sandbox/geraud/fllt.cc (revision 1945) @@ -46,6 +46,7 @@ #include <mln/level/compare.hh> #include <mln/debug/println.hh> #include <mln/labeling/regional_minima.hh> +#include <mln/labeling/regional_maxima.hh> #include <mln/accu/bbox.hh> #include <mln/geom/bbox.hh> #include <mln/pw/all.hh> @@ -55,6 +56,7 @@ #include <mln/literal/colors.hh> #include <mln/util/tree.hh> +#include <mln/util/branch_iter_ind.hh> #include <sstream> @@ -67,9 +69,16 @@ # define fllt_tree(P, V) mln::util::tree< fllt_node_elt<P, V> > # define fllt_node(P, V) mln::util::tree_node< fllt_node_elt<P, V> > +# define fllt_tree_ptr(P, V) mln::util::tree< fllt_node_elt<P, V>* > +# define fllt_node_ptr(P, V) mln::util::tree_node< fllt_node_elt<P, V>* > # define fllt_branch(P, V) mln::util::branch< fllt_node_elt<P, V> > # define fllt_branch_iter_ind(P, V) mln::util::branch_iter_ind< fllt_node_elt<P, V> > + + //Fwd declarations. + template <typename V> struct lower; + template <typename V> struct upper; + template <typename P, typename V> struct fllt_node_elt { @@ -80,10 +89,16 @@ /// if brighter, the node come from the lower level set bool brighter; unsigned npoints; + bool tagged; + + fllt_node_elt() : npoints(0), tagged(false) {} }; - template <typename N_t, typename G> - void update_gN(const N_t& N, G& gN) + template <typename N_t, typename G, typename Set> + void update_gN(const N_t& N, G& gN); + + template <typename N_t, typename G, typename V> + void update_gN(const N_t& N, G& gN, lower<V>) { for (unsigned g = 0; g < 256; ++g) if (N[g]->npoints() != 0) @@ -95,6 +110,21 @@ gN = 255; } + template <typename N_t, typename G, typename V> + void update_gN(const N_t& N, G& gN, upper<V>) + { + for (int g = 255; g >= 0; --g) + { + if (N[g]->npoints() != 0) + { + gN = g; + return; + } + } + // if N is empty, gN is the min value. + gN = 0; + } + template <typename N_t> void print_N(const N_t& N) @@ -142,7 +172,7 @@ image2d<value::int_u8> out = clone(u); - mln_assertion(R_box.is_valid()); + mln_assertion(R_box.npoints() > 0); mln_piter_(box2d) p(R_box); for_all(p) if (is(p) == in_R) @@ -249,42 +279,136 @@ // std::cout << " <<<<<<<exiting blob." << std::endl; } + template <typename P, typename V> + void + move_A_to_R(p_array& A, + image2d<int>& deja_vu, + fllt_node(P, V)* current_cc, + image2d< fllt_node(P, V)* >& smallest_shapes, + int in_R, + int in_N, + const V& g) + { + typedef p_array arr_t; + + mln_piter(arr_t) a(A); + for_all(a) + { + mln_invariant(deja_vu(a) == in_N); + mln_invariant(smallest_shapes(a) != current_cc); + + deja_vu(a) = in_R; + current_cc->elt().npoints++; + if (!smallest_shapes(a)) + { + smallest_shapes(a) = current_cc; + current_cc->elt().points.append(a); + } + else + if (!smallest_shapes(a)->parent()) + if (smallest_shapes(a)->elt().value == g) + { + mln_piter(arr_t) p(smallest_shapes(a)->elt().points); + // Todo optimization here. + for_all(p) + smallest_shapes(p) = 0; - template <typename I, typename Nbh> + delete smallest_shapes(a); + smallest_shapes(a) = current_cc; + current_cc->elt().points.append(a); + } + else + smallest_shapes(a)->set_parent(current_cc); + // N_box is not re-computed so that we save time; + // N_box is always growing while looping from step 3. + } + } + + // LOWER LEVEL SET : region = c4, border = c8 + template <typename V> + struct lower + { + typedef upper<V> opposite; + static bool + compare(const V& u, const V& v) + { + return u < v; + } + + template <typename I, typename N, typename L> + static mln_ch_value(I, L) + regional_extremum(const Image& input, const Neighborhood<N>& nbh, L& nlabels) + { + return labeling::regional_minima(input, nbh, nlabels); + } + + static const bool parent_is_brighter = true; + + static const neighb2d& bdr_nbh() { return c8(); } + static const neighb2d& reg_nbh() { return c4(); } + }; + + + // UPPER LEVEL SET : region = c8, border = c4 + template <typename V> + struct upper + { + typedef lower<V> opposite; + + static bool + compare(const V& u, const V& v) + { return u > v; } + + template <typename I, typename N, typename L> + static mln_ch_value(I, L) + regional_extremum(const Image& input, const Neighborhood<N>& nbh, L& nlabels) + { return labeling::regional_maxima(input, nbh, nlabels); } + + static const bool parent_is_brighter = false; + static const neighb2d& bdr_nbh() { return c4(); } + static const neighb2d& reg_nbh() { return c8(); } + }; + + template <typename I, typename Set> fllt_tree(mln_point(I), mln_value(I))& - fllt(const Image& input_, const Neighborhood<Nbh>& nbh_) + level_set(const Image& input_, + image2d< fllt_node(mln_point(I), mln_value(I))* >& smallest_shapes) { - typedef fllt_node(mln_point(I), mln_value(I)) node_type; - typedef fllt_tree(mln_point(I), mln_value(I)) tree_type; + + typedef mln_point(I) P; + typedef mln_value(I) V; + typedef fllt_node(P, V) node_type; + typedef fllt_tree(P, V) tree_type; const I& input = exact(input_); - const Nbh& nbh = exact(nbh_); + + mln_assertion(input.domain() == smallest_shapes.domain()); unsigned l = 0, l_max; - mln_ch_value(I, unsigned) reg_min = labeling::regional_minima(input, nbh, l_max); + mln_ch_value(I, unsigned) reg_min = Set::regional_extremum(input, Set::reg_nbh(), l_max); std::vector<bool> tag(l_max + 1, false); tag[0] = true; // Variables. I u = mln::clone(input); - mln_point(I) x0; - mln_value(I) g, gN; + P x0; + V g, gN; mln_fwd_piter(I) p(input.domain()); p.start(); - - + level::fill(smallest_shapes, 0); node_type* current_cc; + unsigned in_N = 1, in_R = 2; image2d<int> deja_vu(input.domain().to_larger(1)); level::fill(deja_vu, 0); - typedef p_array<mln_point(I)> arr_t; + typedef p_array arr_t; arr_t* A = new arr_t(); arr_t* N[256]; for (unsigned i = 0; i < 256; ++i) N[i] = new arr_t(); - accu::bbox<mln_point(I)> N_box; + accu::bbox N_box; unsigned n_step_1 = 0, n_step_3 = 0, n_comps = 0, n_holes = 0; @@ -302,6 +426,8 @@ x0 = p; g = input(x0); current_cc = new node_type(); + current_cc->elt().value = g; + } // Step 2. @@ -319,15 +445,20 @@ N_box.take(x0); deja_vu(x0) = in_R; + smallest_shapes(x0) = current_cc; + current_cc->elt().points.append(x0); + current_cc->elt().npoints++; + } // Step 3. step_3: { +// save_u(u, deja_vu, N_box, in_R); ++n_step_3; mln_piter(arr_t) a(*A); - mln_niter(Nbh) x(nbh, a); + mln_niter(neighb2d) x(Set::reg_nbh(), a); // R <- R U A @@ -349,7 +480,7 @@ } } // gN = min u(x) for all x in N - update_gN(N, gN); + update_gN(N, gN, Set()); // FIXME: update the number of CC of the border of R } @@ -358,30 +489,24 @@ step_4: { // a) - if (g < gN) + if (Set::compare(g, gN)) { g = gN; ++n_comps; blob(deja_vu, N, in_N, current_cc); - node_type* parent = current_cc; + node_type* child = current_cc; + current_cc = new node_type(); + current_cc->elt().value = g; + child->set_parent(current_cc); n_holes += current_cc->elt().holes.npoints(); - save_u(u, deja_vu, N_box, in_R); arr_t* tmp = A; A = N[g]; N[g] = tmp; N[g]->clear(); - mln_piter(arr_t) a(*A); - for_all(a) - { - mln_invariant(deja_vu(a) == in_N); - deja_vu(a) = in_R; - current_cc->elt().points.append(a); - // N_box is not re-computed so that we save time; - // N_box is always growing while looping from step 3. - } + move_A_to_R(*A, deja_vu, current_cc, smallest_shapes, in_R, in_N, g); goto step_3; } // b) @@ -391,13 +516,7 @@ A = N[g]; N[g] = tmp; N[g]->clear(); - mln_piter(arr_t) a(*A); - for_all(a) - { - mln_invariant(deja_vu(a) == in_N); - deja_vu(a) = in_R; - current_cc->elt().points.append(a); - } + move_A_to_R(*A, deja_vu, current_cc, smallest_shapes, in_R, in_N, g); goto step_3; } // c) @@ -420,6 +539,101 @@ return *new tree_type(current_cc); } + + template <typename I> + fllt_tree_ptr(mln_point(I), mln_value(I)) + merge_trees(fllt_tree(mln_point(I), mln_value(I))& lower_tree, + fllt_tree(mln_point(I), mln_value(I))& upper_tree, + const image2d<fllt_node(mln_point(I), mln_value(I))*>& low_reg, + const image2d<fllt_node(mln_point(I), mln_value(I))*>& upp_reg, + const Image& input_) + { + + const I& input = exact(input_); + typedef mln_point(I) P; + typedef mln_value(I) V; + + typedef fllt_node_ptr(P, V) node_ptr_type; + typedef fllt_tree_ptr(P, V) tree_ptr_type; + typedef p_array arr_t; + + fllt_node_ptr(P, V)* root = new node_ptr_type(); + fllt_node_ptr(P, V)* n = root; + root->elt() = &(lower_tree.root()->elt()); + +// fllt_branch_iter_ind(P, V) p(lower_tree.main_branch()); +// for_all(p) +// { +// if (p->elt().tagged) +// continue; + +// p->elt().tagged = true; +// fllt_node_ptr(P, V)* node = new node_ptr_type(); +// node->elt() = &(p->elt()); + +// mln_piter(arr_t) hole(p->elt().holes); +// for_all(hole) +// { + +// } +// } + + return tree_ptr_type(root); + } + + template <typename I> + fllt_tree(mln_point(I), mln_value(I)) + fllt(const Image& input_) + { + typedef mln_point(I) P; + typedef mln_value(I) V; + + const I& input = exact(input_); + + fllt_tree(P, V) upper_tree; + fllt_tree(P, V) lower_tree; + image2d<fllt_node(P, V)*> low_reg(input.domain()); + image2d<fllt_node(P, V)*> upp_reg(input.domain()); + + std::cout << "1/ Compute the lower level set." << std::endl; + lower_tree = level_set<I, lower<V> >(input, low_reg); +// draw_tree(input, lower_tree); + + std::cout << "2/ Compute the upper level set." << std::endl; + upper_tree = level_set<I, upper<V> >(input, upp_reg); +// draw_tree(input, upper_tree); + + fllt_tree_ptr(P, V) result_tree = merge_trees(lower_tree, upper_tree, low_reg, upp_reg, input); + +// draw_tree(lena, tree); + +// return fllt_tree(lower_tree); + } + + template <typename P, typename V> + void + draw_tree(const image2d<V>& ima, + fllt_tree(P, V)& tree) + { + fllt_branch_iter_ind(P, V) p(tree.main_branch()); + for_all(p) + { + std::cout << "region mere : " << (*p).parent() << std::endl; + std::cout << " ^" << std::endl; + std::cout << " |" << std::endl; + std::cout << "region : " << &*p + << " value = " << (*p).elt().value << std::endl + << " npoints = " << (*p).elt().npoints << std::endl + << " holes = " << (*p).elt().holes << std::endl; + + std::cout << std::endl; + + if ((*p).elt().points.npoints() > 0) + debug::println(ima | (*p).elt().points); + std::cout << std::endl; + } + } + } // end of namespace mln::my } // end of namespace mln @@ -427,21 +641,38 @@ int main(int argc, char* argv[]) { + using namespace mln; + using namespace mln::my; + using value::int_u8; + + typedef fllt_tree(point2d, int_u8) tree_type; + if (argc != 2) { std::cerr << "usage: " << argv[0] << " filename" << std::endl; return 1; } - using namespace mln; - using value::int_u8; - image2d<int_u8> lena; io::pgm::load(lena, argv[1]); - my::fllt(lena, c4()); - io::pgm::save(lena, "./out.pgm"); +// int vs[8][9] = { {0,0,0,0,0,0,0,0,0}, +// {0,0,0,0,0,0,0,0,0}, +// {0,1,1,1,1,1,1,1,0}, +// {0,1,0,0,1,3,3,1,0}, +// {0,1,0,0,1,3,3,1,0}, +// {0,1,0,0,1,3,3,1,0}, +// {0,1,1,1,1,1,1,1,0}, +// {0,0,0,0,0,0,0,0,0} }; + +// image2d<int> lena_(make::image2d(vs)); +// image2d<int_u8> lena(lena_.domain()); +// level::fill(lena, lena_); + + tree_type tree = my::fllt(lena); + + io::pgm::save(lena, "./out.pgm"); } // 16891 99970

1 0

1944: Have tools compile again.
by Roland Levillain 14 May '08

14 May '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Roland Levillain <roland(a)lrde.epita.fr> Have tools compile again. * tools/seed2tiling.cc: Don't include mln/draw/mesh.hh. * mln/make/voronoi.hh: s/mln::mesh_p/p_graph/. * mln/core/ops.hh: Add a FIXME about a warning on an uninitialized value. mln/core/ops.hh | 2 ++ mln/make/voronoi.hh | 12 ++++++------ tools/seed2tiling.cc | 1 - 3 files changed, 8 insertions(+), 7 deletions(-) Index: tools/seed2tiling.cc --- tools/seed2tiling.cc (revision 1943) +++ tools/seed2tiling.cc (working copy) @@ -42,7 +42,6 @@ # include <mln/io/pgm/save.hh> # include <mln/labeling/foreground.hh> # include <mln/debug/println.hh> -# include <mln/draw/mesh.hh> # include <mln/geom/seeds2tiling.hh> # include <mln/geom/seeds2tiling_roundness.hh> # include <mln/make/voronoi.hh> Index: mln/make/voronoi.hh --- mln/make/voronoi.hh (revision 1943) +++ mln/make/voronoi.hh (working copy) @@ -1,4 +1,4 @@ -// Copyright (C) 2007 EPITA Research and Development Laboratory +// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory (LRDE) // // This file is part of the Olena Library. This library is free // software; you can redistribute it and/or modify it under the terms @@ -30,14 +30,14 @@ /*! \file mln/make/voronoi.hh * - * \brief Routine to construct a voronoi mln::mesh_p. + * \brief Routine to construct a Voronoi mln::p_graph. */ # include <vector> # include <map> # include <mln/core/concept/neighborhood.hh> -# include <mln/core/mesh_p.hh> +# include <mln/core/p_graph.hh> # include <mln/accu/mean.hh> # include <mln/estim/min_max.hh> @@ -58,7 +58,7 @@ * \return The computed graph. */ template <typename I, typename N> - mesh_p<mln_psite(I)> + p_graph<mln_psite(I)> voronoi (Image& ima_, Image& orig_, const Neighborhood<N>& nbh); @@ -67,7 +67,7 @@ template <typename I, typename N> inline - mesh_p<mln_psite(I)> + p_graph<mln_psite(I)> voronoi (Image& ima_, Image& orig_, const Neighborhood<N>& nbh) @@ -127,7 +127,7 @@ for (; it != m.end (); ++it) gr.add_edge((*it).first.first, (*it).first.second); - mesh_p res(gr, v); + p_graph res(gr, v); return res; } Index: mln/core/ops.hh --- mln/core/ops.hh (revision 1943) +++ mln/core/ops.hh (working copy) @@ -413,6 +413,8 @@ operator-(const Object<O>& rhs) { mlc_converts_to(literal::zero_t, O)::check(); + /* FIXME: The following statements trigger a warning: + ``'p_zero' is used uninitialized in this function.'' */ literal::zero_t* p_zero; return O(*p_zero) - exact(rhs); }

1 0

Re: [Olena-patches] [Olena] #81: Olena 1.0 Reference Documentation
by Olena 14 May '08

14 May '08

#81: Olena 1.0 Reference Documentation -----------------------+---------------------------------------------------- Reporter: levill_r | Owner: Olena Team Type: task | Status: new Priority: minor | Milestone: Olena 1.0 Component: Milena | Version: 1.0 Resolution: | Keywords: doc -----------------------+---------------------------------------------------- Changes (by levill_r): * priority: trivial => minor * keywords: => doc * summary: Olena 1.0ß Reference Documentation => Olena 1.0 Reference Documentation Old description: > Tweak the documentation generated by Doxygen. > > (This is not priority for the Olena 1.0ß milestone.) New description: Tweak the documentation generated by Doxygen. (This is not priority for the Olena 1.0 milestone.) Comment: s/1.0ß/1.0/ -- Ticket URL: <https://trac.lrde.org/olena/ticket/81#comment:1> Olena <http://olena.lrde.epita.fr> Olena, a generic and efficient C++ image processing library.

1 0

Re: [Olena-patches] [Olena] #98: Create demonstration video(s) for Olena 1.0
by Olena 14 May '08

14 May '08

#98: Create demonstration video(s) for Olena 1.0 -----------------------+---------------------------------------------------- Reporter: levill_r | Owner: Olena Team Type: task | Status: new Priority: trivial | Milestone: Olena 1.0 Component: Olena | Version: 1.0 Resolution: | Keywords: -----------------------+---------------------------------------------------- Changes (by levill_r): * type: defect => task * component: Milena => Olena * summary: Create demonstration video(s) for Olena 1.0ß => Create demonstration video(s) for Olena 1.0 Comment: s/1.0ß/1.0/ -- Ticket URL: <https://trac.lrde.org/olena/ticket/98#comment:1> Olena <http://olena.lrde.epita.fr> Olena, a generic and efficient C++ image processing library.

1 0

Re: [Olena-patches] [Olena] #95: Create an Olena 1.0 port for the MacPorts (Mac OS X)
by Olena 14 May '08

14 May '08

#95: Create an Olena 1.0 port for the MacPorts (Mac OS X) -----------------------+---------------------------------------------------- Reporter: levill_r | Owner: Olena Team Type: task | Status: new Priority: minor | Milestone: Olena 1.0 Component: Milena | Version: 1.0 Resolution: | Keywords: -----------------------+---------------------------------------------------- Changes (by levill_r): * summary: Create an Olena 1.0ß port for the MacPorts (Mac OS X) => Create an Olena 1.0 port for the MacPorts (Mac OS X) Old description: > Don't forget to add the right dependencies (and suggested packages, if > available). > > Don't mess with the numbering scheme: think about it twice! For > instance, keep in mind that `1.0beta` is after `1.0` in the lexical > order. Find some documentation about this. > > Ideas: > * Create packages for the different variant of Darwin (at least Tiger > and Leopard). > * Several packages might be useful (not necessarily for the 1.0ß > release, but later). E.g., a package (possibly virtual) for Olena, and > another one for Milena. > > See also #93 and #94. New description: Don't forget to add the right dependencies (and suggested packages, if available). Don't mess with the numbering scheme: think about it twice! For instance, keep in mind that `1.0a` is after `1.0` in the lexical order. Find some documentation about this. Ideas: * Create packages for the different variant of Darwin (at least Tiger and Leopard). * Several packages might be useful (not necessarily for the 1.0 release, but later). E.g., a package (possibly virtual) for Olena, and another one for Milena. See also #93 and #94. Comment: s/1.0ß/1.0/ -- Ticket URL: <https://trac.lrde.org/olena/ticket/95#comment:1> Olena <http://olena.lrde.epita.fr> Olena, a generic and efficient C++ image processing library.

1 0

Re: [Olena-patches] [Olena] #94: Create Olena 1.0 Debian package(s) for amd64
by Olena 14 May '08

14 May '08

#94: Create Olena 1.0 Debian package(s) for amd64 -----------------------+---------------------------------------------------- Reporter: levill_r | Owner: Olena Team Type: task | Status: new Priority: minor | Milestone: Olena 1.0 Component: Olena | Version: 1.0 Resolution: | Keywords: -----------------------+---------------------------------------------------- Changes (by levill_r): * summary: Create Olena 1.0ß Debian package(s) for amd64 => Create Olena 1.0 Debian package(s) for amd64 Old description: > Don't forget to add the right dependencies (and suggested packages, such > as a C++ compiler!). > > Don't mess with the numbering scheme: think about it twice! For > instance, keep in mind that `1.0beta` is after `1.0` in the lexical > order. Find some documentation about this. > > Ideas: > * Create packages for the different versions of Debian (at least > `stable` and `unstable`). > * Create packages for Ubuntu (e.g., if the dependencies differ from > Debian's). > * Several packages might be useful (not necessarily for the 1.0ß > release, but later). E.g., a package (possibly virtual) for Olena, and > another one for Milena. > > See also #93 and #95. New description: Don't forget to add the right dependencies (and suggested packages, such as a C++ compiler!). Don't mess with the numbering scheme: think about it twice! For instance, keep in mind that `1.0a` is after `1.0` in the lexical order. Find some documentation about this. Ideas: * Create packages for the different versions of Debian (at least `stable` and `unstable`). * Create packages for Ubuntu (e.g., if the dependencies differ from Debian's). * Several packages might be useful (not necessarily for the 1.0 release, but later). E.g., a package (possibly virtual) for Olena, and another one for Milena. See also #93 and #95. Comment: s/1.0ß/1.0/ -- Ticket URL: <https://trac.lrde.org/olena/ticket/94#comment:2> Olena <http://olena.lrde.epita.fr> Olena, a generic and efficient C++ image processing library.

1 0

Re: [Olena-patches] [Olena] #93: Create Olena 1.0 Debian package(s) for i386
by Olena 14 May '08

14 May '08

#93: Create Olena 1.0 Debian package(s) for i386 -----------------------+---------------------------------------------------- Reporter: levill_r | Owner: Olena Team Type: task | Status: new Priority: minor | Milestone: Olena 1.0 Component: Olena | Version: 1.0 Resolution: | Keywords: -----------------------+---------------------------------------------------- Changes (by levill_r): * summary: Create Olena 1.0ß Debian package(s) for i386 => Create Olena 1.0 Debian package(s) for i386 Old description: > Don't forget to add the right dependencies (and suggested packages, such > as a C++ compiler!). > > Don't mess with the numbering scheme: think about it twice! For > instance, keep in mind that `1.0beta` is after `1.0` in the lexical > order. Find some documentation about this. > > Ideas: > * Create packages for the different versions of Debian (at least > `stable` and `unstable`). > * Create packages for Ubuntu (e.g., if the dependencies differ from > Debian's). > * Several packages might be useful (not necessarily for the 1.0ß > release, but later). E.g., a package (possibly virtual) for Olena, and > another one for Milena. > > See also #94 and #95. New description: Don't forget to add the right dependencies (and suggested packages, such as a C++ compiler!). Don't mess with the numbering scheme: think about it twice! For instance, keep in mind that `1.0a` is after `1.0` in the lexical order. Find some documentation about this. Ideas: * Create packages for the different versions of Debian (at least `stable` and `unstable`). * Create packages for Ubuntu (e.g., if the dependencies differ from Debian's). * Several packages might be useful (not necessarily for the 1.0 release, but later). E.g., a package (possibly virtual) for Olena, and another one for Milena. See also #94 and #95. Comment: s/1.0ß/1.0/ -- Ticket URL: <https://trac.lrde.org/olena/ticket/93#comment:2> Olena <http://olena.lrde.epita.fr> Olena, a generic and efficient C++ image processing library.

1 0

Re: [Olena-patches] [Olena] #7: Miscellaneous tasks related to the Olena 1.0 release process
by Olena 14 May '08

14 May '08

#7: Miscellaneous tasks related to the Olena 1.0 release process -----------------------+---------------------------------------------------- Reporter: levill_r | Owner: levill_r Type: task | Status: new Priority: major | Milestone: Olena 1.0 Component: other | Version: 1.0 Resolution: | Keywords: release NEWS dates copyright headers sandbox -----------------------+---------------------------------------------------- Old description: > To do before the release: > * use Benoît Sigoure's set of Autoconf macros to check for some Boost > libraries (e.g., the BGL). See also #138; > * move away any information about sandboxes (as they won't be > distributed): > * create a !ChangeLog in source:trunk/milena/sandbox, or several > !ChangeLogs in subdirectories; > * move information from source:trunk/milena/ChangeLog to this (these) > file(s); > * warn the Olena Team about these changes, and have them document > changes in the right place; > * update `NEWS` (see the similar task for Olena 0.10a/0.11: #8); > * update `AUTHORS`; > * ensure all copyright headers are sound (using Benoît Sigoure's > `gplize.sh` script, or a variant of it); > * use Vaucanson' `reheader` script to fix broken headers. > > Split this ticket in several parts if the tasks are either too big or too > numerous. New description: To do before the release: * use Benoît Sigoure's set of Autoconf macros to check for some Boost libraries (e.g., the BGL). See also #138; * ~~move away any information about sandboxes (as they won't be distributed):~~ * ~~create a !ChangeLog in source:trunk/milena/sandbox, or several !ChangeLogs in subdirectories~~ (done in r1923); * ~~move information from source:trunk/milena/ChangeLog to this (these) file(s)~~ (done in r1923); * ~~warn the Olena Team about these changes, and have them document changes in the right place~~; * update `NEWS` (see the similar task for Olena 0.10a/0.11: #8); * update `AUTHORS`; * ensure all copyright headers are sound (using Benoît Sigoure's `gplize.sh` script, or a variant of it); * use Vaucanson' `reheader` script to fix broken headers. Split this ticket in several parts if the tasks are either too big or too numerous. Comment (by levill_r): Strike out sandboxes' !ChangeLog-related items (r1923). -- Ticket URL: <https://trac.lrde.org/olena/ticket/7#comment:10> Olena <http://olena.lrde.epita.fr> Olena, a generic and efficient C++ image processing library.

1 0

Re: [Olena-patches] [Olena] #84: Olena 1.0 Code Examples (as documentation)
by Olena 14 May '08

14 May '08

#84: Olena 1.0 Code Examples (as documentation) -----------------------+---------------------------------------------------- Reporter: levill_r | Owner: Olena Team Type: task | Status: new Priority: major | Milestone: Olena 1.0 Component: Milena | Version: 1.0 Resolution: | Keywords: doc -----------------------+---------------------------------------------------- Old description: > Write code samples (to be installed into ''prefix''`/share/olena` or > something). They shall be highly documented. I (Roland) don't think we > should use Doxygen fo such a task: plain comments are just fine, as one > can read them while having the whole code under one's eyes. > > The description of this ticket might be extended by the kind of code > samples that should be added to Olena 1.0ß, or even a detailed list. If > such a list were to grow significantly, a dedicated Trac Wiki page would > be a better place. New description: Write code samples (to be installed into ''prefix''`/share/olena` or something). They shall be highly documented. I (Roland) don't think we should use Doxygen fo such a task: plain comments are just fine, as one can read them while having the whole code under one's eyes. The description of this ticket might be extended by the kind of code samples that should be added to Olena 1.0, or even a detailed list. If such a list were to grow significantly, a dedicated Trac Wiki page would be a better place. Comment (by levill_r): s/1.0ß/1.0/ -- Ticket URL: <https://trac.lrde.org/olena/ticket/84#comment:2> Olena <http://olena.lrde.epita.fr> Olena, a generic and efficient C++ image processing library.

1 0