[Olena-patches] 3026: Add some experimental code about TUFA and regional minima.

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; Add some experimental code about TUFA and regional minima. * geraud/tufa_2008: New directory. * geraud/tufa_2008/steps.0.cc: New. * geraud/tufa_2008/compute_a.cc: New. * geraud/tufa_2008/steps.1.cc: New. * geraud/tufa_2008/steps.2.cc: New. * geraud/tufa_2008/closing.cc: New. * geraud/tufa_2008/steps.3.cc: New. * geraud/tufa_2008/opening.cc: New. * geraud/tufa_2008/experiment.cc: New. * geraud/tufa_2008/steps.2b.cc: New. * geraud/tufa_2008/filter.cc: New. * geraud/tufa_2008/wst_f_equal_wst_a.cc: New. * geraud/tufa_2008/closed_gradient.cc: New. * geraud/tufa_2008/fz_count.cc: New. * geraud/tufa_2008/regmin_count.cc: New. closed_gradient.cc | 64 +++++++ closing.cc | 64 +++++++ compute_a.cc | 411 +++++++++++++++++++++++++++++++++++++++++++++++++++ experiment.cc | 138 +++++++++++++++++ filter.cc | 169 ++++++++++++++++++++ fz_count.cc | 159 +++++++++++++++++++ opening.cc | 63 +++++++ regmin_count.cc | 173 +++++++++++++++++++++ steps.0.cc | 170 +++++++++++++++++++++ steps.1.cc | 273 +++++++++++++++++++++++++++++++++ steps.2.cc | 184 ++++++++++++++++++++++ steps.2b.cc | 184 ++++++++++++++++++++++ steps.3.cc | 184 ++++++++++++++++++++++ wst_f_equal_wst_a.cc | 140 +++++++++++++++++ 14 files changed, 2376 insertions(+) Index: geraud/tufa_2008/steps.0.cc --- geraud/tufa_2008/steps.0.cc (revision 0) +++ geraud/tufa_2008/steps.0.cc (revision 0) @@ -0,0 +1,170 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/steps.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/morpho/tree/data.hh> + +#include <mln/accu/count.hh> + +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/morpho/closing_area.hh> +#include <mln/level/fill.hh> + + +namespace mln +{ + + template <typename I> + void println_par(const I& par) + { + int nr = par.nrows(), nc = par.ncols(); + for (int r = 0; r < nr; ++r) + { + for (int c = 0; c < nc; ++c) + if (par.at(r,c) == point2d(r,c)) + std::cout << "( ) "; + else + std::cout << par.at(r,c) << ' '; + std::cout << std::endl; + } + } + + template <typename P> + inline + mln_value(P) find_root__(P& par, const mln_value(P)& x) + { + if (par(x) == x) + return x; + else + return par(x) = find_root__(par, par(x)); + } + + + template <typename I, typename A, typename N> + mln_concrete(I) filtering(const I& f, const A& a, const N& nbh, mln_value(A) lambda) + { + typedef p_array<mln_psite(I)> S; + S s = level::sort_psites_increasing(a); + + // s maps increasing attributes. + + mln_concrete(I) out; + initialize(out, f); + + mln_ch_value(I, mln_site(I)) par; + mln_ch_value(I, bool) deja_vu; + { + initialize(par, f); + mln_piter(A) p(par.domain()); + for_all(p) + par(p) = p; + initialize(deja_vu, f); + level::fill(deja_vu, false); + } + + // First pass. + { + mln_site(I) r; + mln_fwd_piter(S) p(s); + mln_niter(N) n(nbh, p); + for_all(p) + { + for_all(n) + if (a.domain().has(n) && deja_vu(n)) + { + r = find_root__(par, n); + if (r != p) + if (f(r) == f(p) || a(r) < lambda) + par(r) = p; // Union. + } + deja_vu(p) = true; + } + } + + // Second pass. + { + mln_bkd_piter(S) p(s); + for_all(p) + if (par(p) == p) + out(p) = f(p); + else + out(p) = out(par(p)); + } + return out; + } + + +} // mln + + + + +int main() +{ + using namespace mln; + using value::int_u8; + + int_u8 n; + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, "../../../img/fly.pgm"); + debug::println("f =", f); + + debug::println("ref =", morpho::closing_area(f, c4(), 10)); + + typedef p_array<point2d> S; + S s = level::sort_psites_decreasing(f); + + // Children go towards lower levels so leafs are regional minima. + // We get a min-tree so that we can perform morphological closings. + + morpho::tree::data<I,S> t(f, s, c4()); + accu::count< util::pix<I> > attr; + + image2d<unsigned> a = morpho::tree::compute_attribute_image(attr, t); + I g = filtering(f, a, c4(), 10); + debug::println("g =", g); + + { + S s = level::sort_psites_decreasing(g); + morpho::tree::data<I,S> t(g, s, c4()); + image2d<unsigned> a_g = morpho::tree::compute_attribute_image(attr, t); + debug::println("a(f) =", a); + debug::println("a(g) =", a_g); + } + +} Index: geraud/tufa_2008/compute_a.cc --- geraud/tufa_2008/compute_a.cc (revision 0) +++ geraud/tufa_2008/compute_a.cc (revision 0) @@ -0,0 +1,411 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/compute_a.cc + +#include <mln/core/image/image2d.hh> +#include <mln/core/image/image_if.hh> +#include <mln/pw/all.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/level/fill.hh> +#include <mln/level/paste.hh> +#include <mln/level/compare.hh> + +#include <mln/morpho/tree/data.hh> +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/labeling/regional_minima.hh> + +#include <mln/accu/count.hh> + + + +namespace mln +{ + + template <typename I> + void println_par(const I& par) + { + int nr = par.nrows(), nc = par.ncols(); + for (int r = 0; r < nr; ++r) + { + for (int c = 0; c < nc; ++c) + if (par.at(r,c) == point2d(r,c)) + std::cout << "( ) "; + else + std::cout << par.at(r,c) << ' '; + std::cout << std::endl; + } + } + + + template <typename P> + inline + mln_value(P) find_root__(P& par, const mln_value(P)& x) + { + if (par(x) == x) + return x; + else + return par(x) = find_root__(par, par(x)); + } + + + + //------------------------------- compute_a + + + + template <typename I, typename N, typename A> + mln_ch_value(I, mln_result(A)) + compute_a(const I& f, const N& nbh, A, unsigned& n_regmins) + { + typedef p_array<mln_psite(I)> S; + S s = level::sort_psites_increasing(f); + // s maps increasing attributes. + + mln_ch_value(I, mln_site(I)) par; + mln_ch_value(I, mln_site(I)) zpar; + mln_ch_value(I, bool) deja_vu, flag; + mln_ch_value(I, A) attr; + + n_regmins = f.domain().nsites(); + + + // Initialization. + { + // parent + initialize(par, f); + initialize(zpar, f); + + // deja_vu + initialize(deja_vu, f); + level::fill(deja_vu, false); + + // flag + initialize(flag, f); + level::fill(flag, true); + + // attr + initialize(attr, f); + } + + + // First Pass. + { + mln_site(I) r; + mln_fwd_piter(S) p(s); + mln_niter(N) n(nbh, p); + for_all(p) + { + // Make-Set. + par(p) = p; + zpar(p) = p; + attr(p).take_as_init(p); + + for_all(n) + if (f.domain().has(n) && deja_vu(n)) + { + r = find_root__(zpar, n); + if (r != p) + { + // Fully compressed union. + zpar(r) = p; + attr(p).take(attr(r)); + + if (f(r) == f(p)) + { + // Weak-Union; only for flat zones. + par(r) = p; + + if (flag(p) == false && flag(r) == false) + { + // Two non-reg-min components merge (same flat + // zone) so we had an extra invalidation. + ++n_regmins; + } + flag(p) = flag(p) && flag(r); + --n_regmins; // So we get the number of flat zones + // minus the non-reg-min flat zones. + } + else + { + mln_invariant(f(r) < f(p)); + if (flag(p) == true) + --n_regmins; // Invalidation. + flag(p) = false; + } + } + } + deja_vu(p) = true; + } + } // end of First Pass. + + + std::cout << "n reg min = " << n_regmins << std::endl; + + + { + unsigned n = 0; + mln_fwd_piter(S) p(s); + for_all(p) + if (par(p) == p && flag(p)) + ++n; + mln_assertion(n == n_regmins); + } + + + // The attr image is not correct on flat zones since there is + // no back-propagation of the attribute value of the component + // root. For instance with f="v v v" we get attr="1 2 3" + // instead of "3 3 3". So a finalization is required. + + mln_ch_value(I, mln_result(A)) a; + initialize(a, f); + level::paste(attr, a); + + // Finalization. + { + mln_bkd_piter(S) p(s); // Reverse. + + unsigned n_non_compressed_par = 0; + for_all(p) + { + a(p) = a(par(p)); + if (par(par(p)) != par(p)) + ++ n_non_compressed_par; + } + std::cout << "n_non_compressed_par = " << n_non_compressed_par << std::endl; + } + + // TODO: compress at least the reg minima! + + + { + image2d<unsigned> regmin; + initialize(regmin, f); + { + unsigned i_regmin = 0; + mln_bkd_piter(S) p(s); + for_all(p) + { + if (par(p) == p) + { + if (flag(p)) + regmin(p) = ++i_regmin; + else + regmin(p) = 0; + } + else + regmin(p) = regmin(par(p)); + } + } + + debug::println("f", f); + + debug::println("flag", flag); + // We can see that some point are at true for components that + // are not reg min; flag is a candidate to be compressed... + + debug::println("regmin", regmin); + + + + // TODO: + + // On veut tester ici dans quel ordre on voit les + // reg min lorsque a croit. Pour tous les points d'un reg min, + // est-ce que le root est vu en premier ? + + + + image2d<bool> seen; + initialize(seen, f); + level::fill(seen, false); + + s = level::sort_psites_increasing(a); + mln_bkd_piter(S) p(s); + for_all(p) + { + if (regmin(p) == 0) + continue; + // p is in a regional minimum. + if (par(p) != p) // A non-root point. + { + mln_assertion(regmin(par(p)) != 0); // Root in a regional minimum. + mln_assertion(regmin(par(p)) == regmin(p)); // and the same as p. + mln_assertion(seen(par(p))); + } + seen(p) = true; + } + debug::println(seen); + +// if (flag(p)) +// std::cout << a(p) << ' ' << p << ' ' << (par(p) == p) << std::endl; + + } + + return a; + } + + + + //------------------------------- filtering + + + +// template <typename I, typename A, typename N> +// mln_concrete(I) filtering(const I& f, const A& a, const N& nbh, +// unsigned n_regmins, unsigned n_wanted) +// { +// typedef p_array<mln_psite(I)> S; +// S s = level::sort_psites_increasing(a); + +// // s maps increasing attributes. + +// mln_concrete(I) out; +// initialize(out, f); + +// mln_ch_value(I, mln_site(I)) par; +// mln_ch_value(I, bool) deja_vu, flag; + +// // Initialization. +// { +// initialize(par, f); +// mln_piter(A) p(par.domain()); +// for_all(p) +// par(p) = p; +// initialize(deja_vu, f); +// level::fill(deja_vu, false); + +// // flag +// initialize(flag, f); +// level::fill(flag, true); +// } + + +// int counter = 0; +// // We are trying to count the number of merges of regional minima... + + +// // First Pass. +// { +// mln_site(I) r; +// mln_fwd_piter(S) p(s); +// mln_niter(N) n(nbh, p); +// for_all(p) +// { +// for_all(n) +// if (a.domain().has(n) && deja_vu(n)) +// { +// r = find_root__(par, n); +// if (r != p) +// if (a(r) == a(p)) +// { +// par(r) = p; // Union. +// if (flag(r) == true && flag(p) == true) +// --counter; +// flag(p) = flag(p) && flag(r); +// } +// else // a(r) != a(p) +// { +// if (flag(r) == true && flag(p) == true) +// ++counter; +// mln_invariant(a(p) > a(r)); +// flag(p) = false; +// } +// } +// deja_vu(p) = true; +// } +// std::cout << counter << std::endl; +// } + +// // // Second Pass. +// // { +// // mln_bkd_piter(S) p(s); +// // for_all(p) +// // if (par(p) == p) +// // out(p) = f(p); +// // else +// // out(p) = out(par(p)); +// // } + +// return out; +// } + + + + +} // mln + + + + +int main(int, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, argv[1]); + // debug::println(f); + + accu::count<point2d> area; + unsigned n_regmins; + image2d<unsigned> a = compute_a(f, c4(), area, n_regmins); + // debug::println(a); + +// { +// // Test of 'n_regmins'. +// unsigned ref; +// labeling::regional_minima(f, c4(), ref); +// mln_assertion(n_regmins == ref); +// } + +// { +// // Test of 'a'. +// typedef p_array<point2d> S; +// S s = level::sort_psites_decreasing(f); +// morpho::tree::data<I,S> t(f, s, c4()); +// accu::count< util::pix<I> > area; +// image2d<unsigned> ref = morpho::tree::compute_attribute_image(area, t); +// mln_assertion(a == ref); +// } + + + + + +// unsigned n_wanted = 10; +// I g = filtering(f, a, c4(), n_regmins, n_wanted); + +} Index: geraud/tufa_2008/steps.1.cc --- geraud/tufa_2008/steps.1.cc (revision 0) +++ geraud/tufa_2008/steps.1.cc (revision 0) @@ -0,0 +1,273 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/steps.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/morpho/tree/data.hh> + +#include <mln/accu/count.hh> +#include <mln/util/set.hh> + +#include <mln/labeling/regional_minima.hh> +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/morpho/closing_area.hh> +#include <mln/level/fill.hh> + + +namespace mln +{ + + template <typename I> + void println_par(const I& par) + { + int nr = par.nrows(), nc = par.ncols(); + for (int r = 0; r < nr; ++r) + { + for (int c = 0; c < nc; ++c) + if (par.at(r,c) == point2d(r,c)) + std::cout << "( ) "; + else + std::cout << par.at(r,c) << ' '; + std::cout << std::endl; + } + } + + template <typename P> + inline + mln_value(P) find_root__(P& par, const mln_value(P)& x) + { + if (par(x) == x) + return x; + else + return par(x) = find_root__(par, par(x)); + } + + + template <typename I, typename A, typename N> + void run_run(const I& f, const A& a, const N& nbh) + { + typedef p_array<mln_psite(I)> S; + S s = level::sort_psites_increasing(a); + // s maps increasing attributes. + + mln_ch_value(I, mln_site(I)) par; + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, util::set<unsigned>) labels; + unsigned nbassins, current_n; + + + // Initialization. + { + mln_piter(A) p(f.domain()); + + // parent + initialize(par, f); + for_all(p) + par(p) = p; + + // deja_vu + initialize(deja_vu, f); + level::fill(deja_vu, false); + + // labels + mln_ch_value(I, unsigned) regmin = labeling::regional_minima(a, nbh, + nbassins); + debug::println(regmin); + initialize(labels, f); + for_all(p) + if (regmin(p) != 0) // p in a reg min of the attribute image + labels(p).insert(regmin(p)); + } + + current_n = nbassins; + + // First pass. + { + mln_site(I) r; + mln_fwd_piter(S) p(s); + mln_niter(N) n(nbh, p); + for_all(p) + { + for_all(n) + if (a.domain().has(n) && deja_vu(n)) + { + r = find_root__(par, n); + if (r != p) + { + par(r) = p; // Union. + + // logging the different cases + if (labels(p).is_empty() && labels(r).is_empty()) + { + std::cout << "x"; + + // It can happen with: + // M M m + // (m a min, M a flat zone) + // During the pass: + // M M {m} + // then the 1st point is processed + // {} M {m} + // and when processing the 2nd point + // r={} p={} {m} + // the (left) neighbor of p has an empty set such as p. + +// // Extra log. +// std::cout << std::endl +// << "p = " << p << " " +// << "r = " << r << std::endl; +// debug::println(labels); + } + else + if (labels(p).is_empty() || labels(r).is_empty()) + { + std::cout << (labels(p).is_empty() ? 'p' : 'r'); + } + else + if (labels(p) == labels(r)) + { + std::cout << "="; +// // Extra log. +// std::cout << std::endl +// << "p = " << p << " " +// << "r = " << r << std::endl; +// debug::println(labels); + } + else + std::cout << "."; + // end of log + + // The invariants below are erroneous. +// mln_invariant(! (labels(p).is_empty() && labels(r).is_empty())); +// mln_invariant(labels(p) != labels(r)); + + // Either: + // one of the two label sets is empty (and the other is not) + // or: + // both label sets are not empty then they differ. + + // More restrictively: +// mln_invariant(! labels(r).is_empty()); + + if (labels(p).is_empty()) + { + labels(p).insert(labels(r)); + } + else + { + +// std::cout << std::endl +// << "at " << p << std::endl +// << "before:" << std::endl; +// debug::println(labels); + + unsigned + np = labels(p).nelements(), + nr = labels(r).nelements(); + labels(p).insert(labels(r)); + unsigned + n = labels(p).nelements(), + dnp = n - np, + dnr = n - nr, + delta_n = std::min(dnp, dnr); + + // The invariant below is erroneous. +// mln_invariant(delta_n != 0); + + +// std::cout << "delta = " << delta_n << std::endl +// << "after: " << std::endl; +// debug::println(labels); + + + // We can have the three cases below: + if (dnr > dnp) + std::cout << '>'; + else + if (dnr < dnp) + std::cout << '<'; + else + std::cout << '~'; + std::cout << '(' << delta_n << ')'; + } + } + } + deja_vu(p) = true; + } + } + + std::cout << std::endl; +// debug::println(labels); + } + + +} // mln + + + + +int main() +{ + using namespace mln; + using value::int_u8; + + int_u8 n; + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, "../../../img/tiny.pgm"); +// debug::println("f =", f); + + typedef p_array<point2d> S; + S s = level::sort_psites_decreasing(f); + + // Children go towards lower levels so leafs are regional minima. + // We get a min-tree so that we can perform morphological closings. + + morpho::tree::data<I,S> t(f, s, c4()); + accu::count< util::pix<I> > attr; + + image2d<unsigned> a = morpho::tree::compute_attribute_image(attr, t); + run_run(f, a, c4()); + +// { +// S s = level::sort_psites_decreasing(g); +// morpho::tree::data<I,S> t(g, s, c4()); +// image2d<unsigned> a_g = morpho::tree::compute_attribute_image(attr, t); +// debug::println("a(f) =", a); +// debug::println("a(g) =", a_g); +// } + +} Index: geraud/tufa_2008/steps.2.cc --- geraud/tufa_2008/steps.2.cc (revision 0) +++ geraud/tufa_2008/steps.2.cc (revision 0) @@ -0,0 +1,184 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/steps.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/morpho/tree/data.hh> + +#include <mln/accu/count.hh> +#include <mln/util/set.hh> + +#include <mln/labeling/regional_minima.hh> +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/morpho/closing_area.hh> +#include <mln/level/fill.hh> + + +namespace mln +{ + + template <typename I> + void println_par(const I& par) + { + int nr = par.nrows(), nc = par.ncols(); + for (int r = 0; r < nr; ++r) + { + for (int c = 0; c < nc; ++c) + if (par.at(r,c) == point2d(r,c)) + std::cout << "( ) "; + else + std::cout << par.at(r,c) << ' '; + std::cout << std::endl; + } + } + + template <typename P> + inline + mln_value(P) find_root__(P& par, const mln_value(P)& x) + { + if (par(x) == x) + return x; + else + return par(x) = find_root__(par, par(x)); + } + + + template <typename I, typename A, typename N> + void run_run(const I& f, const A& a, const N& nbh) + { + typedef p_array<mln_psite(I)> S; + S s = level::sort_psites_increasing(a); + // s maps increasing attributes. + + mln_ch_value(I, mln_site(I)) par; + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, util::set<unsigned>) labels; + unsigned nbassins, current_n; + + + // Initialization. + { + mln_piter(A) p(f.domain()); + + // parent + initialize(par, f); + for_all(p) + par(p) = p; + + // deja_vu + initialize(deja_vu, f); + level::fill(deja_vu, false); + + // labels + mln_ch_value(I, unsigned) regmin = labeling::regional_minima(a, nbh, + nbassins); + initialize(labels, f); + for_all(p) + if (regmin(p) != 0) // p in a reg min of the attribute image + labels(p).insert(regmin(p)); + } + + current_n = nbassins; + + // First pass. + { + mln_site(I) r; + mln_fwd_piter(S) p(s); + mln_niter(N) n(nbh, p); + for_all(p) + { + for_all(n) + if (a.domain().has(n) && deja_vu(n)) + { + r = find_root__(par, n); + if (r != p) + { + par(r) = p; // Union. + + if (labels(r).is_empty()) + // No-op. + ; + else + if (labels(p).is_empty()) + labels(p) = labels(r); + else + if (labels(p) == labels(r)) + // No-op. + ; + else + { + labels(p).insert(labels(r)); + --current_n; + } + } + } + deja_vu(p) = true; + } + } + + std::cout << std::endl; + std::cout << "end = " << current_n << std::endl; + } + + +} // mln + + + + +int main(int, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + int_u8 n; + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, argv[1]); + + typedef p_array<point2d> S; + S s = level::sort_psites_decreasing(f); + + // Children go towards lower levels so leafs are regional minima. + // We get a min-tree so that we can perform morphological closings. + + morpho::tree::data<I,S> t(f, s, c4()); + accu::count< util::pix<I> > attr; + + image2d<unsigned> a = morpho::tree::compute_attribute_image(attr, t); + run_run(f, a, c4()); + +} Index: geraud/tufa_2008/closing.cc --- geraud/tufa_2008/closing.cc (revision 0) +++ geraud/tufa_2008/closing.cc (revision 0) @@ -0,0 +1,64 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/closing.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/core/alias/neighb2d.hh> + +#include <mln/morpho/elementary/gradient.hh> +#include <mln/morpho/closing_volume.hh> + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm lambda output.pgm" << std::endl; + std::abort(); +} + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + if (argc != 4) + usage(argv); + + image2d<int_u8> f; + io::pgm::load(f, argv[1]); + unsigned lambda = std::atoi(argv[2]); + io::pgm::save(morpho::closing_volume(f, + c4(), + lambda), + argv[3]); +} Index: geraud/tufa_2008/steps.3.cc --- geraud/tufa_2008/steps.3.cc (revision 0) +++ geraud/tufa_2008/steps.3.cc (revision 0) @@ -0,0 +1,184 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/steps.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/morpho/tree/data.hh> + +#include <mln/accu/count.hh> +#include <mln/util/set.hh> + +#include <mln/labeling/regional_minima.hh> +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/morpho/closing_area.hh> +#include <mln/level/fill.hh> + + +namespace mln +{ + + template <typename I> + void println_par(const I& par) + { + int nr = par.nrows(), nc = par.ncols(); + for (int r = 0; r < nr; ++r) + { + for (int c = 0; c < nc; ++c) + if (par.at(r,c) == point2d(r,c)) + std::cout << "( ) "; + else + std::cout << par.at(r,c) << ' '; + std::cout << std::endl; + } + } + + template <typename P> + inline + mln_value(P) find_root__(P& par, const mln_value(P)& x) + { + if (par(x) == x) + return x; + else + return par(x) = find_root__(par, par(x)); + } + + + template <typename I, typename A, typename N> + void run_run(const I& f, const A& a, const N& nbh) + { + typedef p_array<mln_psite(I)> S; + S s = level::sort_psites_increasing(a); + // s maps increasing attributes. + + mln_ch_value(I, mln_site(I)) par; + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, util::set<unsigned>) labels; + unsigned nbassins, current_n; + + + // Initialization. + { + mln_piter(A) p(f.domain()); + + // parent + initialize(par, f); + for_all(p) + par(p) = p; + + // deja_vu + initialize(deja_vu, f); + level::fill(deja_vu, false); + + // labels + mln_ch_value(I, unsigned) regmin = labeling::regional_minima(a, nbh, + nbassins); + initialize(labels, f); + for_all(p) + if (regmin(p) != 0) // p in a reg min of the attribute image + labels(p).insert(regmin(p)); + } + + current_n = nbassins; + + // First pass. + { + mln_site(I) r; + mln_fwd_piter(S) p(s); + mln_niter(N) n(nbh, p); + for_all(p) + { + for_all(n) + if (a.domain().has(n) && deja_vu(n)) + { + r = find_root__(par, n); + if (r != p) + { + par(r) = p; // Union. + + if (labels(r).is_empty()) + // No-op. + ; + else + if (labels(p).is_empty()) + labels(p) = labels(r); + else + if (labels(p) == labels(r)) + // No-op. + ; + else + { + labels(p).insert(labels(r)); + --current_n; + } + } + } + deja_vu(p) = true; + } + } + + std::cout << std::endl; + std::cout << "end = " << current_n << std::endl; + } + + +} // mln + + + + +int main(int, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + int_u8 n; + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, argv[1]); + + typedef p_array<point2d> S; + S s = level::sort_psites_decreasing(f); + + // Children go towards lower levels so leafs are regional minima. + // We get a min-tree so that we can perform morphological closings. + + morpho::tree::data<I,S> t(f, s, c4()); + accu::count< util::pix<I> > attr; + + image2d<unsigned> a = morpho::tree::compute_attribute_image(attr, t); + run_run(f, a, c4()); + +} Index: geraud/tufa_2008/opening.cc --- geraud/tufa_2008/opening.cc (revision 0) +++ geraud/tufa_2008/opening.cc (revision 0) @@ -0,0 +1,63 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/opening.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/core/alias/neighb2d.hh> + +#include <mln/morpho/opening_volume.hh> + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm lambda output.pgm" << std::endl; + std::abort(); +} + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + if (argc != 4) + usage(argv); + + image2d<int_u8> f, g; + io::pgm::load(f, argv[1]); + unsigned lambda = std::atoi(argv[2]); + initialize(g, f); + morpho::opening_volume(f, c4(), lambda, g); + io::pgm::save(g, + argv[3]); +} Index: geraud/tufa_2008/experiment.cc --- geraud/tufa_2008/experiment.cc (revision 0) +++ geraud/tufa_2008/experiment.cc (revision 0) @@ -0,0 +1,138 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/experiment.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> +#include <mln/literal/black.hh> + +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/io/ppm/save.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/morpho/tree/data.hh> + +#include <mln/accu/volume.hh> +#include <mln/win/disk2d.hh> + +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/morpho/meyer_wst.hh> +#include <mln/morpho/opening.hh> +#include <mln/morpho/closing_area.hh> + +#include <mln/level/fill.hh> + +#include <mln/level/transform.hh> +#include <mln/level/stretch.hh> + + + +namespace mln +{ + + struct colorize : Function_v2v< colorize > + { + typedef value::rgb8 result; + colorize(unsigned max) + : lut(max + 1) + { + lut[0] = 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; + }; + +} // mln + + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl; + std::abort(); +} + + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + if (argc != 2) + usage(argv); + + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, argv[1]); + + typedef p_array<point2d> S; + S s = level::sort_psites_decreasing(f); + + // Children go towards lower levels so leafs are regional minima. + // We get a min-tree so that we can perform morphological closings. + + morpho::tree::data<I,S> t(f, s, c4()); + accu::volume<I> attr; + + image2d<unsigned> a = morpho::tree::compute_attribute_image(attr, t); + io::pgm::save(level::stretch(int_u8(), a), + "a.pgm"); + + unsigned n; + image2d<unsigned> wst_a = morpho::meyer_wst(a, c4(), n); + io::ppm::save(level::transform(wst_a, colorize(n)), + "wst_a.ppm"); + std::cout << "n(a) = " << n << std::endl; + + // FIXME: ça n'a pas de sens de faire ce qui est dessous... :-( + + + image2d<unsigned> aa = morpho::closing_area(a, c4(), 100); + io::pgm::save(level::stretch(int_u8(), aa), + "aa.pgm"); + + image2d<unsigned> wst_aa = morpho::meyer_wst(aa, c4(), n); + io::ppm::save(level::transform(wst_aa, colorize(n)), + "wst_aa.ppm"); + std::cout << "n(aa) = " << n << std::endl; +} Index: geraud/tufa_2008/steps.2b.cc --- geraud/tufa_2008/steps.2b.cc (revision 0) +++ geraud/tufa_2008/steps.2b.cc (revision 0) @@ -0,0 +1,184 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/steps.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/morpho/tree/data.hh> + +#include <mln/accu/count.hh> +#include <set> + +#include <mln/labeling/regional_minima.hh> +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/morpho/closing_area.hh> +#include <mln/level/fill.hh> + + +namespace mln +{ + + template <typename I> + void println_par(const I& par) + { + int nr = par.nrows(), nc = par.ncols(); + for (int r = 0; r < nr; ++r) + { + for (int c = 0; c < nc; ++c) + if (par.at(r,c) == point2d(r,c)) + std::cout << "( ) "; + else + std::cout << par.at(r,c) << ' '; + std::cout << std::endl; + } + } + + template <typename P> + inline + mln_value(P) find_root__(P& par, const mln_value(P)& x) + { + if (par(x) == x) + return x; + else + return par(x) = find_root__(par, par(x)); + } + + + template <typename I, typename A, typename N> + void run_run(const I& f, const A& a, const N& nbh) + { + typedef p_array<mln_psite(I)> S; + S s = level::sort_psites_increasing(a); + // s maps increasing attributes. + + mln_ch_value(I, mln_site(I)) par; + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, std::set<unsigned>) labels; + unsigned nbassins, current_n; + + + // Initialization. + { + mln_piter(A) p(f.domain()); + + // parent + initialize(par, f); + for_all(p) + par(p) = p; + + // deja_vu + initialize(deja_vu, f); + level::fill(deja_vu, false); + + // labels + mln_ch_value(I, unsigned) regmin = labeling::regional_minima(a, nbh, + nbassins); + initialize(labels, f); + for_all(p) + if (regmin(p) != 0) // p in a reg min of the attribute image + labels(p).insert(regmin(p)); + } + + current_n = nbassins; + + // First pass. + { + mln_site(I) r; + mln_fwd_piter(S) p(s); + mln_niter(N) n(nbh, p); + for_all(p) + { + for_all(n) + if (a.domain().has(n) && deja_vu(n)) + { + r = find_root__(par, n); + if (r != p) + { + par(r) = p; // Union. + + if (labels(r).empty()) + // No-op. + ; + else + if (labels(p).empty()) + labels(p) = labels(r); + else + if (labels(p) == labels(r)) + // No-op. + ; + else + { + labels(p).insert(labels(r).begin(), labels(r).end()); + --current_n; + } + } + } + deja_vu(p) = true; + } + } + + std::cout << std::endl; + std::cout << "end = " << current_n << std::endl; + } + + +} // mln + + + + +int main(int, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + int_u8 n; + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, argv[1]); + + typedef p_array<point2d> S; + S s = level::sort_psites_decreasing(f); + + // Children go towards lower levels so leafs are regional minima. + // We get a min-tree so that we can perform morphological closings. + + morpho::tree::data<I,S> t(f, s, c4()); + accu::count< util::pix<I> > attr; + + image2d<unsigned> a = morpho::tree::compute_attribute_image(attr, t); + run_run(f, a, c4()); + +} Index: geraud/tufa_2008/filter.cc --- geraud/tufa_2008/filter.cc (revision 0) +++ geraud/tufa_2008/filter.cc (revision 0) @@ -0,0 +1,169 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/filter.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> +#include <mln/literal/black.hh> + +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/io/ppm/save.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/level/sort_psites.hh> +#include <mln/level/fill.hh> + +#include <mln/morpho/tree/data.hh> +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/morpho/closing_volume.hh> + + + +namespace mln +{ + + template <typename P> + inline + mln_value(P) find_root__(P& par, const mln_value(P)& x) + { + if (par(x) == x) + return x; + else + return par(x) = find_root__(par, par(x)); + } + + + template <typename I, typename A, typename N> + mln_concrete(I) filtering(const I& f, const A& a, const N& nbh, mln_value(A) lambda) + { + typedef p_array<mln_psite(I)> S; + S s = level::sort_psites_increasing(a); + + // s maps increasing attributes. + + mln_concrete(I) out; + initialize(out, f); + + mln_ch_value(I, mln_site(I)) par; + mln_ch_value(I, bool) deja_vu; + { + initialize(par, f); + mln_piter(A) p(par.domain()); + for_all(p) + par(p) = p; + initialize(deja_vu, f); + level::fill(deja_vu, false); + } + + // First pass. + { + mln_site(I) r; + mln_fwd_piter(S) p(s); + mln_niter(N) n(nbh, p); + for_all(p) + { + for_all(n) + if (a.domain().has(n) && deja_vu(n)) + { + r = find_root__(par, n); + if (r != p) + if (f(r) == f(p) || a(r) < lambda) + par(r) = p; // Union. + } + deja_vu(p) = true; + } + } + + // Second pass. + { + mln_bkd_piter(S) p(s); + for_all(p) + if (par(p) == p) + out(p) = f(p); + else + out(p) = out(par(p)); + } + return out; + } + + +} // mln + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm lambda output.pgm" << std::endl; + std::abort(); +} + + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + if (argc != 4) + usage(argv); + + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, argv[1]); + + unsigned lambda = std::atoi(argv[2]); + + typedef p_array<point2d> S; + S s = level::sort_psites_decreasing(f); + + // Children go towards lower levels so leafs are regional minima. + // We get a min-tree so that we can perform morphological closings. + + morpho::tree::data<I,S> t(f, s, c4()); + accu::volume<I> attr; + + image2d<unsigned> a = morpho::tree::compute_attribute_image(attr, t); + + I g = filtering(f, a, c4(), lambda); + + { + I ref = morpho::closing_volume(f, c4(), lambda); + if (g != ref) + { + io::pgm::save(ref, "ref.pgm"); + std::cerr << "oops!" << std::endl; + } + } + + io::pgm::save(g, argv[3]); +} Index: geraud/tufa_2008/wst_f_equal_wst_a.cc --- geraud/tufa_2008/wst_f_equal_wst_a.cc (revision 0) +++ geraud/tufa_2008/wst_f_equal_wst_a.cc (revision 0) @@ -0,0 +1,140 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/steps.cc + +#include <mln/core/image/image2d.hh> +#include <mln/pw/all.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/ppm/save.hh> +#include <mln/io/pbm/save.hh> +#include <mln/literal/black.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/level/transform.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/morpho/tree/data.hh> + +#include <mln/accu/volume.hh> +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/morpho/meyer_wst.hh> + +#include <mln/labeling/regional_minima.hh> + +#include <mln/core/var.hh> + + +namespace mln +{ + + struct colorize : Function_v2v< colorize > + { + typedef value::rgb8 result; + colorize(unsigned max) + : lut(max + 1) + { + lut[0] = 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; + }; + +} // mln + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl; + std::abort(); +} + + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + if (argc != 2) + usage(argv); + + unsigned nref, n; + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, argv[1]); + + labeling::regional_minima(f, c4(), nref); + std::cout << nref << std::endl; + + + image2d<unsigned> wst_f = morpho::meyer_wst(f, c4(), n); + mln_assertion(n == nref); + + io::ppm::save(level::transform(wst_f, colorize(n)), + "wst_f.ppm"); + mln_VAR(WST_f, (pw::value(wst_f) == pw::cst(0u)) | f.domain()); + io::pbm::save(WST_f, "wst_f.pbm"); + + + typedef p_array<point2d> S; + S s = level::sort_psites_decreasing(f); + + // Children go towards lower levels so leafs are regional minima. + // We get a min-tree so that we can perform morphological closings. + + morpho::tree::data<I,S> t(f, s, c4()); + accu::volume<I> vol; + image2d<unsigned> a = morpho::tree::compute_attribute_image(vol, t); + + labeling::regional_minima(a, c4(), n); + mln_assertion(n == nref); + + + image2d<unsigned> wst_a = morpho::meyer_wst(a, c4(), n); + mln_assertion(n == nref); + + io::ppm::save(level::transform(wst_a, colorize(n)), + "wst_a.ppm"); + mln_VAR(WST_a, (pw::value(wst_a) == pw::cst(0u)) | f.domain()); + io::pbm::save(WST_a, "wst_a.pbm"); + + + io::pbm::save((pw::value(WST_a) != pw::value(WST_f)) | f.domain(), + "diff.pbm"); +} Index: geraud/tufa_2008/closed_gradient.cc --- geraud/tufa_2008/closed_gradient.cc (revision 0) +++ geraud/tufa_2008/closed_gradient.cc (revision 0) @@ -0,0 +1,64 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/closed_gradient.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/core/alias/neighb2d.hh> + +#include <mln/morpho/elementary/gradient.hh> +#include <mln/morpho/closing_volume.hh> + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm lambda output.pgm" << std::endl; + std::abort(); +} + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + if (argc != 4) + usage(argv); + + image2d<int_u8> f; + io::pgm::load(f, argv[1]); + unsigned lambda = std::atoi(argv[2]); + io::pgm::save(morpho::closing_volume(morpho::elementary::gradient(f, c4()), + c4(), + lambda), + argv[3]); +} Index: geraud/tufa_2008/fz_count.cc --- geraud/tufa_2008/fz_count.cc (revision 0) +++ geraud/tufa_2008/fz_count.cc (revision 0) @@ -0,0 +1,159 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/fz_count.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/core/alias/neighb2d.hh> + +#include <mln/labeling/regional_minima.hh> +#include <mln/labeling/flat_zones.hh> +#include <mln/level/fill.hh> + + +namespace mln +{ + + template <typename I> + void println_par(const I& par) + { + int nr = par.nrows(), nc = par.ncols(); + for (int r = 0; r < nr; ++r) + { + for (int c = 0; c < nc; ++c) + if (par.at(r,c) == point2d(r,c)) + std::cout << "( ) "; + else + std::cout << par.at(r,c) << ' '; + std::cout << std::endl; + } + } + + template <typename P> + inline + mln_value(P) find_root__(P& par, const mln_value(P)& x) + { + if (par(x) == x) + return x; + else + return par(x) = find_root__(par, par(x)); + } + + + template <typename I, typename N> + unsigned fz_count(const I& f, const N& nbh) + { + typedef p_array<mln_psite(I)> S; + S s = level::sort_psites_increasing(f); + // s maps increasing attributes. + + mln_ch_value(I, mln_site(I)) par; + mln_ch_value(I, bool) deja_vu; + + unsigned counter = f.domain().nsites(); + + // Initialization. + { + mln_piter(I) p(f.domain()); + + // parent + initialize(par, f); + for_all(p) + par(p) = p; + + // deja_vu + initialize(deja_vu, f); + level::fill(deja_vu, false); + } + + // First pass. + { + mln_site(I) r; + mln_fwd_piter(S) p(s); + mln_niter(N) n(nbh, p); + for_all(p) + { + for_all(n) + if (f.domain().has(n) && deja_vu(n)) + { + r = find_root__(par, n); + if (r != p) // not already merged + { + if (f(r) == f(p)) + { + + // Moving the line below out of this test + // (either before the test or after the block) + // makes the algorithm fail. + // The erroneous result is less than the ref + // result => we want this current block too many + // times. Since we merge (thru "par(r) = p") + // whatever "f(r) == f(p)" is true or not, we + // have more often "f(r) == f(p)" than expected. + par(r) = p; + + --counter; + } + } + } + + deja_vu(p) = true; + } + } + + return counter; + } + + +} // mln + + + + +int main(int, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, argv[1]); + + unsigned ref, n = fz_count(f, c4()); + labeling::flat_zones(f, c4(), ref); + + if (n == ref) + std::cout << "success: n flat zones = " << n << std::endl; + else + std::cout << "FAILURE: found = " << n << " v. ref = " << ref << std::endl; +} Index: geraud/tufa_2008/regmin_count.cc --- geraud/tufa_2008/regmin_count.cc (revision 0) +++ geraud/tufa_2008/regmin_count.cc (revision 0) @@ -0,0 +1,173 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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, 51 Franklin Street, Fifth Floor, +// 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. + +/// \file sandbox/geraud/tufa/regmin_count.cc + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/io/pgm/load.hh> +#include <mln/debug/println.hh> + +#include <mln/core/site_set/p_array.hh> +#include <mln/level/sort_psites.hh> +#include <mln/core/alias/neighb2d.hh> + +#include <mln/labeling/regional_minima.hh> +#include <mln/labeling/flat_zones.hh> +#include <mln/level/fill.hh> + + +namespace mln +{ + + template <typename I> + void println_par(const I& par) + { + int nr = par.nrows(), nc = par.ncols(); + for (int r = 0; r < nr; ++r) + { + for (int c = 0; c < nc; ++c) + if (par.at(r,c) == point2d(r,c)) + std::cout << "( ) "; + else + std::cout << par.at(r,c) << ' '; + std::cout << std::endl; + } + } + + template <typename P> + inline + mln_value(P) find_root__(P& par, const mln_value(P)& x) + { + if (par(x) == x) + return x; + else + return par(x) = find_root__(par, par(x)); + } + + + template <typename I, typename N> + unsigned regmin_count(const I& f, const N& nbh) + { + typedef p_array<mln_psite(I)> S; + S s = level::sort_psites_increasing(f); + // s maps increasing attributes. + + mln_ch_value(I, mln_site(I)) par; + mln_ch_value(I, bool) deja_vu, flag; + + unsigned counter = f.domain().nsites(); + + + // Initialization. + { + mln_piter(I) p(f.domain()); + + // parent + initialize(par, f); + for_all(p) + par(p) = p; + + // flag + initialize(flag, f); + level::fill(flag, true); + + // deja_vu + initialize(deja_vu, f); + level::fill(deja_vu, false); + } + + // First pass. + { + mln_site(I) r; + mln_fwd_piter(S) p(s); + mln_niter(N) n(nbh, p); + for_all(p) + { + unsigned loc = 0; + for_all(n) + if (f.domain().has(n) && deja_vu(n)) + { + r = find_root__(par, n); + if (r != p) + { + if (f(r) == f(p)) + { + par(r) = p; // Union. + if (flag(p) == false && flag(r) == false) + { + // Two non-reg-min components merge (same flat + // zone) so we had an extra invalidation. + ++counter; + } + flag(p) = flag(p) && flag(r); + --counter; // So we get the number of flat zones + // minus the non-reg-min flat zones. + } + else + { + mln_invariant(f(r) < f(p)); + if (flag(p) == true) + { + ++loc; + --counter; // Invalidation. + } + flag(p) = false; + } + } + } + mln_invariant(loc == 0 || loc == 1); + deja_vu(p) = true; + } + } + + return counter; + } + + +} // mln + + + + +int main(int, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> I; + I f; + io::pgm::load(f, argv[1]); + + unsigned ref, n = regmin_count(f, c4()); + labeling::regional_minima(f, c4(), ref); + + if (n == ref) + std::cout << "success: n regional minima = " << n << std::endl; + else + std::cout << "FAILURE: found = " << n << " v. ref = " << ref << std::endl; +}