https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox The code is probably finished (and compiles!), but is buggy: the program gets a SIGSEGV at run time... Index: ChangeLog from Roland Levillain &lt;roland(a)lrde.epita.fr&gt; Finish the computation of the fusion graph test. * beguin/irm.cc: s/node/vertex/. (my_mln::mk2_graph): Disable. (my_mln::region_data, my_mln::abs_diff_threshold): New structs. (my_mln::make_fg, my_mln::merge_fg): New routines. (my_mln::main): Adjust. Compute the fusion graphs and output a result. * beguin/+irm6.pgm: Rename image as... * beguin/irm6.pgm: ...this. irm.cc | 307 +++++++++++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 253 insertions(+), 54 deletions(-) Index: beguin/irm.cc --- beguin/irm.cc (revision 1995) +++ beguin/irm.cc (working copy) @@ -1,6 +1,11 @@ +// g++ -Wall -Wextra -DNDEBUG -O3 -I../.. -o irm irm.cc +// g++ -Wall -Wextra -ggdb -I../.. -o irm irm.cc + #include <iterator> #include <iostream> #include <algorithm> + +#include <set> #include <map> #include <mln/core/image2d.hh> @@ -12,8 +17,12 @@ #include <mln/io/pgm/load.hh> #include <mln/io/pgm/save.hh> +#include <mln/io/ppm/save.hh> #include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> +#include <mln/literal/colors.hh> + #include <mln/level/transform.hh> #include <mln/convert/to_window.hh> @@ -32,7 +41,7 @@ #include <mln/debug/println.hh> //#include <segm_to_pregraph.hh> -#include <fusion_graph2.hh> +// #include <fusion_graph2.hh> mln::value::int_u8 foo(unsigned u) { @@ -42,9 +51,10 @@ } -namespace mln - +namespace my_mln { + using namespace mln; + template <typename I> void region_10(const I& lbl) @@ -163,12 +173,12 @@ // // // Graph. // util::graph<float> g; -// // Nodes. +// // Vertices. // for (unsigned i = 0; i <= nlabels; ++i) // { // accu::mean_<unsigned> m; // level::take(irm | (pw::value(lbl) == pw::cst(i)), m); -// g.add_node(m.to_result()); +// g.add_vertex(m.to_result()); // } // // Edges. // for (unsigned l1 = 1; l1 <= nlabels; ++l1) @@ -176,11 +186,16 @@ // if (adja[l1][l2]) // g.add_edge(l1, l2); // //g.print_debug (std::cout); -// for( int i = 0 ; i < g.nnodes() ; ++i) -// std::cout << "[" << i << " : " << g.node_data(i) << " ] "; +// for( int i = 0 ; i < g.nvertices() ; ++i) +// std::cout << "[" << i << " : " << g.vertex_data(i) << " ] "; // std::cout << std::endl; // } + + /*------------------------. + | Region adjacent graph. | + `------------------------*/ + /// Build a region adjacency graph. template <typename I, typename J, typename N> util::graph< accu::mean_<mln_value(J)> > @@ -190,19 +205,19 @@ for (unsigned l = 0; l <= nlabels; ++l) adja[l].resize(nlabels + 1, false); - mln_piter(I) p(lbl.domain()); + mln_piter(I) p(label.domain()); mln_niter(N) n(nbh, p); // We compute the adjacency matrix of the RAG. for_all(p) - if (lbl(p) == 0) // wshed + if (label(p) == 0) // wshed { std::set<unsigned> s; - for_all(n) if (lbl.has(n)) + for_all(n) if (label.has(n)) { - if (lbl(n) == 0) + if (label(n) == 0) continue; - s.insert(lbl(n)); + s.insert(label(n)); } if (s.size() < 2) { @@ -227,12 +242,12 @@ // Graph. util::graph< accu::mean_<mln_value(J)> > g; - // Nodes. + // Vertices. for (unsigned i = 0; i <= nlabels; ++i) { accu::mean_<mln_value(J)> m; - level::take(irm | (pw::value(lbl) == pw::cst(i)), m); - g.add_node(m); + level::take(input | (pw::value(label) == pw::cst(i)), m); + g.add_vertex(m); } // Edges. for (unsigned l1 = 1; l1 <= nlabels; ++l1) @@ -240,57 +255,210 @@ if (adja[l1][l2]) g.add_edge(l1, l2); //g.print_debug (std::cout); - for( int i = 0 ; i < g.nnodes() ; ++i) - std::cout << "[" << i << " : " << g.node_data(i) << " ] "; + for( unsigned i = 0 ; i < g.nvertices() ; ++i) + std::cout << "[" << i << " : " << g.vertex_data(i) << " ] "; std::cout << std::endl; return g; } - template <typename I, typename J, typename N> - void mk2_graph(const I& lbl , const J& irm , const N& - nbh, int nbasins) +// template <typename I, typename J, typename N> +// void mk2_graph(const I& label , const J& input , const N& +// nbh, int nbasins) +// { +// typedef mln::p_set<mln_point(I)> R; +// std::map<int, R > regions; + +// mln_piter(I) p(label.domain()); +// for_all(p) +// { +// if (unsigned int i = label(p)) +// regions[i].insert(p); +// } + +// std::cout << "Boite englobante de la r�gion 10 : " << regions[10].bbox() << std::endl; + +// /* +// for (unsigned i = 0 ; i<nbasins ; ++i) +// { + +// } +// for (unsigned i = 0 ; i<regions.size() ; ++i) +// g.add_vertex(regions[]); +// */ +// } + + + /*---------------. + | Fusion graph. | + `---------------*/ + + /// The data associated to the vertex of a fusion graph. + template <typename A> + struct region_data { - typedef mln::p_set<mln_point(I)> R; - std::map<int, R > regions; + // The subregions composing the region. + std::set<util::vertex_id> subregions; + // An accumulator. + A accumulator; + }; - mln_piter(I) p(lbl.domain()); - for_all(p) + + /// Build a fusion graph from a RAG. + template <typename A> + util::graph< region_data<A> > + make_fg(const util::graph<A>& rag) + { + typedef util::graph<A> region_adjacency_graph; + typedef util::graph< region_data<A> > fusion_graph; + + fusion_graph fg; + + // Vertices. + for (util::vertex_id i = 0; i < rag.nvertices(); ++i) + { + // The singleton subregion: { I }. + std::set<util::vertex_id> sr; + sr.insert(i); + // The accumulator. + A a = rag.vertex_data(i); + // The region data + region_data<A> data = { sr, a }; + + // Add the vertex at index I. + fg.add_vertex(data); + } + + // Edges. + for (typename region_adjacency_graph::edges_t::const_iterator e = + rag.edges().begin(); e != rag.edges().end(); ++e) + { + util::vertex_id v1 = (*e)->v1(); + util::vertex_id v2 = (*e)->v2(); + fg.add_edge(v1, v2); + } + + return fg; + } + + /// Merge regions of the fusion graph. + template <typename A, typename C> + util::graph< region_data<A> > + merge_fg(const util::graph< region_data<A> >& input, const C& criterion) + { + typedef util::graph< region_data<A> > fusion_graph; + + // A disjoint set of vertices represented as a forest of rooted + // trees (disjoint set forest). + std::vector<util::vertex_id> parent(input.nvertices()); + for (util::vertex_id i = 0; i < parent.size(); ++i) + // � make_set �. + parent[i] = i; + + // Parent. + for (typename fusion_graph::edges_t::const_iterator e = + input.edges().begin(); e != input.edges().end(); ++e) { - if (unsigned int i = lbl(p)) - regions[i].insert(p); + util::vertex_id v1 = (*e)->v1(); + util::vertex_id v2 = (*e)->v2(); + A accu1 = input.vertex_data(v1).accumulator; + A accu2 = input.vertex_data(v2).accumulator; + if (criterion(accu1.to_result(), accu2.to_result())) + // � union �. + parent[v1] = v2; } - std::cout << "Boite englobante de la région 10 : " << regions[10].bbox() << std::endl; + // Children. + typedef std::multimap<util::vertex_id, util::vertex_id> children_t; + children_t children; + for (typename util::vertex_id i = 0; i < parent.size(); ++i) + { + // � find �. + util::vertex_id r = i; + while (r != parent[r]) + r = parent[r]; + children.insert(std::make_pair(r, i)); + } + + /// Data associated to merged regions. + typedef std::map< util::vertex_id, region_data<A> > merged_region_data_t; + merged_region_data_t merged_region_data; + for (typename children_t::const_iterator i = children.begin(); + i != children.end(); ++i) + { + util::vertex_id root = i->first; + util::vertex_id child = i->second; + + region_data<A>& root_data = merged_region_data[root]; + const region_data<A>& child_data = input.vertex_data(child); + + root_data.subregions.insert(child_data.subregions.begin(), + child_data.subregions.end()); + root_data.accumulator.take(child_data.accumulator); + } - /* - for (unsigned i = 0 ; i<nbasins ; ++i) + // Fusion graph. + fusion_graph output; + std::map<util::vertex_id, util::vertex_id> old_to_new_root; + // Vertices. + for (typename merged_region_data_t::const_iterator i = + merged_region_data.begin(); i != merged_region_data.end(); ++i) + { + util::vertex_id v = output.add_vertex(i->second); + old_to_new_root[i->first] = v; + } + // Edges. + for (typename fusion_graph::edges_t::const_iterator e = + input.edges().begin(); e != input.edges().end(); ++e) { + // The root of the first vertex connected to E. + util::vertex_id r1 = (*e)->v1(); + while (r1 != parent[r1]) + r1 = parent[r1]; + // The root of the second vertex connected to E. + util::vertex_id r2 = (*e)->v2(); + while (r2 != parent[r2]) + r2 = parent[r2]; + if (r1 != r2) + { + // The regions represented by R1 and R2 are connected. + output.add_edge(old_to_new_root[r1], old_to_new_root[r2]); + } } - for (unsigned i = 0 ; i<regions.size() ; ++i) - g.add_node(regions[]); - */ + return output; } +} // end of namespace my_mln -} // end of namespace mln -template <typename A> -struct region_data +struct abs_diff_threshold { - // The subregions composing the region. - std::set<util::vertex_id> subregions; - // An accumulator. - A a; + abs_diff_threshold(int threshold) + : threshold_(threshold) + { + } + + template <typename T> + bool operator() (const T& a, const T& b) const + { + return std::abs(a - b) < threshold_; } +private: + int threshold_; +}; + + + int main() { using namespace mln; + using namespace my_mln; using value::int_u8; + using value::rgb8; image2d<int_u8> irm; - io::pgm::load(irm, "./+irm6.pgm"); + io::pgm::load(irm, "./irm6.pgm"); // io::pgm::save( morpho::gradient(irm, win::rectangle2d(3,3)), // "tmp_grad_3x3.pgm" ); @@ -316,22 +484,53 @@ //typedef fusion_graph<unsigned int, void> _fg; //_fg fg(wshed , irm , c8() , nbasins); -// fusion_graph fg(lbl&, irm&, c4(), nbasins); -// std::cout << "nnodes = " << fg.nnodes() << std::endl +// fusion_graph fg(label&, irm&, c4(), nbasins); +// std::cout << "nvertices = " << fg.nvertices() << std::endl // << "nregions = " << fg.nn() << std::endl; - // Region adjacency graph. - util::graph< accu::mean_<int_u8> > rag = make_rag(wshed, irm, c4(), nbasins); - - // Fusion graph. - - typedef util::graph< region_data< accu::mean<int_u8> > fusion_graph; - - fusion_graph output; + // Region adjacency graph (RAG). + typedef util::graph< accu::mean_<int_u8> > region_adjacency_graph; + region_adjacency_graph rag = make_rag(wshed, irm, c4(), nbasins); + + // Fusion graph (FG). + typedef util::graph< region_data< accu::mean_<int_u8> > > fusion_graph; + fusion_graph fg = make_fg(rag); fusion_graph work; do - work = output; - output = merge(work); - while (work.nvertices() != output.nvertices()); - + { + work = fg; + fg = merge_fg (work, abs_diff_threshold(20)); + std::cerr << "work.nvertices() = " << work.nvertices() << std::endl; + std::cerr << "fg.nvertices() = " << fg.nvertices() << std::endl; + } + while (work.nvertices() != fg.nvertices()); + + // Compute the values (means) for the output image. + std::map<util::vertex_id, int_u8> region_value; + for (util::vertex_id v = 0; v < fg.nvertices(); ++v) + { + // FIXME: Implicit cast. We should rather change the type of + // the accumulator instead. + int_u8 mean = fg.vertex_data(v).accumulator.to_result(); + for (std::set<util::vertex_id>::const_iterator sr = + fg.vertex_data(v).subregions.begin(); + sr != fg.vertex_data(v).subregions.end(); + ++sr) + region_value[*sr] = mean; + } + + // Output. + image2d<rgb8> output(wshed.domain()); + mln_piter_(image2d<rgb8>) p(output.domain()); + for_all (p) + { + if (wshed(p) == 0) + output(p) = literal::red; + else + { + int_u8 mean = region_value[wshed(p)]; + output(p) = rgb8(mean, mean, mean); + } + } + io::ppm::save(output, "out.ppm"); }