https://svn.lrde.epita.fr/svn/oln/branches/cleanup-2008/milena Index: ChangeLog from Thierry Geraud <thierry.geraud@lrde.epita.fr> Add some morphological reconstruction algorithms in sandbox. * sandbox/geraud/Rd: New directory. * sandbox/geraud/Rd/hybrid.hh, * sandbox/geraud/Rd/queue_based.hh, * sandbox/geraud/Rd/parallel.cc, * sandbox/geraud/Rd/parallel_wo.cc, * sandbox/geraud/Rd/union_find.hh, * sandbox/geraud/Rd/parallel.hh, * sandbox/geraud/Rd/diff.cc, * sandbox/geraud/Rd/sequential_bench.cc, * sandbox/geraud/Rd/sequential.cc, * sandbox/geraud/Rd/utils.hh, * sandbox/geraud/Rd/deco.cc, * sandbox/geraud/Rd/hybrid.cc, * sandbox/geraud/Rd/queue_based.cc, * sandbox/geraud/Rd/sequential_bench.hh, * sandbox/geraud/Rd/min.cc, * sandbox/geraud/Rd/sequential.hh, * sandbox/geraud/Rd/debase.union_find.hh, * sandbox/geraud/Rd/union_find.cc, * sandbox/geraud/Rd/svg.queue_based.hh: New files. debase.union_find.hh | 165 +++++++++++++++++++++++++++++++ deco.cc | 71 +++++++++++++ diff.cc | 28 +++++ hybrid.cc | 51 +++++++++ hybrid.hh | 123 +++++++++++++++++++++++ min.cc | 42 +++++++ parallel.cc | 51 +++++++++ parallel.hh | 93 +++++++++++++++++ parallel_wo.cc | 46 ++++++++ queue_based.cc | 49 +++++++++ queue_based.hh | 129 ++++++++++++++++++++++++ sequential.cc | 51 +++++++++ sequential.hh | 96 ++++++++++++++++++ sequential_bench.cc | 51 +++++++++ sequential_bench.hh | 100 +++++++++++++++++++ svg.queue_based.hh | 118 ++++++++++++++++++++++ union_find.cc | 50 +++++++++ union_find.hh | 173 ++++++++++++++++++++++++++++++++ utils.hh | 268 +++++++++++++++++++++++++++++++++++++++++++++++++++ 19 files changed, 1755 insertions(+) Index: sandbox/geraud/Rd/hybrid.hh --- sandbox/geraud/Rd/hybrid.hh (revision 0) +++ sandbox/geraud/Rd/hybrid.hh (revision 0) @@ -0,0 +1,123 @@ +// 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_MORPHO_RD_HYBRID_HH +# define MLN_MORPHO_RD_HYBRID_HH + +# include <queue> +# include "utils.hh" + + +namespace mln +{ + + namespace morpho + { + + namespace Rd + { + + + template <typename I, typename N> + I hybrid(const I& f, const I& g, const N& nbh, + bool echo = false) + { + typedef mln_point(I) point; + std::queue<point> q; + + f.name_it("f"); + g.name_it("g"); + + // initialisation + I o(f.domain()); + o.name_it("o"); + level::paste(f, o); + // WAS: I o = clone(f); + + unsigned n_init_pushs = 0, n_body_pushs = 0, n_pops = 0; + + // sequence + { + mln_bkd_piter(I) p(f.domain()); + for_all(p) + o(p) = min( max_Nminus(o, p,nbh), g(p) ); + } + { + mln_fwd_piter(I) p(f.domain()); + mln_niter(N) n(nbh, p); + for_all(p) + { + o(p) = min( max_Nplus(o, p,nbh), g(p) ); + for_all(n) if (f.has(n) and n < p) // N+ + if (o(n) < o(p) and o(n) < g(n)) + { + q.push(p); + ++n_init_pushs; + } + } + } + + // propagation + { + point p; + mln_niter(N) n(nbh, p); + while (not q.empty()) + { + p = q.front(); + if (echo) std::cout << std::endl << "pop " << p << " :"; + q.pop(); + ++n_pops; + for_all(n) if (f.has(n)) + if (o(n) < o(p) and o(n) != g(n)) + { + o(n) = min(o(p), g(n)); + if (echo) std::cout << " push " << n; + q.push(n); + ++n_body_pushs; + } + } + if (echo) std::cout << std::endl; + } + + + std::cout << "n_init_pushs = " << n_init_pushs << std::endl + << "n_body_pushs = " << n_body_pushs << std::endl + << "n_pops = " << n_pops << std::endl; + + print_counts(); + + return o; + } + + } // end of namespace mln::morpho::Rd + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_RD_HYBRID_HH Index: sandbox/geraud/Rd/queue_based.hh --- sandbox/geraud/Rd/queue_based.hh (revision 0) +++ sandbox/geraud/Rd/queue_based.hh (revision 0) @@ -0,0 +1,129 @@ +// 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_MORPHO_RD_QUEUE_BASED_HH +# define MLN_MORPHO_RD_QUEUE_BASED_HH + +# include <queue> +# include "utils.hh" + + +namespace mln +{ + + namespace morpho + { + + namespace Rd + { + + template <typename I, typename N> + I queue_based(const I& f, const I& g, const N& nbh, + bool echo = false) + { + + mln_ch_value(I, bool) que(f.domain()); + level::fill(que, false); + + + if (echo) std::cout << std::endl; + + f.name_it("f"); + g.name_it("g"); + + typedef mln_point(I) point; + std::queue<point> q; + I o; + o.name_it("o"); + + unsigned n_init_pushs = 0, n_body_pushs = 0, n_pops = 0; + + // initialisation + { + o = regional_maxima(f, nbh); + // p in M <=> o(p) != 0 + if (echo) debug::println(o); + + mln_piter(I) p(f.domain()); + mln_niter(N) n(nbh, p); + + for_all(p) if (o(p) != 0) // p in M + for_all(n) if (f.has(n) and o(n) == 0) // n not in M + { + q.push(p); + que(p) = true; + ++n_init_pushs; + break; + } + } + + // propagation + { + point p; + mln_niter(N) n(nbh, p); + while (not q.empty()) + { + p = q.front(); + if (echo) std::cout << std::endl << "pop " << p << " :"; + q.pop(); + que(p) = false; + ++n_pops; + for_all(n) if (f.has(n)) + { + if (o(n) < o(p) and o(n) != g(n)) + { + o(n) = min(o(p), g(n)); + if (echo) std::cout << " push " << n; + if (que(n) == false) + { + q.push(n); + que(n) = true; + ++n_body_pushs; + } + } + } + } + if (echo) std::cout << std::endl; + } + + std::cout << "n_init_pushs = " << n_init_pushs << std::endl + << "n_body_pushs = " << n_body_pushs << std::endl + << "n_pops = " << n_pops << std::endl; + + print_counts(); + + return o; + } + + } // end of namespace mln::morpho::Rd + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_RD_QUEUE_BASED_HH Index: sandbox/geraud/Rd/parallel.cc --- sandbox/geraud/Rd/parallel.cc (revision 0) +++ sandbox/geraud/Rd/parallel.cc (revision 0) @@ -0,0 +1,51 @@ +#include <iostream> + +#include <mln/core/image2d.hh> +#include <mln/core/neighb2d.hh> +#include <mln/value/int_u8.hh> + +#include <mln/debug/println.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "parallel.hh" + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " f.pgm g.pgm c output.pgm" << std::endl + << "reconstruction by dilation (parallel version; may 2007)" << std::endl + << "f = marker (to be dilated)" << std::endl + << "g = mask (constraint >= f)" << std::endl + << "c: 4 or 8" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 5) + usage(argv); + + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::pgm::load<int_u8>(argv[1]); + I g = io::pgm::load<int_u8>(argv[2]); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::pgm::save(morpho::Rd::parallel(f, g, + (c == 4 ? c4() : c8())), + argv[4]); +} Index: sandbox/geraud/Rd/parallel_wo.cc --- sandbox/geraud/Rd/parallel_wo.cc (revision 0) +++ sandbox/geraud/Rd/parallel_wo.cc (revision 0) @@ -0,0 +1,46 @@ +#include <iostream> + +#include <mln/core/image2d.hh> +#include <mln/core/neighb2d.hh> +#include <mln/value/int_u8.hh> + +#include <mln/debug/println.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "parallel.hh" + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " f.pgm g.pgm c output.pgm" << std::endl + << "reconstruction by dilation (parallel version; may 2007)" << std::endl + << "f = marker (to be dilated)" << std::endl + << "g = mask (constraint >= f)" << std::endl + << "c: 4 or 8" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 5) + usage(argv); + + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::pgm::load<int_u8>(argv[1]); + I g = io::pgm::load<int_u8>(argv[2]); + + io::pgm::save(morpho::Rd::parallel(f, g, + (c == 4 ? c4() : c8()), + false), + argv[4]); +} Index: sandbox/geraud/Rd/union_find.hh --- sandbox/geraud/Rd/union_find.hh (revision 0) +++ sandbox/geraud/Rd/union_find.hh (revision 0) @@ -0,0 +1,173 @@ +// 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_MORPHO_RD_UNION_FIND_HH +# define MLN_MORPHO_RD_UNION_FIND_HH + +# include "utils.hh" + + +namespace mln +{ + + namespace morpho + { + + namespace Rd + { + + + template <typename I, typename N> + struct union_find_t + { + typedef mln_point(I) point; + typedef mln_value(I) value; + + // in: + const I f, g; + N nbh; + + // out: + I o; + + // aux: + std::vector<point> S; + // was: I data; + mln_ch_value(I, bool) isproc; + mln_ch_value(I, point) parent; + + union_find_t(const I& f, const I& g, const N& nbh) + : f(f), g(g), nbh(nbh) + { + f.name_it("f"); + g.name_it("g"); + initialize(o, f); + o.name_it("o"); + initialize(parent, f); + parent.name_it("parent"); + initialize(isproc, f); + isproc.name_it("isproc"); + // was: initialize(data, f); + + // init + + level::fill(isproc, false); + S = histo_reverse_sort(g); + level::paste(f, o); // new: replace make_set(p) { data(p) = f(p) } + + // first pass + + for (unsigned i = 0; i < S.size(); ++i) + { + point p = S[i]; + + make_set(p); + mln_niter(N) n(nbh, p); + for_all(n) + { + if (f.has(n)) + mln_invariant(isproc(n) == is_proc(n, p)); + if (f.has(n) and isproc(n)) + do_union(n, p); + } + isproc(p) = true; + } + + // second pass + + for (int i = S.size() - 1; i >= 0; --i) + { + point p = S[i]; + if (parent(p) == p) + if (o(p) == mln_max(value)) + o(p) = g(p); + else + o(p) = o(parent(p)); + } + + print_counts(); + + } + + bool is_proc(const point& n, const point& p) const + { + return g(n) > g(p) or (g(n) == g(p) and n < p); + } + + void make_set(const point& p) + { + parent(p) = p; + // was: data(p) = f(p); + // now: in "initialization" + } + + point find_root(const point& x) + { + if (parent(x) == x) + return x; + else + return parent(x) = find_root(parent(x)); + } + + void do_union(const point& n, const point& p) + { + point r = find_root(n); + if (r != p) + { + // NEW: o replaces data + + if (g(r) == g(p) or g(p) >= o(r)) // equiv test + { + parent(r) = p; + if (o(r) > o(p)) + o(p) = o(r); // increasing criterion + } + else + o(p) = mln_max(value); + } + } + + }; + + + template <typename I, typename N> + I union_find(const I& f, const I& g, const N& nbh) + { + mln_precondition(f <= g); + union_find_t<I, N> run(f, g, nbh); + return run.o; + } + + + } // end of namespace mln::morpho::Rd + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_RD_UNION_FIND_HH Index: sandbox/geraud/Rd/parallel.hh --- sandbox/geraud/Rd/parallel.hh (revision 0) +++ sandbox/geraud/Rd/parallel.hh (revision 0) @@ -0,0 +1,93 @@ +// 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_MORPHO_RD_PARALLEL_HH +# define MLN_MORPHO_RD_PARALLEL_HH + +# include "utils.hh" + + +namespace mln +{ + + namespace morpho + { + + namespace Rd + { + + + template <typename I, typename N> + I parallel(const I& f, const I& g, const N& nbh, + bool test = true) + { +// if (test) +// mln_precondition(f <= g); + + f.name_it("f"); + g.name_it("g"); + + I o_(f.domain()); + o_.name_it("o_"); + mln_piter(I) p(f.domain()); + + // initialisation + I o(f.domain()); + o.name_it("o"); + level::paste(f, o); + // WAS: I o = clone(f); + + bool stability; + do + { + level::paste(o, o_); // memorisation + + // opere + for_all(p) + o(p) = max_N(o_, p, nbh); + // conditionne + for_all(p) + o(p) = min(o(p), g(p)); + + stability = (o == o_); + } + while (not stability); + + print_counts(); + + // mln_postcondition(o <= g); + return o; + } + + } // end of namespace mln::morpho::Rd + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_RD_PARALLEL_HH Index: sandbox/geraud/Rd/diff.cc --- sandbox/geraud/Rd/diff.cc (revision 0) +++ sandbox/geraud/Rd/diff.cc (revision 0) @@ -0,0 +1,28 @@ +#include <oln/core/2d/image2d.hh> +#include <oln/io/load_pgm.hh> +#include <oln/level/compare.hh> + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " 1.pgm 2.pgm" << std::endl + << "(may 2007)" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 3) + usage(argv); + + using namespace oln; + typedef image2d<unsigned char> I; + + I ima1 = io::load_pgm(argv[1]); + I ima2 = io::load_pgm(argv[2]); + + if (ima1 != ima2) + std::cout << "images differ" << std::endl; + return 0; +} Index: sandbox/geraud/Rd/sequential_bench.cc --- sandbox/geraud/Rd/sequential_bench.cc (revision 0) +++ sandbox/geraud/Rd/sequential_bench.cc (revision 0) @@ -0,0 +1,51 @@ +#include <mln/core/image2d.hh> +#include <mln/core/neighb2d.hh> + +#include <mln/debug/println.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "sequential_bench.hh" + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " f.pgm g.pgm c output.pgm" << std::endl + << "reconstruction by dilation (sequential version; may 2007)" << std::endl + << "f = marker (to be dilated)" << std::endl + << "g = mask (constraint >= f)" << std::endl + << "c: 4 or 8" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 5) + usage(argv); + + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::pgm::load<int_u8>(argv[1]); + f.name_it("main.f"); + + I g = io::pgm::load<int_u8>(argv[2]); + g.name_it("main.g"); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::pgm::save(morpho::Rd::sequential(f, g, + (c == 4 ? c4() : c8())), + argv[4]); +} Index: sandbox/geraud/Rd/sequential.cc --- sandbox/geraud/Rd/sequential.cc (revision 0) +++ sandbox/geraud/Rd/sequential.cc (revision 0) @@ -0,0 +1,51 @@ +#include <mln/core/image2d.hh> +#include <mln/core/neighb2d.hh> + +#include <mln/debug/println.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "sequential.hh" + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " f.pgm g.pgm c output.pgm" << std::endl + << "reconstruction by dilation (sequential version; may 2007)" << std::endl + << "f = marker (to be dilated)" << std::endl + << "g = mask (constraint >= f)" << std::endl + << "c: 4 or 8" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 5) + usage(argv); + + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::pgm::load<int_u8>(argv[1]); + f.name_it("main.f"); + + I g = io::pgm::load<int_u8>(argv[2]); + g.name_it("main.g"); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::pgm::save(morpho::Rd::sequential(f, g, + (c == 4 ? c4() : c8())), + argv[4]); +} Index: sandbox/geraud/Rd/utils.hh --- sandbox/geraud/Rd/utils.hh (revision 0) +++ sandbox/geraud/Rd/utils.hh (revision 0) @@ -0,0 +1,268 @@ +// 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_MORPHO_RD_UTILS_HH +# define MLN_MORPHO_RD_UTILS_HH + +# include <vector> + +# include <mln/core/concept/image.hh> +# include <mln/core/clone.hh> + +# include <mln/level/fill.hh> +# include <mln/level/paste.hh> +# include <mln/level/compare.hh> + + + +namespace mln +{ + + namespace morpho + { + + namespace Rd + { + + + template <typename T> + T min(T v1, T v2) + { + return v1 < v2 ? v1 : v2; + } + + + template <typename I> + I minimun(const I& ima1, const I& ima2) + { + mln_precondition(ima1.has_data() and ima2.has_data()); + mln_precondition(ima1.domain() == ima2.domain()); + I out(ima1.domain()); + mln_piter(I) p(ima1.domain()); + for_all(p) + out(p) = ima1(p) < ima2(p) ? ima1(p) : ima2(p); + return out; + } + + + template <typename I, typename P, typename N> + mln_value(I) max_N(const I& ima, const P& p, const N& nbh) + { + mln_value(I) v = ima(p); + mln_niter(N) n(nbh, p); + for_all(n) + if (ima.has(n) and ima(n) > v) + v = ima(n); + return v; + } + + + template <typename I, typename P, typename N> + mln_value(I) max_Nminus(const I& ima, const P& p, const N& nbh) + { + mln_value(I) v = ima(p); + mln_niter(N) n(nbh, p); + for_all(n) + if (ima.has(n) and n > p and ima(n) > v) + v = ima(n); + return v; + } + + template <typename I, typename P, typename N> + mln_value(I) max_Nplus(const I& ima, const P& p, const N& nbh) + { + mln_value(I) v = ima(p); + mln_niter(N) n(nbh, p); + for_all(n) + if (ima.has(n) and n < p and ima(n) > v) + v = ima(n); + return v; + } + + + template <typename I> + std::vector<unsigned> compute_histo(const I& ima) + { + std::vector<unsigned> h(256, 0); + mln_piter(I) p(ima.domain()); + for_all(p) + ++h[ima(p)]; + return h; + } + + +// template <typename I> +// std::vector<mln_point(I)> histo_sort(const I& ima) +// { +// std::vector<unsigned> h = compute_histo(ima); +// // preparing output data +// std::vector<int> loc(256); +// loc[0] = 0; +// for (int l = 1; l < 256; ++l) +// loc[l] = loc[l-1] + h[l-1]; +// std::vector<mln_point(I)> vec(ima.points().npoints()); +// // storing output data +// mln_piter(I) p(ima.points()); +// for_all(p) +// vec[loc[ima(p)]++] = p; +// return vec; +// } + + + template <typename I> + std::vector<mln_point(I)> histo_reverse_sort(const I& ima) + { + std::vector<unsigned> h = compute_histo(ima); + // preparing output data + std::vector<int> loc(256); + loc[255] = 0; + for (int l = 254; l >= 0; --l) + loc[l] = loc[l+1] + h[l+1]; + std::vector<mln_point(I)> vec(ima.domain().npoints()); + // storing output data + mln_piter(I) p(ima.domain()); + for_all(p) + vec[loc[ima(p)]++] = p; + return vec; + } + + + + template <typename I, typename N> + struct regional_maxima_t + { + typedef mln_point(I) point; + typedef mln_ch_value(I, bool) image_bool; + typedef mln_ch_value(I, point) image_point; + + // in: + I f; + N nbh; + + // out: + I o; + + // aux: + std::vector<point> S; + image_bool is_proc; + image_bool attr; + image_point parent; + + regional_maxima_t(const I& f, const N& nbh) + : f(f), nbh(nbh) + { + f.name_it("rm.f"); + initialize(o, f); + o.name_it("rm.o"); + initialize(parent, f); + parent.name_it("rm.parent"); + initialize(attr, f); + attr.name_it("rm.attr"); + initialize(is_proc, f); + is_proc.name_it("rm.is_proc"); + + // init + + level::fill(is_proc, false); + S = histo_reverse_sort(f); + + // first pass + + for (unsigned i = 0; i < S.size(); ++i) + { + point p = S[i]; + + make_set(p); + mln_niter(N) n(nbh, p); + for_all(n) + if (f.has(n) and is_proc(n)) + { + if (f(n) == f(p)) + do_union(n, p); + else // f(n) > f(p) + attr(p) = false; + } + is_proc(p) = true; + } + + // second pass + + const mln_value(I) zero = 0; + for (int i = S.size() - 1; i >= 0; --i) + { + point p = S[i]; + if (parent(p) == p) + o(p) = attr(p) ? f(p) : zero; + else + o(p) = o(parent(p)); + } + } + + void make_set(const point& p) + { + parent(p) = p; + attr(p) = true; + } + + point find_root(const point& x) + { + if (parent(x) == x) + return x; + else + return parent(x) = find_root(parent(x)); + } + + void do_union(const point& n, const point& p) + { + point r = find_root(n); + if (r != p) + { + parent(r) = p; + attr(p) = attr(p) and attr(r); + } + } + + }; + + + template <typename I, typename N> + I + regional_maxima(const I& f, const N& nbh) + { + regional_maxima_t<I, N> run(f, nbh); + return run.o; + } + + + } // end of namespace mln::morpho::Rd + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_RD_UTILS_HH Index: sandbox/geraud/Rd/deco.cc --- sandbox/geraud/Rd/deco.cc (revision 0) +++ sandbox/geraud/Rd/deco.cc (revision 0) @@ -0,0 +1,71 @@ +// Copyright (C) 2007 EPITA Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +/*! \file tests/decorated_image.cc + * + * \brief Tests on mln::decorated_image. + * \todo Make this test work. + */ + +#include <mln/core/image2d.hh> +#include <mln/core/decorated_image.hh> + + +unsigned count_read = 0, count_write = 0; + +template <typename I> +struct counter +{ + void reading(const I&, const mln_psite(I)&) const + { + ++count_read; + } + void writing(I&, const mln_psite(I)&, const mln_value(I)&) + { + ++count_write; + } +}; + + +int main() +{ + using namespace mln; + + typedef image2d<int> I; + I ima(1, 1); + point2d p = make::point2d(0, 0); + + decorated_image< I, counter<I> > ima_ = decorate(ima, counter<I>()); + ima_(p) = ima_(p) = 51; + + std::cout << count_read << ' ' << count_write << std::endl; + + mln_assertion(count_read == 1 && count_write == 2); + + const I& imac = ima; + decorated_image< const I, counter<I> > cima_ = decorate(imac, counter<I>()); +} Index: sandbox/geraud/Rd/hybrid.cc --- sandbox/geraud/Rd/hybrid.cc (revision 0) +++ sandbox/geraud/Rd/hybrid.cc (revision 0) @@ -0,0 +1,51 @@ +#include <iostream> + +#include <mln/core/image2d.hh> +#include <mln/core/neighb2d.hh> +#include <mln/value/int_u8.hh> + +#include <mln/debug/println.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "hybrid.hh" + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " f.pgm g.pgm c output.pgm" << std::endl + << "reconstruction by dilation (hybrid version; may 2007)" << std::endl + << "f = marker (to be dilated)" << std::endl + << "g = mask (constraint >= f)" << std::endl + << "c: 4 or 8" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 5) + usage(argv); + + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::pgm::load<int_u8>(argv[1]); + I g = io::pgm::load<int_u8>(argv[2]); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::pgm::save(morpho::Rd::hybrid(f, g, + (c == 4 ? c4() : c8())), + argv[4]); +} Index: sandbox/geraud/Rd/queue_based.cc --- sandbox/geraud/Rd/queue_based.cc (revision 0) +++ sandbox/geraud/Rd/queue_based.cc (revision 0) @@ -0,0 +1,49 @@ +#include <mln/core/image2d.hh> +#include <mln/core/neighb2d.hh> +#include <mln/value/int_u8.hh> + +#include <mln/debug/println.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "queue_based.hh" + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " f.pgm g.pgm c output.pgm" << std::endl + << "reconstruction by dilation (queue_based version; may 2007)" << std::endl + << "f = marker (to be dilated)" << std::endl + << "g = mask (constraint >= f)" << std::endl + << "c: 4 or 8" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 5) + usage(argv); + + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::pgm::load<int_u8>(argv[1]); + I g = io::pgm::load<int_u8>(argv[2]); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::pgm::save(morpho::Rd::queue_based(f, g, + (c == 4 ? c4() : c8())), + argv[4]); +} Index: sandbox/geraud/Rd/sequential_bench.hh --- sandbox/geraud/Rd/sequential_bench.hh (revision 0) +++ sandbox/geraud/Rd/sequential_bench.hh (revision 0) @@ -0,0 +1,100 @@ +// 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_MORPHO_RD_SEQUENTIAL_HH +# define MLN_MORPHO_RD_SEQUENTIAL_HH + +# include "utils.hh" + + +namespace mln +{ + + namespace morpho + { + + namespace Rd + { + + template <typename I, typename N> + I sequential(const I& f, const I& g, const N& nbh) + { + mln_precondition(f <= g); + + f.name_it("f"); + g.name_it("g"); + + I o_(f.domain()); + o_.name_it("o_"); + + unsigned nloops = 0; + + // initialisation + I o = clone(f); + o.name_it("o"); + + bool stability; + do + { + ++nloops; + + level::paste(o, o_); // memorisation + + // passe 1 + { + mln_bkd_piter(I) p(f.domain()); + for_all(p) + o(p) = min( max_Nminus(o, p, nbh), g(p) ); + } + + // passe 2 + { + mln_fwd_piter(I) p(f.domain()); + for_all(p) + o(p) = min( max_Nplus(o, p, nbh), g(p) ); + } + + stability = (o == o_); + } + while (not stability); + + std::cout << "nloops = " << nloops << std::endl; + + print_counts(); + + mln_postcondition(o <= g); + return o; + } + + } // end of namespace mln::morpho::Rd + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_RD_SEQUENTIAL_HH Index: sandbox/geraud/Rd/min.cc --- sandbox/geraud/Rd/min.cc (revision 0) +++ sandbox/geraud/Rd/min.cc (revision 0) @@ -0,0 +1,42 @@ +#include <mln/core/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " 1.pgm 2.pgm out.pgm" << std::endl + << "(2008 Feb)" << std::endl; + exit(1); +} + + +template <typename I> +I min(const I& ima1, const I& ima2) +{ + mln_precondition(ima1.has_data() and ima2.has_data()); + mln_precondition(ima1.domain() == ima2.domain()); + I out(ima1.domain()); + mln_piter(I) p(ima1.domain()); + for_all(p) + out(p) = ima1(p) < ima2(p) ? ima1(p) : ima2(p); + return out; +} + + +int main(int argc, char* argv[]) +{ + if (argc != 4) + usage(argv); + + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + + I ima1 = io::pgm::load<int_u8>(argv[1]); + I ima2 = io::pgm::load<int_u8>(argv[2]); + + io::pgm::save(min(ima1, ima2), argv[3]); +} Index: sandbox/geraud/Rd/sequential.hh --- sandbox/geraud/Rd/sequential.hh (revision 0) +++ sandbox/geraud/Rd/sequential.hh (revision 0) @@ -0,0 +1,96 @@ +// 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_MORPHO_RD_SEQUENTIAL_HH +# define MLN_MORPHO_RD_SEQUENTIAL_HH + +# include "utils.hh" + + +namespace mln +{ + + namespace morpho + { + + namespace Rd + { + + template <typename I, typename N> + I sequential(const I& f, const I& g, const N& nbh) + { + mln_precondition(f <= g); + + f.name_it("f"); + g.name_it("g"); + + I o_(f.domain()); + o_.name_it("o_"); + + // initialisation + I o(f.domain()); + o.name_it("o"); + level::paste(f, o); + // WAS: I o = clone(f); + + bool stability; + do + { + level::paste(o, o_); // memorisation + + // passe 1 + { + mln_bkd_piter(I) p(f.domain()); + for_all(p) + o(p) = min( max_Nminus(o, p, nbh), g(p) ); + } + + // passe 2 + { + mln_fwd_piter(I) p(f.domain()); + for_all(p) + o(p) = min( max_Nplus(o, p, nbh), g(p) ); + } + + stability = (o == o_); + } + while (not stability); + + print_counts(); + + mln_postcondition(o <= g); + return o; + } + + } // end of namespace mln::morpho::Rd + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_RD_SEQUENTIAL_HH Index: sandbox/geraud/Rd/debase.union_find.hh --- sandbox/geraud/Rd/debase.union_find.hh (revision 0) +++ sandbox/geraud/Rd/debase.union_find.hh (revision 0) @@ -0,0 +1,165 @@ +// 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 OLN_MORPHO_RD_UNION_FIND_HH +# define OLN_MORPHO_RD_UNION_FIND_HH + +# include <oln/morpho/Rd/utils.hh> + + +namespace oln +{ + + namespace morpho + { + + namespace Rd + { + + + template <typename I, typename N> + struct union_find_t + { + typedef oln_point(I) point; + + // in: + I f, g; + N nbh; + + // out: + I o; + + // aux: + std::vector<point> S; + I data; + oln_plain_value(I, bool) is_proc; + oln_plain_value(I, point) parent; + + union_find_t(const I& f, const I& g, const N& nbh) + : f(f), g(g), nbh(nbh) + { + prepare(o, with, f); + prepare(parent, with, f); + prepare(is_proc, with, f); + prepare(data, with, f); + + // init + + std::cout << "0 "; + level::fill(inplace(is_proc), false); + S = histo_reverse_sort(g); + + // first pass + + std::cout << "1 "; + for (unsigned i = 0; i < S.size(); ++i) + { + point p = S[i]; + + make_set(p); + oln_niter(N) n(nbh, p); + for_all(n) + if (f.has(n) and is_proc(n)) + do_union(n, p); + is_proc(p) = true; + } + + // second pass + + std::cout << "2 "; + level::fill(inplace(is_proc), false); + for (int i = S.size() - 1; i >= 0; --i) + { + point p = S[i]; + assert(is_proc(p) == false); + if (parent(p) == p) + o(p) = data(p) == 255 ? g(p) : data(p); + else + { + assert(is_proc(parent(p)) == true); + o(p) = o(parent(p)); + } + is_proc(p) = true; + } + + } + + void make_set(const point& p) + { + parent(p) = p; + data(p) = f(p); + } + + point find_root(const point& x) + { + if (parent(x) == x) + return x; + else + return parent(x) = find_root(parent(x)); + } + + bool equiv(const point& r, const point& p) + { + return g(r) == g(p) or g(p) >= data(r); + } + + void do_union(const point& n, const point& p) + { + point r = find_root(n); + if (r != p) + { + if (equiv(r, p)) + { + parent(r) = p; + if (data(r) > data(p)) + data(p) = data(r); // increasing criterion + } + else + data(p) = 255; + } + } + + }; + + + template <typename I, typename N> + I union_find(const I& f, const I& g, const N& nbh) + { + precondition(f <= g); + union_find_t<I, N> run(f, g, nbh); + return run.o; + } + + + } // end of namespace oln::morpho::Rd + + } // end of namespace oln::morpho + +} // end of namespace oln + + +#endif // ! OLN_MORPHO_RD_UNION_FIND_HH Index: sandbox/geraud/Rd/union_find.cc --- sandbox/geraud/Rd/union_find.cc (revision 0) +++ sandbox/geraud/Rd/union_find.cc (revision 0) @@ -0,0 +1,50 @@ +#include <mln/core/image2d.hh> +#include <mln/core/neighb2d.hh> +#include <mln/value/int_u8.hh> + +#include <mln/debug/println.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "union_find.hh" + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " f.pgm g.pgm c output.pgm" << std::endl + << "reconstruction by dilation (union_find version; may 2007)" << std::endl + << "f = marker (to be dilated)" << std::endl + << "g = mask (constraint >= f)" << std::endl + << "c: 4 or 8" << std::endl; + exit(1); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 5) + usage(argv); + + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::pgm::load<int_u8>(argv[1]); + I g = io::pgm::load<int_u8>(argv[2]); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::pgm::save(morpho::Rd::union_find(f, g, + (c == 4 ? c4() : c8())), + argv[4]); +} Index: sandbox/geraud/Rd/svg.queue_based.hh --- sandbox/geraud/Rd/svg.queue_based.hh (revision 0) +++ sandbox/geraud/Rd/svg.queue_based.hh (revision 0) @@ -0,0 +1,118 @@ +// 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_MORPHO_RD_QUEUE_BASED_HH +# define MLN_MORPHO_RD_QUEUE_BASED_HH + +# include <queue> +# include "utils.hh" + + +namespace mln +{ + + namespace morpho + { + + namespace Rd + { + + template <typename I, typename N> + I queue_based(const I& f, const I& g, const N& nbh, + bool echo = false) + { + if (echo) std::cout << std::endl; + + f.name_it("f"); + g.name_it("g"); + + typedef mln_point(I) point; + std::queue<point> q; + I o; + o.name_it("o"); + + unsigned n_init_pushs = 0, n_body_pushs = 0, n_pops = 0; + + // initialisation + { + o = regional_maxima(f, nbh); + // p in M <=> o(p) != 0 + if (echo) debug::println(o); + + mln_piter(I) p(f.domain()); + mln_niter(N) n(nbh, p); + + for_all(p) if (o(p) != 0) // p in M + for_all(n) if (f.has(n) and o(n) == 0) // n not in M + { + q.push(p); + ++n_init_pushs; + break; + } + } + + // propagation + { + point p; + mln_niter(N) n(nbh, p); + while (not q.empty()) + { + p = q.front(); + if (echo) std::cout << std::endl << "pop " << p << " :"; + q.pop(); + ++n_pops; + for_all(n) if (f.has(n)) + { + if (o(n) < o(p) and o(n) != g(n)) + { + o(n) = min(o(p), g(n)); + if (echo) std::cout << " push " << n; + q.push(n); + ++n_body_pushs; + } + } + } + if (echo) std::cout << std::endl; + } + + std::cout << "n_init_pushs = " << n_init_pushs << std::endl + << "n_body_pushs = " << n_body_pushs << std::endl + << "n_pops = " << n_pops << std::endl; + + print_counts(); + + return o; + } + + } // end of namespace mln::morpho::Rd + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_RD_QUEUE_BASED_HH