[Olena-patches] 978: Add raw code for morphological reconstruction by d.

https://svn.lrde.epita.fr/svn/oln/trunk/olena Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; Add raw code for morphological reconstruction by d. * oln/morpho/Rd: New directory. * oln/morpho/Rd/hybrid.hh, * oln/morpho/Rd/queue_based.hh, * oln/morpho/Rd/parallel.cc, * oln/morpho/Rd/hybrid.cc, * oln/morpho/Rd/union_find.hh, * oln/morpho/Rd/queue_based.cc, * oln/morpho/Rd/sequential.hh, * oln/morpho/Rd/debase.union_find.hh, * oln/morpho/Rd/parallel.hh, * oln/morpho/Rd/diff.cc, * oln/morpho/Rd/union_find.cc, * oln/morpho/Rd/sequential.cc, * oln/morpho/Rd/utils.hh, * oln/io/save_pgm.hh, * oln/io/load_pgm.hh: New. * oln/debug/format.hh: Update. * oln/core/gen/neighb.hh (niter): New associated type. * oln/level/fill.hh (iterator): Add include. core/gen/neighb.hh | 4 debug/format.hh | 15 +- io/load_pgm.hh | 169 +++++++++++++++++++++++++++ io/save_pgm.hh | 68 +++++++++++ level/fill.hh | 1 morpho/Rd/debase.union_find.hh | 165 +++++++++++++++++++++++++++ morpho/Rd/diff.cc | 28 ++++ morpho/Rd/hybrid.cc | 46 +++++++ morpho/Rd/hybrid.hh | 103 ++++++++++++++++ morpho/Rd/parallel.cc | 46 +++++++ morpho/Rd/parallel.hh | 80 +++++++++++++ morpho/Rd/queue_based.cc | 46 +++++++ morpho/Rd/queue_based.hh | 106 +++++++++++++++++ morpho/Rd/sequential.cc | 46 +++++++ morpho/Rd/sequential.hh | 85 +++++++++++++ morpho/Rd/union_find.cc | 46 +++++++ morpho/Rd/union_find.hh | 165 +++++++++++++++++++++++++++ morpho/Rd/utils.hh | 250 +++++++++++++++++++++++++++++++++++++++++ 18 files changed, 1466 insertions(+), 3 deletions(-) Index: oln/debug/format.hh --- oln/debug/format.hh (revision 977) +++ oln/debug/format.hh (working copy) @@ -28,6 +28,9 @@ #ifndef OLN_DEBUG_FORMAT_HH # define OLN_DEBUG_FORMAT_HH +# include <string> +# include <sstream> + namespace oln { @@ -41,7 +44,7 @@ const T& format(const T& value); - unsigned + std::string format(const unsigned char& value); char @@ -56,9 +59,15 @@ return value; } - unsigned format(const unsigned char& value) + std::string format(const unsigned char& value) { - return value; + std::ostringstream ostr; + if (value < 10) + ostr << "00"; + else if (value < 100) + ostr << '0'; + ostr << unsigned(value); + return ostr.str(); } char format(bool value) Index: oln/core/gen/neighb.hh --- oln/core/gen/neighb.hh (revision 977) +++ oln/core/gen/neighb.hh (working copy) @@ -31,6 +31,7 @@ # include <oln/core/internal/dpoints_impl.hh> # include <oln/core/internal/neighborhood_base.hh> +# include <oln/core/gen/dpoints_piter.hh> namespace oln @@ -65,6 +66,9 @@ { public: + typedef stc_type(Dp, point) P; + typedef dpoints_fwd_piter_<P> niter; + neighb_(); neighb_<Dp>& take(const Dp& dp); Index: oln/morpho/Rd/hybrid.hh --- oln/morpho/Rd/hybrid.hh (revision 0) +++ oln/morpho/Rd/hybrid.hh (revision 0) @@ -0,0 +1,103 @@ +// 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_HYBRID_HH +# define OLN_MORPHO_RD_HYBRID_HH + +# include <queue> +# include <oln/morpho/Rd/utils.hh> + + +namespace oln +{ + + namespace morpho + { + + namespace Rd + { + + + template <typename I, typename N> + I hybrid(const I& f, const I& g, const N& nbh, + bool echo = false) + { + typedef oln_point(I) point; + std::queue<point> q; + + // initialisation + I o = level::clone(f); + + // sequence + { + oln_bkd_piter(I) p(f.points()); + for_all(p) + o(p) = min( max_Nminus(o, p,nbh), g(p) ); + } + { + oln_fwd_piter(I) p(f.points()); + oln_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); + } + } + + // propagation + { + point p; + oln_niter(N) n(nbh, p); + while (not q.empty()) + { + p = q.front(); + if (echo) std::cout << std::endl << "pop " << p << " :"; + q.pop(); + 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); + } + } + if (echo) std::cout << std::endl; + } + + return o; + } + + } // end of namespace oln::morpho::Rd + + } // end of namespace oln::morpho + +} // end of namespace oln + + +#endif // ! OLN_MORPHO_RD_HYBRID_HH Index: oln/morpho/Rd/queue_based.hh --- oln/morpho/Rd/queue_based.hh (revision 0) +++ oln/morpho/Rd/queue_based.hh (revision 0) @@ -0,0 +1,106 @@ +// 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_QUEUE_BASED_HH +# define OLN_MORPHO_RD_QUEUE_BASED_HH + +# include <queue> +# include <oln/morpho/Rd/utils.hh> + + +namespace oln +{ + + 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; + + typedef oln_point(I) point; + std::queue<point> q; + I o; + + // initialisation + { + o = regional_maxima(f, nbh); + if (echo) debug::println(o); + + oln_piter(I) p(f.points()); + oln_niter(N) n(nbh, p); + + for_all(p) + if (o(p) != 0) + { + bool ok = false; + for_all(n) if (f.has(n)) + if (o(n) = 0) + ok = true; + if (ok) + q.push(p); + } + } + + // propagation + { + point p; + oln_niter(N) n(nbh, p); + while (not q.empty()) + { + p = q.front(); + if (echo) std::cout << std::endl << "pop " << p << " :"; + q.pop(); + 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); + } + } + } + if (echo) std::cout << std::endl; + } + + return o; + } + + } // end of namespace oln::morpho::Rd + + } // end of namespace oln::morpho + +} // end of namespace oln + + +#endif // ! OLN_MORPHO_RD_QUEUE_BASED_HH Index: oln/morpho/Rd/parallel.cc --- oln/morpho/Rd/parallel.cc (revision 0) +++ oln/morpho/Rd/parallel.cc (revision 0) @@ -0,0 +1,46 @@ +#include <oln/core/2d/image2d.hh> +#include <oln/core/2d/neighb2d.hh> + +#include <oln/debug/println.hh> +#include <oln/io/load_pgm.hh> +#include <oln/io/save_pgm.hh> +#include <oln/morpho/Rd/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 oln; + typedef image2d<unsigned char> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::load_pgm(argv[1]); + I g = io::load_pgm(argv[2]); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::save_pgm(morpho::Rd::parallel(f, g, + (c = 4 ? c4 : c8)), + argv[4]); + return 0; +} Index: oln/morpho/Rd/hybrid.cc --- oln/morpho/Rd/hybrid.cc (revision 0) +++ oln/morpho/Rd/hybrid.cc (revision 0) @@ -0,0 +1,46 @@ +#include <oln/core/2d/image2d.hh> +#include <oln/core/2d/neighb2d.hh> + +#include <oln/debug/println.hh> +#include <oln/io/load_pgm.hh> +#include <oln/io/save_pgm.hh> +#include <oln/morpho/Rd/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 oln; + typedef image2d<unsigned char> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::load_pgm(argv[1]); + I g = io::load_pgm(argv[2]); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::save_pgm(morpho::Rd::hybrid(f, g, + (c = 4 ? c4 : c8)), + argv[4]); + return 0; +} Index: oln/morpho/Rd/union_find.hh --- oln/morpho/Rd/union_find.hh (revision 0) +++ oln/morpho/Rd/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: oln/morpho/Rd/queue_based.cc --- oln/morpho/Rd/queue_based.cc (revision 0) +++ oln/morpho/Rd/queue_based.cc (revision 0) @@ -0,0 +1,46 @@ +#include <oln/core/2d/image2d.hh> +#include <oln/core/2d/neighb2d.hh> + +#include <oln/debug/println.hh> +#include <oln/io/load_pgm.hh> +#include <oln/io/save_pgm.hh> +#include <oln/morpho/Rd/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 oln; + typedef image2d<unsigned char> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::load_pgm(argv[1]); + I g = io::load_pgm(argv[2]); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::save_pgm(morpho::Rd::queue_based(f, g, + (c = 4 ? c4 : c8)), + argv[4]); + return 0; +} Index: oln/morpho/Rd/sequential.hh --- oln/morpho/Rd/sequential.hh (revision 0) +++ oln/morpho/Rd/sequential.hh (revision 0) @@ -0,0 +1,85 @@ +// 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_SEQUENTIAL_HH +# define OLN_MORPHO_RD_SEQUENTIAL_HH + +# include <oln/morpho/Rd/utils.hh> + + +namespace oln +{ + + namespace morpho + { + + namespace Rd + { + + template <typename I, typename N> + I sequential(const I& f, const I& g, const N& nbh) + { + precondition(f <= g); + I o_(f.points()); + + // initialisation + I o = level::clone(f); + + bool stability; + do + { + level::paste(o, inplace(o_)); // memorisation + + // passe 1 + { + oln_bkd_piter(I) p(f.points()); + for_all(p) + o(p) = min( max_Nminus(o, p,nbh), g(p) ); + } + + // passe 2 + { + oln_fwd_piter(I) p(f.points()); + for_all(p) + o(p) = min( max_Nplus(o, p,nbh), g(p) ); + } + + stability = (o = o_); + } + while (not stability); + postcondition(o <= g); + return o; + } + + } // end of namespace oln::morpho::Rd + + } // end of namespace oln::morpho + +} // end of namespace oln + + +#endif // ! OLN_MORPHO_RD_SEQUENTIAL_HH Index: oln/morpho/Rd/debase.union_find.hh --- oln/morpho/Rd/debase.union_find.hh (revision 0) +++ oln/morpho/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: oln/morpho/Rd/parallel.hh --- oln/morpho/Rd/parallel.hh (revision 0) +++ oln/morpho/Rd/parallel.hh (revision 0) @@ -0,0 +1,80 @@ +// 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_PARALLEL_HH +# define OLN_MORPHO_RD_PARALLEL_HH + +# include <oln/morpho/Rd/utils.hh> + + +namespace oln +{ + + namespace morpho + { + + namespace Rd + { + + template <typename I, typename N> + I parallel(const I& f, const I& g, const N& nbh) + { + precondition(f <= g); + I o_(f.points()); + oln_piter(I) p(f.points()); + + // initialisation + I o = level::clone(f); + + bool stability; + do + { + level::paste(o, inplace(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); + + postcondition(o <= g); + return o; + } + + } // end of namespace oln::morpho::Rd + + } // end of namespace oln::morpho + +} // end of namespace oln + + +#endif // ! OLN_MORPHO_RD_PARALLEL_HH Index: oln/morpho/Rd/diff.cc --- oln/morpho/Rd/diff.cc (revision 0) +++ oln/morpho/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: oln/morpho/Rd/union_find.cc --- oln/morpho/Rd/union_find.cc (revision 0) +++ oln/morpho/Rd/union_find.cc (revision 0) @@ -0,0 +1,46 @@ +#include <oln/core/2d/image2d.hh> +#include <oln/core/2d/neighb2d.hh> + +#include <oln/debug/println.hh> +#include <oln/io/load_pgm.hh> +#include <oln/io/save_pgm.hh> +#include <oln/morpho/Rd/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 oln; + typedef image2d<unsigned char> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::load_pgm(argv[1]); + I g = io::load_pgm(argv[2]); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::save_pgm(morpho::Rd::union_find(f, g, + (c = 4 ? c4 : c8)), + argv[4]); + return 0; +} Index: oln/morpho/Rd/sequential.cc --- oln/morpho/Rd/sequential.cc (revision 0) +++ oln/morpho/Rd/sequential.cc (revision 0) @@ -0,0 +1,46 @@ +#include <oln/core/2d/image2d.hh> +#include <oln/core/2d/neighb2d.hh> + +#include <oln/debug/println.hh> +#include <oln/io/load_pgm.hh> +#include <oln/io/save_pgm.hh> +#include <oln/morpho/Rd/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 oln; + typedef image2d<unsigned char> I; + + int c = atoi(argv[3]); + if (c != 4 and c != 8) + usage(argv); + + I f = io::load_pgm(argv[1]); + I g = io::load_pgm(argv[2]); + + if (not (f <= g)) + { + std::cerr << "pb" << std::endl; + return 1; + } + + io::save_pgm(morpho::Rd::sequential(f, g, + (c = 4 ? c4 : c8)), + argv[4]); + return 0; +} Index: oln/morpho/Rd/utils.hh --- oln/morpho/Rd/utils.hh (revision 0) +++ oln/morpho/Rd/utils.hh (revision 0) @@ -0,0 +1,250 @@ +// 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_UTILS_HH +# define OLN_MORPHO_RD_UTILS_HH + +# include <vector> + +# include <oln/core/concept/image.hh> +# include <oln/core/internal/f_ch_value.hh> + +# include <oln/level/fill.hh> +# include <oln/level/paste.hh> +# include <oln/level/clone.hh> +# include <oln/level/compare.hh> + + + +namespace oln +{ + + namespace morpho + { + + namespace Rd + { + + + template <typename T> + T min(T v1, T v2) + { + return v1 < v2 ? v1 : v2; + } + + + template <typename I> + std::vector<unsigned> compute_histo(const I& ima) + { + std::vector<unsigned> h(256, 0); + oln_piter(I) p(ima.points()); + for_all(p) + ++h[ima(p)]; + return h; + } + + + template <typename I> + std::vector<oln_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<oln_point(I)> vec(ima.points().npoints()); + // storing output data + oln_piter(I) p(ima.points()); + for_all(p) + vec[loc[ima(p)]++] = p; + return vec; + } + + + template <typename I> + std::vector<oln_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<oln_point(I)> vec(ima.points().npoints()); + // storing output data + oln_piter(I) p(ima.points()); + for_all(p) + vec[loc[ima(p)]++] = p; + return vec; + } + + + template <typename I, typename P, typename N> + oln_value(I) max_Nplus(const I& ima, const P& p, const N& nbh) + { + oln_value(I) v = ima(p); + oln_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> + oln_value(I) max_Nminus(const I& ima, const P& p, const N& nbh) + { + oln_value(I) v = ima(p); + oln_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> + oln_value(I) max_N(const I& ima, const P& p, const N& nbh) + { + oln_value(I) v = ima(p); + oln_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 N> + struct regional_maxima_t + { + typedef oln_point(I) point; + typedef oln_plain_value(I, bool) image_bool; + typedef oln_plain_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) + { + prepare(o, with, f); + prepare(parent, with, f); + prepare(attr, with, f); + prepare(is_proc, with, f); + + // init + + level::fill(inplace(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); + oln_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 + + for (int i = S.size() - 1; i >= 0; --i) + { + point p = S[i]; + if (parent(p) = p) + o(p) = attr(p) ? f(p) : 0; + 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 oln::morpho::Rd + + } // end of namespace oln::morpho + +} // end of namespace oln + + +#endif // ! OLN_MORPHO_RD_UTILS_HH Index: oln/level/fill.hh --- oln/level/fill.hh (revision 977) +++ oln/level/fill.hh (working copy) @@ -33,6 +33,7 @@ # include <oln/core/concept/image.hh> # include <oln/core/concept/point.hh> +# include <oln/core/concept/iterator.hh> # include <oln/core/gen/fun.hh> # include <oln/core/gen/inplace.hh> Index: oln/io/save_pgm.hh --- oln/io/save_pgm.hh (revision 0) +++ oln/io/save_pgm.hh (revision 0) @@ -0,0 +1,68 @@ +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 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_IO_SAVE_PGM_HH +# define OLN_IO_SAVE_PGM_HH + +# include <iostream> +# include <fstream> + +# include <oln/core/2d/image2d.hh> + + +namespace oln +{ + + namespace io + { + + void save_pgm(const image2d<unsigned char>& ima, const std::string& filename) + { + std::ofstream file(filename.c_str()); + if (not file) + { + std::cerr << "error: cannot open file '" << filename + << "'!"; + abort(); + } + file << "P5" << std::endl; + file << "# olena" << std::endl; + file << ima.ncols() << ' ' << ima.nrows() << std::endl; + file << "255" << std::endl; + int col = ima.min_col(); + size_t len = ima.ncols() * sizeof(unsigned char); + for (int row = ima.min_row(); row <= ima.max_row(); ++row) + file.write((char*)(&(ima.at(row, col))), len); + } + + } // end of namespace oln::io + +} // end of namespace oln + + +#endif // ! OLN_IO_SAVE_PGM_HH Index: oln/io/load_pgm.hh --- oln/io/load_pgm.hh (revision 0) +++ oln/io/load_pgm.hh (revision 0) @@ -0,0 +1,169 @@ +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 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_IO_LOAD_PGM_HH +# define OLN_IO_LOAD_PGM_HH + +# include <iostream> +# include <fstream> +# include <string> + +# include <oln/core/2d/image2d.hh> + + +namespace oln { + + namespace io { + + namespace internal { + + void abort() + { + std::cerr << " aborting." << std::endl; + exit(0); + } + + bool read_pnm_header(std::istream& istr, + char& type, + int& nrows, int& ncols, + bool test = false) + { + // check magic + if (istr.get() != 'P' ) + goto err; + type = istr.get(); + if (type < '1' or type > '6') + goto err; + if (istr.get() != '\n') + goto err; + + // skip comments + while (istr.peek() = '#') + { + std::string line; + std::getline(istr, line); + } + + // get size + istr >> ncols >> nrows; + if (nrows <= 0 or ncols <= 0) + goto err; + + // skip maxvalue + if (istr.get() != '\n') + goto err; + if (type != '1' and type != '4') + { + std::string line; + std::getline(istr, line); + } + return true; + + err: + if (not test) + { + std::cerr << "error: badly formed header!"; + abort(); + } + return false; + } + + void read_pnm_header(char ascii, char raw, + std::istream& istr, + char& type, + int& nrows, int& ncols) + { + read_pnm_header(istr, type, nrows, ncols); + if (not (type = ascii or type = raw)) + { + std::cerr << "error: bad pnm type; " + << "expected P" << ascii + << " or P" << raw + << ", get P" << type << "!"; + abort(); + } + } + + + /// load_ascii. + template <typename I> + void load_pgm_ascii(std::ifstream& file, I& ima) + { + oln_fwd_piter(I) p(ima.points()); + for_all(p) + { + unsigned value; + file >> value; + ima(p) = value; + // FIXME: Test alt code below. + // file >> ima(p); + } + } + + + /// load_raw_2d. + template <typename I> + void load_pgm_raw_2d(std::ifstream& file, I& ima) + { + int col = ima.min_col(); + size_t len = ima.ncols() * sizeof(oln_value(I)); + for (int row = ima.min_row(); row <= ima.max_row(); ++row) + file.read((char*)(&(ima.at(row, col))), len); + } + + + } // end of namespace oln::io::internal + + + image2d<unsigned char> load_pgm(const std::string& filename) + { + std::ifstream file(filename.c_str()); + if (not file) + { + std::cerr << "error: file '" << filename + << "' not found!"; + abort(); + } + char type; + int nrows, ncols; + internal::read_pnm_header('2', '5', file, type, nrows, ncols); + image2d<unsigned char> ima(nrows, ncols); + if (type = '5') + internal::load_pgm_raw_2d(file, ima); + else + if (type = '2') + internal::load_pgm_ascii(file, ima); + return ima; + } + + } // end of namespace oln::io + +} // end of namespace oln + + +#endif // ! OLN_IO_LOAD_PGM_HH