https://svn.lrde.epita.fr/svn/oln/trunk/olena
Index: ChangeLog
from Thierry Geraud <thierry.geraud(a)lrde.epita.fr>
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
