[Olena-patches] proto-1.0 198: Add shortest path watershed transform

https://svn.lrde.epita.fr/svn/oln/prototypes/proto-1.0 ChangeLog | 23 +++ oln/core/abstract/niter.hh | 5 oln/core/gen/image_with_nbh.hh | 8 + oln/core/gen/internal/value_box.hh | 14 + oln/level/extrema_components.hh | 35 +++- oln/level/level_components.hh | 9 - oln/makefile.src | 3 oln/morpho/lower_completion.hh | 9 - oln/morpho/shortest_path_watershed.hh | 207 +++++++++++++++++++++++++++++ tests/core/tests/niter2d | 41 +++++ tests/morpho/tests/shortest_path_watershed | 30 ++++ 11 files changed, 365 insertions(+), 19 deletions(-) Index: olena/ChangeLog from Roland Levillain &lt;roland(a)lrde.epita.fr&gt; Add shortest path watershed transform. * oln/morpho/shortest_path_watershed.hh: New file. * tests/morpho/tests/shortest_path_watershed: New test. * oln/core/gen/image_with_nbh.hh (image_with_nbh(const oln_type_of(I, size)&, const abstract::neighborhood<N>&)): New ctor. * oln/level/level_components.hh (level_components), * oln/morpho/lower_completion.hh (lower_completion): Use it. * oln/core/abstract/niter.hh (nbh): New method. * tests/core/tests/niter2d: New test. * oln/level/extrema_components.hh (extrema_components): Use unsigned to compute the level components. * oln/core/gen/internal/value_box.hh (operator>(const value_box<II>&), operator>(const V&)): New operators. * oln/makefile.src (OLN_DEP): Add level/extrema_components.hh, level/level_components.hh and morpho/shortest_path_watershed.hh. Index: olena/tests/morpho/tests/shortest_path_watershed --- olena/tests/morpho/tests/shortest_path_watershed (revision 0) +++ olena/tests/morpho/tests/shortest_path_watershed (revision 0) @@ -0,0 +1,30 @@ + // -*- C++ -*- +#include "data.hh" + +#include <ntg/int.hh> +#include <oln/core/2d/image2d.hh> +#include <oln/core/gen/image_with_nbh.hh> +#include <oln/io/read_image.hh> +#include <oln/morpho/shortest_path_watershed.hh> +#include <oln/utils/md5.hh> + +using namespace oln; + +bool check() +{ + // MD5 sum of object.pbm's watershed transform result. + utils::key::value_type data_key[16] = + { 0x47, 0x69, 0xe9, 0x5c, 0x82, 0xe6, 0x2e, 0xc4, + 0x0, 0x72, 0x4a, 0xf5, 0x5f, 0x69, 0xf, 0x2b }; + utils::key key(data_key); + + image2d<ntg::int_u8> raw_input; + raw_input = io::read(rdata("squares.pgm")); + typedef image_with_nbh<image2d<ntg::int_u8>, neighborhood2d> image_type; + image_type input = join(raw_input, neighb_c4()); + + if (utils::md5(morpho::shortest_path_watershed<ntg::int_u8>(input)) != key) + return true; + + return false; +} Index: olena/tests/core/tests/niter2d --- olena/tests/core/tests/niter2d (revision 0) +++ olena/tests/core/tests/niter2d (revision 0) @@ -0,0 +1,41 @@ + // -*- C++ -*- +#include <iostream> + +#include "ntg/all.hh" +#include "oln/basics2d.hh" +#include "oln/core/gen/image_with_nbh.hh" + +#include "check.hh" +#include "data.hh" + +bool check() +{ + oln::image2d<int> ima_without_nbh(10, 10); + oln::image_with_nbh<oln::image2d<int>, oln::neighborhood2d> ima = + oln::join(ima_without_nbh, oln::neighb_c4()); + int cpt = 0; + + oln::fwd_piter2d fwd_it(ima.size()); + for_all_p (fwd_it) + { + oln::fwd_niter2d fwd_nit(ima); + for_all_n_of_p (fwd_nit, fwd_it) + cpt++; + } + + if (cpt != 400) + return true; + + cpt = 0; + for_all_p (fwd_it) + { + oln::bkd_niter2d bkd_nit(ima); + for_all_n_of_p (bkd_nit, fwd_it) + cpt++; + } + + if (cpt != 400) + return true; + + return false; +} Index: olena/oln/core/abstract/niter.hh --- olena/oln/core/abstract/niter.hh (revision 197) +++ olena/oln/core/abstract/niter.hh (working copy) @@ -119,6 +119,11 @@ this->invalidate(); } + neighb_type nbh() + { + return nbh_; + } + protected: const neighb_type nbh_; // copy is safe Index: olena/oln/core/gen/image_with_nbh.hh --- olena/oln/core/gen/image_with_nbh.hh (revision 197) +++ olena/oln/core/gen/image_with_nbh.hh (working copy) @@ -80,6 +80,14 @@ this->image_ = tmp; } + image_with_nbh(const oln_type_of(I, size)& size, + const abstract::neighborhood<N>& nbh) : + nbh_(nbh.exact()) + { + I tmp(size); // FIXME: hack + this->image_ = tmp; + } + image_with_nbh(abstract::image<I>& image, const abstract::neighborhood<N>& nbh) : super_type(image), Index: olena/oln/core/gen/internal/value_box.hh --- olena/oln/core/gen/internal/value_box.hh (revision 197) +++ olena/oln/core/gen/internal/value_box.hh (working copy) @@ -114,6 +114,20 @@ return this->value() < rhs; } + /// Operator > (rhs is a value_box). + template <typename II> + bool operator>(const value_box<II>& rhs) const + { + return this->value() > rhs.value(); + } + + /// Operator > (rhs is a value). + template <typename V> + bool operator>(const V& rhs) const + { + return this->value() > rhs; + } + /*! \brief op= ** FIXME:... Index: olena/oln/makefile.src --- olena/oln/makefile.src (revision 197) +++ olena/oln/makefile.src (working copy) @@ -150,7 +150,9 @@ \ level/arith.hh \ level/compare.hh \ + level/extrema_components.hh \ level/fill.hh \ + level/level_components.hh \ level/logic.hh \ \ morpho/closing.hh \ @@ -166,6 +168,7 @@ morpho/reconstruction_by_erosion.hh \ morpho/reconstruction_selfdual.hh \ morpho/reconstruction.hh \ + morpho/shortest_path_watershed.hh \ morpho/stat.hh \ morpho/tags.hh \ \ Index: olena/oln/morpho/lower_completion.hh --- olena/oln/morpho/lower_completion.hh (revision 197) +++ olena/oln/morpho/lower_completion.hh (working copy) @@ -50,10 +50,12 @@ typedef std::queue<point_type> queue_type; queue_type q; - typename ch_value_type<I, bool>::ret processed(input.size()); + typename ch_value_type<I, bool>::ret processed(input.size(), + input.nbh_get ()); level::fill (processed, false); - typename ch_value_type<I, unsigned>::ret distance(input.size()); + typename ch_value_type<I, unsigned>::ret distance(input.size(), + input.nbh_get()); unsigned cur_dist = 1; oln_type_of(I, piter) p(input.size()); @@ -108,7 +110,8 @@ | Lower completion. | `-------------------*/ - typename ch_value_type<I, DestValue>::ret output(input.size()); + typename ch_value_type<I, DestValue>::ret output(input.size(), + input.nbh_get()); for_all_p (p) if (distance[p] == ntg_zero_val(DestValue)) output[p] = cur_dist * input[p].value(); Index: olena/oln/morpho/shortest_path_watershed.hh --- olena/oln/morpho/shortest_path_watershed.hh (revision 0) +++ olena/oln/morpho/shortest_path_watershed.hh (revision 0) @@ -0,0 +1,207 @@ +// Copyright (C) 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_SHORTEST_PATH_WATERSHED_HH +# define OLENA_MORPHO_SHORTEST_PATH_WATERSHED_HH + +# include <queue> + +# include "oln/level/fill.hh" +# include "oln/level/extrema_components.hh" +# include "oln/morpho/lower_completion.hh" + +namespace oln { + + namespace morpho { + + namespace internal { + + /// Functor used in ordered queues of points. + template <typename I> + class greater_point_type + { + public: + typedef oln_type_of(I, point) point_type; + + greater_point_type(const abstract::image<I>& im) : + im_ (im) + { + } + + /// Is \a x greater than \a y? + bool operator()(const point_type& x, const point_type& y) + { + return im_[x] > im_[y]; + } + + private: + const abstract::image<I>& im_; + }; + + /// Facade to build a oln::level::greater_point_type. + template <typename I> + greater_point_type<I> + greater_point(const abstract::image<I>& im) + { + return greater_point_type<I>(im); + } + + + // FIXME: To be rewritten. Moreover, the distance d(p, q) is not + // taken into account in this version of cost (but it should not + // bother us as long as we are using first-order neighborhoods). + template <typename I> + int + cost (const abstract::image_with_nbh<I>& im, + const oln_type_of(I, point)& p, const oln_type_of(I, point)& q) + { + mlc_is_a(I, abstract::scalar_valued_image)::ensure(); + + oln_type_of(I, niter) n(im); + + int ls_p = 0; + for_all_n_of_p (n, p) + if (im.hold(n) and int(im[p] - im[n]) > ls_p) + ls_p = im[p] - im[n]; + + int ls_q = 0; + for_all_n_of_p (n, q) + if (im.hold(n) and int(im[q] - im[n]) > ls_q) + ls_q = im[q] - im[n]; + + if (im[p] > im[q]) + return ls_p; + else if (im[q] > im[p]) + return ls_q; + else + return (ls_p + ls_q) / 2; + } + + } // End of namespace internal. + + + template <typename DestValue, typename I> + typename ch_value_type<I, DestValue>::ret + shortest_path_watershed(const abstract::image_with_nbh<I>& input) + { + mlc_is_a(I, abstract::scalar_valued_image)::ensure(); + + const DestValue wshed = ntg_min_val(DestValue); + + // Get a lower complete image of the input. + typename ch_value_type<I, unsigned>::ret lc = + lower_completion<unsigned>(input); + + // We keep a track of already processed points. + typename ch_value_type<I, bool>::ret processed (lc.size(), + lc.nbh_get()); + level::fill (processed, false); + + // Initialise output with the minima components. + typename ch_value_type<I, DestValue>::ret output = + level::minima_components<DestValue>(lc); + + // Distance. + typedef typename ch_value_type<I, unsigned>::ret dist_type; + dist_type dist (lc.size(), lc.nbh_get()); + level::fill(dist, ntg_max_val(DestValue)); + // Initialize distance with values of minima, and mark these + // points as processed (remember that points of LC who have a + // value greater than ntg_min_val(DestValue) belongs to a + // minimum). + oln_type_of(I, piter) p(lc.size()); + for_all_p (p) + if (output[p].value() > ntg_min_val(DestValue)) + dist[p] = lc[p]; + + // Ordered queue. + typedef oln_type_of(I, point) point_type; + typedef + std::priority_queue<point_type, + std::vector<point_type>, + internal::greater_point_type<dist_type> > + ordered_queue_type; + ordered_queue_type q(internal::greater_point(dist)); + // Fill the ordered queue with the points of the border of the + // minima of the (lower complete) input. + for_all_p (p) + if (output[p].value() > ntg_min_val(DestValue)) + { + bool border_p = false; + oln_type_of(I, niter) n(lc); + for_all_n_of_p (n, p) + if (lc.hold(n) and output[n].value() == ntg_min_val(DestValue)) + { + // P is both in a minima and adjacent to a non minima: + // it's on the border of a minima. + border_p = true; + break; + }; + if (border_p) + q.push(p); + else + // Inner points of minima have already been processed. + processed[p] = true; + } + + while (!q.empty()) + { + point_type p = q.top(); + q.pop(); + if (processed[p]) + continue; + processed[p] = true; + + oln_type_of(I, niter) n(lc); + for_all_n_of_p (n, p) + { + if (lc.hold(n)) + if (dist[p] + internal::cost(lc, p, n) < dist[n]) + { + dist[n] = dist[p] + internal::cost(lc, p, n); + output[n] = output[p]; + q.push(n); + + } + else if (output[n] != wshed and + dist[p] + internal::cost(lc, p, n) == dist[n] and + output[n] != output[p]) + { + output[n] = wshed; + q.push(n); + } + } + } + + return output; + } + + } // end of namespace oln::morpho + +} // end of namespace oln + +#endif // ! OLENA_MORPHO_SHORTEST_PATH_WATERSHED_HH Index: olena/oln/level/extrema_components.hh --- olena/oln/level/extrema_components.hh (revision 197) +++ olena/oln/level/extrema_components.hh (working copy) @@ -50,16 +50,17 @@ Cmp cmp; // Compute level components. - typename ch_value_type<I, DestValue>::ret comps = - level_components<DestValue>(input); - std::set<DestValue> extrema; - std::set<DestValue> non_extrema; + typedef unsigned comp_type; + typename ch_value_type<I, comp_type>::ret comps = + level_components<comp_type>(input); + std::set<comp_type> extrema; + std::set<comp_type> non_extrema; // Search extrema components. oln_type_of(I, piter) p(input.size()); for_all_p (p) { - DestValue comp = comps[p]; + comp_type comp = comps[p]; if (non_extrema.find(comp) == non_extrema.end ()) { // A new level is a (potential) extrema by default. @@ -71,23 +72,27 @@ { extrema.erase(comp); non_extrema.insert(comp); + break; } } } // Re-label the extrema. label_map hold the assigned labels. - std::map<DestValue, DestValue> label_map; + std::map<comp_type, DestValue> label_map; { DestValue cur_label = ntg_min_val(DestValue) + 1; - for (typename std::set<DestValue>::const_iterator i = + for (typename std::set<comp_type>::const_iterator i = extrema.begin(); i != extrema.end(); ++i, ++cur_label) + { label_map[*i] = cur_label; } - typename ch_value_type<I, DestValue>::ret output (input.size()); + } + typename ch_value_type<I, DestValue>::ret output (input.size(), + input.nbh_get()); level::fill (output, ntg_min_val(DestValue)); for_all_p (p) { - DestValue comp = comps[p]; + comp_type comp = comps[p]; if (label_map.find(comp) != label_map.end()) output[p] = label_map[comp]; } @@ -103,8 +108,10 @@ minima_components(const abstract::image_with_nbh<I>& input) { mlc_is_a(I, abstract::scalar_valued_image)::ensure(); - return internal::extrema_components< DestValue, - std::less<DestValue> >(input); + typedef oln_type_of(I, value) input_value_type; + return + internal::extrema_components< DestValue, + std::less<input_value_type> >(input); } /// Find the maxima level components of \a input. @@ -113,8 +120,10 @@ maxima_components(const abstract::image_with_nbh<I>& input) { mlc_is_a(I, abstract::scalar_valued_image)::ensure(); - return internal::extrema_components< DestValue, - std::greater<DestValue> >(input); + typedef oln_type_of(I, value) input_value_type; + return + internal::extrema_components< DestValue, + std::greater<input_value_type> >(input); } } // end of namespace oln::level. Index: olena/oln/level/level_components.hh --- olena/oln/level/level_components.hh (revision 197) +++ olena/oln/level/level_components.hh (working copy) @@ -43,8 +43,10 @@ { mlc_is_a(I, abstract::scalar_valued_image)::ensure(); - typename ch_value_type<I, DestValue>::ret labels(input.size()); - typename ch_value_type<I, bool>::ret processed(input.size()); + typename ch_value_type<I, DestValue>::ret labels(input.size(), + input.nbh_get()); + typename ch_value_type<I, bool>::ret processed(input.size(), + input.nbh_get()); level::fill (processed, false); DestValue cur_label = ntg_min_val(DestValue); @@ -57,6 +59,7 @@ if (!processed[p]) { labels[p] = cur_label; + processed[p] = true; q.push(p); while (!q.empty()) { @@ -64,7 +67,7 @@ q.pop(); oln_type_of(I, niter) n(input); for_all_n_of_p (n, s) - if (input.hold(n) && input[n] == input[s] && !processed[n]) + if (input.hold(n) && !processed[n] && input[n] == input[s]) { labels[n] = cur_label; processed[n] = true;