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

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-03-06 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; 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