* edwin/TreeAlgorithmsComp/Makefile: New. * edwin/TreeAlgorithmsComp/bench.hh, * edwin/TreeAlgorithmsComp/nnodes.hh: Shared headers. * edwin/TreeAlgorithmsComp/parallel_union_find.hh, * edwin/TreeAlgorithmsComp/salembier.hh, * edwin/TreeAlgorithmsComp/max_tree.cc, * edwin/TreeAlgorithmsComp/max_tree_with_rank.cc, * edwin/TreeAlgorithmsComp/salembier_.cc, * edwin/TreeAlgorithmsComp/tbb_union_find_8.cc: Routines to build components trees (Minor updates). * edwin/TreeAlgorithmsComp/compute_histo_.cc: Tool to output the histogram of an image whatever quantized. --- milena/sandbox/ChangeLog | 17 ++ milena/sandbox/edwin/TreeAlgorithmsComp/Makefile | 87 ++++++++++ .../edwin/TreeAlgorithmsComp/bench.hh} | 48 +++++- .../edwin/TreeAlgorithmsComp/compute_histo_.cc | 48 ++++++ .../sandbox/edwin/TreeAlgorithmsComp/max_tree.cc | 61 ++++---- .../edwin/TreeAlgorithmsComp/max_tree_with_rank.cc | 131 +++++++++------- milena/sandbox/edwin/TreeAlgorithmsComp/nnodes.hh | 24 +++ .../TreeAlgorithmsComp/parallel_union_find.hh | 19 ++- .../sandbox/edwin/TreeAlgorithmsComp/salembier.hh | 169 ++++++++++++++++++++ .../sandbox/edwin/TreeAlgorithmsComp/salembier_.cc | 51 ++++++ .../edwin/TreeAlgorithmsComp/tbb_union_find_8.cc | 48 ++---- 11 files changed, 570 insertions(+), 133 deletions(-) create mode 100644 milena/sandbox/edwin/TreeAlgorithmsComp/Makefile copy milena/{mln/fun/x2p/essential.hh => sandbox/edwin/TreeAlgorithmsComp/bench.hh} (58%) create mode 100644 milena/sandbox/edwin/TreeAlgorithmsComp/compute_histo_.cc create mode 100644 milena/sandbox/edwin/TreeAlgorithmsComp/nnodes.hh create mode 100644 milena/sandbox/edwin/TreeAlgorithmsComp/salembier.hh create mode 100644 milena/sandbox/edwin/TreeAlgorithmsComp/salembier_.cc diff --git a/milena/sandbox/ChangeLog b/milena/sandbox/ChangeLog index 15f7fa8..5d3a05f 100644 --- a/milena/sandbox/ChangeLog +++ b/milena/sandbox/ChangeLog @@ -1,3 +1,20 @@ +2010-06-02 edwin carlinet &lt;carlinet(a)lrde.epita.fr&gt; + + * edwin/TreeAlgorithmsComp/Makefile: New. + * edwin/TreeAlgorithmsComp/bench.hh, + * edwin/TreeAlgorithmsComp/nnodes.hh: Shared headers. + + * edwin/TreeAlgorithmsComp/parallel_union_find.hh, + * edwin/TreeAlgorithmsComp/salembier.hh, + * edwin/TreeAlgorithmsComp/max_tree.cc, + * edwin/TreeAlgorithmsComp/max_tree_with_rank.cc, + * edwin/TreeAlgorithmsComp/salembier_.cc, + * edwin/TreeAlgorithmsComp/tbb_union_find_8.cc: + Routines to build components trees (Minor updates). + + * edwin/TreeAlgorithmsComp/compute_histo_.cc: + Tool to output the histogram of an image whatever quantized. + 2010-05-04 edwin carlinet &lt;carlinet(a)lrde.epita.fr&gt; Comparison of tree build algorithms. diff --git a/milena/sandbox/edwin/TreeAlgorithmsComp/Makefile b/milena/sandbox/edwin/TreeAlgorithmsComp/Makefile new file mode 100644 index 0000000..fea59aa --- /dev/null +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/Makefile @@ -0,0 +1,87 @@ +-include makefile.rules +#DEBUG = 1 #comment +WALL_CLOCK=1 + +MLN_DIR = ~/git/olena/milena +TBB_DIR = /home/edwin/Telechargement/tbb30_20100406oss +TBB_INCLUDE_DIR = $(TBB_DIR)/include +TBB_LIB_DIR = $(TBB_DIR)/lib/intel64/cc3.4.3_libc2.3.4_kernel2.6.9 +CXXFLAGS += -W -Wall +CXXFLAGS += $(if $(DEBUG), -g -ggdb, -DNDEBUG $(if $(RELEASE), -O3, -O1)) +CXXFLAGS += -I $(MLN_DIR) +CXXFLAGS += $(if $(WALL_CLOCK), -DWALL_CLOCK) +LDFLAGS = $(if $(WALL_CLOCK), -ltbb) + +CC = g++ + +tbb_CXXFLAGS = $(CXXFLAGS) -I $(TBB_INCLUDE_DIR) -DWALL_CLOCK +tbb_LDFLAGS = -L $(TBB_LIB_DIR) -ltbb + +TARGETS = \ + tbb_union_find_8 \ + salembier_08 \ + salembier_12 \ + salembier_16 \ + max_tree_08 \ + max_tree_12 \ + max_tree_16 \ + max_tree_with_rank_08 \ + max_tree_with_rank_12 \ + max_tree_with_rank_16 \ + convert_8_to_12 \ + convert_8_to_16 \ + compute_histo_8 \ + compute_histo_12 \ + compute_histo_16 + + +all: $(TARGETS) + +tbb_union_find_8: tbb_union_find_8.cc + $(CXX) $(tbb_CXXFLAGS) $(tbb_LDFLAGS) $< -o $@ + +clean: + rm -f *.o + rm -f $(TARGETS) + +convert_8_to_12: convert_8_to_.cc + $(CXX) -DQ=int_u12 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +convert_8_to_16: convert_8_to_.cc + $(CXX) -DQ=int_u16 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +compute_histo_8: compute_histo_.cc + $(CXX) -DQ=int_u8 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +compute_histo_12: compute_histo_.cc + $(CXX) -DQ=int_u12 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +compute_histo_16: compute_histo_.cc + $(CXX) -DQ=int_u16 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +max_tree_with_rank_08: max_tree_with_rank.cc + $(CXX) -DQ=int_u8 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +max_tree_with_rank_12: max_tree_with_rank.cc + $(CXX) -DQ=int_u12 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +max_tree_with_rank_16: max_tree_with_rank.cc + $(CXX) -DQ=int_u16 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +salembier_08: salembier_.cc + $(CXX) -DQ=int_u8 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +salembier_12: salembier_.cc + $(CXX) -DQ=int_u12 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +salembier_16: salembier_.cc + $(CXX) -DQ=int_u16 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +max_tree_08: max_tree.cc + $(CXX) -DQ=int_u8 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +max_tree_12: max_tree.cc + $(CXX) -DQ=int_u12 $(CXXFLAGS) $(LDFLAGS) $< -o $@ + +max_tree_16: max_tree.cc + $(CXX) -DQ=int_u16 $(CXXFLAGS) $(LDFLAGS) $< -o $@ \ No newline at end of file diff --git a/milena/mln/fun/x2p/essential.hh b/milena/sandbox/edwin/TreeAlgorithmsComp/bench.hh similarity index 58% copy from milena/mln/fun/x2p/essential.hh copy to milena/sandbox/edwin/TreeAlgorithmsComp/bench.hh index fb8fd10..a68ab71 100644 --- a/milena/mln/fun/x2p/essential.hh +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/bench.hh @@ -23,14 +23,46 @@ // exception does not however invalidate any other reasons why the // executable file might be covered by the GNU General Public License. -#ifndef MLN_FUN_X2P_ESSENTIAL_HH -# define MLN_FUN_X2P_ESSENTIAL_HH +#ifndef BENCH_HH +# define BENCH_HH -/*! \file - * - * \brief File that includes essential functions from point to value. - */ +# ifndef WALL_CLOCK -# include <mln/fun/x2p/all.hh> +# include <ctime> -#endif // ! MLN_FUN_X2P_ESSENTIAL_HH +# define START_BENCH(N_ITERATION) \ + { \ + clock_t start__ = std::clock(); \ + const int n_iteration__ = N_ITERATION; \ + for (int i__ = 0; i__ < n_iteration__; ++i__) { + + +# define END_BENCH(MSG) \ + } \ + std::cout << MSG \ + << (float)(std::clock() - start__) * 1000.0 / \ + (CLOCKS_PER_SEC * n_iteration__) \ + << " ms." << std::endl; \ + } + +# else // !WALL_CLOCK + +# include <tbb/tick_count.h> + +# define START_BENCH(N_ITERATION) \ + { \ + tbb::tick_count t0__ = tbb::tick_count::now(); \ + const int n_iteration__ = N_ITERATION; \ + for (int i__ = 0; i__ < n_iteration__; ++i__) { + +# define END_BENCH(MSG) \ + } \ + tbb::tick_count t1__ = tbb::tick_count::now(); \ + std::cout << MSG \ + << (t1__ - t0__).seconds() * 1000 / n_iteration__ \ + << " ms." << std::endl; \ + } + +# endif + +#endif // !BENCH_HH diff --git a/milena/sandbox/edwin/TreeAlgorithmsComp/compute_histo_.cc b/milena/sandbox/edwin/TreeAlgorithmsComp/compute_histo_.cc new file mode 100644 index 0000000..a56e809 --- /dev/null +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/compute_histo_.cc @@ -0,0 +1,48 @@ +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/int_u12.hh> +#include <mln/value/int_u16.hh> +#include <mln/histo/compute.hh> +#include <mln/io/pgm/load.hh> +#include <fstream> + +#ifndef Q +# define Q int_u12 +#endif + +void usage(char** argv) +{ + std::cerr << "Usage: " << argv[0] << " ima.pgm" << std::endl; + abort(); +} + +int main(int argc, char** argv) +{ + using namespace mln; + using value::int_u8; + using value::int_u12; + using value::int_u16; + typedef Q V; + + if (argc < 2) + usage(argv); + const char* finput = argv[1]; + typedef image2d<V> I; + + I input; + io::pgm::load(input, finput); + + histo::array<V> h = histo::compute(input); + + std::ofstream fo("histo.cvs"); + int nvalues = mln_card(V); + int c = 0; + for (int i = 0; i < nvalues; ++i) + { + if (h[i] != 0) + ++c; + fo << i << ", " << h[i] << std::endl; + } + + fo.close(); +} diff --git a/milena/sandbox/edwin/TreeAlgorithmsComp/max_tree.cc b/milena/sandbox/edwin/TreeAlgorithmsComp/max_tree.cc index 95d5dc1..2052485 100644 --- a/milena/sandbox/edwin/TreeAlgorithmsComp/max_tree.cc +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/max_tree.cc @@ -1,29 +1,21 @@ #include <mln/core/image/image2d.hh> #include <mln/core/alias/neighb2d.hh> -#include <mln/io/pgm/load.hh> -#include <mln/io/pgm/save.hh> - #include <mln/value/int_u8.hh> +#include <mln/value/int_u12.hh> +#include <mln/value/int_u16.hh> + +#include <mln/io/pgm/load.hh> #include <mln/data/sort_offsets.hh> #include <mln/extension/adjust.hh> #include <mln/extension/fill.hh> -#include <mln/debug/println.hh> -#include <mln/util/timer.hh> - -static mln::util::timer tm; -static const int ITERATION = 4; +#include "bench.hh" -#define START_BENCH() \ - tm.restart(); \ - for (int i = 0; i < ITERATION; ++i) { - -#define END_BENCH() \ - } \ - std::cout << ((float)(tm.stop() * 1000.0) / ITERATION) << std::endl; - +#ifndef Q +# define Q int_u8 +#endif namespace mln { @@ -73,7 +65,9 @@ namespace mln radix_decreasing_sort(const I& f, util::array< util::array<unsigned> >& s) { typedef mln_value(I) T; - const unsigned n = unsigned(mln_max(T)) + 1; + const unsigned n = mln_card(T); + // std::cout << "Sort 0" << std::endl; + s.resize(n); // h @@ -82,12 +76,19 @@ namespace mln for_all(pxl) ++h[pxl.val()]; + // std::cout << "Sort 0" << std::endl; + for (unsigned v = 0; v < n; ++v) s[v].reserve(h[v]); + // std::cout << "Sort 0" << std::endl; + + // computing output data for_all(pxl) s[pxl.val()].append(pxl.offset()); + + // std::cout << "Sort 0" << std::endl; } @@ -99,7 +100,7 @@ namespace mln bool run_tree_canonization = true) { typedef mln_value(I) T; - const unsigned n_values = unsigned(mln_max(T)) + 1; + const unsigned n_values = mln_card(T); extension::adjust(f, nbh); extension::fill(f, mln_max(T)); @@ -199,7 +200,7 @@ namespace mln void usage(char* argv[]) { - std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl; + std::cerr << "usage: " << argv[0] << " input.pgm nBench" << std::endl; std::cerr << " max-tree with full optimizations" << std::endl; std::abort(); } @@ -208,25 +209,25 @@ void usage(char* argv[]) int main(int argc, char* argv[]) { - if (argc != 2) + if (argc != 3) usage(argv); using namespace mln; + using value::int_u8; + using value::int_u12; + using value::int_u16; + typedef Q V; neighb2d nbh = c4(); - - image2d<value::int_u8> f; + image2d<V> f; io::pgm::load(f, argv[1]); // go go go image2d<unsigned> parent = compute_max_tree(f, nbh); - - tm.start(); - - START_BENCH(); - compute_max_tree(f, nbh); - END_BENCH(); - - // debug::println(parent); std::cout << "nnodes = " << nnodes(f, parent) << std::endl; + + int nb_bench = atoi(argv[2]); + START_BENCH(nb_bench); + parent = compute_max_tree(f, nbh); + END_BENCH("union find: "); } diff --git a/milena/sandbox/edwin/TreeAlgorithmsComp/max_tree_with_rank.cc b/milena/sandbox/edwin/TreeAlgorithmsComp/max_tree_with_rank.cc index 1e80b12..932cc58 100644 --- a/milena/sandbox/edwin/TreeAlgorithmsComp/max_tree_with_rank.cc +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/max_tree_with_rank.cc @@ -1,18 +1,21 @@ #include <mln/core/image/image2d.hh> #include <mln/core/alias/neighb2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/int_u12.hh> +#include <mln/value/int_u16.hh> + #include <mln/io/pgm/load.hh> -#include <mln/io/pgm/save.hh> -#include <mln/io/dump/save.hh> #include <mln/data/sort_offsets.hh> #include <mln/extension/adjust.hh> #include <mln/extension/fill.hh> -#include <mln/util/timer.hh> -#include <mln/debug/println.hh> - +#include "bench.hh" +#ifndef Q +# define Q int_u8 +#endif namespace mln { @@ -35,7 +38,7 @@ namespace mln void sort_decreasing(const I& input, std::vector<unsigned>& vec) { - const unsigned n = 256; + const unsigned n = mln_card(mln_value(I)); // h std::vector<unsigned> h(n, 0); @@ -59,7 +62,7 @@ namespace mln template <typename I, typename N> mln_ch_value(I, unsigned) compute_max_tree(const I& f, const N& nbh, - bool run_tree_canonization = false) + bool run_tree_canonization = true) { const unsigned npixels = f.domain().nsites(); @@ -74,58 +77,69 @@ namespace mln util::array<int> dp = offsets_wrt(f, nbh); const unsigned n_nbhs = dp.nelements(); - mln_ch_value(I, unsigned) parent, zpar, rnk, last; + mln_ch_value(I, unsigned) parent, zpar, rnk, zpar_to_par; // initialization. initialize(parent, f); initialize(rnk, f); - data::fill(rnk, 0); initialize(zpar, f); + initialize(zpar_to_par, f); + data::fill(rnk, 0); data::fill(zpar, 0); - initialize(last, f); - // Tree construction. for (unsigned i = 0; i < npixels; ++i) { unsigned p = s[i]; - + parent.element(p) = p; zpar.element(p) = p; - last.element(p) = i; + zpar_to_par.element(p) = p; + //last.element(p) = i; for (unsigned j = 0; j < n_nbhs; ++j) { unsigned n = p + dp[j]; - + if (zpar.element(n) == 0) // not deja-vu continue; - - unsigned r_ = zfind_root(zpar, n); - unsigned p_ = zfind_root(zpar, p); - if (r_ != p_) + + unsigned r = zfind_root(zpar, n); + unsigned q = zfind_root(zpar, p); + if (r != q) { - unsigned r = s[last.element(r_)]; - parent.element(r) = p; + unsigned r_ = zpar_to_par.element(r); + parent.element(r_) = p; // Parenthood goes towards root, e.g., to lower levels: - mln_invariant(f.element(p) <= f.element(r)); - - // Link. - if (rnk.element(r_) <= rnk.element(p_)) - { - zpar.element(r_) = p_; - ++rnk.element(p_); - if (last.element(r_) > last.element(p_)) - last.element(p_) = last.element(r_); - } - else - { - zpar.element(p_) = r_; - ++rnk.element(r_); - if (last.element(p_) > last.element(r_)) - last.element(r_) = last.element(p_); - } + // mln_invariant(f.element(p) <= f.element(r)); + if (rnk.element(q) < rnk.element(r)) { + zpar.element(q) = r; + zpar_to_par.element(r) = p; + } else if (rnk.element(q) > rnk.element(r)) { + zpar.element(r) = q; + zpar_to_par.element(q) = p; + } else { + zpar.element(r) = q; + zpar_to_par.element(q) = p; + ++rnk.element(q); + } + + // // Link. + // if (rnk.element(r_) <= rnk.element(p_)) + // { + // zpar.element(r_) = p_; + // ++rnk.element(p_); + // if (last.element(r_) > last.element(p_)) + // last.element(p_) = last.element(r_); + // } + // else + // { + // zpar.element(p_) = r_; + // ++rnk.element(r_); + // if (last.element(p_) > last.element(r_)) + // last.element(r_) = last.element(p_); + // } } } // j @@ -133,16 +147,16 @@ namespace mln } // i -// if (run_tree_canonization) -// // Tree canonization. -// for (int i = npixels - 1; i >= 0; --i) -// { -// unsigned -// p = s[i], -// q = parent.element(p); -// if (f.element(parent.element(q)) == f.element(q)) -// parent.element(p) = parent.element(q); -// } + if (run_tree_canonization) + // Tree canonization. + for (int i = npixels - 1; i >= 0; --i) + { + unsigned + p = s[i], + q = parent.element(p); + if (f.element(parent.element(q)) == f.element(q)) + parent.element(p) = parent.element(q); + } return parent; } @@ -169,7 +183,7 @@ namespace mln void usage(char* argv[]) { - std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl; + std::cerr << "usage: " << argv[0] << " input.pgm nBench" << std::endl; std::cerr << " max-tree with full optimizations" << std::endl; std::abort(); } @@ -178,24 +192,25 @@ void usage(char* argv[]) int main(int argc, char* argv[]) { - if (argc != 2) + if (argc != 3) usage(argv); using namespace mln; + using value::int_u8; + using value::int_u12; + using value::int_u16; + typedef Q V; neighb2d nbh = c4(); - image2d<value::int_u8> f; + image2d<V> f; io::pgm::load(f, argv[1]); - util::timer t; - t.start(); - - // go go go image2d<unsigned> parent = compute_max_tree(f, nbh); - - float ts = t.stop(); - std::cout << ts * 1000 << std::endl; - std::cout << "nnodes = " << nnodes(f, parent) << std::endl; + + int nb_bench = atoi(argv[2]); + START_BENCH(nb_bench); + parent = compute_max_tree(f, nbh); + END_BENCH("Max Tree with Rank: "); } diff --git a/milena/sandbox/edwin/TreeAlgorithmsComp/nnodes.hh b/milena/sandbox/edwin/TreeAlgorithmsComp/nnodes.hh new file mode 100644 index 0000000..c900f9e --- /dev/null +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/nnodes.hh @@ -0,0 +1,24 @@ +#ifndef NNODES_HH +# define NNODES_HH + +namespace mln +{ + + template <typename I, typename J> + unsigned nnodes(const I& f, const J& parent) + { + mlc_equal(mln_value(J), unsigned)::check(); + unsigned n = 0; + mln_piter(I) p(parent.domain()); + for_all(p) + { + unsigned o = parent.index_of_point(p); + if (parent(p) == o || f.element(parent(p)) != f(p)) + ++n; + } + return n; + } + +} + +#endif // !NNODES_HH diff --git a/milena/sandbox/edwin/TreeAlgorithmsComp/parallel_union_find.hh b/milena/sandbox/edwin/TreeAlgorithmsComp/parallel_union_find.hh index a11253b..c01be7d 100644 --- a/milena/sandbox/edwin/TreeAlgorithmsComp/parallel_union_find.hh +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/parallel_union_find.hh @@ -48,7 +48,14 @@ std::cout << "End task: " << task__.str() << " in " << (t1-t0).seconds() << std::endl; \ } - +#ifdef TRACE +# include <tbb/compat/thread> +# define TRACE_START(TASK) BENCH_START(TASK) +# define TRACE_END() BENCH_END() +#else +# define TRACE_START(TASK) ; +# define TRACE_END() ; +#endif namespace mln @@ -248,7 +255,7 @@ ll_union_find<I, N>::operator() (const tbb::blocked_range<int>& r) subdomain.pmax()[0] = r.end() - 1; this->build(subdomain); - std::cout << "Warning: Not parallel merge. " << std::endl; + //std::cout << "Warning: Not parallel merge. " << std::endl; this->merge(subdomain); this->domain.pmax() = subdomain.pmax(); } @@ -259,7 +266,7 @@ template <typename I, typename N> void ll_union_find<I, N>::build(const box2d& subdomain) { - BENCH_START("Build: " << subdomain) + TRACE_START("Build: " << subdomain << " T: " << std::this_thread::get_id()) mln_precondition(subdomain.is_valid()); const unsigned n_values = unsigned(mln_max(T)) + 1; @@ -307,7 +314,7 @@ ll_union_find<I, N>::build(const box2d& subdomain) zpar.element(p) = zpar.element(zpar.element(p)); } } - BENCH_END() + TRACE_END() } @@ -318,7 +325,7 @@ ll_union_find<I, N>::merge(const box2d& d2) const int stride = f.delta_index(dpoint2d(1, 0)); box2d& d1 = this->domain; - BENCH_START("Merge: " << d1 << "&" << d2) + TRACE_START("Merge: " << d1 << "&" << d2 << " T: " << std::this_thread::get_id()) // Split along x axis assert(d1.pmin()[1] == d2.pmin()[1]); @@ -335,7 +342,7 @@ ll_union_find<I, N>::merge(const box2d& d2) ++x2; } - BENCH_END() + TRACE_END() } /// Merge two trees. diff --git a/milena/sandbox/edwin/TreeAlgorithmsComp/salembier.hh b/milena/sandbox/edwin/TreeAlgorithmsComp/salembier.hh new file mode 100644 index 0000000..d58a8dd --- /dev/null +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/salembier.hh @@ -0,0 +1,169 @@ +// Copyright (C) 2008, 2009 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of Olena. +// +// Olena is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free +// Software Foundation, version 2 of the License. +// +// Olena 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 Olena. If not, see <http://www.gnu.org/licenses/>. +// +// As a special exception, you may use this file as part of a free +// software project without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to produce +// an executable, this file does not by itself cause the resulting +// executable to be covered by the GNU General Public License. This +// exception does not however invalidate any other reasons why the +// executable file might be covered by the GNU General Public License. + +#ifndef SALEMBIER_HH +# define SALEMBIER_HH + +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/extension/fill.hh> +#include <mln/extension/adjust.hh> +#include <mln/data/fill.hh> +#include <mln/util/hqueues.hh> + +namespace mln +{ + + template <typename T> + struct salembier + { + enum { nvalues = mln_card(T) }; + + template <typename N> + salembier(const image2d<T>& f, const N& n); + + private: + unsigned flood(unsigned h); + + public: + const image2d<T>& f; + image2d<unsigned> parent; + + private: + image2d<bool> deja_vu; + unsigned levroot[nvalues]; + bool is_node_at_level[nvalues]; + util::hqueues<unsigned, T> hqueues; + + // Aux + unsigned hmin; + util::array<int> dp; + unsigned n_nbhs; + }; + +# ifndef MLN_INCLUDE_ONLY + + template <typename T> + template <typename N> + salembier<T>::salembier(const image2d<T>& f_, const N& nbh) + : f(f_) + { + // Retrieve pmin, hmin + unsigned pmin; + T vmin; + { + util::array< unsigned > h(nvalues, 0); + + pmin = f.index_of_point(f.domain().pmin()); + vmin = f.element(pmin); + + mln_fwd_pixter(const image2d<T>) pxl(f); + for_all(pxl) + { + ++h[pxl.val()]; + if (pxl.val() < vmin) + { + vmin = pxl.val(); + pmin = pxl.offset(); + } + } + + for (unsigned i = 0; i < nvalues; ++i) + hqueues[i].reserve(h[i]); + } + + // Initialize aux data + { + dp = offsets_wrt(f, nbh); + n_nbhs = dp.nelements(); + + extension::adjust(f, nbh); + initialize(parent, f); + initialize(deja_vu, f); + + data::fill(deja_vu, false); + extension::fill(deja_vu, true); + + //std::fill(levroot, levroot + 256, 0); + std::fill(is_node_at_level, is_node_at_level + nvalues, false); + } + + // Start flooding + hmin = vmin; + hqueues[hmin].push(pmin); + levroot[hmin] = pmin; + is_node_at_level[hmin] = true; + deja_vu.element(pmin) = true; + flood(hmin); + } + + + template <typename T> + unsigned + salembier<T>::flood(unsigned h) + { + while (!hqueues[h].empty()) + { + unsigned p = hqueues[h].pop_front(); + parent.element(p) = levroot[h]; + + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p + dp[j]; + // offset 0 does not belong to domain. + // parent.element(n) == 0 imples: + // n has not been processed and n belongs to domain. + if (!deja_vu.element(n)) + { + unsigned l = f.element(n); + hqueues[l].push(n); + deja_vu.element(n) = true; + + if (!is_node_at_level[l]) + { + levroot[l] = n; + is_node_at_level[l] = true; + } + + while (h < l) + l = flood(l); + } + } + } + + // Attach to parent. + unsigned x = levroot[h]; + is_node_at_level[h] = false; + while (h > hmin && !is_node_at_level[h]) + --h; + parent.element(x) = levroot[h]; + return h; + } + +# endif // ! MLN_INCLUDE_ONLY + +} + +#endif // !SALEMBIER_HH diff --git a/milena/sandbox/edwin/TreeAlgorithmsComp/salembier_.cc b/milena/sandbox/edwin/TreeAlgorithmsComp/salembier_.cc new file mode 100644 index 0000000..a638d14 --- /dev/null +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/salembier_.cc @@ -0,0 +1,51 @@ +#include <mln/core/image/image2d.hh> + +#include <mln/value/int_u8.hh> +#include <mln/value/int_u12.hh> +#include <mln/value/int_u16.hh> + +#include <mln/io/pgm/load.hh> + +#include "salembier.hh" +#include "nnodes.hh" +#include "bench.hh" + +#ifndef Q +# define Q int_u8 +#endif + +void usage(char** argv) +{ + std::cout << "Usage: " << argv[0] << " input.pgm nBench" << std::endl; + abort(); +} + +int main(int argc, char** argv) +{ + if (argc != 3) + usage(argv); + + using namespace mln; + using value::int_u8; + using value::int_u12; + using value::int_u16; + typedef Q V; + + neighb2d nbh = c4(); + image2d<V> input; + io::pgm::load(input, argv[1]); + + salembier<V> s(input, nbh); + image2d<unsigned> parent = s.parent; + std::cout << "nnodes = " << nnodes(input, parent) << std::endl; + + int nb_bench = atoi(argv[2]); + START_BENCH(nb_bench); + salembier<V> s(input, nbh); + parent = s.parent; + END_BENCH("Salembier: "); +} + + + + diff --git a/milena/sandbox/edwin/TreeAlgorithmsComp/tbb_union_find_8.cc b/milena/sandbox/edwin/TreeAlgorithmsComp/tbb_union_find_8.cc index 4de9132..1b71abd 100644 --- a/milena/sandbox/edwin/TreeAlgorithmsComp/tbb_union_find_8.cc +++ b/milena/sandbox/edwin/TreeAlgorithmsComp/tbb_union_find_8.cc @@ -1,4 +1,7 @@ +#undef MLN_WO_GLOBAL_VARS + #include <tbb/task_scheduler_init.h> + #include "parallel_union_find.hh" #include <mln/core/image/image2d.hh> @@ -6,54 +9,37 @@ #include <mln/value/int_u8.hh> #include <mln/io/pgm/load.hh> -#include <mln/debug/println.hh> - - namespace mln - { - - template <typename I, typename J> - unsigned nnodes(const I& f, const J& parent) - { - mlc_equal(mln_value(J), unsigned)::check(); - unsigned n = 0; - mln_piter(I) p(parent.domain()); - for_all(p) - { - unsigned o = parent.index_of_point(p); - if (parent(p) == o || f.element(parent(p)) != f(p)) - ++n; - } - return n; - } - - } +#include "nnodes.hh" + +#undef BENCH_START +#undef BENCH_END +#include "bench.hh" void usage(char* argv[]) { - std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl; + std::cerr << "usage: " << argv[0] << " input.pgm nbExec [nbThread]" << std::endl; std::cerr << " max-tree parallel" << std::endl; std::abort(); } int main(int argc, char* argv[]) { - if (argc != 2) + if (argc < 3) usage(argv); - tbb::task_scheduler_init init; - using namespace mln; + int nthread = (argc == 4) ? std::atoi(argv[3]) : tbb::task_scheduler_init::automatic; + tbb::task_scheduler_init init(nthread); + using namespace mln; neighb2d nbh = c4(); image2d<value::int_u8> f; io::pgm::load(f, argv[1]); image2d<unsigned> parent = compute_max_tree(f, nbh); - // debug::println(parent); - - // tm.start(); - // START_BENCH(1); - // //compute_max_tree(f, nbh); - // END_BENCH(); std::cout << "nnodes = " << nnodes(f, parent) << std::endl; + + START_BENCH(std::atoi(argv[2])); + parent = compute_max_tree(f, nbh); + END_BENCH("Parallel union find: "); } -- 1.5.6.5