[Olena-patches] r3496: Add watershed implementation for IRM images

6 Mar 2009

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

ChangeLog:
2009-03-06  Fabien Freling  <fabien.freling@lrde.epita.fr>

	Add watershed implementation for IRM images.
	* fabien/bin/Makefile: Update.
	* fabien/bin/dicom2dump.cc: Update.
	* fabien/bin/dicom_mask.cc: Update.
	* fabien/igr/Makefile: Update.
	* fabien/igr/dumps: Remove.
	* fabien/igr/watershed.cc: Watershed implementation.

---
 bin/Makefile      |    3 
 bin/dicom2dump.cc |   19 ---
 bin/dicom_mask.cc |   26 ++++
 igr/Makefile      |    9 -
 igr/watershed.cc  |  320 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 5 files changed, 353 insertions(+), 24 deletions(-)

Index: trunk/milena/sandbox/fabien/igr/watershed.cc
===================================================================

--- trunk/milena/sandbox/fabien/igr/watershed.cc	(revision 0)
+++ trunk/milena/sandbox/fabien/igr/watershed.cc	(revision 3496)
@@ -0,0 +1,320 @@
+#include <iostream>
+#include <mln/core/image/image2d.hh>
+
+#include <mln/core/alias/neighb2d.hh>
+#include <mln/core/alias/window2d.hh>
+#include <mln/core/image/image_if.hh>
+
+#include <mln/io/ppm/save.hh>
+#include <mln/io/ppm/load.hh>
+#include <mln/io/pgm/load.hh>
+#include <mln/io/dicom/load.hh>
+#include <mln/io/pgm/save.hh>
+#include <mln/io/dump/save.hh>
+
+#include <mln/value/rgb8.hh>
+#include <mln/value/int_u8.hh>
+#include <mln/value/int_u12.hh>
+#include <mln/value/label_16.hh>
+
+#include <mln/level/transform.hh>
+#include <mln/level/stretch.hh>
+
+#include <mln/labeling/mean_values.hh>
+
+#include <mln/convert/to_fun.hh>
+
+#include <mln/make/graph.hh>
+
+#include <mln/morpho/gradient.hh>
+#include <mln/morpho/closing/area.hh>
+#include <mln/morpho/meyer_wst.hh>
+
+#include <mln/fun/l2l/wrap.hh>
+
+#include <mln/core/var.hh>
+#include <mln/morpho/elementary/dilation.hh>
+
+#include <mln/core/routine/extend.hh>
+
+#include <mln/util/graph.hh>
+
+#include <mln/essential/2d.hh>
+#include <mln/core/alias/vec3d.hh>
+#include <mln/debug/draw_graph.hh>
+#include <mln/util/graph.hh>
+#include <mln/accu/center.hh>
+#include <mln/io/dump/all.hh>
+#include <mln/value/label_8.hh>
+#include <mln/value/rgb16.hh>
+#include <mln/accu/compute.hh>
+#include <mln/core/alias/dpoint2d.hh>
+#include <mln/draw/box.hh>
+#include <mln/level/stretch.hh>
+#include <mln/fun/v2v/id.hh>
+#include <mln/core/image/line_graph_elt_neighborhood.hh>
+#include <mln/morpho/elementary/dilation.hh>
+#include <mln/labeling/mean_values.hh>
+#include <mln/extension/adjust_fill.hh>
+#include <mln/extract/all.hh>
+#include <mln/make/region_adjacency_graph.hh>
+
+
+
+// Given a color image and a wshed image, computes the component graph.
+// Vertex values are computed thanks to a RGB image.
+
+namespace mln
+{
+
+  template <typename V>
+  value::int_u8 dist(const V& c1, const V& c2)
+  {
+    unsigned d = math::diff_abs(c1.red(), c2.red());
+    unsigned d_;
+    d_ = math::diff_abs(c1.green(), c2.green());
+
+    if (d_ > d)
+      d = d_;
+
+    d_ = math::diff_abs(c1.blue(), c2.blue());
+
+    if (d_ > d)
+      d = d_;
+    return d;
+  }
+
+
+  // ima_v, image on graph vertices; value = mean color per vertex (watershed basin)
+  template <typename I>
+  inline
+  pw::image<fun::i2v::array<value::int_u8>,
+	    p_edges<util::graph, fun::i2v::array<mln_site(I)> > >
+  make_edge_graph_image(const I& ima_v, const util::graph& g)
+  {
+    // edge sites.
+    typedef fun::i2v::array<mln_site(I)> edge_site_t;
+    edge_site_t edge_site(g.e_nmax(), literal::origin);
+    typedef p_edges<util::graph, edge_site_t > pe_t;
+    pe_t pe(g, edge_site);
+
+    // edge values
+    typedef fun::i2v::array<value::int_u8> edge_values_t;
+    edge_values_t edge_values(g.e_nmax());
+
+    // image on graph edges
+    typedef pw::image<edge_values_t, pe_t> ima_e_t;
+    ima_e_t ima_e = (edge_values | pe);
+
+    mln_piter(ima_e_t) e(ima_e.domain());
+    for_all(e) // in ima_e
+      ima_e(e) = math::diff_abs(ima_v.function()(e.element().v1()), ima_v.function()(e.element().v2()));
+
+    return ima_e;
+  }
+
+
+  template <typename I, typename J>
+  inline
+  util::array<mln_value(I)>
+  mean_color_values(const I& input, const J& w, mln_value(J) nbasins)
+  {
+    // Cf. sandbox/theo/color/segment_rgb_pixels.cc
+
+    util::array<vec3d_f> m_3f = labeling::compute(accu::mean<mln_value(I)>(),
+						  input, // input color image
+						  w, // watershed labeling
+						  nbasins);
+    m_3f[0] = literal::zero;
+
+    util::array<mln_value(I)> m;
+    convert::from_to(m_3f, m);
+    m[0] = literal::yellow;
+
+    io::ppm::save(level::transform(w,
+	  convert::to< fun::i2v::array<mln_value(I)> >(m)),
+	"wst_rag_wshd_color.ppm");
+
+    return m;
+  }
+
+
+  template <typename I, typename J>
+  pw::image<fun::i2v::array<mln_value(I)>, p_vertices<util::graph, fun::i2v::array<mln_site(I)> > >
+  make_vertex_graph_image(const util::graph& g, const I&input, const J& w, const mln_value(J)& nbasins)
+  {
+    typedef util::array<mln_site(I)> vertex_sites_t;
+    vertex_sites_t site_values;
+    convert::from_to(labeling::compute(accu::center<mln_site(I)>(), w, nbasins),
+				       site_values);
+
+    typedef fun::i2v::array<mln_site(J)> f_sites_t;
+    f_sites_t sites;
+    convert::from_to(site_values, sites);
+
+    // p_vertices
+    typedef p_vertices<util::graph, f_sites_t> S;
+    S pv(g, sites);
+
+
+    typedef fun::i2v::array<mln_value(I)> vertex_values_t;
+    vertex_values_t vertex_values;
+    // convert::from_to(mean_color_values(input, w, nbasins), vertex_values);
+
+    mln_VAR(ima_v, (vertex_values | pv));
+
+    return ima_v;
+  }
+
+
+  template <typename I, typename V>
+  struct edge_to_color : Function_p2v< edge_to_color<I,V> >
+  {
+    typedef V result;
+
+    edge_to_color(const I& ima)
+      : ima_(ima)
+    {
+    }
+
+    V
+    operator()(const unsigned& e) const
+    {
+      return convert::to<V>(ima_.function()(e));
+    }
+
+    const I& ima_;
+  };
+
+  template <typename V>
+  inline
+  unsigned
+  find_root(util::array<V>& parent, const unsigned& x)
+  {
+    if (parent[x] == x)
+      return x;
+    else
+      return parent[x] = find_root(parent, parent[x]);
+  }
+
+
+  template <typename I, typename V, typename E>
+  inline
+  image2d<mln_value(I)>
+  make_debug_graph_image(const I& input,
+			 const V& ima_v, const E& ima_e,
+			 unsigned box_size, const value::rgb8& bg)
+  {
+    image2d<mln_value(I)> ima;
+    initialize(ima, input);
+
+    data::fill(ima, bg);
+    debug::draw_graph(ima, ima_v.domain(),
+		      pw::cst(mln_value(I)(literal::green)),
+		      edge_to_color<E, mln_value(I)>(ima_e));
+
+    dpoint2d tl(-box_size,-box_size);
+    dpoint2d br(box_size,box_size);
+    mln_piter(V) p(ima_v.domain());
+    for_all(p)
+    {
+      box2d b(p + tl, p + br);
+      b.crop_wrt(ima.domain());
+      data::fill((ima | b).rw(), convert::to<mln_value(I)>(ima_v(p)));
+    }
+
+    return ima;
+  }
+
+}
+
+
+///////////////////
+//               //
+// Main Function //
+//               //
+///////////////////
+
+
+int main(int argc, char *argv[])
+{
+  using namespace mln;
+  using value::int_u8;
+  using value::int_u12;
+  using value::rgb8;
+  using value::label_16;
+
+  if (argc < 4)
+  {
+    std::cout << "Usage: " << argv[0] << " <ima.dcm> <closure_lambda> <box_size>"
+	      << std::endl;
+    return 1;
+  }
+
+  unsigned box_size = atoi(argv[3]);
+
+  image2d<int_u12> dcm;
+  io::dicom::load(dcm, argv[1]);
+
+  image2d<int_u12> grad = morpho::gradient(dcm, win_c4p());
+  image2d<int_u12> clo = morpho::closing::area(grad, c4(), atoi(argv[2]));
+
+  label_16 nbasins;
+  image2d<label_16> wshed = morpho::meyer_wst(clo, c4(), nbasins);
+
+  mln_VAR(vol2_, morpho::elementary::dilation(extend(wshed | (pw::value(wshed) == 0u), wshed), c8()));
+
+  data::fill((wshed | (pw::value(wshed) == 0u)).rw(), vol2_);
+
+  /// Build graph
+  util::graph g = make::graph(wshed, c4(), nbasins);
+  // Build graph images and compute distance values with a RGB image.
+  mln_VAR(ima_v, make_vertex_graph_image(g, dcm, wshed, nbasins));
+  mln_VAR(ima_e, make_edge_graph_image(ima_v, g));
+
+  /// Try to merge vertices.
+  mln_piter_(ima_e_t) e(ima_e.domain());
+  util::array<label_16> parent(g.v_nmax());
+  for (unsigned i = 0; i < parent.nelements(); ++i)
+    parent[i] = i;
+
+  for_all(e)
+  {
+    unsigned v1 = e.element().v1();
+    unsigned v2 = e.element().v2();
+    if (ima_e(e) <= (unsigned)(atoi(argv[5]))
+	&& find_root(parent, v1) != find_root(parent, v2))
+      parent[find_root(parent, v1)] = find_root(parent, v2);
+  }
+
+  fun::i2v::array<label_16> f(parent.nelements());
+  std::vector<unsigned> new_label(parent.nelements(), 0);
+  unsigned nbasins2 = 0;
+  for (unsigned i = 0; i < parent.nelements(); ++i)
+  {
+    unsigned p = find_root(parent, i);
+    mln_assertion(parent[p] == find_root(parent, i));
+    if (new_label[p] == 0)
+      new_label[p] = nbasins2++;
+    f(i) = new_label[p];
+  }
+  mln_invariant(f(0) == 0u);
+  --nbasins2; // nbasins2 does not count the basin with label 0.
+  image2d<label_16> wsd2 = level::transform(wshed, f);
+
+
+  /// Reconstruct a graph from the simplified image.
+  util::graph g2 = make::graph(wsd2, c4(), nbasins2);
+
+  // Compute distance values with a RGB image.
+  mln_VAR(ima_v2, make_vertex_graph_image(g2, dcm, wsd2, nbasins2));
+  mln_VAR(ima_e2, make_edge_graph_image(ima_v2, g2));
+
+  mln_VAR(wsd2_, morpho::elementary::dilation(extend(wsd2 | (pw::value(wsd2) == 0u), wsd2), c8()));
+
+  data::fill((wsd2 | (pw::value(wsd2) == 0u)).rw(), wsd2_);
+
+  //io::ppm::save(labeling::mean_values(dcm, wsd2, nbasins2), "wst_rag_mean_colors.ppm");
+  //io::ppm::save(make_debug_graph_image(dcm, ima_v2, ima_e2, box_size, literal::white), "wst_rag_graph_image2_white.ppm");
+  //io::ppm::save(make_debug_graph_image(dcm, ima_v2, ima_e2, box_size, literal::black), "wst_rag_graph_image2_black.ppm");
+}
Index: trunk/milena/sandbox/fabien/igr/Makefile
===================================================================
--- trunk/milena/sandbox/fabien/igr/Makefile	(revision 3495)
+++ trunk/milena/sandbox/fabien/igr/Makefile	(revision 3496)
@@ -8,7 +8,7 @@
 	    -framework CoreFoundation
 CXXFLAGS = -DNDEBUG -O1
 
