[Olena-patches] proto-1.0 107: Add hybrid and sequential versions of reconstruction

Index: ChangeLog from Damien Thivolle &lt;damien(a)lrde.epita.fr&gt; * oln/convert/ng_to_se.hh: Convert a neighborhood to a structuring element. * oln/core/value_box.hh: Fix operator==. * oln/core/properties.hh: Add window property. * oln/core/abstract/struct_elt.hh: Add operator-. * oln/core/abstract/witer.hh (nth): Add method. * oln/core/abstract/dpoint.hh (nth): Add method. * oln/core/abstract/neighborhood.hh: Add a window property. * oln/core/2d/dpoint2d.hh (nth): Add method. * oln/core/2d/neighborhood2d.hh: Add a window property. * oln/core/2d/fwd_witer2d.hh: Replace unknown type in constructor. * oln/utils/clone.hh: Use concrete_type. * oln/makefile.src: Add new files. * oln/morpho/reconstruction.hh: New. Sequential and hybrid versions. * oln/morpho/splitse.hh: New. Split a structuring element. * oln/morpho/stat.hh: Typo fix. * oln/morpho/erosion.hh: Fix erroneus function name . * oln/io/write_image_2d_pnm.hh: Invert binary values (0=>1, 1=>0). * oln/io/read_image_2d_pnm.hh: Invert binary values. convert/ng_to_se.hh | 88 +++++ core/2d/dpoint2d.hh | 21 + core/2d/fwd_witer2d.hh | 4 core/2d/neighborhood2d.hh | 18 - core/abstract/dpoint.hh | 5 core/abstract/neighborhood.hh | 22 - core/abstract/struct_elt.hh | 38 +- core/abstract/witer.hh | 13 core/properties.hh | 1 core/value_box.hh | 18 - io/read_image_2d_pnm.hh | 4 io/write_image_2d_pnm.hh | 2 makefile.src | 4 morpho/erosion.hh | 6 morpho/reconstruction.hh | 646 ++++++++++++++++++++++++++++++++++++++++++ morpho/splitse.hh | 208 +++++++++++++ morpho/stat.hh | 2 utils/clone.hh | 7 18 files changed, 1062 insertions, 45 deletions Index: oln/convert/ng_to_se.hh --- oln/convert/ng_to_se.hh (revision 0) +++ oln/convert/ng_to_se.hh (revision 0) @@ -0,0 +1,88 @@ +// Copyright (C) 2002, 2004, 2005 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, 59 Temple Place - Suite 330, 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 OLENA_CONVERT_NG_TO_SE_HH +# define OLENA_CONVERT_NG_TO_SE_HH + +# include <oln/basics.hh> + +// because of the internal function in this file +//# include <oln/basics1d.hh> +# include <oln/basics2d.hh> +//# include <oln/basics3d.hh> +//# include <oln/core/1d/neighborhood1d.hh> +# include <oln/core/2d/neighborhood2d.hh> +//# include <oln/core/3d/neighborhood3d.hh> + +namespace oln { + namespace convert { + /*! Convert a neighborhood to a window. + ** + ** \see ng_to_cse + */ + template<class N> + oln_type_of(N, window) + ng_to_se(const oln::abstract::neighborhood<N>& ng) + { + oln_type_of(N, window) output; + + for (unsigned i = 0; i < ng.card(); i++) + output.add(ng[i]); + return output; + } + + void dpoint_zero(dpoint2d& dp) + { + dp.row() = 0; + dp.col() = 0; + } + + /*! Convert a neighborhood to a window and add the center. + ** + ** \see ng_to_cs + */ + template<class N> + oln_type_of(N, window) + ng_to_cse(const oln::abstract::neighborhood<N>& ng) + { + oln_type_of(N, window) output; + + for (unsigned i = 0; i < ng.card(); i++) + output.add(ng[i]); + + oln_type_of(N, dpoint) zero; + dpoint_zero(zero); + output.add(zero); + return output; + } + + + } // convert +} // oln + + +#endif // OLENA_CONVERT_NG_TO_SE_HH Index: oln/core/value_box.hh --- oln/core/value_box.hh (revision 106) +++ oln/core/value_box.hh (working copy) @@ -81,11 +81,17 @@ /// op== - bool operator==(const value_box<const I>& rhs) const + bool operator==(const value_box<I>& rhs) const { return this->value() == rhs.value(); } + template <typename II> + bool operator==(const value_box<II>& rhs) const + { + return this->value() < rhs.value(); + } + template <typename V> bool operator==(const V& rhs) const { @@ -113,7 +119,7 @@ template <typename II> bool operator<(const value_box<II>& rhs) const { - return this->value() < value.value(); + return this->value() < rhs.value(); } /*! \brief specialized version for value_box<I> @@ -122,7 +128,7 @@ */ bool operator<(const value_box<I>& rhs) const { - return this->value() < value.value(); + return this->value() < rhs.value(); } /*! \brief op= @@ -277,6 +283,12 @@ return this->value() == rhs.value(); } + template <typename II> + bool operator==(const value_box<const II>& rhs) const + { + return this->value() < rhs.value(); + } + template <typename V> bool operator==(const V& rhs) const { Index: oln/core/properties.hh --- oln/core/properties.hh (revision 106) +++ oln/core/properties.hh (working copy) @@ -68,6 +68,7 @@ struct delegated_type; struct size_type; struct se_type; + struct window_type; struct image_constness_type; struct image_dimension_type; Index: oln/core/abstract/struct_elt.hh --- oln/core/abstract/struct_elt.hh (revision 106) +++ oln/core/abstract/struct_elt.hh (working copy) @@ -47,7 +47,7 @@ // category template <typename E> - struct set_category< abstract::struct_elt<E> > + struct set_category< abstract::struct_elt<E> > { typedef category::struct_elt ret; }; @@ -71,7 +71,7 @@ << " fwd_witer_type = " << typeid(fwd_witer_type).name() << "," << std::endl << " dpoint_type = " << typeid(dpoint_type).name() << " }" << std::endl; } - + }; mlc_register_prop(category::struct_elt, dpoint_type); @@ -137,7 +137,7 @@ get_delta() const { return this->exact().impl_get_delta(); - } + } coord_t delta_update(const dpoint_type& dp) @@ -145,14 +145,36 @@ return this->exact().impl_delta_update(dp); } + exact_type + operator-() const + { + exact_type se(this->exact()); + + se.sym(); + return se; + } + + void + sym() + { + this->exact().impl_sym(); + } + protected: + void + impl_sym() + { + for (unsigned i = 0; i < this->card(); ++i) + dp_[i] = - dp_[i]; + } + bool impl_has(const dpoint_type& dp) const { return std::find(dp_.begin(), dp_.end(), dp) != dp_.end(); - } - + } + exact_type& impl_add(const dpoint_type& dp) { @@ -172,14 +194,14 @@ impl_card() const { return dp_.size(); - } + } const dpoint_type impl_at(unsigned i) const { precondition(i < this->card()); return dp_[i]; - } + } struct_elt() : dp_(), delta_(0) {}; @@ -208,6 +230,6 @@ o << "]"; return o; } - + #endif // ! OLENA_CORE_STRUCT_ELT_HH Index: oln/core/abstract/witer.hh --- oln/core/abstract/witer.hh (revision 106) +++ oln/core/abstract/witer.hh (working copy) @@ -64,7 +64,7 @@ } }; - mlc_register_prop(category::witer, se_type); + mlc_register_prop(category::witer, se_type); namespace abstract { @@ -104,7 +104,12 @@ this->exact().impl_invalidate(); postcondition(! this->is_valid()); } - + + coord_t nth(unsigned i) + { + return this->se_[this->pos_].nth(i); + } + protected: void impl_start() @@ -127,8 +132,8 @@ pos_ = se_.card(); } - witer(const se_type& se) - : se_(se), pos_(0) + witer(const se_type& se) + : se_(se), pos_(0) {} const se_type& se_; Index: oln/core/abstract/dpoint.hh --- oln/core/abstract/dpoint.hh (revision 106) +++ oln/core/abstract/dpoint.hh (working copy) @@ -65,6 +65,11 @@ return not this->operator==(rhs); } + coord_t nth(unsigned i) const + { + return this->exact().impl_nth(i); + } + protected: dpoint() {} Index: oln/core/abstract/neighborhood.hh --- oln/core/abstract/neighborhood.hh (revision 106) +++ oln/core/abstract/neighborhood.hh (working copy) @@ -47,7 +47,7 @@ // category template <typename E> - struct set_category< abstract::neighborhood<E> > + struct set_category< abstract::neighborhood<E> > { typedef category::neighborhood ret; }; @@ -61,6 +61,7 @@ mlc_decl_prop(category::neighborhood, dpoint_type); mlc_decl_prop(category::neighborhood, size_type); + mlc_decl_prop(category::neighborhood, window_type); static void echo(std::ostream& ostr) { @@ -69,11 +70,12 @@ << " size_type = " << typeid(size_type).name() << "," << std::endl << " dpoint_type = " << typeid(dpoint_type).name() << " }" << std::endl; } - + }; mlc_register_prop(category::neighborhood, dpoint_type); mlc_register_prop(category::neighborhood, size_type); + mlc_register_prop(category::neighborhood, window_type); namespace abstract { @@ -134,7 +136,7 @@ get_delta() const { return this->exact().impl_get_delta(); - } + } coord_t delta_update(const dpoint_type& dp) @@ -148,8 +150,8 @@ impl_has(const dpoint_type& dp) const { return std::find(dp_.begin(), dp_.end(), dp) != dp_.end(); - } - + } + exact_type& impl_add(const dpoint_type& dp) { @@ -161,25 +163,25 @@ this->delta_update(dp); return this->exact(); } - + coord_t impl_get_delta() const { return delta_; } - + unsigned impl_card() const { return dp_.size(); - } + } const dpoint_type impl_at(unsigned i) const { precondition(i < this->card()); return dp_[i]; - } + } neighborhood() : dp_(), delta_(0) {}; @@ -208,6 +210,6 @@ o << "]"; return o; } - + #endif // ! OLENA_CORE_NEIGHBORHOOD_HH Index: oln/core/2d/dpoint2d.hh --- oln/core/2d/dpoint2d.hh (revision 106) +++ oln/core/2d/dpoint2d.hh (working copy) @@ -94,10 +94,31 @@ const coord_t row() const { return row_; } const coord_t col() const { return col_; } + //FIXME : name it impl_nth when dpoint2d derives from abstract::dpoint + const coord_t nth(unsigned i) const + { + assert(i < 2); + if (i == 0) + return row_; + else + return col_; + } + + coord_t& row() { return row_; } coord_t& col() { return col_; } + //FIXME : name it impl_nth when dpoint2d derives from abstract::dpoint + coord_t& nth(unsigned i) + { + assert(i < 2); + if (i == 0) + return row_; + else + return col_; + } + protected: coord_t row_, col_; }; Index: oln/core/2d/neighborhood2d.hh --- oln/core/2d/neighborhood2d.hh (revision 106) +++ oln/core/2d/neighborhood2d.hh (working copy) @@ -30,6 +30,7 @@ # include <oln/core/2d/dpoint2d.hh> # include <oln/core/2d/size2d.hh> +# include <oln/core/2d/window2d.hh> # include <oln/core/coord.hh> # include <oln/core/abstract/neighborhood.hh> @@ -39,7 +40,7 @@ // category template <> - struct set_category< neighborhood2d > + struct set_category< neighborhood2d > { typedef category::neighborhood ret; }; // super_type @@ -50,18 +51,19 @@ }; template <> - struct set_props< category::neighborhood, neighborhood2d > + struct set_props< category::neighborhood, neighborhood2d > : public props_of<category::neighborhood> { typedef dpoint2d dpoint_type; typedef size2d size_type; + typedef window2d window_type; }; class neighborhood2d : public abstract::neighborhood<neighborhood2d> { public: - + typedef abstract::neighborhood<neighborhood2d> super_type; /*! @@ -86,13 +88,13 @@ { for (unsigned i = 0; i < 2 * n; i += 2) this->add(dpoint_type(crd[i], crd[i+1])); - } + } neighborhood2d& add(const dpoint_type& dp) { return this->exact().impl_add(dp); - } + } neighborhood2d& add(coord_t row, coord_t col) @@ -100,8 +102,8 @@ dpoint_type dp(row, col); return add(dp); } - + /// Return the name of the type. static std::string name() @@ -116,7 +118,7 @@ delta_(abs(dp.col())); return delta_; } - + }; inline const neighborhood2d& @@ -134,7 +136,7 @@ static const neighborhood2d neighb(4, crd); return neighb; } - + } // end of oln #endif // OLENA_CORE_NEIGHBORHOOD2D_HH Index: oln/core/2d/fwd_witer2d.hh --- oln/core/2d/fwd_witer2d.hh (revision 106) +++ oln/core/2d/fwd_witer2d.hh (working copy) @@ -48,13 +48,13 @@ { typedef window2d se_type; }; - + struct fwd_witer2d : public abstract::witer< fwd_witer2d > { typedef abstract::witer<fwd_witer2d> super_type; - fwd_witer2d(const se_type& se) : + fwd_witer2d(const window2d& se) : super_type(se) { this->exact_ptr = this; Index: oln/utils/clone.hh --- oln/utils/clone.hh (revision 106) +++ oln/utils/clone.hh (working copy) @@ -60,10 +60,11 @@ namespace impl { template <typename I> - struct clone_type : public abstract::image_unary_operator<I, I, clone_type<I> > + struct clone_type + : public abstract::image_unary_operator<oln_type_of(I, concrete), I, clone_type<I> > // FIXME: use concrete_type; Cf. erosion.hh { - typedef abstract::image_unary_operator<I, I, clone_type<I> > super_type; + typedef abstract::image_unary_operator<oln_type_of(I, concrete), I, clone_type<I> > super_type; clone_type(const abstract::image<I>& input) : super_type(input) @@ -72,7 +73,7 @@ void impl_run() { - I tmp(this->input.size()); // FIXME: trick + oln_type_of(I, concrete) tmp(this->input.size()); // FIXME: trick this->output = tmp; oln_type_of(I, fwd_piter) p(this->input.size()); Index: oln/makefile.src --- oln/makefile.src (revision 106) +++ oln/makefile.src (working copy) @@ -13,6 +13,7 @@ basics3d.hh \ config/pconf.hh \ config/system.hh \ + convert/ng_to_se.hh \ convert/value_to_point.hh \ core/1d/array1d.hh \ core/1d/dpoint1d.hh \ @@ -87,6 +88,9 @@ io/write_image_2d_pnm.hh \ level/compare.hh \ level/fill.hh \ + morpho/dilation.hh \ morpho/erosion.hh \ + morpho/reconstruction.hh \ + morpho/splitse.hh \ morpho/stat.hh \ utils/clone.hh Index: oln/morpho/reconstruction.hh --- oln/morpho/reconstruction.hh (revision 0) +++ oln/morpho/reconstruction.hh (revision 0) @@ -0,0 +1,646 @@ +// Copyright (C) 2001, 2002, 2004, 2005 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, 59 Temple Place - Suite 330, Boston, +// MA 02111-1307, USA. +// +// As a special exception, you may use this filek 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 OLENA_MORPHO_RECONSTRUCTION_HH +# define OLENA_MORPHO_RECONSTRUCTION_HH + +# include <queue> + +# include <mlc/contract.hh> +# include <oln/convert/ng_to_se.hh> +# include <oln/morpho/splitse.hh> +# include <oln/level/compare.hh> +# include <oln/core/abstract/neighborhood.hh> +# include <oln/core/properties.hh> +# include <oln/core/abstract/image_operator.hh> +# include <oln/utils/clone.hh> + +// FIXME: ADD TESTS !!!! + +namespace oln { + + namespace morpho { + + // fwd decl + template <typename I, typename E> struct reconstruction_ret; + + } + + // category + template <typename I, typename E> + struct set_category< morpho::reconstruction_ret<I,E> > { typedef category::image ret; }; + + // super_type + template <typename I, typename E> + struct set_super_type< morpho::reconstruction_ret<I,E> > + { + typedef abstract::image_binary_operator<I, I, I, morpho::reconstruction_ret<I, E> > ret; + // FIXME: see below + }; + + + namespace morpho { + + template <typename I, typename N> + struct reconstruction_ret : public abstract::image_binary_operator<I, I, I, reconstruction_ret<I, N> > + // FIXME: abstract::image_binary_operator<oln_type_of(I, concrete), ... + { + typedef abstract::image_binary_operator<I, I, I, reconstruction_ret<I, N> > super_type; + + const N ng; + + reconstruction_ret(const abstract::non_vectorial_image<I>& input1, + const abstract::non_vectorial_image<I>& input2, + const abstract::neighborhood<N>& ng) : + super_type(input1.exact(), input2.exact()), + ng(ng.exact()) + {} + + }; + + namespace sequential { + /*! + ** \brief Perform a geodesic reconstruction dilation. + ** + ** Compute the reconstruction by dilation of marker with respect + ** to the mask image using se as structuring element. Soille + ** p.160. The algorithm used is the one defined as sequential in + ** Vincent(1993), Morphological grayscale reconstruction in + ** image analysis: applications and efficient algorithms, itip, + ** 2(2), 176--201. + ** + ** \pre Mask must be greater or equal than marker. + ** + ** \param I1 Exact type of image marker. + ** \param I2 Exact type of image mask. + ** \param N Exact type of neighborhood. + ** + ** \arg marker Image to work on. + ** \arg mask Image used for geodesic dilation. + ** \arg Ng Neighborhood to use. + ** + ** \code + ** #include <oln/basics2d.hh> + ** #include <oln/morpho/opening.hh> + ** #include <oln/morpho/reconstruction.hh> + ** #include <oln/level/compare.hh> + ** #include <ntg/all.hh> + ** int main() + ** { + ** typedef oln::image2d<ntg::int_u8> im_type; + ** + ** im_type im1(oln::load(IMG_IN "lena128.pgm")); + ** im_type im2 = oln::morpho::opening(im1, oln::win_c4p()); + ** + ** oln::save(oln::morpho::sequential::geodesic_reconstruction_dilation(im2, + ** im1, + ** oln::neighb_c4()), + ** IMG_OUT "oln_morpho_sequential_geodesic_reconstruction_dilation.pbm"); + ** return 0; + ** } + ** \endcode + ** + ** \image html lena128_pgm.png + ** \image latex lena128_pgm.png + ** => + ** \image html oln_morpho_sequential_geodesic_reconstruction_dilation.png + ** \image latex oln_morpho_sequential_geodesic_reconstruction_dilation.png + ** + */ + + namespace impl { + + template <typename I, typename N> + struct reconstruction_dilation_ret : public reconstruction_ret<I, N> + { + typedef reconstruction_ret<I, N> super_type; + + reconstruction_dilation_ret(const abstract::non_vectorial_image<I>& input1, //marker + const abstract::non_vectorial_image<I>& input2, //mask + const abstract::neighborhood<N>& ng) + + : super_type(input1, input2, ng) + {} + + void impl_run() + { + mlc::is_true<mlc::type::eq<oln_type_of(I, size), + oln_type_of(N, size)>::ret>::ensure(); + precondition(this->input1.size() == this->input2.size()); + precondition(level::is_greater_or_equal(this->input2, this->input1)); + + // Conversion of neighborhood into a SE. + oln_type_of(N, window) se_plus = get_plus_se_p(convert::ng_to_cse(this->ng)); + oln_type_of(N, window) se_minus = get_minus_se_p(convert::ng_to_cse(this->ng)); + + I output; + output = utils::clone(this->input1); + bool non_stability = true; + oln_type_of(I, fwd_piter) fwd_p(output.size()); + oln_type_of(I, fwd_piter) bkd_p(output.size()); + while (non_stability) + { + I work; + work = utils::clone(output); + for_all (fwd_p) + work[fwd_p] = ntg::min(morpho::max(work, fwd_p, se_plus), this->input2[fwd_p].value()); + for_all (bkd_p) + work[bkd_p] = ntg::min(morpho::max(work, bkd_p, se_minus), this->input2[bkd_p].value()); + non_stability = !(level::is_equal(work, output)); + output = work; + } + this->output = output; + } + }; + + } + + template<class I, class N> + reconstruction_ret<I, N> + geodesic_reconstruction_dilation(const abstract::non_vectorial_image<I> & marker, + const abstract::non_vectorial_image<I> & mask, + const abstract::neighborhood<N>& ng) + { + impl::reconstruction_dilation_ret<I, N> tmp(marker, mask, ng); + tmp.run(); + return tmp; + } + }// sequential + + + namespace hybrid { + + namespace internal { + + /*! + ** \brief Check if it exists initialization for dilation. + ** + ** \arg p Point to consider. + ** \arg marker Image to work on. + ** \arg mask Image used as mask. + ** \arg Ng Neighborhood to use. + */ + template<class P, class I1, class I2, class E> inline + static bool + exist_init_dilation(const abstract::point<P>& p, + const abstract::non_vectorial_image<I1>& marker, + const abstract::non_vectorial_image<I2>& mask, + const abstract::struct_elt<E>& se) + { + oln_type_of(E, fwd_witer) dp(se.exact()); + for_all (dp) + { + P q = (oln_type_of(E, dpoint))dp + p.exact(); + if (marker.hold(q) && (marker[q] < marker[p.exact()]) && + (marker[q] < mask[q])) + return true; + } + return false; + } + + } //internal + + /*! + ** \brief Perform a geodesic reconstruction dilation. + ** + ** Compute the reconstruction by dilation of marker with + ** respect to the mask image using se as structuring + ** element. Soille p.160. The algorithm used is the one defined + ** as hybrid in Vincent(1993), Morphological grayscale + ** reconstruction in image analysis: applications and efficient + ** algorithms, itip, 2(2), 176--201. + ** + ** \pre Mask must be greater or equal than marker. + ** + ** \param I1 Exact type of image marker. + ** \param I2 Exact type of image mask. + ** \param N Exact type of neighborhood. + ** + ** \arg marker Image to work on. + ** \arg mask Image used for geodesic dilation. + ** \arg Ng Neighborhood to use. + ** + ** \code + ** #include <oln/basics2d.hh> + ** #include <oln/morpho/opening.hh> + ** #include <oln/morpho/reconstruction.hh> + ** #include <oln/level/compare.hh> + ** #include <ntg/all.hh> + ** int main() + ** { + ** typedef oln::image2d<ntg::int_u8> im_type; + ** + ** im_type im1(oln::load(IMG_IN "lena128.pgm")); + ** im_type im2 = oln::morpho::opening(im1, oln::win_c4p()); + ** + ** oln::save(oln::morpho::hybrid::geodesic_reconstruction_dilation(im2, + ** im1, + ** oln::neighb_c4()), + ** IMG_OUT "oln_morpho_hybrid_geodesic_reconstruction_dilation.pbm"); + ** return 0; + ** } + ** \endcode + ** + ** \image html lena128_pgm.png + ** \image latex lena128_pgm.png + ** => + ** \image html oln_morpho_hybrid_geodesic_reconstruction_dilation.png + ** \image latex oln_morpho_hybrid_geodesic_reconstruction_dilation.png + ** + */ + + namespace impl { + + template <typename I, typename N> + struct reconstruction_dilation_ret : public reconstruction_ret<I, N> + { + typedef reconstruction_ret<I, N> super_type; + + reconstruction_dilation_ret(const abstract::non_vectorial_image<I>& input1, //marker + const abstract::non_vectorial_image<I>& input2, //mask + const abstract::neighborhood<N>& ng) + + : super_type(input1, input2, ng) + {} + + void impl_run() + { + mlc::is_true<mlc::type::eq<oln_type_of(I, size), + oln_type_of(N, size)>::ret>::ensure(); + precondition(this->input1.size() == this->input2.size()); + precondition(level::is_greater_or_equal(this->input2, this->input1)); + + oln_type_of(I, concrete) output; + output = utils::clone(this->input1); + { + oln_type_of(N, window) se_plus = get_plus_se_p(convert::ng_to_cse(this->ng)); + oln_type_of(N, window) se_minus = get_minus_se_p(convert::ng_to_cse(this->ng)); + oln_type_of(I, fwd_piter) fwd_p(output.size()); + oln_type_of(I, fwd_piter) bkd_p(output.size()); + + for_all (fwd_p) + output[fwd_p] = ntg::min(morpho::max(output, fwd_p, se_plus), + this->input2[fwd_p].value()); + + std::queue<oln_type_of(I, point) > fifo; + for_all (bkd_p) + { + output[bkd_p] = ntg::min(morpho::max(output, bkd_p, se_minus), + this->input2[bkd_p].value()); + if (internal::exist_init_dilation((oln_type_of(I, point))bkd_p, output, + this->input2, se_minus)) + fifo.push(bkd_p); + } + // Propagation Step + while (!fifo.empty()) + { + oln_type_of(I, point) p = fifo.front(); + fifo.pop(); + typedef oln_type_of(N, window) window_type; + window_type w = convert::ng_to_se(this->ng); + oln_type_of(window_type, fwd_witer) dp(w); + + for_all (dp) + { + oln_type_of(I, point) q = (oln_type_of(window_type, dpoint))dp + p; + if (output.hold(q)) + { + if ((output[q] < output[p]) && + (this->input2[q] != output[q])) + { + output[q] = ntg::min(output[p].value(), + this->input2[q].value()); + fifo.push(q); + } + } + } + } + } + this->output = output; + } + }; + + } + + template<class I, class N> + reconstruction_ret<I, N> + geodesic_reconstruction_dilation(const abstract::non_vectorial_image<I> & marker, + const abstract::non_vectorial_image<I> & mask, + const abstract::neighborhood<N>& ng) + { + impl::reconstruction_dilation_ret<I, N> tmp(marker, mask, ng); + tmp.run(); + return tmp; + } + }// hybrid + + + namespace sequential { + + /*! + ** \brief Perform a geodesic reconstruction erosion. + ** + ** Compute the reconstruction by erosion of marker with respect + ** to the mask image using se as structuring element. Soille + ** p.160. The algorithm used is the one defined as sequential + ** in Vincent(1993), Morphological grayscale reconstruction in + ** image analysis: applications and efficient algorithms, itip, + ** 2(2), 176--201. + ** + ** \pre Marker must be greater or equal than mask. + ** + ** \param I1 Exact type of image marker. + ** \param I2 Exact type of image mask. + ** \param N Exact type of neighborhood. + ** + ** \arg marker Image to work on. + ** \arg mask Image used for geodesic erosion. + ** \arg Ng Neighborhood to use. + ** + ** \code + ** #include <oln/basics2d.hh> + ** #include <oln/morpho/opening.hh> + ** #include <oln/morpho/reconstruction.hh> + ** #include <oln/level/compare.hh> + ** #include <ntg/all.hh> + ** int main() + ** { + ** typedef oln::image2d<ntg::int_u8> im_type; + ** + ** im_type im1(oln::load(IMG_IN "lena128.pgm")); + ** im_type im2 = oln::morpho::opening(im1, oln::win_c4p()); + ** + ** oln::save(oln::morpho::sequential::geodesic_reconstruction_erosion(im1, + ** im2, + ** oln::neighb_c4()), + ** IMG_OUT "oln_morpho_sequential_geodesic_reconstruction_erosion.pbm"); + ** return 0; + ** } + ** \endcode + ** + ** \image html lena128_pgm.png + ** \image latex lena128_pgm.png + ** => + ** \image html oln_morpho_sequential_geodesic_reconstruction_erosion.png + ** \image latex oln_morpho_sequential_geodesic_reconstruction_erosion.png + ** + */ + + + namespace impl { + + template <typename I, typename N> + struct reconstruction_erosion_ret : public reconstruction_ret<I, N> + { + typedef reconstruction_ret<I, N> super_type; + + reconstruction_erosion_ret(const abstract::non_vectorial_image<I>& input1, //marker + const abstract::non_vectorial_image<I>& input2, //mask + const abstract::neighborhood<N>& ng) + + : super_type(input1, input2, ng) + {} + + void impl_run() + { + mlc::is_true<mlc::type::eq<oln_type_of(I, size), + oln_type_of(N, size)>::ret>::ensure(); + precondition(this->input1.size() == this->input2.size()); + precondition(level::is_greater_or_equal(this->input2, this->input1)); + + // Conversion of neighborhood into a SE. + oln_type_of(N, window) se_plus = get_plus_se_p(convert::ng_to_cse(this->ng)); + oln_type_of(N, window) se_minus = get_minus_se_p(convert::ng_to_cse(this->ng)); + + I output; + output = utils::clone(this->input1); + bool non_stability = true; + oln_type_of(I, fwd_piter) fwd_p(output.size()); + oln_type_of(I, fwd_piter) bkd_p(output.size()); + while (non_stability) + { + I work; + work = utils::clone(output); + for_all (fwd_p) + work[fwd_p] = ntg::max(morpho::min(work, fwd_p, se_plus), + this->input2[fwd_p].value()); + for_all (bkd_p) + work[bkd_p] = ntg::max(morpho::min(work, bkd_p, se_minus), + this->input2[bkd_p].value()); + non_stability = !(level::is_equal(work, output)); + output = work; + } + this->output = output; + } + }; + + } + + template<class I, class N> + reconstruction_ret<I, N> + geodesic_reconstruction_erosion(const abstract::non_vectorial_image<I> & marker, + const abstract::non_vectorial_image<I> & mask, + const abstract::neighborhood<N>& ng) + { + impl::reconstruction_erosion_ret<I, N> tmp(marker, mask, ng); + tmp.run(); + return tmp; + } + + } // sequential + + namespace hybrid { + + namespace internal { + + /*! + ** \brief Check if it exists initialization for erosion. + ** + ** \arg p Point to consider. + ** \arg marker Image to work on. + ** \arg mask Image used as mask. + ** \arg Ng Neighborhood to use. + */ + template<class P, class I1, class I2, class E> inline + static bool + exist_init_erosion(const abstract::point<P>& p, + const abstract::non_vectorial_image<I1>& marker, + const abstract::non_vectorial_image<I2>& mask, + const abstract::struct_elt<E>& se) + { + oln_type_of(E, fwd_witer) dp(se.exact()); + for_all (dp) + { + P q = (oln_type_of(E, dpoint))dp + p.exact(); + if (marker.hold(q) && (marker[q] > marker[p.exact()]) && + (marker[q] > mask[q])) + return true; + } + return false; + } + + } //internal + + /*! + ** \brief Perform a geodesic reconstruction erosion. + ** + ** Compute the reconstruction by erosion of marker with respect + ** to the mask mask image using se as structuring + ** element. Soille p.160. The algorithm used is the one defined + ** as hybrid in Vincent(1993), Morphological grayscale + ** reconstruction in image analysis: applications and efficient + ** algorithms, itip, 2(2), 176--201. + ** + ** \pre Marker must be greater or equal than mask. + ** + ** \param I1 Exact type of image marker. + ** \param I2 Exact type of image mask. + ** \param N Exact type of neighborhood. + ** + ** \arg marker Image to work on. + ** \arg mask Image used for geodesic erosion. + ** \arg Ng Neighborhood to use. + ** + ** \code + ** #include <oln/basics2d.hh> + ** #include <oln/morpho/opening.hh> + ** #include <oln/morpho/reconstruction.hh> + ** #include <oln/level/compare.hh> + ** #include <ntg/all.hh> + ** int main() + ** { + ** typedef oln::image2d<ntg::int_u8> im_type; + ** + ** im_type im1(oln::load(IMG_IN "lena128.pgm")); + ** im_type im2 = oln::morpho::opening(im1, oln::win_c4p()); + ** + ** oln::save(oln::morpho::hybrid::geodesic_reconstruction_erosion(im1, + ** im2, + ** oln::neighb_c4()), + ** IMG_OUT "oln_morpho_hybrid_geodesic_reconstruction_erosion.pbm"); + ** return 0; + ** } + ** \endcode + ** + ** \image html lena128_pgm.png + ** \image latex lena128_pgm.png + ** => + ** \image html oln_morpho_hybrid_geodesic_reconstruction_erosion.png + ** \image latex oln_morpho_hybrid_geodesic_reconstruction_erosion.png + ** + */ + + namespace impl { + + template <typename I, typename N> + struct reconstruction_erosion_ret : public reconstruction_ret<I, N> + { + typedef reconstruction_ret<I, N> super_type; + + reconstruction_erosion_ret(const abstract::non_vectorial_image<I>& input1, //marker + const abstract::non_vectorial_image<I>& input2, //mask + const abstract::neighborhood<N>& ng) + + : super_type(input1, input2, ng) + {} + + void impl_run() + { + mlc::is_true<mlc::type::eq<oln_type_of(I, size), + oln_type_of(N, size)>::ret>::ensure(); + precondition(this->input1.size() == this->input2.size()); + precondition(level::is_greater_or_equal(this->input2, this->input1)); + + oln_type_of(I, concrete) output; + output = utils::clone(this->input1); + { + oln_type_of(N, window) se_plus = get_plus_se_p(convert::ng_to_cse(this->ng)); + oln_type_of(N, window) se_minus = get_minus_se_p(convert::ng_to_cse(this->ng)); + oln_type_of(I, fwd_piter) fwd_p(output.size()); + oln_type_of(I, fwd_piter) bkd_p(output.size()); + + for_all (fwd_p) + output[fwd_p] = ntg::max(morpho::min(output, fwd_p, se_plus), + this->input2[fwd_p].value()); + + std::queue<oln_type_of(I, point) > fifo; + for_all (bkd_p) + { + output[bkd_p] = ntg::max(morpho::min(output, bkd_p, se_minus), + this->input2[bkd_p].value()); + if (internal::exist_init_erosion((oln_type_of(I, point))bkd_p, output, + this->input2, se_minus)) + fifo.push(bkd_p); + } + // Propagation Step + while (!fifo.empty()) + { + oln_type_of(I, point) p = fifo.front(); + fifo.pop(); + typedef oln_type_of(N, window) window_type; + window_type w = convert::ng_to_se(this->ng); + oln_type_of(window_type, fwd_witer) dp(w); + + for_all (dp) + { + oln_type_of(I, point) q = (oln_type_of(window_type, dpoint))dp + p; + if (output.hold(q)) + { + if ((output[q] > output[p]) && + (this->input2[q] != output[q])) + { + output[q] = ntg::max(output[p].value(), + this->input2[q].value()); + fifo.push(q); + } + } + } + } + } + this->output = output; + } + }; + + } + + template<class I, class N> + reconstruction_ret<I, N> + geodesic_reconstruction_erosion(const abstract::non_vectorial_image<I> & marker, + const abstract::non_vectorial_image<I> & mask, + const abstract::neighborhood<N>& ng) + { + impl::reconstruction_erosion_ret<I, N> tmp(marker, mask, ng); + tmp.run(); + return tmp; + } + }// hybrid + + } + +} + +#endif // ! OLENA_MORPHO_RECONSTRUCTION_HH Index: oln/morpho/splitse.hh --- oln/morpho/splitse.hh (revision 0) +++ oln/morpho/splitse.hh (revision 0) @@ -0,0 +1,208 @@ +// Copyright (C) 2001, 2002, 2004, 2005 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, 59 Temple Place - Suite 330, 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 OLENA_MORPHO_SPLITSE_HH +# define OLENA_MORPHO_SPLITSE_HH + +# include <oln/basics.hh> + +# include <oln/core/1d/size1d.hh> +# include <oln/core/2d/size2d.hh> +# include <oln/core/3d/size3d.hh> + +namespace oln { + namespace morpho { + + template <class S> + struct dim_traits + { + }; + + template <> + struct dim_traits<size1d> + { + enum {dim = 1}; + }; + + template <> + struct dim_traits<size2d> + { + enum {dim = 2}; + }; + + template <> + struct dim_traits<size3d> + { + enum {dim = 3}; + }; + + /*! + ** \brief Get a sub part of a structuring element. + ** + ** \param E Exact type of the structuring element. + ** + ** \arg se The structuring element. + ** + ** A point p take part of the new structuring element if it exists + ** a i that belongs to [[0..dim-1]] like p(i) < 0 and for all j + ** that belongs to [[0..i-1]] p(j) = 0. + ** + */ + + + template<class E> + E + get_plus_se_only(const abstract::struct_elt<E>& se) + { + oln_type_of(E, fwd_witer) dp(se.exact()); + E out; + + for_all (dp) + { + unsigned n; + for (n = 0; n < dim_traits<oln_type_of(E, size)>::dim; ++n) + if (dp.nth(n) < 0) { + out.add(dp); + break; + } else if (dp.nth(n) > 0) { + break; + } + } + return out; + } + + /*! + ** \brief Get a sub part of a structuring element. + ** + ** \param E Exact type of the structuring element. + ** + ** \arg se The structuring element. + ** + ** A point p take part of the new structuring element if it exists + ** a i that belongs to [[0..dim-1]] like p(i) < 0 and for all j + ** that belongs to [[0..i-1]] p(j) = 0 or if for all i that + ** belongs to [[0..dim-1]] p(i) = 0. + ** + */ + template<class E> + E + get_plus_se_p(const abstract::struct_elt<E>& se) + { + oln_type_of(E, fwd_witer) dp(se.exact()); + E out; + + for_all (dp) + { + unsigned n; + for (n = 0; n < dim_traits<oln_type_of(E, size)>::dim; ++n) + if (dp.nth(n) < 0) { + out.add(dp); + break; + } else if (dp.nth(n) > 0) { + break; + } + // All p.nth(n) are 0. + if (n == dim_traits<oln_type_of(E, size)>::dim) + out.add(dp); + } + return out; + } + + /*! + ** \brief Get a sub part of a structuring element. + ** + ** \param E Exact type of the structuring element. + ** + ** \arg se The structuring element. + ** + ** A point p take part of the new structuring element if it exists + ** a i that belongs to [[0..dim-1]] like p(i) > 0 and for all j + ** that belongs to [[0..i-1]] p(j) = 0. + ** + */ + template<class E> + E + get_minus_se_only(const abstract::struct_elt<E>& se) + { + oln_type_of(E, fwd_witer) dp(se.exact()); + E out; + + for_all (dp) + { + unsigned n; + for (n = 0; n < dim_traits<oln_type_of(E, size)>::dim; ++n) + if (dp.nth(n) > 0) { + out.add(dp); + break; + } else if (dp.nth(n) < 0) { + break; + } + } + return out; + } + + /*! + ** \brief Get a sub part of a structuring element. + ** + ** \param E Exact type of the structuring element. + ** + ** \arg se The structuring element. + ** + ** A point p take part of the new structuring element if it exists + ** a i that belongs to [[0..dim-1]] like p(i) > 0 and for all j + ** that belongs to [[0..i-1]] p(j) = 0 or if for all i that + ** belongs to [[0..dim-1]] p(i) = 0. + ** + */ + template<class E> + E + get_minus_se_p(const abstract::struct_elt<E>& se) + { + oln_type_of(E, fwd_witer) dp(se.exact()); + E out; + + for_all (dp) + { + unsigned n; + for (n = 0; n < dim_traits<oln_type_of(E, size)>::dim; ++n) + if (dp.nth(n) > 0) { + out.add(dp); + break; + } else if (dp.nth(n) < 0) { + break; + } + // All p.nth(n) are 0. + if (n == dim_traits<oln_type_of(E, size)>::dim) + out.add(dp); + } + return out; + } + + } // morpho +} // oln + +#endif // OLENA_MORPHO_SPLITSE_HH Index: oln/morpho/stat.hh --- oln/morpho/stat.hh (revision 106) +++ oln/morpho/stat.hh (working copy) @@ -47,7 +47,7 @@ ** \param E Structuring element type. ** \param V Associated value type. */ - template <class I, class E, class V =oln_type_of(I, value)> + template <class I, class E, class V = oln_type_of(I, value)> struct stat_ { /*! Index: oln/morpho/erosion.hh --- oln/morpho/erosion.hh (revision 106) +++ oln/morpho/erosion.hh (working copy) @@ -43,7 +43,7 @@ /// Erosion as a procedure (do not use it; prefer morpho::erosion). template<typename I, typename S> - oln_type_of(I, concrete) dilation(const abstract::image<I>& input, + oln_type_of(I, concrete) erosion(const abstract::image<I>& input, const abstract::struct_elt<S>& se) { oln_type_of(I, concrete) output(input.size()); @@ -131,14 +131,14 @@ typedef typename super_type::output_type output_type; const E se; - + erosion_ret(const abstract::image<I>& input, const abstract::struct_elt<E>& se) : super_type(input), se(se.exact()) { } - + }; Index: oln/io/write_image_2d_pnm.hh --- oln/io/write_image_2d_pnm.hh (revision 106) +++ oln/io/write_image_2d_pnm.hh (working copy) @@ -122,7 +122,7 @@ bin_v = 0; ret = true; } - if (c == value_type(1)) + if (c == value_type(0)) bin_v |= 1 << bin_offset; bin_offset--; return ret; Index: oln/io/read_image_2d_pnm.hh --- oln/io/read_image_2d_pnm.hh (revision 106) +++ oln/io/read_image_2d_pnm.hh (working copy) @@ -172,9 +172,9 @@ offset = 7; } if ((int)(v & (1<<offset--)) == 0) + c = 0; + else c = 1; - else - c = 0; } };