4721: Add some bench + canvas + subsampling + browsing code. - Olena-patches

17 Nov 2009

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 image2dvalue::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);
+
+}

    