-all: 2d 3d
+all: 2d 3d wsd
 
 2d: seg_vol_irm.hh seg2d.cc
 	g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} seg2d.cc -o seg2d
@@ -16,11 +16,8 @@
 3d: seg_vol_irm.hh seg3d.cc
 	g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} seg3d.cc -o seg3d
 
-grad: grad_clo_and_wshd.cc
-	g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} $^ -o grad_clo
-
-wst: wst_rag.cc
-	g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} $^ -o wst_rag
+wsd: watershed.cc
+	g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} $^ -o wsd
 
 clean:
 	rm -rf *.dump *.p?m *.plot *.log *.csv
Index: trunk/milena/sandbox/fabien/bin/dicom2dump.cc
===================================================================
--- trunk/milena/sandbox/fabien/bin/dicom2dump.cc	(revision 3495)
+++ trunk/milena/sandbox/fabien/bin/dicom2dump.cc	(revision 3496)
@@ -2,7 +2,6 @@
 #include <mln/core/image/image2d.hh>
 #include <mln/core/image/image3d.hh>
 
-#include <mln/value/int_u8.hh>
 #include <mln/value/int_u12.hh>
 #include <mln/io/dicom/load.hh>
 #include <mln/io/dump/save.hh>
@@ -17,31 +16,15 @@
 }
 
 
