[Olena-patches] 4331: Add some experimental code.

27 Jul 2009

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox


Index: ChangeLog
from  Thierry Geraud  <thierry.geraud@lrde.epita.fr>

	Add some experimental code.

	* icdar/2009/pscomp: New directory.
	* icdar/2009/pscomp/iz.cc: New.
	Compute an IZ image along with distances between objects. 
	* icdar/2009/pscomp/rect_filter.cc: New.
	Algebraic filtering with a combinaison of bounding box width
	and height.

	* theo/experimental/depeche: New directory.
	* theo/experimental/depeche/homogen.cc: New.
	* theo/experimental/depeche/col.cc: New.
	* theo/experimental/depeche/row_mm.cc: New.
	* theo/experimental/depeche/row.cc: New.

 icdar/2009/pscomp/iz.cc              |  174 +++++++++++++++++++++
 icdar/2009/pscomp/rect_filter.cc     |  165 ++++++++++++++++++++
 theo/experimental/depeche/col.cc     |  286 +++++++++++++++++++++++++++++++++++
 theo/experimental/depeche/homogen.cc |   43 +++++
 theo/experimental/depeche/row.cc     |  117 ++++++++++++++
 theo/experimental/depeche/row_mm.cc  |   96 +++++++++++
 6 files changed, 881 insertions(+)

Index: icdar/2009/pscomp/iz.cc

