https://svn.lrde.epita.fr/svn/oln/branches/cleanup-2008/milena Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; Add a 2D simplicial complex (toy code) + a double-window type. * sandbox/geraud/cs2d: New directory. * sandbox/geraud/cs2d/cs2d.cc: New. * sandbox/geraud/cs2d/dbl_neighb.cc: New. * sandbox/geraud/cs2d/cs2d.hh: New. * sandbox/geraud/cs2d/cs2d_utils.hh: New. * sandbox/geraud/cs2d/dbl_neighb.hh: New. * sandbox/geraud/cs2d/cs2d_morpho.hh: New. cs2d.cc | 110 ++++++++++++++++ cs2d.hh | 172 ++++++++++++++++++++++++++ cs2d_morpho.hh | 372 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ cs2d_utils.hh | 251 ++++++++++++++++++++++++++++++++++++++ dbl_neighb.cc | 112 +++++++++++++++++ dbl_neighb.hh | 167 +++++++++++++++++++++++++ 6 files changed, 1184 insertions(+) Index: sandbox/geraud/cs2d/cs2d.cc --- sandbox/geraud/cs2d/cs2d.cc (revision 0) +++ sandbox/geraud/cs2d/cs2d.cc (revision 0) @@ -0,0 +1,110 @@ +# include "cs2d.hh" +# include "cs2d_morpho.hh" +# include "dbl_neighb.hh" + +# include <mln/debug/println.hh> + +# include <mln/io/pgm/load.hh> +# include <mln/io/pgm/save.hh> +# include <mln/io/ppm/save.hh> + +# include <mln/level/paste.hh> +# include <mln/level/fill.hh> + +# include <mln/morpho/dilation.hh> +# include <mln/morpho/meyer_wst.hh> +# include <mln/morpho/closing_area.hh> + +# include <mln/logical/not.hh> +# include <mln/labeling/blobs.hh> +# include <mln/accu/mean.hh> + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm lambda output.pgm" << std::endl; + abort(); +} + + +struct is_row_odd_t +{ + bool operator()(const mln::point2d& p) const + { + return p.row() % 2; + } +} is_row_odd; + + + +int main(int argc, char* argv[]) +{ + if (argc != 4) + usage(argv); + + + using namespace mln; + using value::int_u8; + + // trace::quiet = false; + + cmorpho::dbl_neighb2d nbh_e; + nbh_e + .insert_even(dpoint2d(-1,-1)) + .insert_even(dpoint2d(-1,+1)) + .insert_even(dpoint2d( 0,-2)) + .insert_odd(dpoint2d(-2, 0)) + .insert_odd(dpoint2d(-1,-1)) + .insert_odd(dpoint2d(-1,+1)); + + // display_masks(make::box2d(2,3)); + + image2d<int_u8> lena; + io::pgm::load(lena, argv[1]); + + cs2d<int_u8, int_u8, int_u8> cs(lena.domain()); + + fill_cells_with(cs, lena); + level::fill(cs.vertex, 0); + + cmorpho::dbl_neighb2d nbh_g; + nbh_g.insert_odd(up).insert_even(left); + level::paste(cmorpho::gradient(cs.edge, nbh_g), + cs.edge); + // io::pgm::save(cs.impl, "grad.pgm"); + + int lambda = atoi(argv[2]); + if (lambda < 0) + usage(argv); + + if (lambda != 0) + level::paste(morpho::closing_area(cs.edge, nbh_e, lambda), cs.edge); + + cs2d<int_u8, bool, bool> cs_(lena.domain()); + level::fill(cs_.impl, false); + fill_cells_with(cs_, lena); + + unsigned l; + level::fill(cs_.edge, + pw::value(cmorpho::meyer_wst(cs.edge, nbh_e, l)) == pw::cst(0)); + std::cout << "l = " << l << std::endl; + + level::paste(morpho::dilation(cs_.vertex, + convert::to_window(c4())), + cs_.vertex); + + unsigned L; + // debug::println(cs_.impl); + image2d<unsigned> label = labeling::blobs(logical::not_(cs_.impl), c4(), L); + // debug::println(label); + mln_invariant(L == l); + + std::vector<int_u8> v(l + 1); + accu::compute<accu::mean>(cs_.cell, label, v); + + level::fill(cs_.cell, level::transform(label, v)); + io::pgm::save(cell_image(cs_), argv[3]); + +// io::ppm::save(pack(cs_, colorize, 7), "log.ppm"); +} Index: sandbox/geraud/cs2d/dbl_neighb.cc --- sandbox/geraud/cs2d/dbl_neighb.cc (revision 0) +++ sandbox/geraud/cs2d/dbl_neighb.cc (revision 0) @@ -0,0 +1,112 @@ +# include <cstdlib> + +# include <mln/value/int_u8.hh> +# include <mln/value/rgb8.hh> + +# include <mln/core/image2d.hh> +# include <mln/core/neighb2d.hh> +# include <mln/core/window2d.hh> + +# include <mln/convert/to_p_set.hh> +# include <mln/convert/to_image.hh> +# include <mln/debug/println.hh> + +# include "dbl_neighb.hh" + +# include <mln/io/pgm/load.hh> +# include <mln/io/pgm/save.hh> +# include <mln/io/pbm/save.hh> +# include <mln/io/ppm/save.hh> + +# include <mln/level/transform.hh> +# include <mln/literal/black.hh> +# include <mln/pw/all.hh> + +# include <mln/morpho/closing_area.hh> +# include <mln/morpho/gradient.hh> +# include <mln/morpho/meyer_wst.hh> + + + + +struct is_chess_t +{ + bool operator()(const mln::point2d& p) const + { + return p.col() % 2 == p.row() % 2; + } +}; + + +struct colorize : mln::Function_v2v< colorize > +{ + typedef mln::value::rgb8 result; + colorize(unsigned max) + : lut(max + 1) + { + lut[0] = mln::literal::black; + for (unsigned i = 1; i <= max; ++i) + lut[i] = result(100 + std::rand() % 150, + 100 + std::rand() % 150, + 100 + std::rand() % 150); + } + result operator()(unsigned i) const + { + return lut[i]; + } + std::vector<result> lut; +}; + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm lambda output.ppm" << std::endl; + abort(); +} + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + if (argc != 4) + usage(argv); + + image2d<int_u8> lena; + io::pgm::load(lena, argv[1]); + + bool vert[] = { 1, 1, 0, + 1, 0, 1, + 0, 1, 1 }; + + bool hori[] = { 0, 1, 1, + 1, 0, 1, + 1, 1, 0 }; + + dbl_neighb_<dpoint2d, is_chess_t> nbh; + nbh + .when_true (make::neighb2d(vert)) + .when_false(make::neighb2d(hori)); + + image2d<int_u8> clo, grad = morpho::gradient(lena, nbh); + // io::pgm::save(grad, "grad.pgm"); + + int lambda = atoi(argv[2]); + if (lambda > 1) + { + clo = morpho::closing_area(grad, nbh, lambda); + // io::pgm::save(clo, "clo.pgm"); + } + else + clo = grad; + + unsigned l; + image2d<unsigned> wst = morpho::meyer_wst(clo, nbh, l); + std::cout << "l = " << l << std::endl; + + io::ppm::save(level::transform(wst, colorize(l)), argv[3]); + +// io::pbm::save((pw::value(wst) == pw::cst(0)) | lena.domain(), +// argv[3]); +} Index: sandbox/geraud/cs2d/cs2d.hh --- sandbox/geraud/cs2d/cs2d.hh (revision 0) +++ sandbox/geraud/cs2d/cs2d.hh (revision 0) @@ -0,0 +1,172 @@ +# include <mln/core/image2d.hh> +# include <mln/core/neighb2d.hh> +# include <mln/core/window2d.hh> +# include <mln/convert/to_window.hh> + +# include <mln/core/image_if.hh> +# include <mln/debug/println.hh> + + +# include "cs2d_utils.hh" + + +namespace mln +{ + + + template <typename C, typename E, typename X> + struct cs2d + { + typedef C cell_value; + typedef E edge_value; + typedef X vertex_value; + + cs2d(unsigned nrows, unsigned ncols) + { + init(nrows, ncols); + } + + cs2d(const box2d& b) + { + init(b.nrows(), b.ncols()); + } + + void init(unsigned nrows, unsigned ncols) + { + mln_precondition(nrows && ncols); + + // Raw data (impl). + b = make::box2d(2 * nrows - 1, 2 * ncols - 1); + cimpl = image2d<C>(b, 0); + eimpl = image2d<E>(b, 0); + ximpl = image2d<X>(b, 0); + + // Images. + cell = cimpl | is_cell; + edge = eimpl | is_edge; + vertex = ximpl | is_vertex; + } + + unsigned nrows() const { return (b.nrows() + 1) / 2; } + unsigned ncols() const { return (b.ncols() + 1) / 2; } + + box2d b; + image2d<C> cimpl; + image2d<E> eimpl; + image2d<X> ximpl; + + typedef image_if< image2d<C>, is_cell_t > cell_t; + cell_t cell; + + typedef image_if< image2d<E>, is_edge_t > edge_t; + edge_t edge; + + typedef image_if< image2d<X>, is_vertex_t > vertex_t; + vertex_t vertex; + }; + + + + + template <typename T> + struct cs2d< T, T, T> + { + typedef T cell_value; + typedef T edge_value; + typedef T vertex_value; + + cs2d(unsigned nrows, unsigned ncols) + { + init(nrows, ncols); + } + + cs2d(const box2d& b) + { + init(b.nrows(), b.ncols()); + } + + void init(unsigned nrows, unsigned ncols) + { + mln_precondition(nrows && ncols); + + // Raw data (impl). + b = make::box2d(2 * nrows - 1, 2 * ncols - 1); + impl = image2d<T>(b, 0); + + // Images. + cell = impl | is_cell; + edge = impl | is_edge; + vertex = impl | is_vertex; + } + + unsigned nrows() const { return (b.nrows() + 1) / 2; } + unsigned ncols() const { return (b.ncols() + 1) / 2; } + + box2d b; + image2d<T> impl; + + typedef image_if< image2d<T>, is_cell_t > cell_t; + cell_t cell; + + typedef image_if< image2d<T>, is_edge_t > edge_t; + edge_t edge; + + typedef image_if< image2d<T>, is_vertex_t > vertex_t; + vertex_t vertex; + }; + + + + + + template <typename C, typename T> + struct cs2d< C, T, T> + { + typedef C cell_value; + typedef T edge_value; + typedef T vertex_value; + + cs2d(unsigned nrows, unsigned ncols) + { + init(nrows, ncols); + } + + cs2d(const box2d& b) + { + init(b.nrows(), b.ncols()); + } + + void init(unsigned nrows, unsigned ncols) + { + mln_precondition(nrows && ncols); + + // Raw data (impl). + b = make::box2d(2 * nrows - 1, 2 * ncols - 1); + cimpl = image2d<C>(b, 0); + impl = image2d<T>(b, 0); + + // Images. + cell = cimpl | is_cell; + edge = impl | is_edge; + vertex = impl | is_vertex; + } + + unsigned nrows() const { return (b.nrows() + 1) / 2; } + unsigned ncols() const { return (b.ncols() + 1) / 2; } + + box2d b; + image2d<C> cimpl; + image2d<T> impl; + + typedef image_if< image2d<C>, is_cell_t > cell_t; + cell_t cell; + + typedef image_if< image2d<T>, is_edge_t > edge_t; + edge_t edge; + + typedef image_if< image2d<T>, is_vertex_t > vertex_t; + vertex_t vertex; + }; + + +} // mln Index: sandbox/geraud/cs2d/cs2d_utils.hh --- sandbox/geraud/cs2d/cs2d_utils.hh (revision 0) +++ sandbox/geraud/cs2d/cs2d_utils.hh (revision 0) @@ -0,0 +1,251 @@ +# include <mln/value/rgb8.hh> +# include <mln/value/int_u8.hh> +# include <mln/literal/black.hh> +# include <mln/literal/white.hh> +# include <mln/literal/colors.hh> +# include <mln/win/rectangle2d.hh> +# include <mln/level/fill.hh> + + +namespace mln +{ + + // 0 1 2 + // 0 c c c + // 1 c c c + + // 0 1 2 3 4 + // 0 c e c e c + // 1 e x e x e + // 2 c e c e c + + + struct is_cell_t : Function_p2b<is_cell_t> + { + typedef bool result; + bool operator()(const point2d& p) const + { + return p.row() % 2 == 0 && p.col() % 2 == 0; + } + } + is_cell; + + + struct is_edge_t : Function_p2b<is_edge_t> + { + typedef bool result; + bool operator()(const point2d& p) const + { + return p.row() % 2 + p.col() % 2 == 1; + } + } + is_edge; + + + struct is_vertex_t : Function_p2b<is_vertex_t> + { + typedef bool result; + bool operator()(const point2d& p) const + { + return p.row() % 2 && p.col() % 2; + } + } + is_vertex; + + + + void display_masks(box2d b) + { + box2d b_(2 * b.nrows() - 1, 2 * b.ncols() - 1); + image2d<char> tmp(b_); + mln_piter_(box2d) p(b_); + for_all(p) + if (is_cell(p)) + { + mln_invariant(! is_edge(p) && ! is_vertex(p)); + tmp(p) = 'O'; + } + else if (is_edge(p)) + { + mln_invariant(! is_cell(p) && ! is_vertex(p)); + tmp(p) = p.to_point().row() % 2 ? '-' : '|'; + } + else if (is_vertex(p)) + { + mln_invariant(! is_edge(p) && ! is_cell(p)); + tmp(p) = '+'; + } + else + mln_invariant(0); + debug::println(tmp); + } + + + struct colorize_t + { + typedef value::rgb8 result; + result operator()(value::int_u8 i) const + { + return result(i,i,i); + } + result operator()(bool b) const + { + if (b) + return literal::red; + else + return literal::black; + } + } + colorize; + + + template <typename Cs, typename F> + image2d<mln_result(F)> pack(const Cs& cs, F f) + { + image2d<mln_result(F)> output(cs.b); + mln_piter(box2d) p(cs.b); + for_all(p) + if (is_cell(p)) + output(p) = f(cs.cell(p)); + else if (is_edge(p)) + output(p) = f(cs.edge(p)); + else + output(p) = f(cs.vertex(p)); + return output; + } + + + template <typename Cs, typename F> + image2d<mln_result(F)> pack(const Cs& cs, F f, unsigned clen) + { + typedef mln_result(F) R; + + mln_precondition(clen % 2); // oddity + + typedef win::rectangle2d Sqr; + Sqr sqr(clen, clen); + + typedef win::hline2d H; + H hline(clen); + + typedef win::vline2d V; + V vline(clen); + + // 0 1 2 3 4 + // 0 c e c e c + // 1 e x e x e + // 2 c e c e c + + // 0 1 2 3 4 5 6 7 8 9 10 + + // 0 c c c . e . c c c . e . c c c + // 1 c c c . e . c c c . e . c c c + // 2 c c c . e . c c c . e . c c c + // 3 . . . . . . . . . . . . . . . + // 4 e e e . x . e e e . x . e e e + // 5 . . . . . . . . . . . . . . . + // 6 c c c . e . c c c . e . c c c + // 7 c c c . e . c c c . e . c c c + // 8 c c c . e . c c c . e . c c c + + box2d b(cs.nrows() * clen + 3 * (cs.nrows() - 1), + cs.ncols() * clen + 3 * (cs.ncols() - 1)); + image2d<R> output(b); + level::fill(output, literal::black); + + // Cells. + mln_piter(Cs::cell_t) c(cs.cell.domain()); + for_all(c) + { + R v = f(cs.cell(c)); + point2d c_ = c.to_point(), p_; + p_.row() = clen / 2 + (clen + 3) * (c_.row() / 2); + p_.col() = clen / 2 + (clen + 3) * (c_.col() / 2); + mln_qiter(Sqr) q(sqr, p_); + for_all(q) + output(q) = v; + } + + // Points. + mln_piter(Cs::vertex_t) x(cs.vertex.domain()); + for_all(x) + { + point2d x_ = x, q; + q.row() = clen + 1 + (clen + 3) * (x_.row() / 2); + q.col() = clen + 1 + (clen + 3) * (x_.col() / 2); + output(q) = f(cs.vertex(x)); + } + + + // 0 1 2 3 4 + // 0 c e c e c + // 1 e x e x e + // 2 c e c e c + + // 0 1 2 3 4 5 6 7 8 9 10 + + // 0 c c c . e . c c c . e . c c c + // 1 c c c . e . c c c . e . c c c + // 2 c c c . e . c c c . e . c c c + // 3 . . . . . . . . . . . . . . . + // 4 e e e . x . e e e . x . e e e + // 5 . . . . . . . . . . . . . . . + // 6 c c c . e . c c c . e . c c c + // 7 c c c . e . c c c . e . c c c + // 8 c c c . e . c c c . e . c c c + + + // Edges. + mln_piter(Cs::edge_t) e(cs.edge.domain()); + for_all(e) + { + R v = f(cs.edge(e)); + point2d e_ = e, p_; + if (e_.row() % 2) // Odd => horizontal. + { + p_.row() = clen + 1 + (clen + 3) * (e_.row() / 2); // vertex-like + p_.col() = clen / 2 + (clen + 3) * (e_.col() / 2); // cell-like + mln_qiter(V) q(hline, p_); + for_all(q) if (output.has(q)) + output(q) = v; + } + else // Even => vertical. + { + p_.row() = clen / 2 + (clen + 3) * (e_.row() / 2); // cell-like + p_.col() = clen + 1 + (clen + 3) * (e_.col() / 2); // vertex-like + mln_qiter(H) q(vline, p_); + for_all(q) if (output.has(q)) + output(q) = v; + } + } + + return output; + } + + + + template <typename Cs, typename T> + void fill_cells_with(Cs& cs, const image2d<T>& input) + { + mlc_equal(T, typename Cs::cell_value)::check(); + mln_piter(box2d) p(input.domain()); + mln_piter(Cs::cell_t) c(cs.cell.domain()); + for_all_2(p, c) + cs.cell(c) = input(p); + } + + template <typename Cs> + image2d<typename Cs::cell_value> + cell_image(const Cs& cs) + { + image2d<typename Cs::cell_value> output(cs.nrows(), cs.ncols()); + mln_piter(box2d) p(output.domain()); + mln_piter(Cs::cell_t) c(cs.cell.domain()); + for_all_2(p, c) + output(p) = cs.cell(c); + return output; + } + +} // mln + + Index: sandbox/geraud/cs2d/dbl_neighb.hh --- sandbox/geraud/cs2d/dbl_neighb.hh (revision 0) +++ sandbox/geraud/cs2d/dbl_neighb.hh (revision 0) @@ -0,0 +1,167 @@ +# include <mln/core/concept/neighborhood.hh> + + +namespace mln +{ + + + // Fwd decl. + template <typename D, typename F> class dbl_niter_; + + + template <typename D, typename F> + struct dbl_neighb_ : public Neighborhood< dbl_neighb_<D,F> > + { + typedef D dpoint; + typedef mln_point(D) point; + + typedef dbl_niter_<D,F> fwd_niter; + typedef dbl_niter_<D,F> bkd_niter; + typedef dbl_niter_<D,F> niter; + + dbl_neighb_(F f) + : f(f) + { + } + + dbl_neighb_() + : f() + { + } + + bool is_empty() const + { + return trues.size() == 0 && falses.size() == 0; + } + + dbl_neighb_& insert_true(const D& dp) + { + trues.push_back(dp); + trues.push_back(-dp); + return *this; + } + + dbl_neighb_& insert_false(const D& dp) + { + falses.push_back(dp); + falses.push_back(-dp); + return *this; + } + + template <typename N> + dbl_neighb_& when_true(const Neighborhood<N>& nbh_) + { + const N& nbh = exact(nbh_); + for (unsigned i = 0; i < nbh.ndpoints(); ++i) + trues.push_back(nbh.dp(i)); + return *this; + } + + template <typename N> + dbl_neighb_& when_false(const Neighborhood<N>& nbh_) + { + const N& nbh = exact(nbh_); + for (unsigned i = 0; i < nbh.ndpoints(); ++i) + falses.push_back(nbh.dp(i)); + return *this; + } + + std::vector<D> trues, falses; + F f; + }; + + + + template <typename D, typename F> + class dbl_niter_ + : public internal::point_iterator_base_< mln_point(D), dbl_niter_<D,F> > + { + public: + + typedef mln_point(D) point; + + template <typename Pref> + dbl_niter_(const Neighborhood< dbl_neighb_<D,F> >& nbh, + const Point_Site<Pref>& p_ref) : + dps_true(exact(nbh).trues), + dps_false(exact(nbh).falses), + p_ref_(exact(p_ref).to_point()), + f(exact(nbh).f) + { + invalidate(); + } + + /// Convertion to point. + operator point () const + { + mln_precondition(is_valid()); + return p_; + } + + /// Reference to the corresponding point. + const point& to_point() const + { + return p_; + } + + /// Test the iterator validity. + bool is_valid() const + { + return i_ != dps_true.size(); + } + + /// Invalidate the iterator. + void invalidate() + { + i_ = dps_true.size(); + } + + /// Start an iteration. + void start() + { + i_ = 0; + update(); + } + + /// Go to the next point. + void next_() + { + ++i_; + update(); + } + + /// Give the i-th coordinate. + int operator[](unsigned i) const + { + mln_precondition(is_valid()); + return p_[i]; + } + + /// The point around which this iterator moves. + const point& center_point() const + { + return p_ref_; + } + + /// Force this iterator to update its location to take into + /// account that its center point may have moved. + void update() + { + if (is_valid()) + p_ = p_ref_ + (f(p_ref_) ? dps_true[i_] : dps_false[i_]); + } + + protected: + + const std::vector<D>& dps_true, dps_false; + const point& p_ref_; // reference point (or "center point") + F f; + + unsigned i_; + point p_; // location of this iterator; p_ makes this iterator be + // itself a potential center point. + }; + + + +} // mln Index: sandbox/geraud/cs2d/cs2d_morpho.hh --- sandbox/geraud/cs2d/cs2d_morpho.hh (revision 0) +++ sandbox/geraud/cs2d/cs2d_morpho.hh (revision 0) @@ -0,0 +1,372 @@ +# include <mln/accu/max.hh> +# include <mln/accu/min_max.hh> +# include <mln/labeling/regional_minima.hh> + +# include <mln/core/p_queue_fast.hh> + + +namespace mln +{ + + namespace cmorpho + { + + template <typename T> + T abs_diff(T t1, T t2) + { + return t1 > t2 ? t1 - t2 : t2 - t1; + } + + template <typename I> + void gradient(I& input) + { + mln_piter(I) p(input.domain()); + for_all(p) + { + mln_point(I) p_ = p; + input(p) = + p_.row() % 2 ? + // Odd so horizontal edge and vertical cells. + abs_diff(input(p_ + left), input(p_ + right)) : + // Otherwise + abs_diff(input(p_ + up), input(p_ + down)); + } + } + + + template <typename I, typename N> + inline + mln_concrete(I) gradient(const Image<I>& input_, + const Neighborhood<N>& nbh_) + { + const I& input = exact(input_); + const N& nbh = exact(nbh_); + + mln_concrete(I) output; + initialize(output, input); + + accu::min_max_<mln_value(I)> m; + + mln_piter(I) p(input.domain()); + mln_niter(N) n(nbh, p); + for_all(p) + { + m.init(); + for_all(n) if (input.owns_(n)) + m.take(input(n)); + output(p) = m.second() - m.first(); + } + return output; + } + + + class dbl_niter; + + struct dbl_neighb2d : public Neighborhood< dbl_neighb2d > + { + typedef dpoint2d dpoint; + typedef point2d point; + + typedef dbl_niter fwd_niter; + typedef dbl_niter bkd_niter; + typedef dbl_niter niter; + + dbl_neighb2d() + { + } + + dbl_neighb2d& insert_odd(const dpoint2d& dp) + { + odds.push_back(dp); + odds.push_back(-dp); + return *this; + } + dbl_neighb2d& insert_even(const dpoint2d& dp) + { + evens.push_back(dp); + evens.push_back(-dp); + return *this; + } + + std::vector<dpoint2d> odds, evens; + }; + + + class dbl_niter + : public internal::point_iterator_base_< point2d, dbl_niter > + { + public: + + typedef dpoint2d D; + + template <typename N, typename Pref> + dbl_niter(const N& nbh, const Point_Site<Pref>& p_ref) : + dps_odd(exact(nbh).odds), + dps_even(exact(nbh).evens), + p_ref_(exact(p_ref).to_point()) + { + + invalidate(); + } + + /// Convertion to point. + operator point2d () const + { + mln_precondition(is_valid()); + return p_; + } + + /// Reference to the corresponding point. + const point2d& to_point() const + { + return p_; + } + + /// Test the iterator validity. + bool is_valid() const + { + return i_ != dps_odd.size(); + } + + /// Invalidate the iterator. + void invalidate() + { + i_ = dps_odd.size(); + } + + /// Start an iteration. + void start() + { + i_ = 0; + update(); + } + + /// Go to the next point. + void next_() + { + ++i_; + update(); + } + + /// Give the i-th coordinate. + int operator[](unsigned i) const + { + mln_precondition(is_valid()); + return p_[i]; + } + + /// The point around which this iterator moves. + const point2d& center_point() const + { + return p_ref_; + } + + /// Force this iterator to update its location to take into + /// account that its center point may have moved. + void update() + { + if (is_valid()) + p_ = p_ref_ + (p_ref_.row() % 2 ? dps_odd[i_] : dps_even[i_]); + } + + protected: + + const std::vector<dpoint2d>& dps_odd, dps_even; + const point2d& p_ref_; // reference point (or "center point") + + unsigned i_; + point2d p_; // location of this iterator; p_ makes this iterator be + // itself a potential center point. + }; + + + template <typename L, typename I, typename N> + mln_ch_value(I, L) + meyer_wst(const Image<I>& input_, const Neighborhood<N>& nbh_, + L& nbasins) + { + trace::entering("morpho::meyer_wst"); + + const I& input = exact(input_); + const N& nbh = exact(nbh_); + + /* FIXME: Ensure the input image has scalar values. */ + + const L unmarked = 0; + + // Initialize the output with the markers (minima components). + mln_ch_value(I, L) output; + output = labeling::regional_minima(input, nbh, nbasins); + + // Ordered queue. + std::vector< p_queue_fast<mln_psite(I)> > qu(256); + unsigned qu_size = 0; + + // 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.has(n) && output(n) != unmarked) + { + qu[input(p)].push_force(p); + ++qu_size; + break; + } + } + + /* Until the queue is empty, extract a psite p from the + hierarchical queue, at the highest priority level, that is, + the lowest level. */ + mln_psite(I) p; + mln_niter(N) n(nbh, p); + + while (qu_size != 0) + { + + for (unsigned l = 0; l < 256; ++l) + if (! qu[l].is_empty()) + { + p = qu[l].pop_front(); + // qu[l].pop(); + break; + } + --qu_size; + + // Last seen marker adjacent to P. + L adjacent_marker = unmarked; + // Has P a single adjacent marker? + bool single_adjacent_marker_p = true; + for_all(n) + if (output.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.has(n) && output(n) == unmarked) + { + qu[input(n)].push_force(n); + ++qu_size; + } + } + } + + trace::exiting("morpho::meyer_wst"); + return output; + } + + } // mln::cmorpho + + + namespace accu + { + + template <typename A_, typename I, typename L> + inline + std::vector<mln_accu_with(A_, mln_value(I))::result> + compute(const Image<I>& input_, + const Image<L>& label_, + mln_value(L) n) + { + mlc_is_a(A_, Meta_Accumulator)::check(); + + trace::entering("accu::compute"); + + const I& input = exact(input_); + const L& label = exact(label_); + + typedef mln_accu_with(A_, mln_value(I)) A; + std::vector<A> a(n + 1); + + mln_piter(I) p(input.domain()); + for_all(p) + a[label(p)].take(input(p)); + + typedef typename A::result R; + std::vector<R> v(n + 1); + for (unsigned l = 1; l <= n; ++l) + v[l] = a[l].to_result(); + + trace::exiting("accu::compute"); + return v; + } + + + template <typename A_, typename I, typename L, typename R> + inline + void + compute(const Image<I>& input_, + const Image<L>& label_, + std::vector<R>& v) + { + mlc_is_a(A_, Meta_Accumulator)::check(); + + trace::entering("accu::compute"); + + const I& input = exact(input_); + const L& label = exact(label_); + + typedef mln_accu_with(A_, mln_value(I)) A; + const unsigned n = v.size(); + std::vector<A> a(n); + + mln_piter(I) p(input.domain()); + for_all(p) + a[label(p)].take(input(p)); + + for (unsigned l = 1; l < n; ++l) + v[l] = a[l].to_result(); + + trace::exiting("accu::compute"); + } + + } // mln::accu + + + namespace level + { + + template<typename I, typename R> + mln_ch_value(I, R) + transform(const Image<I>& input_, const std::vector<R>& v) + { + trace::entering("level::transform"); + const I& input = exact(input_); + + mln_ch_value(I, R) output; + initialize(output, input); + + mln_piter(I) p(input.domain()); + for_all(p) + output(p) = v[input(p)]; + + trace::exiting("level::transform"); + return output; + } + + } // mln::level + + +} // mln