-/*int main(int argc, char* argv[])
-{
-  using namespace mln;
-  using value::int_u12;
-
-  if (argc != 3)
-    return usage(argv);
-
-  image3d<int_u12> ima;
-  io::dicom::load(ima, argv[1]);
-  io::dump::save(ima, argv[2]);
-
-  return 0;
-}*/
-
 int main(int argc, char* argv[])
 {
   using namespace mln;
-  using value::int_u8;
   using value::int_u12;
 
   if (argc != 3)
     return usage(argv);
 
-  image2d<int_u8> ima;
+  image3d<int_u12> ima;
   io::dicom::load(ima, argv[1]);
   io::dump::save(ima, argv[2]);
 
Index: trunk/milena/sandbox/fabien/bin/dicom_mask.cc
===================================================================
--- trunk/milena/sandbox/fabien/bin/dicom_mask.cc	(revision 3495)
+++ trunk/milena/sandbox/fabien/bin/dicom_mask.cc	(revision 3496)
@@ -1,6 +1,9 @@
 #include <mln/core/concept/image.hh>
 #include <mln/core/image/image2d.hh>
 
+#include <mln/geom/ncols.hh>
+#include <mln/geom/nrows.hh>
+
 #include <mln/value/int_u8.hh>
 #include <mln/io/dicom/load.hh>
 #include <mln/io/pbm/save.hh>
@@ -33,6 +36,29 @@
   image2d<int_u8> input;
   io::dicom::load(input, argv[1]);
 
+  util::array<unsigned> xproj(geom::nrows(input), 0);
+  util::array<unsigned> yproj(geom::ncols(input), 0);
+
+  mln_piter_(image2d<int_u8>) p(input.domain());
+  for_all(p)
+  {
+    xproj[p.row()] += input(p);
+    yproj[p.col()] += input(p);
+  }
+
+  // Plot files
+  std::ofstream fout_x("x.plot");
+  for (unsigned int i = 0; i < xproj.nelements(); ++i)
+  {
+    fout_x << i << " " << xproj[i] << std::endl;
+  }
+
+  std::ofstream fout_y("y.plot");
+  for (unsigned int i = 0; i < yproj.nelements(); ++i)
+  {
+    fout_y << i << " " << yproj[i] << std::endl;
+  }
+
   image2d<bool> ima = level::transform(input, fun::v2b::threshold<int_u8>(1));
   io::pbm::save(ima, argv[2]);
 }
Index: trunk/milena/sandbox/fabien/bin/Makefile
===================================================================
--- trunk/milena/sandbox/fabien/bin/Makefile	(revision 3495)
+++ trunk/milena/sandbox/fabien/bin/Makefile	(revision 3496)
@@ -21,3 +21,6 @@
 
 dump: dicom2dump.cc
 	g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} $^ -o dicom2dump
+
+pgm: dicom2pgm.cc
+	g++ -I../../../ ${DICOM_INC} ${DICOM_LIB} ${CXXFLAGS} $^ -o dicom2pgm

    