https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud <thierry.geraud@lrde.epita.fr> Add some bench + canvas + subsampling + browsing code. * bench: New directory. * bench/input_iz.pgm.gz: New. * bench/fastest_statistical_tour_browsing.cc: New. * bench/fastest_forall_p_browsing.cc: New. * bench/iz.cc: New. * bench/fastest_statistical_tour_nbh_browsing.cc: New. * bench/z_sub_browsing: New. * bench/z_sub_browsing/fast.cc: New. * bench/z_sub_browsing/debase.hh: New. * bench/z_sub_browsing/integral.hh: New. * bench/z_sub_browsing/+inc: New. * bench/z_sub_browsing/in.pgm.gz: New. * bench/z_sub_browsing/debase.cc: New. * bench/z_sub_browsing/integral.cc: New. * bench/z_sub_browsing/README: New. * theo/mln/morpho/canvas/lena_blurred.pgm.gz: New. * theo/mln/morpho/canvas/lena.pgm.gz: New. * theo/mln/morpho/canvas/g.pbm.gz: New. * theo/mln/morpho/canvas/lena_min.pgm.gz: New. * theo/mln/morpho/canvas/reconstruction_on_set.hh: Layout. * theo/mln/morpho/canvas/f_and_g.pbm.gz: New. * theo/mln/morpho/canvas/regminid.pbm.gz: New. * theo/mln/morpho/canvas/one_domain.cc: New. * theo/mln/subsampling: New directory. * theo/mln/subsampling/sizes.cc: New. * theo/mln/subsampling/debase.hh: New. * theo/mln/subsampling/integral.hh: New. * theo/mln/subsampling/in.pgm.gz: New. * theo/mln/subsampling/debase.cc: New. * theo/mln/subsampling/integral.cc: New. * theo/mln/browsing: New directory. * theo/mln/browsing/window_sliding.cc: New. bench/fastest_forall_p_browsing.cc | 384 ++++++++++++++++ bench/fastest_statistical_tour_browsing.cc | 168 +++++++ bench/fastest_statistical_tour_nbh_browsing.cc | 186 ++++++++ bench/iz.cc | 386 ++++++++++++++++ bench/z_sub_browsing/README | 5 bench/z_sub_browsing/debase.cc | 23 + bench/z_sub_browsing/debase.hh | 351 +++++++++++++++ bench/z_sub_browsing/fast.cc | 139 ++++++ bench/z_sub_browsing/integral.cc | 23 + bench/z_sub_browsing/integral.hh | 165 +++++++ theo/mln/browsing/window_sliding.cc | 52 ++ theo/mln/morpho/canvas/one_domain.cc | 569 +++++++++++++++++++++++++ theo/mln/subsampling/debase.cc | 23 + theo/mln/subsampling/debase.hh | 351 +++++++++++++++ theo/mln/subsampling/integral.cc | 38 + theo/mln/subsampling/integral.hh | 295 ++++++++++++ theo/mln/subsampling/sizes.cc | 45 + 17 files changed, 3203 insertions(+) Index: bench/input_iz.pgm.gz Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: bench/input_iz.pgm.gz ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Index: bench/fastest_statistical_tour_browsing.cc --- bench/fastest_statistical_tour_browsing.cc (revision 0) +++ bench/fastest_statistical_tour_browsing.cc (revision 0) @@ -0,0 +1,168 @@ +// Copyright (C) 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. + +#include <mln/core/image/image2d.hh> +#include <mln/pw/all.hh> +#include <mln/data/compare.hh> +#include <mln/debug/println.hh> +#include <mln/util/timer.hh> + + + +# define loop(n) for (unsigned i = 0; i < n; ++i) + + + +namespace mln +{ + + + // test_ + + template <typename I> + void test_(const std::string& routine, float t, const I& ima, mln_value(I) n) + { + std::cout << routine << ' ' << t << std::endl; + if (ima != (pw::cst(n) | ima.domain())) + std::cerr << "bug in " << routine << '!' << std::endl; + } + + + + // piter + + template <typename I> + void piter_based(I& ima) + { + mln_piter(I) p(ima.domain()); + for_all(p) + if ((p.row() + p.col()) % 2 == 0) + ++ima(p); + for_all(p) + if ((p.row() + p.col()) % 2 == 1) + ++ima(p); + } + + + // row_col + + template <typename I> + void row_col_based(I& ima) + { + const unsigned nrows = ima.nrows(), ncols = ima.ncols(); + for (int row = 0; row < nrows; ++row) + for (int col = row % 2; col < ncols; col += 2) + ++ima.at_(row, col); + for (int row = 0; row < nrows; ++row) + for (int col = ! (row % 2); col < ncols; col += 2) + ++ima.at_(row, col); + } + + + // run_ptr + + template <typename I> + void run_ptr_based(I& ima) + { + mln_box_runstart_piter(I) s(ima.domain()); + const unsigned n = s.run_length(); + + // half-tour 0 + unsigned skip = 0; + for_all(s) + { + mln_value(I)* ptr = &ima(s) + skip; + for (unsigned i = 0; i < n; i += 2) + { + ++*ptr; + ptr += 2; + } + skip = ! skip; + } + + // half-tour 1 + skip = 1; + for_all(s) + { + mln_value(I)* ptr = &ima(s) + skip; + for (unsigned i = 0; i < n; i += 2) + { + ++*ptr; + ptr += 2; + } + skip = ! skip; + } + } + + +} // end of namespace mln + + + + + +int main() +{ + using namespace mln; + + const unsigned n = 10, size = 1000; + image2d<char> ima(size, size); + util::timer t; + std::string routine; + + + // piter + + { + routine = "piter"; + data::fill(ima, 0); + t.start(); + loop(n) piter_based(ima); + test_(routine, t.stop(), ima, n); + } + + + // row_col + + { + routine = "row_col"; + data::fill(ima, 0); + t.start(); + loop(n) row_col_based(ima); + test_(routine, t.stop(), ima, n); + } + + + // run_ptr + + { + routine = "run_ptr"; + data::fill(ima, 0); + t.start(); + loop(n) run_ptr_based(ima); + test_(routine, t.stop(), ima, n); + } + +} Index: bench/fastest_forall_p_browsing.cc --- bench/fastest_forall_p_browsing.cc (revision 0) +++ bench/fastest_forall_p_browsing.cc (revision 0) @@ -0,0 +1,384 @@ +// Copyright (C) 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. + +#include <mln/core/image/image2d.hh> +#include <mln/util/timer.hh> + + + +# define loop(n) for (unsigned i = 0; i < n; ++i) + + + +namespace mln +{ + + + // (r, c) + + template <typename I> + void bench_00(I& ima) + { + const unsigned nr = ima.nrows(), nc = ima.ncols(); + for (int r = 0; r < nr; ++r) + for (int c = 0; c < nc; ++c) + ima.at_(r, c) = 0; + } + + + // ptr + + template <typename I> + void bench_0(I& ima) + { + const unsigned n = ima.nelements(); + mln_value(I)* ptr = ima.buffer(); + for (unsigned i = 0; i < n; ++i) + *ptr++ = 0; + } + + template <typename I> + void bench_0(I& ima1, const I& ima2) + { + const unsigned n = ima1.nelements(); + mln_value(I)* ptr1 = ima1.buffer(); + const mln_value(I)* ptr2 = ima2.buffer(); + for (unsigned i = 0; i < n; ++i) + *ptr1++ = *ptr2++; + } + + + // piter + + template <typename I> + void bench_1(I& ima) + { + mln_piter(I) p(ima.domain()); + for_all(p) + ima(p) = 0; + } + + + // pixter + + template <typename I> + void bench_2(I& ima) + { + mln_fwd_pixter(I) pxl(ima); + for_all(pxl) + pxl.val() = 0; + } + + template <typename I> + void bench_2(I& ima1, const I& ima2) + { + mln_fwd_pixter(I) pxl1(ima1); + mln_fwd_pixter(const I) pxl2(ima2); + for_all_2(pxl1, pxl2) + pxl1.val() = pxl2.val(); + } + + + // pixter -> offset + + template <typename I> + void bench_3(I& ima) + { + mln_fwd_pixter(I) pxl(ima); + for_all(pxl) + { + unsigned p = pxl.offset(); + ima.element(p) = 0; + } + } + + template <typename I> + void bench_3(I& ima1, const I& ima2) + { + // border::equalize(ima1, ima2); + mln_fwd_pixter(I) pxl(ima1); + for_all(pxl) + { + unsigned p = pxl.offset(); + ima1.element(p) = ima2.element(p); + } + } + + + // pixter -> point + + template <typename I> + void bench_4(I& ima) + { + mln_fwd_pixter(I) pxl(ima); + for_all(pxl) + { + unsigned o = pxl.offset(); + mln_psite(I) p = ima.point_at_index(o); + ima(p) = 0; + } + } + + + + // box_runstart_piter -> ptr + + template <typename I> + void bench_5(I& ima) + { + mln_box_runstart_piter(I) s(ima.domain()); + const unsigned n = s.run_length(); + for_all(s) + { + mln_value(I)* ptr = & ima(s); + for (unsigned i = 0; i < n; ++i) + { + *ptr++ = 0; + } + } + } + + + // box_runstart_piter -> (ptr, point) + + template <typename I> + void bench_6(I& ima) + { + typedef mln_site(I) P; + const unsigned dim_1 = P::dim - 1; + mln_box_runstart_piter(I) s(ima.domain()); + const unsigned n = s.run_length(); + for_all(s) + { + mln_value(I)* ptr = & ima(s); + P p = s; + for (unsigned i = 0; i < n; ++i) + { + *ptr++ = 0; + ++p[dim_1]; + } + } + } + + template <typename I> + void bench_6(I& ima1, const I& ima2) + { + typedef mln_site(I) P; + const unsigned dim_1 = P::dim - 1; + mln_box_runstart_piter(I) s(ima1.domain()); + const unsigned n = s.run_length(); + for_all(s) + { + mln_value(I)* ptr1 = & ima1(s); + const mln_value(I)* ptr2 = & ima2(s); + P p = s; + for (unsigned i = 0; i < n; ++i) + { + *ptr1++ = *ptr2++; + ++p[dim_1]; + } + } + } + + +// template <typename I> +// void bench_6_alt(I& ima1, const I& ima2) +// { +// typedef mln_site(I) P; +// const unsigned dim_1 = P::dim - 1; +// mln_fwd_by_run_piter(I) p(ima1.domain()); +// for_all(p) +// { +// if (p.at_runstart()) +// { +// p.update(ptr1, ima1); +// p.update(ptr2, ima2); +// } +// *ptr1++ = *ptr2++; +// p; // converts to point2d +// } +// } + + + + // box_runstart_piter -> (offset, point) + + template <typename I> + void bench_7(I& ima) + { + typedef mln_site(I) P; + const unsigned dim_1 = P::dim - 1; + mln_box_runstart_piter(I) s(ima.domain()); + const unsigned n = s.run_length(); + for_all(s) + { + unsigned offset = ima.index_of_point(s); + P p = s; + for (unsigned i = 0; i < n; ++i) + { + ima.element(offset) = 0; + ++offset; + ++p[dim_1]; + } + } + } + + + +} // end of namespace mln + + + + + +int main() +{ + using namespace mln; + + const unsigned n = 100, size = 1000; + image2d<int> ima(size, size), ima2(size, size); + util::timer t; + + + // (r, c) + + { + t.start(); + loop(n) bench_00(ima); + float ts = t.stop(); + std::cout << "(r, c) " << ts << std::endl; + } + + + // ptr + + { + t.start(); + loop(n) bench_0(ima); + float ts = t.stop(); + std::cout << "ptr " << ts << std::endl; + } + +// { +// t.start(); +// loop(n) bench_0(ima, ima2); +// float ts = t.stop(); +// std::cout << "2 ptrs " << ts << std::endl; +// } + + + // piter + + { + t.start(); + loop(n) bench_1(ima); + std::cout << "piter " << t.stop() << std::endl; + } + + + // pixter + + { + t.start(); + loop(n) bench_2(ima); + float ts = t.stop(); + std::cout << "pixter " << ts << std::endl; + } + +// { +// t.start(); +// loop(n) bench_2(ima, ima2); +// float ts = t.stop(); +// std::cout << "2 pixters " << ts << std::endl; +// } + + + // pixter -> offset + + { + t.start(); + loop(n) bench_3(ima); + std::cout << "pixter -> offset " << t.stop() << std::endl; + } + +// { +// t.start(); +// loop(n) bench_3(ima, ima2); +// std::cout << t.stop() << std::endl; +// } + + + // pixter -> point + + { + t.start(); + loop(n) bench_4(ima); + std::cout << "pixter -> point " << t.stop() << std::endl; + } + + + // box_runstart_piter -> ptr + + { + t.start(); + loop(n) bench_5(ima); + std::cout << "box_runstart_piter -> ptr " << t.stop() << std::endl; + } + + + // box_runstart_piter -> (ptr, point) + + { + t.start(); + loop(n) bench_6(ima); + float ts = t.stop(); + std::cout << "runstart -> (ptr, point) " << ts << std::endl; + } + +// { +// t.start(); +// loop(n) bench_6(ima, ima2); +// float ts = t.stop(); +// std::cout << "runstart -> (2 ptrs, point) " << ts << std::endl; +// } + + + // box_runstart_piter -> (offset, point) + + { + t.start(); + loop(n) bench_7(ima); + float ts = t.stop(); + std::cout << "runstart -> (offset, point) " << ts << std::endl; + } + +// { +// t.start(); +// loop(n) bench_7(ima, ima2); +// float ts = t.stop(); +// std::cout << "runstart -> (2 offsets, point) " << ts << std::endl; +// } + +} Index: bench/iz.cc --- bench/iz.cc (revision 0) +++ bench/iz.cc (revision 0) @@ -0,0 +1,386 @@ +// Copyright (C) 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. + +#include <queue> +#include <mln/core/routine/duplicate.hh> +#include <mln/data/fill.hh> +#include <mln/extension/fill.hh> + +#include <mln/transform/influence_zone_geodesic.hh> + +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/util/timer.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + + + +# define loop(n) for (unsigned i = 0; i < n; ++i) + + + +namespace mln +{ + + + template <typename I, typename N> + mln_concrete(I) do_iz(const I& input, const N& nbh) + { + mln_ch_value(I, unsigned) dmap; + std::queue<unsigned> q; + mln_concrete(I) output; + + // Initialization. + { + output = duplicate(input); // <-- init + extension::fill(output, 1); // For the extension to be ignored. + + mln_pixter(const I) p(input); + mln_nixter(const I, N) n(p, nbh); + for_all(p) + if (p.val() != 0) // <-- inqueue_p_wrt_input_p + { + ; // <-- init_p + // FIXME: dmap.element(p.offset()) = 0; + for_all(n) + if (n.val() == 0) // <-- inqueue_p_wrt_input_n + { + q.push(p.offset()); + break; + } + } + } + + // Propagation. + { + unsigned p; + + util::array<int> dp = offsets_wrt(input, nbh); + const unsigned n_nbhs = dp.nelements(); + + while (! q.empty()) + { + p = q.front(); + q.pop(); + // pas de saturation + if (output.element(p) == 0) + std::cerr << "oops!" << std::endl; + for (unsigned i = 0; i < n_nbhs; ++i) + { + unsigned n = p + dp[i]; + if (output.element(n) == 0) + { + output.element(n) = output.element(p); // <- process + q.push(n); + } + } + } + } + + return output; + } + + + + template <typename I, typename N> + mln_concrete(I) do_iz__bis(const I& input, const N& nbh) + { + std::queue<unsigned> q; + mln_concrete(I) output; + + util::array<int> dp = offsets_wrt(input, nbh); + const unsigned n_nbhs = dp.nelements(); + + // Initialization. + { + output = duplicate(input); // <-- init + extension::fill(output, 1); // For the extension to be ignored. + + mln_box_runstart_piter(I) s(input.domain()); + const unsigned len = s.run_length(); + for_all(s) + { + unsigned p = input.index_of_point(s); + for (unsigned i = 0; i < len; ++i, ++p) + if (input.element(p) != 0) + { + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p + dp[j]; + if (input.element(n) == 0) + { + q.push(p); + break; + } + } + } + } + } + + // Propagation. + { + unsigned p; + + while (! q.empty()) + { + p = q.front(); + q.pop(); + mln_invariant(output.element(p) != 0); + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p + dp[j]; + if (output.element(n) == 0) + { + output.element(n) = output.element(p); // <- process + q.push(n); + } + } + } + } + + return output; + } + + + + + + template <typename I, typename N> + mln_concrete(I) do_iz__ter(const I& input, const N& nbh) + { + std::queue<unsigned> q; + mln_concrete(I) output; + + util::array<int> dp = offsets_wrt(input, nbh); + const unsigned n_nbhs = dp.nelements(); + const unsigned + skip_border = 2 * input.border(), + nrows = input.nrows(), + ncols = input.ncols(); + + // Initialization. + { + output = duplicate(input); // <-- init + extension::fill(output, 1); // For the extension to be ignored. + + unsigned p = input.index_of_point(point2d(0,0)); + for (unsigned row = 0; row < nrows; ++row, p += skip_border) + { + for (unsigned i = 0; i < ncols; ++i, ++p) + if (input.element(p) != 0) + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p + dp[j]; + if (input.element(n) == 0) + { + q.push(p); + break; + } + } + } + } + + // Propagation. + { + unsigned p; + + while (! q.empty()) + { + p = q.front(); + q.pop(); + mln_invariant(output.element(p) != 0); + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p + dp[j]; + if (output.element(n) == 0) + { + output.element(n) = output.element(p); // <- process + q.push(n); + } + } + } + } + + return output; + } + + + + + template <typename I, typename N> + mln_concrete(I) do_iz__ptr(const I& input, const N& nbh) + { + std::queue<mln_value(I)*> q; + mln_concrete(I) output; + + util::array<int> dp = offsets_wrt(input, nbh); + const unsigned n_nbhs = dp.nelements(); + const unsigned + skip_border = 2 * input.border(), + nrows = input.nrows(), + ncols = input.ncols(); + + // Initialization. + { + output = duplicate(input); + // For the extension to be ignored: + extension::fill(input, 0); // in initialization + extension::fill(output, 1); // in propagation + + const unsigned nelts = input.nelements(); + const mln_value(I)* p_i = & input.at_(0, 0); + mln_value(I)* p_o = & output.at_(0, 0); + for (unsigned i = 0; i < nelts; ++i, ++p_i, ++p_o) + { + if (*p_i == 0) + continue; + for (unsigned j = 0; j < n_nbhs; ++j) + { + const mln_value(I)* n_i = p_i + dp[j]; + if (*n_i == 0) + { + q.push(p_o); + break; + } + } + } + + } + + // Propagation. + { + mln_value(I)* ptr; + + while (! q.empty()) + { + ptr = q.front(); + q.pop(); + mln_invariant(*ptr != 0); + for (unsigned j = 0; j < n_nbhs; ++j) + { + mln_value(I)* ntr = ptr + dp[j]; + if (*ntr == 0) + { + *ntr = *ptr; + q.push(ntr); + } + } + } + } + + return output; + } + + + +} // end of namespace mln + + + + + +int main() +{ + using namespace mln; + using value::int_u8; + + image2d<int_u8> in; + io::pgm::load(in, "input_iz.pgm"); + + util::timer t; + neighb2d nbh = c4(); + +// { +// t.start(); +// image2d<int_u8> ref = transform::influence_zone_geodesic(in, nbh, mln_max(unsigned)); +// float ts = t.stop(); +// std::cout << ts << std::endl; +// io::pgm::save(ref, "ref.pgm"); +// } + +// { +// t.start(); +// image2d<int_u8> out = do_iz__bis(in, nbh); +// float ts = t.stop(); +// std::cout << ts << std::endl; +// io::pgm::save(out, "iz_bis.pgm"); +// } + +// { +// t.start(); +// image2d<int_u8> out = do_iz__ter(in, nbh); +// float ts = t.stop(); +// std::cout << ts << std::endl; +// io::pgm::save(out, "iz_ter.pgm"); +// } + +// { +// t.start(); +// image2d<int_u8> out = do_iz__ptr(in, nbh); +// float ts = t.stop(); +// std::cout << ts << std::endl; +// io::pgm::save(out, "iz_ptr.pgm"); +// } + + + const unsigned n_times = 10; + + { + t.start(); + loop(n_times) transform::influence_zone_geodesic(in, nbh, mln_max(unsigned)); + float ts = t.stop(); + std::cout << ts << std::endl; + } + + { + t.start(); + loop(n_times) do_iz(in, nbh); + float ts = t.stop(); + std::cout << ts << std::endl; + } + + { + t.start(); + loop(n_times) do_iz__bis(in, nbh); + float ts = t.stop(); + std::cout << ts << std::endl; + } + + { + t.start(); + loop(n_times) do_iz__ter(in, nbh); + float ts = t.stop(); + std::cout << ts << std::endl; + } + + { + t.start(); + loop(n_times) do_iz__ptr(in, nbh); + float ts = t.stop(); + std::cout << ts << std::endl; + } + +} Index: bench/fastest_statistical_tour_nbh_browsing.cc --- bench/fastest_statistical_tour_nbh_browsing.cc (revision 0) +++ bench/fastest_statistical_tour_nbh_browsing.cc (revision 0) @@ -0,0 +1,186 @@ +// Copyright (C) 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. + +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/pw/all.hh> +#include <mln/data/compare.hh> +#include <mln/debug/println.hh> +#include <mln/util/timer.hh> + + + +# define loop(n) for (unsigned i = 0; i < n; ++i) + + + +namespace mln +{ + + + // test_ + + template <typename I> + void test_(const std::string& routine, float t, const I& ima, unsigned n) + { + dpoint2d dp(1, 1); + box2d + b = ima.domain(), + inner(b.pmin() + dp, b.pmax() - dp); + std::cout << routine << ' ' << t << std::endl; + if ((ima | inner) != (pw::cst(4 * n) | inner)) + std::cerr << "bug in " << routine << '!' << std::endl; + // FIXME: other tests... + } + + + + // piter + + template <typename I, typename N> + void piter_based(I& ima, const N& nbh) + { + mln_piter(I) p(ima.domain()); + mln_niter(N) n(nbh, p); + for_all(p) + if ((p.row() + p.col()) % 2 == 0) + for_all(n) + ++ima(n); + for_all(p) + if ((p.row() + p.col()) % 2 == 1) + for_all(n) + ++ima(n); + } + + + // row_col + + template <typename I, typename N> + void row_col_based(I& ima, const N& nbh) + { + const unsigned nrows = ima.nrows(), ncols = ima.ncols(); + unsigned n_nbhs = nbh.size(); + for (int row = 0; row < nrows; ++row) + for (int col = row % 2; col < ncols; col += 2) + for (unsigned j = 0; j < n_nbhs; ++j) + ++ima.at_(row + nbh.dp(j).row(), col + nbh.dp(j).col()); + for (int row = 0; row < nrows; ++row) + for (int col = ! (row % 2); col < ncols; col += 2) + for (unsigned j = 0; j < n_nbhs; ++j) + ++ima.at_(row + nbh.dp(j).row(), col + nbh.dp(j).col()); + } + + + // run_ptr + + template <typename I, typename N> + void run_ptr_based(I& ima, const N& nbh) + { + mln_box_runstart_piter(I) s(ima.domain()); + const unsigned n = s.run_length(); + + util::array<int> dp = offsets_wrt(ima, nbh); + const unsigned n_nbhs = dp.nelements(); + + // half-tour 0 + unsigned skip = 0; + for_all(s) + { + mln_value(I)* ptr = &ima(s) + skip; + for (unsigned i = 0; i < n; i += 2) + { + for (unsigned j = 0; j < n_nbhs; ++j) + ++*(ptr + dp(j)); + ptr += 2; + } + skip = ! skip; + } + + // half-tour 1 + skip = 1; + for_all(s) + { + mln_value(I)* ptr = &ima(s) + skip; + for (unsigned i = 0; i < n; i += 2) + { + for (unsigned j = 0; j < n_nbhs; ++j) + ++*(ptr + dp(j)); + ptr += 2; + } + skip = ! skip; + } + } + + +} // end of namespace mln + + + + + +int main() +{ + using namespace mln; + + const unsigned n = 10, size = 1000; + image2d<char> ima(size, size); + neighb2d nbh = c4(); + util::timer t; + std::string routine; + + + // piter + + { + routine = "piter"; + data::fill(ima, 0); + t.start(); + loop(n) piter_based(ima, nbh); + test_(routine, t.stop(), ima, n); + } + + + // row_col + + { + routine = "row_col"; + data::fill(ima, 0); + t.start(); + loop(n) row_col_based(ima, nbh); + test_(routine, t.stop(), ima, n); + } + + + // run_ptr + + { + routine = "run_ptr"; + data::fill(ima, 0); + t.start(); + loop(n) run_ptr_based(ima, nbh); + test_(routine, t.stop(), ima, n); + } + +} Index: bench/z_sub_browsing/fast.cc --- bench/z_sub_browsing/fast.cc (revision 0) +++ bench/z_sub_browsing/fast.cc (revision 0) @@ -0,0 +1,139 @@ +#include <mln/core/image/image2d.hh> +#include <mln/util/array.hh> +#include <mln/util/timer.hh> +#include <mln/data/fill.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + + +namespace mln +{ + +// template <typename I> +// void +// test1(const I& input, I& out) +// { +// unsigned s = 2; +// util::timer timer_; +// timer_.start(); + +// mln_box_runstart_piter(I) sp(out.domain()); + +// const unsigned n = sp.run_length(); + +// unsigned row_offset = input.domain().ncols() +// + 2 * input.border(); + +// util::array< mln_value(I)*> ptr(s); +// util::array<const mln_value(I)*> ptr_input(s); + +// sp.start(); +// while (sp.is_valid()) +// { +// mln_site(I) +// site = sp; // Site at scale 1. + +// unsigned nptr = 0; +// for (; sp.is_valid() && nptr < ptr.size(); ++nptr) +// { +// ptr(nptr) = & out(sp); +// ptr_input(nptr) = & input(sp); +// sp.next(); +// } +// ptr.resize(nptr); + +// for (unsigned col = 0; col < n; col += s, site[1] += s) +// { +// for (unsigned i = 0; i < ptr.size() ; ++i) +// for (unsigned j = 0; j < s; ++j) +// *ptr(i)++ = *ptr_input(i)++; +// } +// } + +// float t_ = timer_; +// std::cout << "Test 1 - " << t_ << std::endl; +// } + + + template <typename I> + void + test2(const I& input, I& out, unsigned s) + { + const unsigned s_2 = s * s, round_it = s_2 / 2; + util::timer timer_; + timer_.start(); + + typedef mln_value(I) V; + typedef mln_sum(V) S; + + mln_box_runstart_piter(I) sp(input.domain()); + + const unsigned n = sp.run_length(); + + int row_offset = input.domain().ncols() + 2 * input.border(); + int row_offset_s = row_offset - s; + int row_offset_next = - (s - 1) * row_offset; + + mln_value(I)* ptr = out.buffer(); + const mln_value(I)* ptr_input; + + sp.start(); + while (sp.is_valid()) + { + mln_site(I) + site = sp; // Site at scale 1. + + ptr_input = & input(sp); + + for (unsigned col = 0; col < n; col += s, site[1] += s) + { + + // tile + S sum = 0; + for (unsigned i = 0; i < s;) + { + for (unsigned j = 0; j < s; ++j) + sum += *ptr_input++; + + ++i; + if (i < s) + ptr_input += row_offset_s; + } + + ptr_input += row_offset_next; + *ptr++ = sum / s_2; + } + + for (unsigned i = 0; i < s; ++i) + sp.next(); + } + + float t_ = timer_; + std::cout << t_ << std::endl; + } + + + +} // end of namespace mln + +int main() +{ + using namespace mln; + using value::int_u8; + + def::coord + r = 3144, + c = 2252; + + image2d<int_u8> input; + io::pgm::load(input, "in.pgm"); + + unsigned s = 3; + box2d b = make::box2d((input.nrows() + s - 1) / s, + (input.ncols() + s - 1) / s); + image2d<int_u8> output(b, 0); + test2(input, output, s); + + io::pgm::save(output, "out.pgm"); + +} Index: bench/z_sub_browsing/debase.hh --- bench/z_sub_browsing/debase.hh (revision 0) +++ bench/z_sub_browsing/debase.hh (revision 0) @@ -0,0 +1,351 @@ +// Copyright (C) 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 MLN_SUBSAMPLING_DEBASE_HH +# define MLN_SUBSAMPLING_DEBASE_HH + +/// \file +/// +/// Debase subsampling. + +#include <mln/core/concept/image.hh> + +#define MLN_FLOAT double +#include <sandbox/theo/exec/gaussian_directional_2d.hh> + + + +namespace mln +{ + + namespace subsampling + { + + /// FIXME: Doc. + + template <typename I> + mln_concrete(I) + debase(const Image<I>& input, + unsigned gap, + const mln_deduce(I, site, delta)& shift); // FIXME: Add round_up_size. + + + template <typename I> + mln_concrete(I) + debase(const Image<I>& input, unsigned gap); + + + +# ifndef MLN_INCLUDE_ONLY + + + // Implementation. + + namespace impl + { + + + template <typename I> + inline + mln_concrete(I) + debase(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase"); + + const I& input = exact(input_); + + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + + box<P> b(pmin, pmax); + mln_concrete(I) output(b); + mln_piter(I) p(output.domain()); + for_all(p) + output(p) = input(p * gap + shift); + + trace::exiting("subsampling::impl::debase"); + return output; + } + + + template <typename I> + inline + mln_concrete(I) + debase_fastest(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase_fastest"); + + const I& input = exact(input_); + + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + + box<P> b(pmin, pmax); + + typedef mln_concrete(I) O; + O output(b); + + mln_box_runstart_piter(O) s(output.domain()); + const unsigned n = s.run_length(); + for_all(s) + { + const mln_value(I)* pi = & input(s * gap + shift); + mln_value(O)* po = & output(s); + for (unsigned i = 0; i < n; ++i) + { + *po++ = *pi; + pi += gap; + } + } + + trace::exiting("subsampling::impl::debase_fastest"); + return output; + } + + + template <typename I> + inline + mln_concrete(I) + debase_gaussian_antialiased_fastest(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase_gaussian_antialiased_fastest"); + + const I& input = exact(input_); + + image2d<MLN_FLOAT> temp(input.domain()); + data::fill(temp, input); + temp = linear::gaussian_directional_2d(temp, 1, 1.5, 0); + temp = linear::gaussian_directional_2d(temp, 0, 1.5, 0); + + typedef mln_value(I) V; + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + + box<P> b(pmin, pmax); + + typedef mln_concrete(I) O; + O output(b); + + mln_box_runstart_piter(O) s(output.domain()); + const unsigned n = s.run_length(); + for_all(s) + { + const MLN_FLOAT* pi = & temp(s * gap + shift); + mln_value(O)* po = & output(s); + for (unsigned i = 0; i < n; ++i) + { + *po++ = V(*pi + 0.49); + pi += gap; + } + } + + trace::entering("subsampling::impl::debase_gaussian_antialiased_fastest"); + return output; + } + + + + + template <typename I> + inline + mln_concrete(I) + debase_2d_antialias_pointer_fastest(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase_2d_antialias_pointer_fastest"); + + const I& input = exact(input_); + + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + + box<P> b(pmin, pmax); + + typedef mln_concrete(I) O; + O output(b); + + const unsigned nrows = input.nrows(); + const unsigned ncols = input.ncols(); + + typedef mln_value(I) V; + typedef mln_sum(V) S; + for (unsigned row = 0; row < nrows; row += gap) + { + const V* ptr1 = & input.at_(row, 0); + const V* ptr2 = & input.at_(row + 1, 0); + const V* ptr3 = & input.at_(row + 2, 0); + for (unsigned col = 0; col < ncols; col += gap) + { + S s; + + // mean + s = *ptr1++ + *ptr1++ + *ptr1++; + s += *ptr2++ + *ptr2++ + *ptr2++; + s += *ptr3++ + *ptr3++ + *ptr3++; + output.at_(row / gap, col / gap) = (s + 4) / 9; +// // cut +// s = 1 * *ptr1++ + 2 * *ptr1++ + 1 * *ptr1++; +// s += 2 * *ptr2++ + 3 * *ptr2++ + 2 * *ptr2++; +// s += 1 * *ptr3++ + 2 * *ptr3++ + 1 * *ptr3++; +// output.at_(row / gap, col / gap) = (s + 7) / 15; + } + } + + trace::exiting("subsampling::impl::debase_2d_antialias_pointer_fastest"); + return output; + } + + + + template <typename I> + inline + mln_concrete(I) + debase_2d_antialias_offset_fastest(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase_2d_antialias_offset_fastest"); + + const I& input = exact(input_); + + typedef mln_value(I) V; + typedef mln_sum(V) S; + + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + box<P> b(pmin, pmax); + typedef mln_concrete(I) O; + O output(b); + + const V* ptr1 = & input.at_(0, 0); + const V* ptr2 = & input.at_(1, 0); + const V* ptr3 = & input.at_(2, 0); + + mln_box_runstart_piter(O) s(output.domain()); + const unsigned n = s.run_length(); + + unsigned offset = input.delta_index(point2d(3,0) - point2d(0,3*n)); + + for_all(s) + { + mln_value(O)* po = & output(s); + for (unsigned i = 0; i < n; ++i) + { + S s; + s = *ptr1++ + *ptr1++ + *ptr1++; + s += *ptr2++ + *ptr2++ + *ptr2++; + s += *ptr3++ + *ptr3++ + *ptr3++; + *po++ = (s + 4) / 9; + } + ptr1 += offset; + ptr2 += offset; + ptr3 += offset; + } + + trace::exiting("subsampling::impl::debase_2d_antialias_offset_fastest"); + return output; + } + + + } // end of namespace mln::subsampling::impl + + + + // Dispatch. + + namespace internal + { + + template <typename I> + inline + mln_concrete(I) + debase_dispatch(const Image<I>& input, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + return impl::debase_2d_antialias_pointer_fastest /* debase_2d_antialias_offset_fastest */(input, gap, shift); + } + + + } // end of namespace mln::subsampling::internal + + + + // Facades. + + template <typename I> + inline + mln_concrete(I) + debase(const Image<I>& input, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::debase"); + + typedef mln_site(I) P; + + mln_precondition(exact(input).is_valid()); + mln_precondition(exact(input).domain().pmin() == literal::origin); + mln_precondition(gap > 1); + for (unsigned i = 0; i < P::dim; ++i) + mln_precondition(shift[i] < gap); + + mln_concrete(I) output; + output = internal::debase_dispatch(input, gap, shift); + + trace::exiting("subsampling::debase"); + return output; + } + + + template <typename I> + inline + mln_concrete(I) + debase(const Image<I>& input, unsigned gap) + { + return debase(input, gap, literal::zero); + } + + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::subsampling + +} // end of namespace mln + + +#endif // ! MLN_SUBSAMPLING_DEBASE_HH Index: bench/z_sub_browsing/integral.hh --- bench/z_sub_browsing/integral.hh (revision 0) +++ bench/z_sub_browsing/integral.hh (revision 0) @@ -0,0 +1,165 @@ +// Copyright (C) 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 MLN_SUBSAMPLING_INTEGRAL_HH +# define MLN_SUBSAMPLING_INTEGRAL_HH + +/// \file +/// +/// Both subsampling and integral image computation. + +#include <mln/core/concept/image.hh> + + + +namespace mln +{ + + namespace subsampling + { + + /// FIXME: Doc. + + template <typename I> + mln_concrete(I) + integral(const Image<I>& input, unsigned scale); + + + +# ifndef MLN_INCLUDE_ONLY + + + // Implementation. + + namespace impl + { + + + template <unsigned scale, typename I> + inline + mln_concrete(I) + integral_(const Image<I>& input_) + { + trace::entering("subsampling::impl::integral_"); + + const I& input = exact(input_); + const unsigned scale2 = scale * scale, round_it = scale2 / 2; + + typedef mln_value(I) V; + typedef mln_sum(V) S; + typedef mln_site(I) P; + + box<P> b = make::box2d(input.nrows() / scale, input.ncols() / scale); + const unsigned nrows = b.nrows(), ncols = b.ncols(); + + typedef mln_concrete(I) O; + O output(b, 0); // no external border in 'output' + V* po = output.buffer(); + + unsigned + tile_cr = input.delta_index(dpoint2d(scale, - ncols * scale)), + cr = input.delta_index(dpoint2d(1, - (int)scale)); + + const V* pi = & input.at_(0, 0); // first point of tiles + for (unsigned row = 0; row < nrows; ++row) + { + S hsum = 0, hsum2 = 0; // horizontal integral + for (unsigned col = 0; col < ncols; ++col) + { + S sum = 0, sum2 = 0; + const V* ptr = pi; + + for (unsigned r = 0; r < scale; ++r) + { + for (unsigned c = 0; c < scale; ++c) + { + sum += *ptr; + // sum2 += *ptr * *ptr; + ++ptr; + } + ptr += cr; + } + + // ... + pi += scale; + *po++ = (round_it + sum) / scale2; + } + pi += tile_cr; + } + + trace::exiting("subsampling::impl::integral_"); + return output; + } + + + template <typename I> + inline + mln_concrete(I) + integral(const Image<I>& input, unsigned scale) + { + // mln_precondition(input.nrows() % scale == 0); + // mln_precondition(input.ncols() % scale == 0); + if (scale == 3) + return integral_<3>(input); + else + std::cerr << "NYI!" << std::endl; + } + + } // end of namespace mln::subsampling::impl + + + + + // Facades. + + template <typename I> + inline + mln_concrete(I) + integral(const Image<I>& input_, unsigned scale) + { + trace::entering("subsampling::integral"); + + const I& input = exact(input_); + + mln_precondition(input.is_valid()); + mln_precondition(input.domain().pmin() == literal::origin); + mln_precondition(scale > 1); + + mln_concrete(I) output; + output = impl::integral(input, scale); + + trace::exiting("subsampling::integral"); + return output; + } + + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::subsampling + +} // end of namespace mln + + +#endif // ! MLN_SUBSAMPLING_INTEGRAL_HH Index: bench/z_sub_browsing/in.pgm.gz Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: bench/z_sub_browsing/in.pgm.gz ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Index: bench/z_sub_browsing/debase.cc --- bench/z_sub_browsing/debase.cc (revision 0) +++ bench/z_sub_browsing/debase.cc (revision 0) @@ -0,0 +1,23 @@ +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/util/timer.hh> + +#include "debase.hh" + + +int main() +{ + using namespace mln; + using value::int_u8; + + image2d<int_u8> in; + io::pgm::load(in, "in.pgm"); + + util::timer t; + t.start(); + image2d<int_u8> out = subsampling::debase(in, 3, dpoint2d(1,1)); + float t_ = t.stop(); + std::cout << t_ << std::endl; + + io::pgm::save(out, "out.pgm"); +} Index: bench/z_sub_browsing/integral.cc --- bench/z_sub_browsing/integral.cc (revision 0) +++ bench/z_sub_browsing/integral.cc (revision 0) @@ -0,0 +1,23 @@ +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/util/timer.hh> + +#include "integral.hh" + + +int main() +{ + using namespace mln; + using value::int_u8; + + image2d<int_u8> in; + io::pgm::load(in, "in.pgm"); + + util::timer t; + t.start(); + image2d<int_u8> out = subsampling::integral(in, 3); + float t_ = t.stop(); + std::cout << t_ << std::endl; + + io::pgm::save(out, "out.pgm"); +} Index: bench/z_sub_browsing/README --- bench/z_sub_browsing/README (revision 0) +++ bench/z_sub_browsing/README (revision 0) @@ -0,0 +1,5 @@ +debase en 03 DNDEBUG: 0.08 +integral: 0.16 +fast: ??? + +les 3 codes font une moyenne de tiles pour le sub-sampling Index: theo/mln/morpho/canvas/lena_blurred.pgm.gz Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: theo/mln/morpho/canvas/lena_blurred.pgm.gz ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Index: theo/mln/morpho/canvas/lena.pgm.gz Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: theo/mln/morpho/canvas/lena.pgm.gz ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Index: theo/mln/morpho/canvas/g.pbm.gz Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: theo/mln/morpho/canvas/g.pbm.gz ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Index: theo/mln/morpho/canvas/lena_min.pgm.gz Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: theo/mln/morpho/canvas/lena_min.pgm.gz ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Index: theo/mln/morpho/canvas/reconstruction_on_set.hh Index: theo/mln/morpho/canvas/f_and_g.pbm.gz Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: theo/mln/morpho/canvas/f_and_g.pbm.gz ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Index: theo/mln/morpho/canvas/regminid.pbm.gz Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: theo/mln/morpho/canvas/regminid.pbm.gz ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Index: theo/mln/morpho/canvas/one_domain.cc --- theo/mln/morpho/canvas/one_domain.cc (revision 0) +++ theo/mln/morpho/canvas/one_domain.cc (revision 0) @@ -0,0 +1,569 @@ + +# include <mln/core/concept/image.hh> +# include <mln/core/concept/neighborhood.hh> +# include <mln/morpho/canvas/internal/find_root.hh> + +# include <mln/data/fill.hh> +# include <mln/data/sort_psites.hh> +# include <mln/border/equalize.hh> +# include <mln/border/fill.hh> + + +// cc + +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/debug/println.hh> +#include <mln/io/pbm/load.hh> +#include <mln/io/pgm/load.hh> +#include <mln/io/pbm/save.hh> +#include <mln/io/pgm/save.hh> +#include <mln/data/compare.hh> +#include <mln/util/timer.hh> + +#include <mln/labeling/regional_maxima.hh> +#include <mln/labeling/regional_minima.hh> +#include <mln/pw/all.hh> + + + +namespace mln +{ + + namespace morpho + { + + namespace canvas + { + + + template <typename I, typename N, typename F> + inline + mln_result(F) + one_domain_union_find(const Image<I>& f_, const Neighborhood<N>& nbh_, + F& functor) + { + trace::entering("morpho::canvas::one_domain_union_find"); + + const I& f = exact(f_); + const N& nbh = exact(nbh_); + + typedef typename F::S S; + typedef mln_psite(I) P; + typedef mln_value(I) V; + + // Auxiliary data. + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, P) parent; + mln_result(F) output; + + // Initialization. + { + initialize(output, f); + functor.output = output; + + initialize(parent, f); + initialize(deja_vu, f); + data::fill(deja_vu, false); + + functor.set_default(); // <-- set_default + } + + // First pass. + { + mln_bkd_piter(S) p(functor.s()); + mln_niter(N) n(nbh, p); + for_all(p) + if (functor.domain(p)) // <-- domain + { + // Make-Set. + parent(p) = p; + functor.init(p); // <-- init + + for_all(n) if (f.domain().has(n)) + { + // See below (*) + if (functor.domain(n)) // <-- domain + { + if (deja_vu(n)) + { + if (functor.equiv(p, n)) // <-- equiv + { + // Do-Union. + P r = internal::find_root(parent, n); + if (r != p) + { +// if (do_merge(p, r)) + { + parent(r) = p; + functor.merge(p, r); // <-- merge + } +// else +// no_merge(p, r); + } + } + else + functor.do_no_union(p, n); // <-- do_no_union + } + } + else + { + mln_invariant(deja_vu(n) == false); + functor.border(p, n); // <-- border + } + } + deja_vu(p) = true; + } + } + + // Second pass. + { + mln_fwd_piter(S) p(functor.s()); + for_all(p) + if (functor.domain(p)) // <-- domain + if (parent(p) == p) + functor.root(p); // <-- root + else + output(p) = output(parent(p)); + } + + trace::exiting("morpho::canvas::one_domain_union_find"); + return output; + } + + } // end of namespace mln::morpho::canvas + + } // end of namespace mln::morpho + + + + // base functor + + struct functor_base + { + void set_default() + { + } + template <typename P> + void init(const P&) + { + } + template <typename P> + void root(const P&) + { + } + template <typename P> + bool domain(const P&) + { + return true; + } + template <typename P1, typename P2> + void do_no_union(const P1&, const P2&) + { + } + template <typename P1, typename P2> + void border(const P1&, const P2&) + { + } + template <typename P1, typename P2> + bool equiv(const P1&, const P2&) + { + return true; + } + template <typename P1, typename P2> + void merge(const P1&, const P2&) + { + } + // no default for: + // typedef S + // .s() + }; + + + // "binary reconstruction by dilation" functor + + template <typename F, typename G> + struct rd_functor : functor_base + { + typedef mln_concrete(F) result; + + const F& f; + const G& g; + + typedef mln_domain(F) S; + + rd_functor(const F& f, const G& g) + : f(f), g(g) + { + } + + result output; + typedef mln_psite(F) P; + + void set_default() + { + mln::data::fill(output, false); + } + + bool domain(const P& p) const + { + return g(p) == true; + } + + void init(const P& p) + { + output(p) = f(p); + } + + void merge(const P& p, const P& r) + { + output(p) = output(p) || output(r); + } + + S s() const + { + return f.domain(); + } + + }; + + + // "labeling" functor + + template <typename I, typename L> + struct lab_functor : functor_base + { + typedef mln_ch_value(I, L) result; + + const I& input; + mln_value(I) val; + + L nlabels; + bool status; + + lab_functor(const I& input, const mln_value(I)& val) + : input(input), val(val) + { + } + + result output; + typedef mln_psite(I) P; + + typedef mln_domain(I) S; + + S s() const + { + return input.domain(); + } + + void set_default() + { + data::fill(output, L(literal::zero)); + nlabels = 0; + status = true; + } + + bool domain(const P& p) const + { + return input(p) == val; + } + + bool equiv(const P& p, const P& n) + { + return input(n) == val; + } + + void root(const P& p) + { + if (status == false) + return; + + if (nlabels == mln_max(L)) + { + status = false; + trace::warning("labeling aborted!"); + } + else + output(p) = ++nlabels; + } + + }; + + + // "regional maxima labeling" functor + + template <typename I, typename L> + struct regmaxlab_functor : functor_base + { + typedef mln_ch_value(I, L) result; + + const I& input; + + L nlabels; + bool status; + mln_ch_value(I, bool) attr; + + typedef mln_psite(I) P; + typedef p_array<P> S; + S s_; + + result output; + + regmaxlab_functor(const I& input) + : input(input) + { + s_ = data::sort_psites_increasing(input); + } + + const S& s() const + { + return s_; + } + + void set_default() + { + initialize(attr, input); + data::fill(attr, true); + data::fill(output, L(literal::zero)); + nlabels = 0; + status = true; + } + + bool equiv(const P& p, const P& n) + { + return input(n) == input(p); + } + + void merge(const P& p, const P& r) + { + attr(p) = attr(p) && attr(r); + } + + void do_no_union(const P& p, const P& n) + { + mln_invariant(input(n) > input(p)); + attr(p) = false; + (void)n; + } + + void root(const P& p) + { + if (attr(p) == false || status == false) + return; + + if (nlabels == mln_max(L)) + { + status = false; + trace::warning("labeling aborted!"); + } + else + output(p) = ++nlabels; + } + + }; + + + // "regional minima identification" functor + + template <typename I> + struct regminid_functor : functor_base + { + typedef mln_ch_value(I, bool) result; + + const I& input; + + typedef mln_psite(I) P; + typedef p_array<P> S; + S s_; + + result output; + + regminid_functor(const I& input) + : input(input) + { + s_ = data::sort_psites_decreasing(input); + } + + const S& s() const + { + return s_; + } + + void init(const P& p) + { + output(p) = true; + } + + bool equiv(const P& p, const P& n) + { + return input(n) == input(p); + } + + void merge(const P& p, const P& r) + { + output(p) = output(p) && output(r); + } + + void do_no_union(const P& p, const P& n) + { + if (input(n) < input(p)) + output(p) = false; + } + + }; + + + + // "reconstruction by dilation on function" functor + + template <typename F, typename G> + struct rdf_functor : functor_base + { + typedef mln_concrete(F) result; + + const F& f; + const G& g; + + typedef mln_value(F) V; + typedef mln_psite(F) P; + typedef p_array<P> S; + S s_; + + result output; + + rdf_functor(const F& f, const G& g) + : f(f), g(g) + { + s_ = data::sort_psites_decreasing(g); + } + + const S& s() const + { + return f.domain(); + } + + void set_default() + { + mln::data::fill(output, f); + } + + bool equiv(const P& p, const P& n) + { + return g(r) == g(p) || g(p) >= output(r); + } + + void merge(const P& p, const P& r) + { + if (output(r) > output(p)) + output(p) = output(r); + } + + void do_no_union(const P& p, const P& n) + { + output(p) = mln_max(V); + } + + void root(const P& p) + { + if (output(p) == mln_max(V)) + output(p) = g(p); + } + + }; + + + +} // end of namespace mln + + + +int main() +{ + using namespace mln; + + typedef value::int_u8 L; + neighb2d nbh = c4(); + + util::timer t; + + + // On sets. + +// { +// typedef image2d<bool> I; + +// I f, g, ref; + +// io::pbm::load(f, "f_and_g.pbm"); +// io::pbm::load(g, "g.pbm"); + +// { +// t.start(); +// rd_functor<I,I> fun(f, g); +// I rd = morpho::canvas::one_domain_union_find(f, nbh, fun); +// std::cout << "rd: " << t.stop() << std::endl; +// io::pbm::save(rd, "rd.pbm"); +// } + +// { +// t.start(); +// lab_functor<I,L> fun(f, true); +// image2d<L> lab = morpho::canvas::one_domain_union_find(f, nbh, fun); +// std::cout << "lab: " << t.stop() << std::endl; +// io::pgm::save(lab, "lab.pgm"); +// } +// } + + + // On functions. + + { + typedef image2d<value::int_u8> I; + I f; + io::pgm::load(f, "./lena_blurred.pgm"); + +// { +// // regional maxima labeling +// t.start(); +// regmaxlab_functor<I,L> fun(f); +// image2d<L> lab = morpho::canvas::one_domain_union_find(f, nbh, fun); +// std::cout << "regmax lab: " << t.stop() << std::endl; +// io::pgm::save(lab, "regmaxlab.pgm"); + +// L nlabels_ref; +// if (lab != labeling::regional_maxima(f, nbh, nlabels_ref) +// || fun.nlabels != nlabels_ref) +// std::cerr << "oops" << std::endl; +// } + +// { +// // regional minima identification +// t.start(); +// regminid_functor<I> fun(f); +// image2d<bool> id = morpho::canvas::one_domain_union_find(f, nbh, fun); +// std::cout << "regmin id: " << t.stop() << std::endl; +// io::pbm::save(id, "regminid.pbm"); + +// L forget_it; +// if (id != ((pw::value(labeling::regional_minima(f, nbh, forget_it)) != pw::cst(0)) | id.domain())) +// std::cerr << "oops" << std::endl; +// } + + { + // reconstruction by dilation + I f, g; + io::pgm::load(f, "lena.pgm"); + io::pgm::load(g, "lena_min.pgm"); + if (! (f <= g)) + std::cerr << "OK" << std::endl; + + rdf_functor<I,I> fun(f, g); + I rdf = morpho::canvas::one_domain_union_find(f, nbh, fun); + std::cout << "rdf: " << t.stop() << std::endl; + io::pgm::save(rdf, "rdf.pbm"); + + } + + } + +} Index: theo/mln/subsampling/sizes.cc --- theo/mln/subsampling/sizes.cc (revision 0) +++ theo/mln/subsampling/sizes.cc (revision 0) @@ -0,0 +1,45 @@ +#include <iostream> +#include <cstdlib> + + +void usage(const char* argv[]) +{ + std::cerr << argv[0] << " n_scales s q nr nc" << std::endl; + std::abort(); +} + + +unsigned sub(unsigned nbr, unsigned down_scaling) +{ + return (nbr + down_scaling - 1) / down_scaling; +} + + +int main(int argc, const char* argv[]) +{ + unsigned n_scales = std::atoi(argv[1]); + unsigned s = std::atoi(argv[2]); + unsigned q = std::atoi(argv[3]); + + unsigned nr_1 = std::atoi(argv[4]); + unsigned nc_1 = std::atoi(argv[5]); + + unsigned* nc; + nc = new unsigned[n_scales + 1]; + + nc[1] = nc_1; + nc[2] = sub(nc[1], s); + for (unsigned i = 3; i <= n_scales; ++i) + nc[i] = sub(nc[i - 1], q); + + for (unsigned i = 1; i <= n_scales; ++i) + std::cout << nc[i] << ' '; + std::cout << std::endl; + + for (unsigned i = n_scales; i > 2; --i) + nc[i - 1] = q * nc[i]; + + for (unsigned i = 1; i <= n_scales; ++i) + std::cout << nc[i] << ' '; + std::cout << std::endl; +} Index: theo/mln/subsampling/debase.hh --- theo/mln/subsampling/debase.hh (revision 0) +++ theo/mln/subsampling/debase.hh (revision 0) @@ -0,0 +1,351 @@ +// Copyright (C) 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 MLN_SUBSAMPLING_DEBASE_HH +# define MLN_SUBSAMPLING_DEBASE_HH + +/// \file +/// +/// Debase subsampling. + +#include <mln/core/concept/image.hh> + +#define MLN_FLOAT double +#include <sandbox/theo/exec/gaussian_directional_2d.hh> + + + +namespace mln +{ + + namespace subsampling + { + + /// FIXME: Doc. + + template <typename I> + mln_concrete(I) + debase(const Image<I>& input, + unsigned gap, + const mln_deduce(I, site, delta)& shift); // FIXME: Add round_up_size. + + + template <typename I> + mln_concrete(I) + debase(const Image<I>& input, unsigned gap); + + + +# ifndef MLN_INCLUDE_ONLY + + + // Implementation. + + namespace impl + { + + + template <typename I> + inline + mln_concrete(I) + debase(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase"); + + const I& input = exact(input_); + + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + + box<P> b(pmin, pmax); + mln_concrete(I) output(b); + mln_piter(I) p(output.domain()); + for_all(p) + output(p) = input(p * gap + shift); + + trace::exiting("subsampling::impl::debase"); + return output; + } + + + template <typename I> + inline + mln_concrete(I) + debase_fastest(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase_fastest"); + + const I& input = exact(input_); + + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + + box<P> b(pmin, pmax); + + typedef mln_concrete(I) O; + O output(b); + + mln_box_runstart_piter(O) s(output.domain()); + const unsigned n = s.run_length(); + for_all(s) + { + const mln_value(I)* pi = & input(s * gap + shift); + mln_value(O)* po = & output(s); + for (unsigned i = 0; i < n; ++i) + { + *po++ = *pi; + pi += gap; + } + } + + trace::exiting("subsampling::impl::debase_fastest"); + return output; + } + + + template <typename I> + inline + mln_concrete(I) + debase_gaussian_antialiased_fastest(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase_gaussian_antialiased_fastest"); + + const I& input = exact(input_); + + image2d<MLN_FLOAT> temp(input.domain()); + data::fill(temp, input); + temp = linear::gaussian_directional_2d(temp, 1, 1.5, 0); + temp = linear::gaussian_directional_2d(temp, 0, 1.5, 0); + + typedef mln_value(I) V; + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + + box<P> b(pmin, pmax); + + typedef mln_concrete(I) O; + O output(b); + + mln_box_runstart_piter(O) s(output.domain()); + const unsigned n = s.run_length(); + for_all(s) + { + const MLN_FLOAT* pi = & temp(s * gap + shift); + mln_value(O)* po = & output(s); + for (unsigned i = 0; i < n; ++i) + { + *po++ = V(*pi + 0.49); + pi += gap; + } + } + + trace::entering("subsampling::impl::debase_gaussian_antialiased_fastest"); + return output; + } + + + + + template <typename I> + inline + mln_concrete(I) + debase_2d_antialias_pointer_fastest(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase_2d_antialias_pointer_fastest"); + + const I& input = exact(input_); + + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + + box<P> b(pmin, pmax); + + typedef mln_concrete(I) O; + O output(b); + + const unsigned nrows = input.nrows(); + const unsigned ncols = input.ncols(); + + typedef mln_value(I) V; + typedef mln_sum(V) S; + for (unsigned row = 0; row < nrows; row += gap) + { + const V* ptr1 = & input.at_(row, 0); + const V* ptr2 = & input.at_(row + 1, 0); + const V* ptr3 = & input.at_(row + 2, 0); + for (unsigned col = 0; col < ncols; col += gap) + { + S s; + + // mean + s = *ptr1++ + *ptr1++ + *ptr1++; + s += *ptr2++ + *ptr2++ + *ptr2++; + s += *ptr3++ + *ptr3++ + *ptr3++; + output.at_(row / gap, col / gap) = (s + 4) / 9; +// // cut +// s = 1 * *ptr1++ + 2 * *ptr1++ + 1 * *ptr1++; +// s += 2 * *ptr2++ + 3 * *ptr2++ + 2 * *ptr2++; +// s += 1 * *ptr3++ + 2 * *ptr3++ + 1 * *ptr3++; +// output.at_(row / gap, col / gap) = (s + 7) / 15; + } + } + + trace::exiting("subsampling::impl::debase_2d_antialias_pointer_fastest"); + return output; + } + + + + template <typename I> + inline + mln_concrete(I) + debase_2d_antialias_offset_fastest(const Image<I>& input_, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::impl::debase_2d_antialias_offset_fastest"); + + const I& input = exact(input_); + + typedef mln_value(I) V; + typedef mln_sum(V) S; + + typedef mln_site(I) P; + P pmin = input.domain().pmin() / gap, + pmax = input.domain().pmax() / gap; + box<P> b(pmin, pmax); + typedef mln_concrete(I) O; + O output(b); + + const V* ptr1 = & input.at_(0, 0); + const V* ptr2 = & input.at_(1, 0); + const V* ptr3 = & input.at_(2, 0); + + mln_box_runstart_piter(O) s(output.domain()); + const unsigned n = s.run_length(); + + unsigned offset = input.delta_index(point2d(3,0) - point2d(0,3*n)); + + for_all(s) + { + mln_value(O)* po = & output(s); + for (unsigned i = 0; i < n; ++i) + { + S s; + s = *ptr1++ + *ptr1++ + *ptr1++; + s += *ptr2++ + *ptr2++ + *ptr2++; + s += *ptr3++ + *ptr3++ + *ptr3++; + *po++ = (s + 4) / 9; + } + ptr1 += offset; + ptr2 += offset; + ptr3 += offset; + } + + trace::exiting("subsampling::impl::debase_2d_antialias_offset_fastest"); + return output; + } + + + } // end of namespace mln::subsampling::impl + + + + // Dispatch. + + namespace internal + { + + template <typename I> + inline + mln_concrete(I) + debase_dispatch(const Image<I>& input, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + return impl::debase_2d_antialias_pointer_fastest /* debase_2d_antialias_offset_fastest */(input, gap, shift); + } + + + } // end of namespace mln::subsampling::internal + + + + // Facades. + + template <typename I> + inline + mln_concrete(I) + debase(const Image<I>& input, + unsigned gap, + const mln_deduce(I, site, delta)& shift) + { + trace::entering("subsampling::debase"); + + typedef mln_site(I) P; + + mln_precondition(exact(input).is_valid()); + mln_precondition(exact(input).domain().pmin() == literal::origin); + mln_precondition(gap > 1); + for (unsigned i = 0; i < P::dim; ++i) + mln_precondition(shift[i] < gap); + + mln_concrete(I) output; + output = internal::debase_dispatch(input, gap, shift); + + trace::exiting("subsampling::debase"); + return output; + } + + + template <typename I> + inline + mln_concrete(I) + debase(const Image<I>& input, unsigned gap) + { + return debase(input, gap, literal::zero); + } + + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::subsampling + +} // end of namespace mln + + +#endif // ! MLN_SUBSAMPLING_DEBASE_HH Index: theo/mln/subsampling/integral.hh --- theo/mln/subsampling/integral.hh (revision 0) +++ theo/mln/subsampling/integral.hh (revision 0) @@ -0,0 +1,295 @@ +// Copyright (C) 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 MLN_SUBSAMPLING_INTEGRAL_HH +# define MLN_SUBSAMPLING_INTEGRAL_HH + +/// \file +/// +/// Both subsampling and integral image computation. + +#include <mln/core/concept/image.hh> +#include <mln/debug/println.hh> + + + +namespace mln +{ + + namespace subsampling + { + + /// FIXME: Doc. + + template <typename I> + mln_concrete(I) + integral(const Image<I>& input, unsigned scale); + + + +# ifndef MLN_INCLUDE_ONLY + + + // Implementation. + + namespace impl + { + + + template <unsigned scale, typename I> + inline + mln_concrete(I) + integral__why_is_it_slow(const Image<I>& input_) + { + trace::entering("subsampling::impl::integral__why_is_it_slow"); + + const I& input = exact(input_); + const unsigned scale2 = scale * scale, round_it = scale2 / 2; + + typedef mln_value(I) V; + typedef mln_sum(V) S; + typedef mln_site(I) P; + + box<P> b = make::box2d(input.nrows() / scale, input.ncols() / scale); + const unsigned nrows = b.nrows(), ncols = b.ncols(); + + typedef mln_concrete(I) O; + O output(b, 0); // no external border in 'output' + V* po = output.buffer(); + + unsigned + tile_cr = input.delta_index(dpoint2d(scale, - ncols * scale)), + cr = input.delta_index(dpoint2d(1, - (int)scale)); + + const V* pi = & input.at_(0, 0); // first point of tiles + + // (row, col) are tile coordinates. + for (unsigned row = 0; row < nrows; ++row) + { + S hsum = 0, hsum2 = 0; // horizontal integral + for (unsigned col = 0; col < ncols; ++col) + { + S sum = 0, sum2 = 0; + const V* ptr = pi; + + for (unsigned r = 0; r < scale; ++r) + { + for (unsigned c = 0; c < scale; ++c) + { + sum += *ptr; + ++ptr; + } + ptr += cr; + } + // ... + pi += scale; + *po++ = (round_it + sum) / scale2; + } + pi += tile_cr; + } + + trace::exiting("subsampling::impl::integral__why_is_it_slow"); + return output; + } + + + + template <unsigned scale, typename I> + inline + mln_concrete(I) + integral_(const Image<I>& input_) + { + trace::entering("subsampling::impl::integral_"); + + const I& input = exact(input_); + const unsigned scale2 = scale * scale, round_it = scale2 / 2; + + typedef mln_value(I) V; + typedef mln_sum(V) S; + typedef mln_site(I) P; + + box<P> b = make::box2d((input.nrows() + scale - 1) / scale, + (input.ncols() + scale - 1) / scale); + typedef mln_concrete(I) O; + O output(b, 0); + V* po = output.buffer(); + + const unsigned nrows = input.nrows(); + const unsigned ncols = input.ncols(); + + typedef const V* Ptr; + Ptr ptr[scale]; + for (unsigned row = 0; row < nrows; row += scale) + { + for (unsigned r = 0; r < scale; ++r) + ptr[r] = & input.at_(row + r, 0); + + for (unsigned col = 0; col < ncols; col += scale) + { + S sum = 0; + for (unsigned r = 0; r < scale; ++r) + for (unsigned c = 0; c < scale; ++c) + // g++ does *not* un-roll this inner loop!!! + // icc does... + sum += *ptr[r]++; + + *po++ = (round_it + sum) / scale2; + } + } + + trace::exiting("subsampling::impl::integral_"); + return output; + } + + + + + template <typename I> + inline + mln_concrete(I) + integral_3(const Image<I>& input_) + { + trace::entering("subsampling::impl::integral_3"); + + const unsigned scale = 3; + + const I& input = exact(input_); + const unsigned area = scale * scale; + + typedef mln_value(I) V; + typedef mln_sum(V) S; + typedef mln_site(I) P; + + box<P> b = make::box2d((input.nrows() + scale - 1) / scale, + (input.ncols() + scale - 1) / scale); + const unsigned up = b.ncols(); + + mln_concrete(I) sub(b, 0); + V* p_sub = sub.buffer(); + + mln_ch_value(I, S) integral_sum(b, 0); + S* p_isum = integral_sum.buffer(); + + mln_ch_value(I, S) integral_sum_2(b, 0); + S* p_isum_2 = integral_sum_2.buffer(); + + const unsigned nrows = input.nrows(); + const unsigned ncols = input.ncols(); + + for (unsigned row = 0; row < nrows; row += scale) + { + S h_sum = 0, h_sum_2 = 0; + const V* ptr1 = & input.at_(row, 0); + const V* ptr2 = & input.at_(row + 1, 0); + const V* ptr3 = & input.at_(row + 2, 0); + for (unsigned col = 0; col < ncols; col += scale) + { + S sum; + sum = *ptr1++ + *ptr1++ + *ptr1++; + sum += *ptr2++ + *ptr2++ + *ptr2++; + sum += *ptr3++ + *ptr3++ + *ptr3++; + + S val = sum / area; + *p_sub++ = val; + + h_sum += val; + h_sum_2 += val * val; + + if (row != 0) + { + *p_isum = h_sum + *(p_isum - up); + *p_isum_2 = h_sum_2 + *(p_isum_2 - up); + } + else + { + // exception + *p_isum = h_sum; + *p_isum_2 = h_sum_2; + } + + p_isum += 1; + p_isum_2 += 1; + } + } + + debug::println(sub); + debug::println(integral_sum); + debug::println(integral_sum_2); + + trace::exiting("subsampling::impl::integral_3"); + return sub; + } + + + template <typename I> + inline + mln_concrete(I) + integral(const Image<I>& input, unsigned scale) + { + // mln_precondition(input.nrows() % scale == 0); + // mln_precondition(input.ncols() % scale == 0); + if (scale == 3) + return integral_3(input); + else + std::cerr << "NYI!" << std::endl; + } + + + } // end of namespace mln::subsampling::impl + + + + + // Facades. + + template <typename I> + inline + mln_concrete(I) + integral(const Image<I>& input_, unsigned scale) + { + trace::entering("subsampling::integral"); + + const I& input = exact(input_); + + mln_precondition(input.is_valid()); + mln_precondition(input.domain().pmin() == literal::origin); + mln_precondition(scale > 1); + + mln_concrete(I) output; + output = impl::integral(input, scale); + + trace::exiting("subsampling::integral"); + return output; + } + + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::subsampling + +} // end of namespace mln + + +#endif // ! MLN_SUBSAMPLING_INTEGRAL_HH Index: theo/mln/subsampling/in.pgm.gz Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: theo/mln/subsampling/in.pgm.gz ___________________________________________________________________ Added: svn:mime-type + application/octet-stream Index: theo/mln/subsampling/debase.cc --- theo/mln/subsampling/debase.cc (revision 0) +++ theo/mln/subsampling/debase.cc (revision 0) @@ -0,0 +1,23 @@ +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/util/timer.hh> + +#include "debase.hh" + + +int main() +{ + using namespace mln; + using value::int_u8; + + image2d<int_u8> in; + io::pgm::load(in, "in.pgm"); + + util::timer t; + t.start(); + image2d<int_u8> out = subsampling::debase(in, 3, dpoint2d(1,1)); + float t_ = t.stop(); + std::cout << t_ << std::endl; + + io::pgm::save(out, "out.pgm"); +} Index: theo/mln/subsampling/integral.cc --- theo/mln/subsampling/integral.cc (revision 0) +++ theo/mln/subsampling/integral.cc (revision 0) @@ -0,0 +1,38 @@ +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/util/timer.hh> + +#include "integral.hh" + +#include <mln/debug/iota.hh> +#include <mln/debug/println.hh> + + +int main() +{ + using namespace mln; + using value::int_u8; + + + image2d<int_u8> in(6, 9); + debug::iota(in); + debug::println(in); + + image2d<int_u8> out = subsampling::integral(in, 3); +// debug::println(out); + + +// { +// image2d<int_u8> in; +// io::pgm::load(in, "in.pgm"); + +// util::timer t; +// t.start(); +// image2d<int_u8> out = subsampling::integral(in, 3); +// float t_ = t.stop(); +// std::cout << t_ << std::endl; + +// io::pgm::save(out, "out.pgm"); +// } + +} Index: theo/mln/browsing/window_sliding.cc --- theo/mln/browsing/window_sliding.cc (revision 0) +++ theo/mln/browsing/window_sliding.cc (revision 0) @@ -0,0 +1,52 @@ +#include <iostream> +#include <cassert> + + + + +void window_sliding(unsigned nr, unsigned nc, + unsigned h, unsigned w) +{ + const unsigned h_2 = h / 2; + const int d_row_up = - 1 - h_2, d_row_down = h_2; + + // top part + for (int row = 0; row <= h_2; ++row) + { + assert(row + d_up < 0); + assert(row + d_down >= 0 && row + d_down < nr); + // top left + // ... + + // top middle + // ... + + // top right + // ... + } + // main part + for (int row = h_2 + 1; row <= nr - 1 - h_2; ++row) + { + assert(row + d_up >= 0 && row + d_up < nr); + assert(row + d_down >= 0 && row + d_down < nr); + } + // bottom part + for (int row = nr - h_2; row < nr; ++row) + { + assert(row + d_up >= 0 && row + d_up < nr); + assert(row + d_down >= nr); + } +} + + + +int main() +{ + +// unsigned nr = 50, nc = 100; + + unsigned nr = 10, nc = 10; + unsigned h = 3, w = 5; + window_sliding(nr, nc, h, w); + +}