[Olena-patches] 1728: Have Meyer's Watershed Transform work on line graph-based images.

https://svn.lrde.epita.fr/svn/oln/trunk/milena There's still a lot of work to do here (mainly clean up). But we should do that after we revamp the types associated to images (sites, etc.). Anyway, we now have a watershed transform working on images based on line graphs, and it's very cool! :) Index: ChangeLog from Roland Levillain &lt;roland(a)lrde.epita.fr&gt; Have Meyer's Watershed Transform work on line graph-based images. * mln/morpho/extrema_components.hh, * mln/morpho/level_components.hh: New. Ressurect and adjust those files from a revision 1700, since mln::labeling::regional_minima (needed by Meyer's WST) doesn't work (yet) on image which are not point-wise accessible. * mln/util/greater_psite.hh: New. Clone and modified from mln/util/greater_point.hh. * mln/morpho/meyer_wst.hh (morpho::meyer_wst): s/point/psite/. Use morpho::minima_components instead of level::regional_minima, since the latter doesn't work on non point-wise accessible images. * tests/core/line_graph_image_wst.cc: New test. * tests/core/Makefile.am (check_PROGRAMS): Add line_graph_image_wst. (line_graph_image_wst_SOURCES): New. mln/morpho/extrema_components.hh | 16 ++-- mln/morpho/level_components.hh | 18 +++- mln/morpho/meyer_wst.hh | 28 ++++--- mln/util/greater_psite.hh | 34 ++++---- tests/core/Makefile.am | 2 tests/core/line_graph_image_wst.cc | 141 +++++++++++++++++-------------------- 6 files changed, 126 insertions(+), 113 deletions(-) Index: mln/morpho/extrema_components.hh --- mln/morpho/extrema_components.hh (revision 1719) +++ mln/morpho/extrema_components.hh (working copy) @@ -75,7 +75,7 @@ template <typename DestValue, typename I, typename W> mln_ch_value(I, DestValue) minima_components(const Image<I>& input, const Window<W>& win, - unsigned& nminima); + DestValue& nminima); /// \brief Like the 3-argument version of mln::moprho::minima_components, /// without the out argument \c nminima. @@ -108,7 +108,7 @@ template <typename DestValue, typename I, typename W> mln_ch_value(I, DestValue) maxima_components(const Image<I>& input, const Window<W>& win, - unsigned& nmaxima); + DestValue& nmaxima); /// \brief Like the 3-argument version of mln::moprho::maxima_components, /// without the out argument \c nmaxima. @@ -126,7 +126,7 @@ template <typename DestValue, typename Cmp, typename I, typename W> mln_ch_value(I, DestValue) extrema_components(const Image<I>& input_, const Window<W>& win_, - unsigned& nextrema) + DestValue& nextrema) { /* FIXME: Errors due to a too small DestValue type should be @@ -170,6 +170,7 @@ mln_qiter(W) q(win, p); for_all(q) + { if (input.has(q) and cmp(input(q), input(p))) { extrema.erase(comp); @@ -178,6 +179,7 @@ } } } + } // Update nextrema. nextrema = extrema.size(); @@ -212,7 +214,7 @@ template <typename DestValue, typename I, typename W> mln_ch_value(I, DestValue) minima_components(const Image<I>& input, const Window<W>& win, - unsigned& nminima) + DestValue& nminima) { /* FIXME: Ensure the input image has scalar values. */ typedef std::less< mln_value(I) > cmp_t; @@ -224,7 +226,7 @@ minima_components(const Image<I>& input, const Window<W>& win) { // Dummy value. - unsigned nminima; + DestValue nminima; return minima_components<DestValue>(input, win, nminima); } @@ -232,7 +234,7 @@ template <typename DestValue, typename I, typename W> mln_ch_value(I, DestValue) maxima_components(const Image<I>& input, const Window<W>& win, - unsigned& nmaxima) + DestValue& nmaxima) { /* FIXME: Ensure the input image has scalar values. */ typedef std::greater< mln_value(I) > cmp_t; @@ -244,7 +246,7 @@ maxima_components(const Image<I>& input, const Window<W>& win) { // Dummy value. - unsigned nmaxima; + DestValue nmaxima; return maxima_components<DestValue>(input, win, nmaxima); } Index: mln/morpho/level_components.hh --- mln/morpho/level_components.hh (revision 1719) +++ mln/morpho/level_components.hh (working copy) @@ -99,13 +99,23 @@ const W& win = exact(win_); mln_ch_value(I, DestValue) labels(input.domain()); - mln_ch_value(I, bool) processed(input.domain()); + /* FIXME: Yet another KLUDGE, this time required by the + specialization std::vector<bool>, which prevents forming + (mutable) reference to any of its elements. This creates + errors later with images using std::vector<bool> to store + their data (e.g., line_graph_image<P, V>). + + Alas, we cannot prevent the compiler to use this + specialization. Our workaround is simply... to use integers + instead of booleans. */ +// mln_ch_value(I, bool) processed(input.domain()); + mln_ch_value(I, int) processed(input.domain()); level::fill (processed, false); DestValue cur_label = mln_min(DestValue); - typedef mln_point(I) point_type; - std::queue<point_type> queue; + typedef mln_psite(I) psite_type; + std::queue<psite_type> queue; mln_piter(I) p(input.domain()); for_all(p) @@ -116,7 +126,7 @@ queue.push(p); while (!queue.empty()) { - point_type s = queue.front(); + psite_type s = queue.front(); queue.pop(); mln_qiter(W) q(win, s); for_all(q) Index: mln/util/greater_psite.hh --- mln/util/greater_psite.hh (revision 1719) +++ mln/util/greater_psite.hh (working copy) @@ -25,8 +25,8 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_UTIL_GREATER_POINT_HH -# define MLN_UTIL_GREATER_POINT_HH +#ifndef MLN_UTIL_GREATER_PSITE_HH +# define MLN_UTIL_GREATER_PSITE_HH # include <mln/core/concept/image.hh> @@ -35,55 +35,55 @@ namespace util { - /** \brief A ``greater than'' functor comparing points w.r.t. the + /** \brief A ``greater than'' functor comparing psites w.r.t. the values they refer to in an image. This functor used in useful to implement ordered queues of - points. */ + psites. */ template <typename I> - class greater_point + class greater_psite { public: - typedef mln_point(I) point; + typedef mln_psite(I) psite; - greater_point(const Image<I>& ima); + greater_psite(const Image<I>& ima); /// Is \a x greater than \a y? - bool operator()(const point& x, const point& y); + bool operator()(const psite& x, const psite& y); private: const I& ima_; }; - /// Helper to build a mln::util::greater_point. + /// Helper to build a mln::util::greater_psite. /* FIXME: To be moved into mln/make/? */ template <typename I> - greater_point<I> - make_greater_point(const Image<I>& ima); + greater_psite<I> + make_greater_psite(const Image<I>& ima); # ifndef MLN_INCLUDE_ONLY template <typename I> - greater_point<I>::greater_point(const Image<I>& ima) + greater_psite<I>::greater_psite(const Image<I>& ima) : ima_ (exact(ima)) { } template <typename I> bool - greater_point<I>::operator()(const mln_point(I)& x, const mln_point(I)& y) + greater_psite<I>::operator()(const mln_psite(I)& x, const mln_psite(I)& y) { return ima_(x) > ima_(y); } template <typename I> - greater_point<I> - make_greater_point(const Image<I>& ima) + greater_psite<I> + make_greater_psite(const Image<I>& ima) { - return greater_point<I>(ima); + return greater_psite<I>(ima); } # endif // ! MLN_INCLUDE_ONLY @@ -92,4 +92,4 @@ } // end of namespace mln -#endif // ! MLN_UTIL_GREATER_POINT_HH +#endif // ! MLN_UTIL_GREATER_PSITE_HH Index: mln/morpho/meyer_wst.hh --- mln/morpho/meyer_wst.hh (revision 1727) +++ mln/morpho/meyer_wst.hh (working copy) @@ -43,9 +43,11 @@ # include <mln/trait/ch_value.hh> -# include <mln/util/greater_point.hh> +# include <mln/util/greater_psite.hh> # include <mln/morpho/includes.hh> -# include <mln/labeling/regional_minima.hh> +// FIXME: See below. +// # include <mln/labeling/regional_minima.hh> +# include <mln/morpho/extrema_components.hh> @@ -112,15 +114,23 @@ const marker unmarked = mln_min(marker); // Initialize the output with the markers (minima components). + /* FIXME: labeling::regional_minima doesn't work on non + point-accessible images! We temporarily reuse the old + morpho::minima_components routine to work-around this + limitation. As soon as labeling::regional_minima work, get rid of + - mln/morpho/level_components.hh, and + - mln/morpho/extrema_components.hh. */ +// mln_ch_value(I, marker) markers = +// labeling::regional_minima (input, nbh, nbasins); mln_ch_value(I, marker) markers = - labeling::regional_minima (input, nbh, nbasins); + minima_components(input, convert::to_window(nbh), nbasins); // Ordered queue. - typedef mln_point(I) point; + typedef mln_psite(I) psite; typedef - std::priority_queue< point, std::vector<point>, util::greater_point<I> > + std::priority_queue< psite, std::vector<psite>, util::greater_psite<I> > ordered_queue_type; - ordered_queue_type queue(util::make_greater_point(input)); + ordered_queue_type queue(util::make_greater_psite(input)); // Insert every neighbor P of every marked area in a // hierarchical queue, with a priority level corresponding to @@ -136,12 +146,12 @@ break; } - /* Until the queue is empty, extract a point P from the + /* 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.empty()) { - point p = queue.top(); + psite p = queue.top(); queue.pop(); // Last seen marker adjacent to P. marker adjacent_marker = unmarked; @@ -161,7 +171,7 @@ single_adjacent_marker_p = false; break; } - /* If the neighborhood of P contains only points with the + /* 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. */ Index: tests/core/line_graph_image_wst.cc --- tests/core/line_graph_image_wst.cc (revision 1719) +++ tests/core/line_graph_image_wst.cc (working copy) @@ -25,17 +25,17 @@ // reasons why the executable file might be covered by the GNU General // Public License. -/// \file tests/core/graph_image.cc -/// \brief Tests on mln::graph_image. +/// \file tests/core/line_graph_image.cc +/// \brief Tests on the Watershed Transform on a mln::line_graph_image. #include <vector> #include <mln/core/point2d.hh> #include <mln/core/line_graph_image.hh> -#include <mln/core/line_graph_elt_window.hh> -#include <mln/core/line_graph_window_piter.hh> +#include <mln/core/line_graph_elt_neighborhood.hh> +#include <mln/core/line_graph_neighborhood_piter.hh> -#include <mln/morpho/dilation.hh> +#include <mln/morpho/meyer_wst.hh> #include <mln/draw/graph.hh> #include <mln/debug/iota.hh> @@ -46,48 +46,58 @@ { using namespace mln; - /*--------. - | Graph. | - `--------*/ - - // Points associated to nodes. - std::vector<point2d> points; - points.push_back(make::point2d(0,0)); // Point associated to node 0. - points.push_back(make::point2d(2,2)); // Point associated to node 1. - points.push_back(make::point2d(0,4)); // Point associated to node 2. - points.push_back(make::point2d(4,3)); // Point associated to node 3. - points.push_back(make::point2d(4,4)); // Point associated to node 4. + /*-------------. + | Line graph. | + `-------------*/ + + /* Actually this graph is from Jean Cousty's PhD dissertation, page 76. + + 0 1 2 3 (rows) + ,------------------------ + | 0 10 5 + 0 | o-----o-----o-----o + | | | | | + | 2| 4| 6| 0| + | | | | | + 1 | o-----o-----o-----o + 3 5 2 + (cols) - // Edges. + In G, nodes and egdes are numbered following in the classical + foward order. */ util::graph<point2d> g; - // Populate the graph with nodes. - for (unsigned i = 0; i < points.size(); ++i) - g.add_node (points[i]); - // Populate the graph with edges. - g.add_edge(0, 1); - g.add_edge(1, 2); - g.add_edge(1, 3); - g.add_edge(3, 4); - g.add_edge(4, 2); - - /*-------. - | Graph. | - `-------*/ - p_line_graph<point2d> plg(g); + // Nodes. + std::vector<int> node_values; + // We don't really care about values on nodes, because the algorithm + // won't see them. So all are set to 0. + g.add_node (make::point2d(0, 0)); node_values.push_back (0); // Node 0. + g.add_node (make::point2d(0, 1)); node_values.push_back (0); // Node 1. + g.add_node (make::point2d(0, 2)); node_values.push_back (0); // Node 2. + g.add_node (make::point2d(0, 3)); node_values.push_back (0); // Node 3. + + g.add_node (make::point2d(1, 0)); node_values.push_back (0); // Node 4. + g.add_node (make::point2d(1, 1)); node_values.push_back (0); // Node 5. + g.add_node (make::point2d(1, 2)); node_values.push_back (0); // Node 6. + g.add_node (make::point2d(1, 3)); node_values.push_back (0); // Node 7. - /*-------------------. - | Line graph image. | - `-------------------*/ - - // Values ("empty" vectors). - std::vector<int> node_values(5); - std::vector<int> edge_values(5); - // FIXME: hand-made iota's. - for (unsigned i = 0; i < node_values.size(); ++i) - node_values[i] = i; - for (unsigned i = 0; i < edge_values.size(); ++i) - edge_values[i] = i; + // Edges. + std::vector<int> edge_values; + g.add_edge (0, 1); edge_values.push_back (0); + g.add_edge (1, 2); edge_values.push_back (10); + g.add_edge (2, 3); edge_values.push_back (5); + + g.add_edge (0, 4); edge_values.push_back (2); + g.add_edge (1, 5); edge_values.push_back (4); + g.add_edge (2, 6); edge_values.push_back (6); + g.add_edge (3, 7); edge_values.push_back (0); + + g.add_edge (4, 5); edge_values.push_back (3); + g.add_edge (5, 6); edge_values.push_back (5); + g.add_edge (6, 7); edge_values.push_back (2); + + // Line graph point set. + p_line_graph<point2d> plg(g); // Line graph image. /* FIXME: We probably don't want to build line_graph_images by hand; @@ -103,40 +113,19 @@ mln_piter_(ima_t) p(ima.domain()); for_all (p) std::cout << "ima (" << p << ") = " << ima(p) << std::endl; - - // Manual iterations over the neighborhoods of each point site of IMA. - typedef line_graph_elt_window<point2d> win_t; - win_t win; - mln_qiter_(win_t) q(win, p); - for_all (p) - { - std::cout << "neighbors of " << p << " (" << ima(p) << "), " - << "including the site itself:" << std::endl; - for_all (q) - std::cout << " " << q << " (level = " << ima(q) << ")" - << std::endl; - } std::cout << std::endl; + typedef line_graph_elt_neighborhood<point2d> nbh_t; + nbh_t nbh; -// /* FIXME: When implementing convert::to_line_graph_image, don't -// forget to give a second argument defaulting to something like -// fun::vv2v::max(). This second argument is a functor used to -// compute the values of the edges of the line graph image. */ -// image2d ima; // = ... -// lg_ima_t lg_ima = convert::to_line_graph_image (ima); - -// // Image2d representation. -// image2d<value_t> out = convert::to_image2d (lg_ima); -// io::pgm::save(out, "out.pgm"); - -// /* FIXME: Then, add some real processing on the line graph image -// before converting to an image2d: -// - erosion/dilation, -// - Beucher's gradient, -// - Meyer's WST (on the gradient of LG_IMA?), -// - attribute filters (see my notes on Laurent's remarks about -// attributes), -// - etc. -// Creating seperate tests for all these would be wise. */ + typedef unsigned wst_val_t; + wst_val_t nbasins; + typedef line_graph_image<point2d, wst_val_t> wst_ima_t; + wst_ima_t wshed = morpho::meyer_wst(ima, nbh, nbasins); + std::cout << "nbasins = " << nbasins << std::endl; + + // Manual iteration over the domain of WSHED. + mln_piter_(wst_ima_t) pw(wshed.domain()); + for_all (pw) + std::cout << "wshed (" << pw << ") = " << wshed(pw) << std::endl; } Index: tests/core/Makefile.am --- tests/core/Makefile.am (revision 1727) +++ tests/core/Makefile.am (working copy) @@ -12,6 +12,7 @@ graph_image \ line_graph_elt_window \ line_graph_image \ + line_graph_image_wst \ mono_obased_rle_image \ mono_rle_image \ obased_rle_image \ @@ -35,6 +36,7 @@ graph_image_SOURCES = graph_image.cc line_graph_elt_window_SOURCES = line_graph_elt_window.cc line_graph_image_SOURCES = line_graph_image.cc +line_graph_image_wst_SOURCES = line_graph_image_wst.cc mono_obased_rle_image_SOURCES = mono_obased_rle_image.cc mono_rle_image_SOURCES = mono_rle_image.cc obased_rle_image_SOURCES = obased_rle_image.cc