Careful, some implementations may not work (especially the ones based on complex pointer manipulations). --- milena/apps/bench/Makefile.am | 30 +- .../{static_window.hh => static_neighborhood.hh} | 200 +++---- milena/apps/bench/trait.hh | 31 +- milena/apps/bench/watershed-lena.cc | 116 ++++ milena/apps/bench/watershed.hh | 580 ++++++++++++++++++++ 5 files changed, 834 insertions(+), 123 deletions(-) copy milena/apps/bench/{static_window.hh => static_neighborhood.hh} (58%) create mode 100644 milena/apps/bench/watershed-lena.cc create mode 100644 milena/apps/bench/watershed.hh diff --git a/milena/apps/bench/Makefile.am b/milena/apps/bench/Makefile.am index d2391c1..c75ea60 100644 --- a/milena/apps/bench/Makefile.am +++ b/milena/apps/bench/Makefile.am @@ -1,4 +1,5 @@ -# Copyright (C) 2010, 2011 EPITA Research and Development Laboratory (LRDE). +# Copyright (C) 2010, 2011, 2012 EPITA Research and Development +# Laboratory (LRDE). # # This file is part of Olena. # @@ -14,6 +15,8 @@ # You should have received a copy of the GNU General Public License # along with Olena. If not, see <http://www.gnu.org/licenses/>. +# Benchmarks for the ICIP 2011 submission and the GRETSI 2011 paper. + # Find Milena headers. AM_CPPFLAGS = -I$(top_srcdir)/milena -I$(top_builddir)/milena # Produce fast code. @@ -21,6 +24,10 @@ AM_CXXFLAGS = $(APPS_CXXFLAGS_NODEBUG) EXTRA_DIST = lena1024.pgm lena2048.pgm +## ---------- ## +## Dilation. ## +## ---------- ## + noinst_PROGRAMS = \ dilation-lena \ dilation-lena-bench-nongen \ @@ -124,3 +131,24 @@ test-dilation-lena-bench: test-dilation-lena-bench.in Makefile CLEANFILES = test-dilation-lena-bench TESTS = test-dilation-lena-bench + + +## --------------------- ## +## Watershed transform. ## +## --------------------- ## + +noinst_PROGRAMS += \ + watershed-lena + +watershed_lena_SOURCES = watershed-lena.cc + +MOSTLYCLEANFILES += \ + watershed-lena-out-512-nongen.pgm \ + watershed-lena-out-512-gen.pgm \ + watershed-lena-out-512-fast.pgm \ + watershed-lena-out-1024-nongen.pgm \ + watershed-lena-out-1024-gen.pgm \ + watershed-lena-out-1024-fast.pgm \ + watershed-lena-out-2048-nongen.pgm \ + watershed-lena-out-2048-gen.pgm \ + watershed-lena-out-2048-fast.pgm diff --git a/milena/apps/bench/static_window.hh b/milena/apps/bench/static_neighborhood.hh similarity index 58% copy from milena/apps/bench/static_window.hh copy to milena/apps/bench/static_neighborhood.hh index 5a8ce21..aa64a98 100644 --- a/milena/apps/bench/static_window.hh +++ b/milena/apps/bench/static_neighborhood.hh @@ -1,5 +1,5 @@ -// Copyright (C) 2007, 2008, 2009, 2010, 2011 EPITA Research and Development -// Laboratory (LRDE) +// Copyright (C) 2007, 2008, 2009, 2010, 2011, 2012 EPITA Research and +// Development Laboratory (LRDE) // // This file is part of Olena. // @@ -24,14 +24,14 @@ // exception does not however invalidate any other reasons why the // executable file might be covered by the GNU General Public License. -#ifndef MLN_CORE_STATIC_WINDOW_HH -# define MLN_CORE_STATIC_WINDOW_HH +#ifndef MLN_CORE_STATIC_NEIGHBORHOOD_HH +# define MLN_CORE_STATIC_NEIGHBORHOOD_HH // FIXME: Review and update the documentation. /*! \file * - * \brief Definition of the generic window class mln::window. + * \brief Definition of the generic neighborhood class mln::neighborhood. * * \todo Make naming coherent: we have window (without '_') but * point_, neighb_, etc. @@ -43,7 +43,6 @@ * \todo Rename dps_hook_ as util_set_. */ -# include <mln/core/internal/window_base.hh> # include <mln/core/concept/gdpoint.hh> # include <mln/metal/is_a.hh> @@ -53,120 +52,91 @@ # include <mln/literal/zero.hh> # include "apps/bench/static_array.hh" +# include "apps/bench/static_window.hh" namespace mln { // Forward declarations. - template <typename D, unsigned n> class static_window; + template <typename D, unsigned n> class static_neighborhood; template <typename V> class dpsites_fwd_piter; template <typename V> class dpsites_bkd_piter; - namespace trait - { - - template <typename D, unsigned n> - struct window_< mln::static_window<D, n> > - { - typedef trait::window::size::fixed size; - typedef trait::window::support::regular support; - typedef trait::window::definition::unique definition; - }; - - } // end of namespace mln::trait - - - - /*! \brief Generic window class. + /*! \brief Generic neighborhood class. * - * This type of window is just like a set of delta-points. The + * This type of neighborhood is just like a set of delta-points. The * parameter is \c D, type of delta-point. */ template <typename D, unsigned n> - class static_window : public internal::window_base< D, static_window<D, n> > + class static_neighborhood + : public Neighborhood< static_neighborhood<D, n> > { public: enum { length = n }; - /// Regular window associated type. - typedef static_window<D, n> regular; - + /// Window associated type. + typedef static_window<D, n> window; - static_window(const mln::util::static_array<D, n>& dps); + /// Get the corresponding window. + window win() const; - /*! \brief Test if the window is centered. - * - * \return True if the delta-point 0 belongs to the window. - */ - bool is_centered() const; + static_neighborhood(const mln::util::static_array<D, n>& dps); - /*! Test if the window is symmetric. - * - * \return True if for every dp of this window, -dp is also in - * this window. - */ - bool is_symmetric() const; - - /// Apply a central symmetry to the target window. - void sym(); - - - /*! \brief Site_Iterator type to browse the points of a basic window + /*! \brief Site_Iterator type to browse the points of a basic neighborhood * w.r.t. the ordering of delta-points. */ - typedef dpsites_fwd_piter< static_window<D, n> > fwd_qiter; + typedef dpsites_fwd_piter< static_neighborhood<D, n> > fwd_niter; - /*! \brief Site_Iterator type to browse the points of a basic window + /*! \brief Site_Iterator type to browse the points of a basic neighborhood * w.r.t. the reverse ordering of delta-points. */ - typedef dpsites_bkd_piter< static_window<D, n> > bkd_qiter; + typedef dpsites_bkd_piter< static_neighborhood<D, n> > bkd_niter; - /*! \brief Site_Iterator type to browse the points of a basic window + /*! \brief Site_Iterator type to browse the points of a basic neighborhood * whatever the ordering of delta-points. */ - typedef fwd_qiter qiter; - + typedef fwd_niter niter; - /// Give the window size, i.e., the number of delta-sites. + /// Give the neighborhood size, i.e., the number of delta-sites. unsigned size() const; - /*! \brief Test if the window is empty (null size; no delta-point). + /*! \brief Test if the neighborhood is empty (null size; no delta-point). */ bool is_empty() const; - // /// Clear the window. + // /// Clear the neighborhood. // void clear(); - /*! \brief Give the maximum coordinate gap between the window - center and a window point. + /*! \brief Give the maximum coordinate gap between the neighborhood + center and a neighborhood point. */ unsigned delta() const; /// Give the \p i-th delta-point. const D& dp(unsigned i) const; - /// Test if \p dp is in this window definition. + /// Test if \p dp is in this neighborhood definition. bool has(const D& dp) const; // /// Insert a delta-point \p dp. - // static_window<D, n>& insert(const D& dp); + // static_neighborhood<D, n>& insert(const D& dp); - // /// Insert another window \p win. + // /// Insert another neighborhood \p nbh. // template <typename W> - // static_window<D, n>& insert(const Window<W>& win); + // static_neighborhood<D, n>& insert(const Neighborhood<W>& nbh); // /// \{ Insertion of a delta-point with different numbers of // /// arguments (coordinates) w.r.t. the dimension. - // static_window<D, n>& insert(const mln_coord(D)& dind); // For 1D or index access. + // static_neighborhood<D, n>& insert(const mln_coord(D)& dind); // For 1D or index access. - // static_window<D, n>& insert(const mln_coord(D)& drow, + // static_neighborhood<D, n>& insert(const mln_coord(D)& drow, // const mln_coord(D)& dcol); // For 2D. - // static_window<D, n>& insert(const mln_coord(D)& dsli, + // static_neighborhood<D, n>& insert(const mln_coord(D)& dsli, // const mln_coord(D)& drow, // const mln_coord(D)& dcol); // For 3D. /// \} @@ -175,7 +145,10 @@ namespace mln /// Hook to the set of D. const mln::util::static_array<D, n>& dps_hook_() const; - /// Print the window definition into \p ostr. + /// Return true by default. + bool is_valid() const; + + /// Print the neighborhood definition into \p ostr. void print(std::ostream& ostr) const; private: @@ -188,22 +161,22 @@ namespace mln - /*! \brief Equality comparison between windows \p lhs and \p rhs. + /*! \brief Equality comparison between neighborhoods \p lhs and \p rhs. * - * \relates mln::window + * \relates mln::neighborhood */ template <typename D, unsigned n> - bool operator==(const static_window<D, n>& lhs, const static_window<D, n>& rhs); + bool operator==(const static_neighborhood<D, n>& lhs, const static_neighborhood<D, n>& rhs); # ifndef MLN_INCLUDE_ONLY - // static_window<D, n> + // static_neighborhood<D, n> template <typename D, unsigned n> inline - static_window<D, n>::static_window(const mln::util::static_array<D, n>& dps) + static_neighborhood<D, n>::static_neighborhood(const mln::util::static_array<D, n>& dps) : dps_(dps) { // FIXME HERE: Was: mln::metal::is_a<D, Dpoint>::check(); @@ -212,39 +185,16 @@ namespace mln template <typename D, unsigned n> inline - bool - static_window<D, n>::is_symmetric() const + static_window<D, n> + static_neighborhood<D, n>::win() const { - static_window<D, n> cpy = *this; - cpy.sym(); - return cpy == *this; + return static_window<D, n>(dps_); } template <typename D, unsigned n> inline bool - static_window<D, n>::is_centered() const - { - return this->dps_.has(literal::zero); - } - - template <typename D, unsigned n> - inline - void - static_window<D, n>::sym() - { - util::static_array<D, n> rev; - // FIXME: Can't we use std::copy and reverse_iterators here? - for (std::size_t i = 0; i < n; ++i) - rev[i] = dps_[n - 1 - i]; - static_window<D, n> tmp(rev); - *this = tmp; - } - - template <typename D, unsigned n> - inline - bool - static_window<D, n>::is_empty() const + static_neighborhood<D, n>::is_empty() const { return n == 0; } @@ -252,7 +202,7 @@ namespace mln // template <typename D, unsigned n> // inline // void - // static_window<D, n>::clear() + // static_neighborhood<D, n>::clear() // { // dps_.clear(); // } @@ -260,7 +210,7 @@ namespace mln template <typename D, unsigned n> inline unsigned - static_window<D, n>::delta() const + static_neighborhood<D, n>::delta() const { unsigned d = 0; for (unsigned i = 0; i < n; ++i) @@ -275,7 +225,7 @@ namespace mln template <typename D, unsigned n> inline unsigned - static_window<D, n>::delta_(int i) const + static_neighborhood<D, n>::delta_(int i) const { return i; } @@ -283,7 +233,7 @@ namespace mln template <typename D, unsigned n> inline unsigned - static_window<D, n>::delta_(const Gdpoint<D>& dp) const + static_neighborhood<D, n>::delta_(const Gdpoint<D>& dp) const { return norm::linfty(exact(dp).to_vec()); } @@ -291,7 +241,7 @@ namespace mln template <typename D, unsigned n> inline unsigned - static_window<D, n>::size() const + static_neighborhood<D, n>::size() const { return n; } @@ -299,7 +249,7 @@ namespace mln template <typename D, unsigned n> inline const D& - static_window<D, n>::dp(unsigned i) const + static_neighborhood<D, n>::dp(unsigned i) const { mln_precondition(i < n); return dps_[i]; @@ -308,15 +258,15 @@ namespace mln template <typename D, unsigned n> inline bool - static_window<D, n>::has(const D& dp) const + static_neighborhood<D, n>::has(const D& dp) const { return dps_.has(dp); } // template <typename D, unsigned n> // inline - // static_window<D, n>& - // static_window<D, n>::insert(const D& dp) + // static_neighborhood<D, n>& + // static_neighborhood<D, n>::insert(const D& dp) // { // dps_.insert(dp); // return *this; @@ -325,20 +275,20 @@ namespace mln // template <typename D, unsigned n> // template <typename W> // inline - // static_window<D, n>& - // static_window<D, n>::insert(const Window<W>& win_) + // static_neighborhood<D, n>& + // static_neighborhood<D, n>::insert(const Neighborhood<W>& nbh_) // { - // const W& win = exact(win_); - // const unsigned n = win.size(); + // const W& nbh = exact(nbh_); + // const unsigned n = nbh.size(); // for (unsigned i = 0; i < n; ++i) - // dps_.insert(win.dp(i)); + // dps_.insert(nbh.dp(i)); // return *this; // } // template <typename D, unsigned n> // inline - // static_window<D, n>& - // static_window<D, n>::insert(const mln_coord(D)& dind) + // static_neighborhood<D, n>& + // static_neighborhood<D, n>::insert(const mln_coord(D)& dind) // { // mlc_bool(D::dim == 1)::check(); // D dp(dind); @@ -347,8 +297,8 @@ namespace mln // template <typename D, unsigned n> // inline - // static_window<D, n>& - // static_window<D, n>::insert(const mln_coord(D)& drow, + // static_neighborhood<D, n>& + // static_neighborhood<D, n>::insert(const mln_coord(D)& drow, // const mln_coord(D)& dcol) // { // mlc_bool(D::dim == 2)::check(); @@ -358,8 +308,8 @@ namespace mln // template <typename D, unsigned n> // inline - // static_window<D, n>& - // static_window<D, n>::insert(const mln_coord(D)& dsli, + // static_neighborhood<D, n>& + // static_neighborhood<D, n>::insert(const mln_coord(D)& dsli, // const mln_coord(D)& drow, // const mln_coord(D)& dcol) // { @@ -371,15 +321,23 @@ namespace mln template <typename D, unsigned n> inline const mln::util::static_array<D, n>& - static_window<D, n>::dps_hook_() const + static_neighborhood<D, n>::dps_hook_() const { return dps_; } template <typename D, unsigned n> inline + bool + static_neighborhood<D, n>::is_valid() const + { + return true; + } + + template <typename D, unsigned n> + inline void - static_window<D, n>::print(std::ostream& ostr) const + static_neighborhood<D, n>::print(std::ostream& ostr) const { ostr << dps_; } @@ -389,7 +347,7 @@ namespace mln template <typename D, unsigned n> bool - operator==(const static_window<D, n>& lhs, const static_window<D, n>& rhs) + operator==(const static_neighborhood<D, n>& lhs, const static_neighborhood<D, n>& rhs) { return lhs.dps_hook_() == rhs.dps_hook_(); } @@ -402,4 +360,4 @@ namespace mln # include <mln/core/dpsites_piter.hh> -#endif // ! MLN_CORE_STATIC_WINDOW_HH +#endif // ! MLN_CORE_STATIC_NEIGHBORHOOD_HH diff --git a/milena/apps/bench/trait.hh b/milena/apps/bench/trait.hh index dd2d7e5..c75e842 100644 --- a/milena/apps/bench/trait.hh +++ b/milena/apps/bench/trait.hh @@ -1,4 +1,4 @@ -// Copyright (C) 2011 EPITA Research and Development Laboratory (LRDE) +// Copyright (C) 2011, 2012 EPITA Research and Development Laboratory (LRDE) // // This file is part of Olena. // @@ -30,6 +30,7 @@ # include <mln/metal/none.hh> # include "apps/bench/static_window.hh" +# include "apps/bench/static_neighborhood.hh" # include "apps/bench/static_dpoints_pixter.hh" @@ -66,6 +67,34 @@ namespace mln typedef static_dpoints_bkd_pixter< const image2d<T>, static_window<D, n> > ret; }; + + // nixter + + template <typename T, typename D, unsigned n> + struct fwd_nixter< image2d<T>, static_neighborhood<D, n> > + { + typedef static_dpoints_fwd_pixter< image2d<T>, static_neighborhood<D, n> > ret; + }; + + template <typename T, typename D, unsigned n> + struct fwd_nixter< const image2d<T>, static_neighborhood<D, n> > + { + typedef static_dpoints_fwd_pixter< const image2d<T>, static_neighborhood<D, n> > ret; + }; + + template <typename T, typename D, unsigned n> + struct bkd_nixter< image2d<T>, static_neighborhood<D, n> > + { + typedef static_dpoints_bkd_pixter< image2d<T>, static_neighborhood<D, n> > ret; + }; + + template <typename T, typename D, unsigned n> + struct bkd_nixter< const image2d<T>, static_neighborhood<D, n> > + { + typedef static_dpoints_bkd_pixter< const image2d<T>, static_neighborhood<D, n> > ret; + }; + + // FIXME: Also handle mln::image1d<T> and mln::image3d<T>. diff --git a/milena/apps/bench/watershed-lena.cc b/milena/apps/bench/watershed-lena.cc new file mode 100644 index 0000000..2c4710a --- /dev/null +++ b/milena/apps/bench/watershed-lena.cc @@ -0,0 +1,116 @@ +// Copyright (C) 2010, 2011, 2012 EPITA Research and Development +// Laboratory (LRDE) +// +// This file is part of Olena. +// +// Olena is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free +// Software Foundation, version 2 of the License. +// +// Olena 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 Olena. If not, see <http://www.gnu.org/licenses/>. +// +// As a special exception, you may use this file as part of a free +// software project 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. + +#include <iostream> + +#include <mln/core/alias/neighb2d.hh> + +#include <mln/morpho/gradient.hh> + +#include <mln/value/label.hh> +#include <mln/data/stretch.hh> + +#include "apps/bench/watershed.hh" +#include "apps/data.hh" + + +// Shortcut macros for run. + +#define WATERSHED_WITH_BUILTIN_WINDOW(Namespace, Suffix, Headline) \ + do \ + { \ + t.start(); \ + w = Namespace::watershed(g, nbasins); \ + t.stop(); \ + std::cout << Headline << t.read() << " s" << std::endl; \ + io::pgm::save(data::stretch(int_u8(), w), \ + prefix + '-' + length + '-' + Suffix + ".pgm"); \ + } \ + while (0) + +#define WATERSHED(Namespace, Nbh, Suffix, Headline) \ + do \ + { \ + t.start(); \ + w = Namespace::watershed(g, Nbh, nbasins); \ + t.stop(); \ + std::cout << Headline << t.read() << " s" << std::endl; \ + io::pgm::save(data::stretch(int_u8(), w), \ + prefix + '-' + length + '-' + Suffix + ".pgm"); \ + } \ + while (0) + + +void +run(const std::string& filename, const std::string& length) +{ + using namespace mln; + using value::int_u8; + using value::label; + + border::thickness = 1; + image2d<int_u8> lena; + io::pgm::load(lena, filename); + image2d<int_u8> g = morpho::gradient(lena, win_c4p()); + + typedef label<24> label_t; + label_t nbasins; + image2d<label_t> w; + util::timer t; + + std::string prefix = "watershed-lena-out"; + std::cout << "== " << filename << std::endl; + + WATERSHED_WITH_BUILTIN_WINDOW(nongen, "nongen", "nongen\t\t"); + // FIXME: Does not work yet... + // WATERSHED_WITH_BUILTIN_WINDOW(nongen_3ptr, "nongen_3ptr", "nongen_3ptr\t"); + + WATERSHED(gen, c4(), "gen", "gen\t\t"); + // // FIXME: Introduce a new test case, gen_static, using a static window + // // and static_qiters. + WATERSHED(fast, c4(), "fast", "fast\t\t"); + + // Static window, neighborhood and qixters. + const unsigned n = 4; + mln::dpoint2d dps[n] = { mln::dpoint2d( 0, -1), + mln::dpoint2d(-1, 0), + mln::dpoint2d(+1, 0), + mln::dpoint2d( 0, +1) }; + mln::util::static_array<mln::dpoint2d, n> sa(dps, dps + n); + mln::static_neighborhood<mln::dpoint2d, n> static_nbh_c4 (sa); + + WATERSHED(fast_static, static_nbh_c4, "fast_static", "fast_static\t"); + + std::cout << std::endl; +} + +int +main () +{ + run(MLN_IMG_DIR "/lena.pgm", "512"); + run(MLN_APPS_DIR "/bench/lena1024.pgm", "1024"); + run(MLN_APPS_DIR "/bench/lena2048.pgm", "2048"); +} diff --git a/milena/apps/bench/watershed.hh b/milena/apps/bench/watershed.hh new file mode 100644 index 0000000..d808647 --- /dev/null +++ b/milena/apps/bench/watershed.hh @@ -0,0 +1,580 @@ +// Copyright (C) 2008, 2009, 2012 EPITA Research and Development +// Laboratory (LRDE) +// +// This file is part of Olena. +// +// Olena is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free +// Software Foundation, version 2 of the License. +// +// Olena 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 Olena. If not, see <http://www.gnu.org/licenses/>. +// +// As a special exception, you may use this file as part of a free +// software project 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 APPS_BENCH_WATERSHED_HH +# define APPS_BENCH_WATERSHED_HH + +/// \file +/// \brief Watershed transform by flooding benchmark cases. +/// +/// Generic code extracted from mln/morpho/watershed/flooding.hh and +/// modified. + +# include <mln/core/image/image2d.hh> +# include <mln/core/alias/window2d.hh> + +# include <mln/value/int_u8.hh> +# include <mln/value/label.hh> + +# include <mln/io/pgm/load.hh> +# include <mln/io/pgm/save.hh> + +# include <mln/literal/zero.hh> +# include <mln/labeling/regional_minima.hh> + +# include <mln/core/site_set/p_queue_fast.hh> +# include <mln/core/site_set/p_priority.hh> + +# include <mln/extension/adjust_fill.hh> + +# include <mln/util/timer.hh> + +# include "apps/bench/static_neighborhood.hh" +# include "apps/bench/static_dpoints_pixter.hh" +# include "apps/bench/trait.hh" + + +/* FIXME: We should not take into account the time passed into + labeling::regional_minima, since this function is the same for in + versions of the algorithm. */ + +/*----------------------. +| Non-generic version. | +`----------------------*/ + +namespace nongen +{ + using namespace mln; + + typedef mln::image2d<mln::value::int_u8> image; + typedef mln::value::label<24> label; + typedef mln::image2d<label> labels; + + labels + watershed(const image& input, label& n_basins) + { + // FIXME: Simplify the body of this routine even more? + + typedef label marker; + const marker unmarked = 0; + + typedef mln_value_(image) V; + const V max = mln_max(V); + + // Initialize the output with the markers (minima components). + typedef labels O; + O output = labeling::regional_minima(input, c4(), n_basins); + + // Ordered queue. + typedef p_queue_fast<point2d> Q; + p_priority<V, Q> queue; + + // In_queue structure to avoid processing sites several times. + mln_ch_value_(image, bool) in_queue; + initialize(in_queue, input); + data::fill(in_queue, false); + + + // Insert every neighbor (R, C) of every marked area in a + // hierarchical queue, with a priority level corresponding to + // the grey level input.at_(R, C). + for (unsigned int r = 0; r < input.nrows(); ++r) // Iterate on rows. + for (unsigned int c = 0; c < input.ncols(); ++c) // Iterate on columns. + { + if (output.at_(r, c) == unmarked) + { + if (( r != 0 && output.at_(r-1, c) != unmarked) + || (c != 0 && output.at_(r, c-1) != unmarked) + || (c != input.ncols() - 1 && output.at_(r, c+1) != unmarked) + || (r != input.nrows() - 1 && output.at_(r+1, c) != unmarked)) + { + queue.push(max - input.at_(r, c), point2d(r, c)); + in_queue.at_(r, c) = true; + } + } + } + + /* Until the queue is empty, extract a point (R, C) from the + hierarchical queue, at the highest priority level, that is, + the lowest level. */ + while (! queue.is_empty()) + { + point2d p = queue.front(); + queue.pop(); + unsigned int r = p.row(); + unsigned int c = p.col(); + + // Last seen marker adjacent to (R, C). + marker adjacent_marker = unmarked; + // Has (R, C) a single adjacent marker? + bool single_adjacent_marker_p = true; + + const size_t nneighbs = 4; + point2d neighb[nneighbs] = { point2d(r-1, c), + point2d(r, c-1), + point2d(r, c+1), + point2d(r+1, c) }; + for (size_t n = 0; n < nneighbs; ++n) + if (output.domain().has(neighb[n]) && output(neighb[n]) != unmarked) + { + if (adjacent_marker == unmarked) + { + adjacent_marker = output(neighb[n]); + single_adjacent_marker_p = true; + } + else + if (adjacent_marker != output(neighb[n])) + { + single_adjacent_marker_p = false; + break; + } + } + /* If the neighborhood of (R, C) contains only points with the + same label, then (R, C) is marked with this label, and its + neighbors that are not yet marked are put into the + hierarchical queue. */ + if (single_adjacent_marker_p) + { + output.at_(r, c) = adjacent_marker; + for (size_t n = 0; n < nneighbs; ++n) + if (output.domain().has(neighb[n]) + && output(neighb[n]) == unmarked + && ! in_queue(neighb[n])) + { + queue.push(max - input(neighb[n]), neighb[n]); + in_queue(neighb[n]) = true; + } + } + } + + return output; + } +} + +/*--------------------------------------------------. +| Non-generic version, using up to three pointers. | +`--------------------------------------------------*/ + +namespace nongen_3ptr +{ + using namespace mln; + + typedef mln::image2d<mln::value::int_u8> image; + typedef mln::value::label<24> label; + typedef mln::image2d<label> labels; + + labels + watershed(const image& input, label& n_basins) + { + // FIXME: Simplify the body of this routine even more? + + typedef label marker; + const marker unmarked = 0; + + typedef mln_value_(image) V; + const V max = mln_max(V); + + + // Initialize the output with the markers (minima components). + typedef label L; + typedef labels O; + O output = labeling::regional_minima(input, c4(), n_basins); + + const size_t nneighbs = 4; + // Offsets corresponding to a 4-c neighborhood on OUTPUT. + ptrdiff_t nbh_offset[nneighbs] = + { &output.at_(-1, 0) - &output.at_(0, 0), + &output.at_( 0, -1) - &output.at_(0, 0), + &output.at_( 0, +1) - &output.at_(0, 0), + &output.at_(+1, 0) - &output.at_(0, 0) }; + + // Ordered queue. + typedef p_queue_fast<L*> Q; + p_priority<V, Q> queue; + + // In_queue structure to avoid processing sites several times. + mln_ch_value_(image, bool) in_queue; + initialize(in_queue, input); + data::fill(in_queue, false); + + // Offset between a the pixel located at the same position in + // OUTPUT and INPUT. + const ptrdiff_t input_offset = + (char*)(&input.at_(0, 0)) - (char*)(&output.at_(0, 0)); + // Offset between a the pixel located at the same position in + // OUTPUT and IN_QUEUE. + const ptrdiff_t in_queue_offset = + (char*)(&in_queue.at_(0, 0)) - (char*)(&output.at_(0, 0)); + + // Insert every neighbor of every marked area in a hierarchical + // queue, with a priority level corresponding to the grey level of + // this area. + for (unsigned int r = 0; r < output.nrows(); ++r) // Iterate on rows. + { + const V* pi = &input.at_(r, 0); + L* po = &output.at_(r, 0); + bool* pq = &in_queue.at_(r, 0); + for (; pi < &input.at_(r, 0) + input.ncols(); ++pi, ++po, ++pq) + { + if (*po == unmarked) + { + if ( + // (-1, 0) neighbor. + ( r != 0 + && *(po + nbh_offset[0]) != unmarked) + // (0, -1) neighbor. + || (po != &output.at_(r, 0) + && *(po + nbh_offset[1]) != unmarked) + // (0, +1) neighbor. + || (po != &output.at_(r, 0) + output.ncols() - 1 + && *(po + nbh_offset[2]) != unmarked) + // (+1, 0) neighbor. + || (r != output.nrows() - 1 + && *(po + nbh_offset[3]) != unmarked)) + { + queue.push(max - *pi, po); + *pq = true; + } + } + } + } + + /* Until the queue is empty, extract a point (R, C) from the + hierarchical queue, at the highest priority level, that is, + the lowest level. */ + while (! queue.is_empty()) + { + L* po = queue.front(); + queue.pop(); + + // Last seen marker adjacent to the current point. + marker adjacent_marker = unmarked; + // Has the current point a single adjacent marker? + bool single_adjacent_marker_p = true; + + for (size_t n = 0; n < nneighbs; ++n) + { + L* no = po + nbh_offset[n]; + // In the border, OUTPUT is unmarked so N is ignored. + if (*no != unmarked) + { + if (adjacent_marker == unmarked) + { + adjacent_marker = *no; + single_adjacent_marker_p = true; + } + else + if (adjacent_marker != *no) + { + single_adjacent_marker_p = false; + break; + } + } + } + /* If the neighborhood of the current point contains only + points with the same label, then this point is marked with + this label, and its neighbors that are not yet marked are + put into the hierarchical queue. */ + if (single_adjacent_marker_p) + { + *po = adjacent_marker; + for (size_t n = 0; n < nneighbs; ++n) + { + L* no = po + nbh_offset[n]; + /* Deduce the pointers NI and NQ corresponding to NO + in IN_QUEUE and INPUT using the memory offsets + bewteen each of these images and OUPUT. This works + because the data of qll these images are arranged + similarly. + + The ugly casts to and back from char* are required, + since the ptrdiff_t values stored in INPUT_OFFSET + and IN_QUEUE_OFFSET are expressed as differences + between chars. */ + V* ni = (V*)((char*) no + input_offset); + // FIXME: Remove when OK. + assert(&input.at_(0, 0) < ni); + assert(ni < &input.element(input.nelements() - 1)); + bool* nq = (bool*)((char*) no + in_queue_offset); + // FIXME: Remove when OK. + assert(&in_queue.at_(0, 0) < nq); + assert(nq < &in_queue.element(in_queue.nelements() - 1)); + if (*no == unmarked && ! *nq) + { + queue.push(max - *ni, no); + *nq = true; + } + } + } + } + + return output; + } +} + + +/*----------------------------------------------. +| Non-generic version, using a single pointer. | +`----------------------------------------------*/ + +namespace nongen_2ptr +{ + // FIXME: To do. +} + + +/*------------------. +| Generic version. | +`------------------*/ + +namespace gen +{ + using namespace mln; + + template <typename L, typename I, typename N> + mln_ch_value(I, L) + watershed(const I& input_, const Neighborhood<N>& nbh_, L& n_basins) + { + const I input = exact(input_); + const N nbh = exact(nbh_); + + typedef L marker; + const marker unmarked = literal::zero; + + typedef mln_value(I) V; + const V max = mln_max(V); + + // Initialize the output with the markers (minima components). + typedef mln_ch_value(I, L) O; + O output = labeling::regional_minima(input, nbh, n_basins); + + typedef mln_psite(I) psite; + + // Ordered queue. + typedef p_queue_fast<psite> Q; + p_priority<V, Q> queue; + + // In_queue structure to avoid processing sites several times. + mln_ch_value(I, bool) in_queue; + initialize(in_queue, input); + data::fill(in_queue, false); + + // Insert every neighbor P of every marked area in a + // hierarchical queue, with a priority level corresponding to + // the grey level input(P). + mln_piter(I) p(output.domain()); + mln_niter(N) n(nbh, p); + for_all(p) + if (output(p) == unmarked) + for_all(n) + if (output.domain().has(n) && output(n) != unmarked) + { + queue.push(max - input(p), p); + in_queue(p) = true; + break; + } + + /* Until the queue is empty, extract a psite P from the + hierarchical queue, at the highest priority level, that is, + the lowest level. */ + while (! queue.is_empty()) + { + psite p = queue.front(); + queue.pop(); + + // Last seen marker adjacent to P. + marker adjacent_marker = unmarked; + // Has P a single adjacent marker? + bool single_adjacent_marker_p = true; + mln_niter(N) n(nbh, p); + for_all(n) + if (output.domain().has(n) && output(n) != unmarked) + { + if (adjacent_marker == unmarked) + { + adjacent_marker = output(n); + single_adjacent_marker_p = true; + } + else + if (adjacent_marker != output(n)) + { + single_adjacent_marker_p = false; + break; + } + } + /* If the neighborhood of P contains only psites with the + same label, then P is marked with this label, and its + neighbors that are not yet marked are put into the + hierarchical queue. */ + if (single_adjacent_marker_p) + { + output(p) = adjacent_marker; + for_all(n) + if (output.domain().has(n) && output(n) == unmarked + && ! in_queue(n)) + { + queue.push(max - input(n), n); + in_queue(n) = true; + } + } + } + + return output; + } +} + +/*-----------------------------------------------. +| Fast version (with one pixter), less generic. | +`-----------------------------------------------*/ + +namespace fast +{ + using namespace mln; + + template <typename I, typename N, typename L> + mln_ch_value(I, L) + watershed(const Image<I>& input_, const Neighborhood<N>& nbh_, + L& n_basins) + { + const I input = exact(input_); + const N nbh = exact(nbh_); + + typedef L marker; + const marker unmarked = literal::zero; + + typedef mln_value(I) V; + const V max = mln_max(V); + + extension::adjust_fill(input, nbh, max); + + // Initialize the output with the markers (minima components). + typedef mln_ch_value(I, L) O; + O output = labeling::regional_minima(input, nbh, n_basins); + extension::fill(output, unmarked); + + // Ordered queue. + typedef p_queue_fast<unsigned> Q; + p_priority<V, Q> queue; + + // In_queue structure to avoid processing sites several times. + mln_ch_value(I, bool) in_queue; + initialize(in_queue, input); + data::fill(in_queue, false); + extension::fill(in_queue, true); + + // Insert every neighbor P of every marked area in a + // hierarchical queue, with a priority level corresponding to + // the grey level input(P). + mln_pixter(const O) p_out(output); + mln_nixter(const O, N) n_out(p_out, nbh); + for_all(p_out) + if (p_out.val() == unmarked) + for_all(n_out) + if (n_out.val() != unmarked) + { + unsigned po = p_out.offset(); + queue.push(max - input.element(po), po); + in_queue.element(po) = true; + break; + } + + /* Until the queue is empty, extract a psite P from the + hierarchical queue, at the highest priority level, that is, + the lowest level. */ + util::array<int> dp = offsets_wrt(input, nbh); + const unsigned n_nbhs = dp.nelements(); + while (! queue.is_empty()) + { + unsigned p = queue.front(); + queue.pop(); + + // Last seen marker adjacent to P. + marker adjacent_marker = unmarked; + // Has P a single adjacent marker? + bool single_adjacent_marker_p = true; + for (unsigned i = 0; i < n_nbhs; ++i) + { + unsigned n = p + dp[i]; + // In the border, OUTPUT is unmarked so N is ignored. + if (output.element(n) != unmarked) + { + if (adjacent_marker == unmarked) + { + adjacent_marker = output.element(n); + single_adjacent_marker_p = true; + } + else + if (adjacent_marker != output.element(n)) + { + single_adjacent_marker_p = false; + break; + } + } + } + /* If the neighborhood of P contains only psites with the + same label, then P is marked with this label, and its + neighbors that are not yet marked are put into the + hierarchical queue. */ + if (single_adjacent_marker_p) + { + output.element(p) = adjacent_marker; + for (unsigned i = 0; i < n_nbhs; ++i) + { + unsigned n = p + dp[i]; + if (output.element(n) == unmarked + // In the border, IN_QUEUE is true so N is ignored. + && ! in_queue.element(n)) + { + queue.push(max - input.element(n), n); + in_queue.element(n) = true; + } + } + } + } + return output; + } +} + + +/*------------------------------------------------------------. +| Fast version (with one pixter) with a static neighborhood. | +`------------------------------------------------------------*/ + +namespace fast_static +{ + using namespace mln; + + template <typename I, typename N, typename L> + mln_ch_value(I, L) + watershed(const Image<I>& input_, const Neighborhood<N>& nbh_, + L& n_basins) + { + // `fast_static' has the same implementation as `fast'. + return ::fast::watershed(input_, nbh_, n_basins); + } +} + +#endif // ! APPS_BENCH_WATERSHED_HH -- 1.7.2.5