[Olena-patches] 3376: Update ICIP 2009 code.

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; Update ICIP 2009 code. * theo/publis, * theo/publis/icip2009: New directories. * theo/tufa_2008/regmin_count.cc: Move... * theo/publis/icip2009/regmin_count.cc: ...here. * theo/tufa_2008/compute_a.cc: Move... * theo/publis/icip2009/compute_a.cc: ...here. Update. compute_a.cc | 366 ++++++++++++++++++++++++++++++----------------------------- 1 file changed, 190 insertions(+), 176 deletions(-) Index: theo/publis/icip2009/compute_a.cc --- theo/publis/icip2009/compute_a.cc (revision 0) +++ theo/publis/icip2009/compute_a.cc (working copy) @@ -1,4 +1,5 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// Copyright (C) 2008, 2009 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 @@ -25,7 +26,7 @@ // reasons why the executable file might be covered by the GNU General // Public License. -/// \file sandbox/geraud/tufa/compute_a.cc +/// \file sandbox/theo/publis/icip2009/compute_a.cc #include <mln/core/image/image2d.hh> #include <mln/core/image/image_if.hh> @@ -47,6 +48,7 @@ #include <mln/accu/count.hh> #include <mln/opt/at.hh> +#include <mln/set/card.hh> @@ -80,6 +82,16 @@ } + template <typename P> + p_array<P> reverse(const p_array<P>& arr) + { + p_array<P> out; + mln_bkd_piter(p_array<P>) p(arr); + for_all(p) + out.insert(p); + return out; + } + //------------------------------- compute_a @@ -87,14 +99,16 @@ 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) + compute_a(const I& f, const N& nbh, A, unsigned& n_regmins, + bool echo = false) { - typedef p_array<mln_psite(I)> S; + typedef mln_site(I) P; + typedef p_array<P> 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, P) par; + mln_ch_value(I, P) zpar; mln_ch_value(I, bool) deja_vu, flag; mln_ch_value(I, A) attr; @@ -122,7 +136,7 @@ // First Pass. { - mln_site(I) r; + P r; mln_fwd_piter(S) p(s); mln_niter(N) n(nbh, p); for_all(p) @@ -138,230 +152,211 @@ r = find_root__(zpar, n); if (r != p) { - // Fully compressed union. - zpar(r) = p; - attr(p).take(attr(r)); + par(r) = p; // Min-tree computation. + zpar(r) = p; // Fully compressed union. + + attr(p).take(attr(r)); // Attribute computation. if (f(r) == f(p)) { - // Weak-Union; only for flat zones. - par(r) = p; - + // Flat zones. if (flag(p) == false && flag(r) == false) { // Two non-reg-min components merge (same flat // zone) so we had an extra invalidation. ++n_regmins; } + else 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)); + mln_invariant(f(r) < f(p)); // Increasing browsing. 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. + // Output is the attribute image. mln_ch_value(I, mln_result(A)) a; initialize(a, f); data::paste(attr, a); - // Finalization. + + // Back-propagate flags to non-representant component sites. { 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; + if (f(par(p)) == f(p) && flag(par(p)) != flag(p)) + flag(p) = flag(par(p)); + if (echo) + debug::println("flag", flag); } - // TODO: compress at least the reg minima! + unsigned + n_sites = 0, + n_roots = 0, + n_nodes = 0; + + // Tree canonization { - image2d<unsigned> regmin; - initialize(regmin, f); - { - unsigned i_regmin = 0; - mln_bkd_piter(S) p(s); + mln_bkd_piter(S) p(s); // Reverse. for_all(p) { + P q = par(p); + if (f(par(q)) == f(q)) + par(p) = par(q); + ++n_sites; if (par(p) == p) + ++n_roots; + if (par(p) == p || f(par(p)) != f(p)) + ++n_nodes; + } + + // Test. { - if (flag(p)) - regmin(p) = ++i_regmin; - else - regmin(p) = 0; + mln_piter(S) p(s); + for_all(p) + { + P q = par(p); + if (q == p || par(q) == q) // Either p or par(p) is root. + continue; + if (f(q) == f(p)) + // p and par(p) within a flat zone + // so q=par(p) is the representant + mln_invariant(f(par(q)) != f(q)); } - 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... + std::cout << "%age nodes = " + << (float(n_nodes) / n_sites * 100.f) << std::endl + << std::endl; + - debug::println("regmin", regmin); + { + // 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_bkd_piter(S) p(s); // Reverse. + for_all(p) + if (f(par(p)) == f(p)) + a(p) = a(par(p)); + } - // 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 ? + // Display tree. + if (echo) + { + mln_ch_value(I, char) tree; + initialize(tree, f); + data::fill(tree, ' '); + mln_piter(I) p(f.domain()); + for_all(p) + if (par(p) == p) // Root. + tree(p) = '@'; + else + if (f(par(p)) != f(p)) // Representative. + tree(p) = 'o'; + else + tree(p) = '.'; + debug::println(tree); + } - image2d<bool> seen; - initialize(seen, f); - data::fill(seen, false); + { + mln_ch_value(I, unsigned) nchildren; + initialize(nchildren, f); + data::fill(nchildren, 0); - s = level::sort_psites_increasing(a); - mln_bkd_piter(S) p(s); + mln_fwd_piter(S) p(s); for_all(p) + if (f(par(p)) != f(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; + // Since the tree is canonized, par(p) is a node + // (a component representative): + P q = par(p); + mln_invariant(par(q) == q || f(par(q)) != f(q)); + ++nchildren(par(p)); } - debug::println(seen); -// if (flag(p)) -// std::cout << a(p) << ' ' << p << ' ' << (par(p) == p) << std::endl; + if (echo) + debug::println("nchildren", nchildren); + // Test. + { + typedef morpho::tree::data<I,S> T; + S s_ = reverse(s); + T t(f, s_, nbh); + + typedef typename T::function C; + mln_ch_value(C, unsigned) ref; + initialize(ref, t.f()); + data::fill(ref, 0); + + mln_fwd_piter(T) p(t.domain()); + for_all(p) + if (t.is_a_non_root_node(p)) + { + mln_invariant(t.is_a_node(t.parent(p))); + ++ref(t.parent(p)); } - return a; + mln_invariant((nchildren | t.domain()) == ref); + } } - //------------------------------- filtering + // Tests. + { + { + unsigned n_ref; + mln_ch_value(I, unsigned) ref; + ref = labeling::regional_minima(f, nbh, n_ref); + // debug::println("ref", ref); + + mln_assertion(n_regmins == n_ref); + mln_assertion(((pw::value(ref) != 0) | ref.domain()) == flag); + } + { + unsigned n = 0; + mln_fwd_piter(S) p(s); + for_all(p) + if (f(par(p)) != f(p) && flag(p)) + ++n; + mln_assertion(n == n_regmins); + } + } // end of Tests. -// 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); -// data::fill(deja_vu, false); - -// // flag -// initialize(flag, f); -// data::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; -// } + return a; + } // compute @@ -370,43 +365,62 @@ -int main(int, char* argv[]) +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm echo" << std::endl; + std::abort(); +} + + + +int main(int argc, char* argv[]) { using namespace mln; using value::int_u8; + if (argc != 3) + usage(argv); + + int echo = std::atoi(argv[2]); + if (echo != 0 && echo != 1) + usage(argv); + + + neighb2d nbh = c4(); + + typedef image2d<int_u8> I; I f; io::pgm::load(f, argv[1]); - // debug::println(f); + if (echo) + debug::println("f", 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); -// } + image2d<unsigned> a = compute_a(f, nbh, area, n_regmins, echo); + if (echo) + debug::println("a", a); + + + + + { + // Test of 'a'. + typedef p_array<point2d> S; + S s = level::sort_psites_decreasing(f); + morpho::tree::data<I,S> t(f, s, nbh); + accu::count< util::pix<I> > area; + image2d<unsigned> a_ref = morpho::tree::compute_attribute_image(area, t); + // debug::println("a_ref", a_ref); + mln_assertion(a == a_ref); + } // unsigned n_wanted = 10; -// I g = filtering(f, a, c4(), n_regmins, n_wanted); +// I g = filtering(f, a, nbh, n_regmins, n_wanted); } Property changes on: theo/publis/icip2009/compute_a.cc ___________________________________________________________________ Added: svn:mergeinfo Property changes on: theo/publis/icip2009/regmin_count.cc ___________________________________________________________________ Added: svn:mergeinfo