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

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; 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; +// } + +}