--- icdar/2009/pscomp/iz.cc	(revision 0)
+++ icdar/2009/pscomp/iz.cc	(revision 0)
@@ -0,0 +1,174 @@
+
+// From mln/canvas/distance_geodesic.hh
+
+# include <mln/core/concept/image.hh>
+# include <mln/core/concept/neighborhood.hh>
+# include <mln/core/routine/duplicate.hh>
+# include <mln/core/site_set/p_queue_fast.hh>
+# include <queue>
+# include <mln/data/fill.hh>
+# include <mln/extension/adjust_fill.hh>
+
+
+// Local.
+
+#include <map>
+#include <mln/core/image/image2d.hh>
+#include <mln/core/alias/neighb2d.hh>
+#include <mln/value/int_u8.hh>
+#include <mln/io/pbm/load.hh>
+#include <mln/io/pgm/save.hh>
+#include <mln/debug/println.hh>
+#include <mln/data/saturate.hh>
+
+#include <mln/labeling/blobs.hh>
+#include <mln/labeling/value.hh>
+
+
+
+
+namespace mln
+{
+
+
+  // From  mln/canvas/distance_geodesic.hh
+  //  and  mln/transform/internal/influence_zone_functor.hh
+  //  and  mln/transform/influence_zone_geodesic.hh
+
+
+  template <typename I, typename N>
+  std::map< std::pair<mln_value(I),mln_value(I)>, unsigned >
+  distances(// in:
+	    const I& input, const N& nbh,
+	    // out:
+	    mln_ch_value(I, unsigned)& dmap,
+	    mln_concrete(I)& iz)
+  {
+    trace::entering("canvas::impl::generic::distance_geodesic");
+
+    typedef mln_value(I) L;
+    std::map< std::pair<L,L>, unsigned > dist; // NEW
+
+    const unsigned max = mln_max(unsigned);
+    typedef mln_site(I) P;
+    p_queue_fast<P> q;
+
+    // Initialization.
+    {
+      iz = duplicate(input); // <-- init
+
+      data::fill(dmap, max);
+
+      mln_piter(I) p(input.domain());
+      mln_niter(N) n(nbh, p);
+      for_all(p)
+	if (input(p) != 0) // <-- inqueue_p_wrt_input_p
+	  {
+	    ; // <-- init_p
+	    dmap(p) = 0;
+	    for_all(n)
+	      if (input.domain().has(n) &&
+		  input(n) == 0) // <-- inqueue_p_wrt_input_n
+		{
+		  q.push(p);
+		  break;
+		}
+	  }
+    }
+
+    // Propagation.
+    {
+      P p;
+      mln_niter(N) n(nbh, p);
+      while (! q.is_empty())
+	{
+	  p = q.pop_front();
+	  if (dmap(p) == max)
+	    {
+	      // Saturation so stop.
+	      q.clear();
+	      break;
+	    }
+	  for_all(n)
+	  {
+	    if (! input.domain().has(n))
+	      continue;
+	    if (dmap(n) == max)
+	      {
+		dmap(n) = dmap(p) + 1;
+		iz(n) = iz(p); // <- process
+		q.push(n);
+	      }
+	    else
+	      {
+		if (iz(n) != iz(p))
+		  {
+		    L l1 = iz(n), l2 = iz(p);
+		    if (l1 > l2)
+		      std::swap(l1, l2);
+		    unsigned& d_ = dist[std::make_pair(l1,l2)];
+		    unsigned d12 = dmap(p) + dmap(n) + 1;
+		    if (d_ == 0 || d12 < d_)
+		      d_ = d12;
+		  }
+	      }
+	  }
+	}
+    }
+
+    return dist;
+  }
+
+
+} // mln
+
+
+
+int main(int argc, char* argv[])
+{
+  using namespace mln;
+
+  typedef unsigned L;
+
+  image2d<bool> input;
+  io::pbm::load(input, argv[1]);
+
+  L n_labels;
+  image2d<L> lab = labeling::value(input, true, c8(), n_labels);
+  std::cout << "n labels = " << n_labels << std::endl;
+
+
+  image2d<unsigned> dmap(input.domain());
+  image2d<L> iz(input.domain());
+
+  typedef std::map< std::pair<L,L>, unsigned > map_t; 
+  map_t d = distances(lab, c8(), dmap, iz);
+
+  io::pgm::save(data::saturate(value::int_u8(), dmap), "dmap.pgm");
+
+  // debug::println("lab = ", lab);
+  std::cout << "map size = " << d.size() << std::endl;
+
+
+//   {
+//     // Compute statistics.
+//     std::map<unsigned, unsigned> h;
+//     for (map_t::const_iterator i = d.begin();
+// 	 i != d.end(); ++i)
+//       ++h[i->second];
+
+//     for (std::map<unsigned, unsigned>::const_iterator i = h.begin();
+// 	 i != h.end(); ++i)
+//       std::cout << i->first << " " << i->second << std::endl;
+//   }
+
+
+//   {
+//     // Print: d(label1, label2) = shortest distance
+//     for (map_t::const_iterator i = d.begin();
+// 	 i != d.end(); ++i)
+//       std::cout << "d(" << i->first.first << ", " << i->first.second << ") = "
+// 		<< i->second << std::endl;
+//   }
+
+}
Index: icdar/2009/pscomp/rect_filter.cc
--- icdar/2009/pscomp/rect_filter.cc	(revision 0)
+++ icdar/2009/pscomp/rect_filter.cc	(revision 0)
@@ -0,0 +1,165 @@
+#include <mln/core/image/image2d.hh>
+#include <mln/core/alias/neighb2d.hh>
+#include <mln/value/int_u8.hh>
+#include <mln/io/pbm/load.hh>
+#include <mln/io/pgm/save.hh>
+#include <mln/debug/println.hh>
+
+#include <mln/transform/distance_geodesic.hh>
+#include <mln/accu/shape/bbox.hh>
+#include <mln/morpho/closing/algebraic.hh>
+#include <mln/morpho/watershed/flooding.hh>
+#include <mln/morpho/watershed/superpose.hh>
+
+#include <mln/core/routine/duplicate.hh>
+#include <mln/data/convert.hh>
+#include <mln/pw/all.hh>
+#include <mln/io/pbm/save.hh>
+#include <mln/io/ppm/save.hh>
+
+
+
+
+namespace mln
+{
+
+  struct pair_and
+  {
+    pair_and()
+    {}
+    pair_and(unsigned first, unsigned second)
+      : first(first),
+	second(second)
+    {}
+//     bool operator<(const pair_and& other) const
+//     {
+//       return first < other.first && second < other.second;
+//     }
+    unsigned first, second;
+  };
+
+
+  bool operator<(const pair_and& lhs, const pair_and& rhs)
+  {
+    return lhs.first < rhs.first || lhs.second < rhs.second;
+  }
+
+
+  struct rectangle_and;
+
+
+  namespace trait
+  {
+
+    template <>
+    struct accumulator_< rectangle_and >
+    {
+      typedef accumulator::has_untake::no	has_untake;
+      typedef accumulator::has_set_value::no	has_set_value;
+      typedef accumulator::has_stop::no		has_stop;
+      typedef accumulator::when_pix::use_p	when_pix;
+    };
+
+  } // end of namespace mln::trait
+
+
+  struct rectangle_and : public accu::internal::base< pair_and, rectangle_and >
+  {
+    typedef point2d argument;
+
+    rectangle_and()
+    {
+    }
+
+    void init()
+    {
+      bb_.init();
+    }
+    void take_as_init_(const point2d& p)
+    {
+      bb_.take_as_init_(p);
+    }
+    void take(const point2d& p)
+    {
+      bb_.take(p);
+    }
+    void take(const rectangle_and& other)
+    {
+      bb_.take(other.bb_);
+    }
+
+    pair_and to_result() const
+    {
+      pair_and tmp(bb_.to_result().len(0),
+		   bb_.to_result().len(1));
+      return tmp;
+    }
+
+    bool is_valid() const
+    {
+      return true;
+    }
+
+  protected:
+    accu::shape::bbox<point2d> bb_;
+  };
+
+
+} // mln
+
+
+void usage(char* argv[])
+{
+  std::cerr << "usage: " << argv[0] << " input.pgm w h output.ppm" << std::endl;
+  std::abort();
+}
+
+
+int main(int argc, char* argv[])
+{
+  using namespace mln;
+  using value::int_u8;
+  using value::rgb8;
+
+  if (argc != 5)
+    usage(argv);
+
+  image2d<bool> input;
+  io::pbm::load(input, argv[1]);
+
+  typedef value::int_u8 D;
+
+  image2d<D> dmap;
+  dmap = transform::distance_geodesic(input,
+				      c8(),
+				      mln_max(D));
+  io::pgm::save(dmap, "temp_dmap.pgm");
+
+  int w = std::atoi(argv[2]), h = std::atoi(argv[3]);
+
+  pair_and lambda(w, h);
+  image2d<D> clo = morpho::closing::algebraic(dmap,
+					      c8(),
+					      rectangle_and(),
+					      lambda);
+
+  io::pgm::save(clo, "temp_clo.pgm");
+
+  typedef unsigned L;
+  L n_basins;
+  image2d<L> ws = morpho::watershed::flooding(clo, c8(), n_basins);
+  std::cout << "n basins = " << n_basins<< std::endl;
+
+  {
+    image2d<rgb8> out = data::convert(rgb8(), dmap);
+    out = morpho::watershed::superpose(out, ws, rgb8(255,0,0));
+    io::ppm::save(out, argv[4]);
+  }
+
+//   {
+//     image2d<bool> out(input.domain());
+//     out = duplicate((pw::value(ws) == pw::cst(0)) | input.domain());
+//     io::pbm::save(out, argv[2]);
+//   }
+
+}
Index: theo/experimental/depeche/homogen.cc
--- theo/experimental/depeche/homogen.cc	(revision 0)
+++ theo/experimental/depeche/homogen.cc	(revision 0)
@@ -0,0 +1,43 @@
+#include <mln/core/image/image2d.hh>
+#include <mln/core/alias/neighb2d.hh>
+#include <mln/morpho/attribute/height.hh>
+#include <mln/morpho/closing/leveling.hh>
+
+#include <mln/arith/minus.hh>
+
+#include <mln/value/int_u8.hh>
+#include <mln/io/pgm/load.hh>
+#include <mln/io/pgm/save.hh>
+
+
+
+void usage(char* argv[])
+{
+  std::cerr << "usage: " << argv[0] << " input.pgm h output.pgm" << std::endl;
+  abort();
+}
+
+
+
+int main(int argc, char* argv[])
+{
+  using namespace mln;
+  using value::int_u8;
+
+  if (argc != 4)
+    usage(argv);
+
+  image2d<int_u8> input;
+  io::pgm::load(input, argv[1]);
+
+  int i = std::atoi(argv[2]);
+  int_u8 h = i;
+
+  morpho::attribute::height< image2d<int_u8> > acc;
+  image2d<int_u8> clo = morpho::closing::leveling(input, c4(),
+						  acc, h);
+  io::pgm::save(clo, "tmp.pgm");
+  arith::minus_inplace(clo, input);
+
+  io::pgm::save(clo, argv[3]);
+}
Index: theo/experimental/depeche/col.cc
--- theo/experimental/depeche/col.cc	(revision 0)
+++ theo/experimental/depeche/col.cc	(revision 0)
@@ -0,0 +1,286 @@
+#include <mln/core/image/image2d.hh>
+#include <mln/core/alias/neighb2d.hh>
+
+#include <mln/core/image/image1d.hh>
+#include <mln/core/alias/neighb1d.hh>
+
+#include <mln/fun/v2v/projection.hh>
+#include <mln/core/image/dmorph/unproject_image.hh>
+#include <mln/core/image/dmorph/image_if.hh>
+#include <mln/pw/value.hh>
+
+#include <mln/accu/stat/mean.hh>
+#include <mln/accu/image/init.hh>
+#include <mln/accu/image/take.hh>
+#include <mln/accu/image/to_result.hh>
+
+#include <mln/value/int_u8.hh>
+#include <mln/io/pgm/load.hh>
+#include <mln/io/plot/save.hh>
+
+#include <mln/linear/gaussian_1d.hh>
+#include <mln/labeling/regional_maxima.hh>
+#include <mln/data/convert.hh>
+#include <mln/io/ppm/save.hh>
+
+#include <mln/data/fill.hh>
+#include <mln/debug/println.hh>
+#include <mln/accu/stat/median_h.hh>
+
+#include <mln/morpho/closing/volume.hh>
+
+
+
+namespace mln
+{
+
+  
+  float get_approx_gap(const util::array<unsigned>& line, unsigned delta)
+  {
+    accu::stat::median_h< value::int_u<12> > med;
+    const unsigned n = line.nelements() - delta;
+    for (unsigned i = 0; i < n; ++i)
+      med.take(line[i + delta] - line[i]);
+
+     std::cout << "delta = " << delta << std::endl;
+     std::cout << med.histo() << std::endl;
+     std::cout << "med = " << med.to_result() << std::endl;
+
+    return float(med.to_result()) / delta;
+  }
+
+  
+  void refine_gap(const util::array<unsigned>& line, unsigned delta, float& gap)
+  {
+     std::cout << "gap = " << gap << std::endl;
+    unsigned count = 0;
+
+    accu::stat::mean<float> mean;
+    const unsigned n = line.nelements() - delta;
+    for (unsigned i = 0; i < n; ++i)
+      {
+	float g = line[i + delta] - line[i];
+	if (g / delta > gap - .2 && g / delta < gap + .2)
+	  {
+	    mean.take(g);
+	    ++count;
+	  }
+      }
+     std::cout << "count = " << count << "  total = " << n << std::endl;
+
+    gap = mean.to_result() / delta;
+     std::cout << "gap = " << gap << std::endl;
+  }
+
+
+  util::array<unsigned>
+  select_lines(const util::array<unsigned>& line, float gap)
+  {
+    std::set<unsigned> s;
+    const unsigned n = line.nelements() - 1;
+    for (unsigned i = 0; i < n; ++i)
+      {
+	float g = line[i + 1] - line[i];
+	if (g > gap - 2 && g < gap + 2)
+	  {
+	    s.insert(line[i]);
+	    s.insert(line[i + 1]);
+	  }
+      }
+
+    std::vector<unsigned> v(s.begin(), s.end());
+    util::array<unsigned> out;
+    out.hook_std_vector_() = v;
+
+    return out;
+  }
+
+
+  util::array<unsigned>
+  image_2_line(const image1d<double>& input)
+  {
+    image1d<bool> is_max(input.domain());
+    unsigned n;
+    data::fill(is_max,
+	       pw::value(labeling::regional_maxima(input, c2(), n)) | input.domain());
+    std::cout << "n = " << n << std::endl;
+
+    util::array<unsigned> line;
+    {
+      mln_piter_(box1d) p(input.domain());
+      for_all(p)
+	if (is_max(p))
+	  line.append(p.ind());
+    }
+
+    return line;
+  }
+
+
+  util::array<unsigned>
+  analyze(const util::array<unsigned>& line, float& gap, bool& ok)
+  {
+    // FIXME: The value of delta in 'get_approx_gap' and in
+    // 'refine_gap' should change w.r.t. to the nature of input data.
+    gap = get_approx_gap(line, 10); // FIXME: 10 / 3
+    float rough_gap = gap; 
+    refine_gap(line, 20, gap);      // FIXME: 20 / 4
+    float finer_gap = gap;
+
+    ok = std::abs(finer_gap - rough_gap) < 1;
+
+    return select_lines(line, gap);
+  }
+
+
+  void draw(const image2d<value::int_u8>& input,
+	    const util::array<unsigned>& line,
+	    value::rgb8 color,
+	    const std::string& filename)
+  {
+    image2d<value::rgb8> output = data::convert(value::rgb8(), input);
+    const int nrows = input.nrows();
+
+    for (unsigned i = 0; i < line.nelements(); ++i)
+      {
+	int col = line[i];
+	for (int row = 0; row < nrows; ++row)
+	  output.at_(row, col) = color;
+      }
+  
+    io::ppm::save(output, filename);
+}
+
+
+} // mln
+
+
+
+void usage(char* argv[])
+{
+  std::cerr << "usage: " << argv[0] << " input.pgm lambda h.txt output.ppm" << std::endl;
+  std::cerr << "hint: lambda = 50; delta = 10 / 20 otherwise 3 / 4" << std::endl;
+  abort();
+}
+
+
+
+int main(int argc, char* argv[])
+{
+  using namespace mln;
+  using value::int_u8;
+  using value::rgb8;
+
+  if (argc != 5)
+    usage(argv);
+
+  image2d<int_u8> input;
+  io::pgm::load(input, argv[1]);
+  box2d b = input.domain();
+
+
+  // h computing
+
+  typedef accu::stat::mean<int_u8, unsigned, double> A;
+  
+  image1d<A> hmean(b.ncols());
+  accu::image::init(hmean);
+
+  fun::v2v::projection<point2d, 0> hproj;
+  accu::image::take(unproject(hmean, b, hproj).rw(), input);
+
+  image1d<double> h = accu::image::to_result(hmean);
+  io::plot::save(h, argv[3]);
+
+
+  // h cleaning
+
+  h = linear::gaussian_1d(h, 1.0, 255);
+  io::plot::save(h, "h_g.txt");
+
+  int lambda = std::atoi(argv[2]);
+  h = morpho::closing::volume(h, c2(), lambda);
+  io::plot::save(h, "h_gc.txt");
+
+
+
+  // analyzing
+
+  util::array<unsigned> line = image_2_line(h);
+  draw(input, line, value::rgb8(255,0,0), "temp_h.ppm");
+
+  float gap;
+  bool ok;
+  line = analyze(line, gap, ok);
+
+  draw(input, line, value::rgb8(255,0,0), "temp_h_sel.ppm");
+
+  if (! ok)
+    {
+      std::cerr << "analyze failed: abort!" << std::endl;
+      std::abort();
+    }
+
+
+
+  // getting a math model
+  float o;
+
+  {
+    accu::stat::mean<float> me;
+    // col = r * gap  where  r = o + i (with offset o and integer i)
+    // col = gap * i + orig (orig = o * gap)
+    for (unsigned l = 0; l < line.nelements(); ++l)
+      {
+	float r = float(line[l]) / gap;
+	int i = int(r + 0.49999);
+	float o = r - float(i);
+ 	me.take(o);
+	std::cout << line[l] << " = " << (gap * o) << " + " << i << " * " << gap << std::endl;
+      }
+    o = me.to_result();
+
+    std::cout << "gap = " << gap << "  offset = " << (gap * o) << std::endl;
+
+    // verif
+    util::array<unsigned> line_;
+    for (unsigned l = 0; l < line.nelements(); ++l)
+      {
+ 	int i = int(float(line[l]) / gap - o + 0.49999);
+ 	line_.append(int((o + i) * gap + 0.49999));
+      }
+
+    std::cout << std::endl;
+
+    // over-display
+    draw(input, line, value::rgb8(255,0,0), "temp_h_found.ppm");
+
+  }
+
+
+  // drawing matrix
+
+  {
+    image2d<rgb8> output = data::convert(rgb8(), input);
+    const int nrows = input.nrows();
+    int ncols = input.ncols();
+
+    int n_max = int(float(ncols) / gap - o + 3);
+    for (int i = 0; i < n_max; ++i)
+      {
+	float col_ = (o + i) * gap;
+	int col = int(col_ + 0.49999);
+	if (col < 0)
+	  continue;
+	if (col > ncols - 1)
+	  break;
+	for (int row = 0; row < nrows; ++row)
+	  output.at_(row, col) = rgb8(0,255,0);
+      }
+
+    io::ppm::save(output, argv[4]);
+  }
+
+
+
+}
Index: theo/experimental/depeche/row_mm.cc
--- theo/experimental/depeche/row_mm.cc	(revision 0)
+++ theo/experimental/depeche/row_mm.cc	(revision 0)
@@ -0,0 +1,96 @@
+#include <mln/core/image/image2d.hh>
+#include <mln/core/alias/neighb2d.hh>
+
+#include <mln/core/image/image1d.hh>
+#include <mln/core/alias/neighb1d.hh>
+
+#include <mln/fun/v2v/projection.hh>
+#include <mln/core/image/dmorph/unproject_image.hh>
+#include <mln/core/image/dmorph/image_if.hh>
+#include <mln/pw/value.hh>
+
+#include <mln/accu/stat/mean.hh>
+#include <mln/accu/image/init.hh>
+#include <mln/accu/image/take.hh>
+#include <mln/accu/image/to_result.hh>
+
+#include <mln/value/int_u8.hh>
+#include <mln/io/pgm/load.hh>
+#include <mln/io/pgm/save.hh>
+#include <mln/io/plot/save.hh>
+
+#include <mln/value/rgb8.hh>
+#include <mln/io/ppm/save.hh>
+// #include <mln/draw/line.hh>
+// #include <mln/literal/colors.hh>
+
+#include <mln/morpho/closing/volume.hh>
+#include <mln/morpho/opening/volume.hh>
+
+#include <mln/labeling/regional_minima.hh>
+#include <mln/labeling/regional_maxima.hh>
+#include <mln/data/convert.hh>
+
+
+
+void usage(char* argv[])
+{
+  std::cerr << "usage: " << argv[0] << " input.pgm lambda output.pgm" << std::endl;
+  abort();
+}
+
+
+
+int main(int argc, char* argv[])
+{
+  using namespace mln;
+  using value::int_u8;
+  using value::rgb8;
+
+  if (argc != 4)
+    usage(argv);
+
+  image2d<int_u8> input;
+  io::pgm::load(input, argv[1]);
+  box2d b = input.domain();
+
+  int lambda = std::atoi(argv[2]);
+
+  typedef accu::stat::mean<int_u8, unsigned, int_u8> A;
+  
+  image1d<A> hmean(b.nrows());
+  accu::image::init(hmean);
+
+  fun::v2v::projection<point2d, 1> hproj;
+  accu::image::take(unproject(hmean, b, hproj).rw(), input);
+
+  image1d<int_u8> h = accu::image::to_result(hmean);
+  io::plot::save(h, "temp_h0.txt");
+  
+  h = morpho::closing::volume(h, c2(), lambda);
+  io::plot::save(h, "temp_h1.txt");
+
+  h = morpho::opening::volume(h, c2(), lambda);
+  io::plot::save(h, "temp_h2.txt");
+
+
+  {
+    int_u8 n;
+    image1d<int_u8>
+      lab_bg = labeling::regional_maxima(h, c2(), n),
+      lab_fg = labeling::regional_minima(h, c2(), n);
+
+    image2d<rgb8> cool = data::convert(rgb8(), input);
+    mln_piter_(box2d) p(cool.domain());
+    for_all(p)
+    {
+      if (lab_bg.at_(p.row()) != 0)
+	cool(p).red() = 255; 
+      if (lab_fg.at_(p.row()) != 0)
+	cool(p).green() = 255; 
+    }
+    io::ppm::save(cool, "temp_cool.ppm");
+  }
+
+  // io::pgm::save(output, argv[2]);
+}
Index: theo/experimental/depeche/row.cc
--- theo/experimental/depeche/row.cc	(revision 0)
+++ theo/experimental/depeche/row.cc	(revision 0)
@@ -0,0 +1,117 @@
+#include <mln/core/image/image2d.hh>
+#include <mln/core/alias/neighb2d.hh>
+
+#include <mln/core/image/image1d.hh>
+#include <mln/core/alias/neighb1d.hh>
+
+#include <mln/fun/v2v/projection.hh>
+#include <mln/core/image/dmorph/unproject_image.hh>
+#include <mln/core/image/dmorph/image_if.hh>
+#include <mln/pw/value.hh>
+
+#include <mln/accu/stat/mean.hh>
+#include <mln/accu/image/init.hh>
+#include <mln/accu/image/take.hh>
+#include <mln/accu/image/to_result.hh>
+
+#include <mln/value/int_u8.hh>
+#include <mln/io/pgm/load.hh>
+#include <mln/io/pgm/save.hh>
+#include <mln/io/plot/save.hh>
+
+#include <mln/value/rgb8.hh>
+#include <mln/io/ppm/save.hh>
+
+#include <mln/linear/gaussian_1d.hh>
+#include <mln/data/fill.hh>
+#include <mln/data/convert.hh>
+
+#include <mln/labeling/regional_minima.hh>
+#include <mln/draw/line.hh>
+#include <mln/literal/colors.hh>
+
+#include <mln/accu/stat/median_h.hh>
+
+
+
+void usage(char* argv[])
+{
+  std::cerr << "usage: " << argv[0] << " input.pgm output.ppm" << std::endl;
+  abort();
+}
+
+
+
+int main(int argc, char* argv[])
+{
+  using namespace mln;
+  using value::int_u8;
+  using value::rgb8;
+
+  if (argc != 3)
+    usage(argv);
+
+  image2d<int_u8> input;
+  io::pgm::load(input, argv[1]);
+  box2d b = input.domain();
+
+  typedef accu::stat::mean<int_u8, unsigned, double> A;
+    
+  image1d<A> hmean(b.nrows());
+  accu::image::init(hmean);
+
+  fun::v2v::projection<point2d, 1> hproj;
+  accu::image::take(unproject(hmean, b, hproj).rw(), input);
+
+  image1d<double> hmean_d = accu::image::to_result(hmean);
+  io::plot::save(hmean_d, "temp_h.txt");
+
+  hmean_d = linear::gaussian_1d(hmean_d, 5, 255);
+  io::plot::save(hmean_d, "temp_hg.txt");
+
+  io::plot::save(hmean, argv[2]);
+
+  int_u8 n;
+  image1d<int_u8> lab = labeling::regional_minima(hmean_d, c2(), n);
+  // io::plot::save(lab, "temp_lab.txt");
+
+
+  {
+    image2d<rgb8> cool = data::convert(rgb8(), input);
+    mln_piter_(box1d) p(lab.domain());
+    for_all(p)
+      if (lab(p) != 0)
+	draw::line(cool,
+		   point2d(p.ind(), 0),
+		   point2d(p.ind(), b.ncols() -1),
+		   literal::red);
+    io::ppm::save(cool, argv[2]);
+  }
+
+//   accu::stat::median_h<int_u8> med;
+//   {
+//     double min = 255;
+//     unsigned row_min;
+//     int last = -1;
+//     mln_piter_(box1d) p(lab.domain());
+//     for_all(p)
+//     {
+//       if (lab(p) == 0)
+// 	continue;
+//       if (hmean_d(p) < min)
+// 	{
+// 	  min = hmean_d(p);
+// 	  row_min = p.ind();
+// 	}
+//       if (last == -1)
+// 	last = p.ind();
+//       else
+// 	{
+// 	  med.take(p.ind() - last);
+// 	  last = p.ind();
+// 	}
+//     }
+//     std::cout << med << ' ' << row_min << std::endl;
+//   }
+
+}

    