[Olena-patches] r3896: Create tools for processing chain on edge images

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-05-27 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; Create tools for processing chain on edge images. * fabien/igr/Makefile: Add targets for edges processing chain. * fabien/igr/clo_vol.cc: Small update. * fabien/igr/dist_edges.cc: Compute distance on edges, now with 3 distances. * fabien/igr/norm.cc: Small update. * fabien/igr/seg_fixed.cc: Remove. * fabien/igr/wst_edges.cc: Compute a volume closing and a watershed. --- Makefile | 17 ++-- clo_vol.cc | 2 dist_edges.cc | 214 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ norm.cc | 2 wst_edges.cc | 142 ++++++++++++++++++++++++++++++++++++++ 5 files changed, 367 insertions(+), 10 deletions(-) Index: trunk/milena/sandbox/fabien/igr/seg_fixed.cc (deleted) =================================================================== Index: trunk/milena/sandbox/fabien/igr/clo_vol.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/clo_vol.cc (revision 3895) +++ trunk/milena/sandbox/fabien/igr/clo_vol.cc (revision 3896) @@ -38,7 +38,7 @@ if (argc != 4) { - std::cout << "Usage: " << argv[0] << " <input.dump> <dimensions> <closure_lambda>" + std::cout << "Usage: " << argv[0] << " input.dump dim lambda output.dump" << std::endl; return 1; } Index: trunk/milena/sandbox/fabien/igr/norm.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/norm.cc (revision 3895) +++ trunk/milena/sandbox/fabien/igr/norm.cc (revision 3896) @@ -93,7 +93,7 @@ return 1; } - image3d<int_u12> ima; + image3d<float> ima; io::dump::load(ima, argv[1]); image3d<float> norm = normalize(ima, atoi(argv[2]), atoi(argv[3])); io::dump::save(norm, argv[4]); Index: trunk/milena/sandbox/fabien/igr/wst_edges.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/wst_edges.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/wst_edges.cc (revision 3896) @@ -0,0 +1,142 @@ +#include <algorithm> + +#include <mln/core/image/image1d.hh> +#include <mln/core/alias/window1d.hh> +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/core/image/image3d.hh> +#include <mln/core/image/slice_image.hh> +#include <mln/core/image/image_if.hh> +#include <mln/core/routine/duplicate.hh> +#include <mln/core/routine/extend.hh> +#include <mln/core/var.hh> + +#include <mln/io/dump/all.hh> +#include <mln/io/pgm/save.hh> +#include <mln/io/ppm/save.hh> + +#include <mln/value/int_u8.hh> +#include <mln/value/int_u12.hh> +#include <mln/value/label_16.hh> +#include <mln/value/rgb8.hh> + +#include <mln/accu/sum.hh> +#include <mln/accu/mean.hh> +#include <mln/accu/image/all.hh> +#include <mln/accu/stat/deviation.hh> +#include <mln/arith/div.hh> +#include <mln/data/fill.hh> +#include <mln/data/paste.hh> +#include <mln/debug/quiet.hh> +#include <mln/convert/from_to.hh> +#include <mln/fun/v2v/fit.hh> +#include <mln/labeling/compute.hh> +#include <mln/labeling/wrap.hh> +#include <mln/level/compute.hh> +#include <mln/level/convert.hh> +#include <mln/level/stretch.hh> +#include <mln/make/image2d.hh> +#include <mln/make/w_window1d.hh> +#include <mln/math/diff_abs.hh> +#include <mln/morpho/dilation.hh> +#include <mln/morpho/erosion.hh> +#include <mln/morpho/closing/volume.hh> +#include <mln/morpho/watershed/flooding.hh> +#include <mln/pw/all.hh> +#include <mln/util/array.hh> +#include <mln/win/segment1d.hh> + +#include <mln/world/inter_pixel/display_edge.hh> +#include <mln/world/inter_pixel/compute.hh> +#include <mln/world/inter_pixel/immerse.hh> +#include <mln/world/inter_pixel/neighb2d.hh> +#include <mln/world/inter_pixel/all.hh> + +#include <mln/labeling/colorize.hh> +#include <mln/debug/println.hh> +#include <mln/trace/quiet.hh> + + +using namespace mln; +using value::int_u8; +using value::int_u12; +using value::label_16; + +const float saturation = 1.0; + + + +namespace mln +{ + + template <typename I> + void io_save_edges_int_u12(const I& input, + value::int_u8 bg, + const std::string& filename) + { + mlc_equal(mln_value(I), value::int_u12)::check(); + mln_ch_value(I, value::int_u8) output; + initialize(output, input); + arith::div_cst(input, 16, output); + io::pgm::save(world::inter_pixel::display_edge(output.unmorph_(), + bg, + 3), + filename); + } + +} // end of namespace mln + + + + + + + +int usage(const char* bin) +{ + std::cout << "Usage: " << bin << " input.dump lambda" << std::endl; + return 1; +} + +int main(int argc, char* argv[]) +{ + if (argc != 3) + return usage(argv[0]); + + unsigned lambda = atoi(argv[2]); + + + // Initialization. + typedef int_u12 input_type; + image2d<input_type> input; + io::dump::load(input, argv[1]); + + + // Closing. + mln_VAR(d_clo, morpho::closing::volume(input | world::inter_pixel::is_separator(), world::inter_pixel::e2e(), lambda)); + io_save_edges_int_u12(d_clo, 0, "d_clo.pgm"); + + + // Watershed. + typedef label_16 L; + L nbasins; + mln_VAR(wst, morpho::watershed::flooding(d_clo, world::inter_pixel::e2e(), nbasins)); + + std::cout << "nbasins: " << nbasins << std::endl; + + + mln_VAR(w_all, wst.unmorph_()); + //io::dump::save(w_all, "watershed_edges.dump"); + //data::fill((w | (!world::inter_pixel::is_separator())).rw(), nbasins.next()); + mln_VAR(w_pixels, w_all | world::inter_pixel::is_pixel()); + data::paste(morpho::dilation(extend(w_pixels, pw::value(w_all)), c4().win()), w_all); + // edges -> dots + mln_VAR(w_dots, w_all | world::inter_pixel::dim2::is_dot()); + data::paste(morpho::erosion(extend(w_dots, pw::value(w_all)), c4().win()), w_all); + + //io::ppm::save(labeling::colorize(value::rgb8(), w, nbasins.next()), "result.ppm"); + io::pgm::save(labeling::wrap(int_u8(), w_all), "watershed.pgm"); + + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/Makefile (revision 3895) +++ trunk/milena/sandbox/fabien/igr/Makefile (revision 3896) @@ -57,19 +57,20 @@ label2gif: label2gif.cc ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} $^ -o label2gif -norm: norm.cc - ${CXX} -I../../../ ${CXXFLAGS} $^ -o norm - -##################### - -seg_fixed: seg_fixed.cc - ${CXX} -I../../../ ${CXXFLAGS} $^ -o seg_fixed - ##################### filter: filter.cc ${CXX} -I../../../ ${CXXFLAGS} $^ -o filter +norm: norm.cc + ${CXX} -I../../../ ${CXXFLAGS} $^ -o norm + +dist_edges: dist_edges.cc + ${CXX} -I../../../ ${CXXFLAGS} $^ -o dist_edges + +wst_edges: wst_edges.cc + ${CXX} -I../../../ ${CXXFLAGS} $^ -o wst_edges + ##################### clean: Index: trunk/milena/sandbox/fabien/igr/dist_edges.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/dist_edges.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/dist_edges.cc (revision 3896) @@ -0,0 +1,214 @@ +#include <algorithm> + +#include <mln/core/image/image1d.hh> +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/core/image/image3d.hh> +#include <mln/core/image/slice_image.hh> +#include <mln/core/image/image_if.hh> +#include <mln/core/routine/duplicate.hh> +#include <mln/core/routine/extend.hh> + +#include <mln/io/dump/all.hh> +#include <mln/io/pgm/save.hh> +#include <mln/io/ppm/save.hh> + +#include <mln/value/int_u8.hh> +#include <mln/value/int_u12.hh> + +#include <mln/accu/sum.hh> +#include <mln/arith/div.hh> +#include <mln/convert/from_to.hh> +#include <mln/level/compute.hh> +#include <mln/level/convert.hh> +#include <mln/level/stretch.hh> +#include <mln/make/image2d.hh> +#include <mln/math/diff_abs.hh> +#include <mln/morpho/dilation.hh> +#include <mln/morpho/erosion.hh> +#include <mln/pw/all.hh> +#include <mln/util/array.hh> + +#include <mln/world/inter_pixel/display_edge.hh> +#include <mln/world/inter_pixel/compute.hh> +#include <mln/world/inter_pixel/immerse.hh> +#include <mln/world/inter_pixel/neighb2d.hh> +#include <mln/world/inter_pixel/all.hh> + + +using namespace mln; +using value::int_u8; +using value::int_u12; + + + +namespace mln +{ + + template <typename I> + void io_save_edges_int_u12(const I& input, + value::int_u8 bg, + const std::string& filename) + { + mlc_equal(mln_value(I), value::int_u12)::check(); + mln_ch_value(I, value::int_u8) output; + initialize(output, input); + arith::div_cst(input, 16, output); + io::pgm::save(world::inter_pixel::display_edge(output.unmorph_(), + bg, + 3), + filename); + } + +} // end of namespace mln + + + +// Distance function. +//------------------- + +struct dist_t : Function_vv2v<dist_t> +{ + typedef int_u12 result; + + template <typename V> + int_u12 operator()(const util::array<V>& v1, + const util::array<V>& v2) const + { + float res = 0.f; + + for (unsigned i = 0; i < v1.nelements(); ++i) + res += std::min(v1[i], v2[i]); + + image1d<V> tmp_ima; + image1d<V> tmp_ima2; + accu::sum<V> accu_sum; + + convert::from_to(v1, tmp_ima); + float sum_v1 = level::compute(accu_sum, tmp_ima); + + convert::from_to(v2, tmp_ima2); + float sum_v2 = level::compute(accu_sum, tmp_ima2); + + if (sum_v1 == 0 && sum_v2 == 0) + return 1; + + res /= std::max(sum_v1, sum_v2); + res = 1 - res; + res = res * 4095; + + return (int) res; + } +} dist; + + +struct dist2_t : Function_vv2v<dist2_t> +{ + typedef int_u12 result; + + template <typename V> + int_u12 operator()(const util::array<V>& v1, + const util::array<V>& v2) const + { + float res = 0.f; + + for (unsigned i = 0; i < v1.nelements(); ++i) + { + unsigned maxi = std::max(v1[i], v2[i]); + if (maxi != 0u) + res += float(std::min(v1[i], v2[i])) / float(std::max(v1[i], v2[i])); + else + res += 1; + } + + res = res / v1.nelements(); + res = 1 - res; + res = res * 4095; + + return (int) res; + } +} dist2; + + +struct dist3_t : Function_vv2v<dist3_t> +{ + typedef int_u12 result; + + template <typename V> + int_u12 operator()(const util::array<V>& v1, + const util::array<V>& v2) const + { + float res = 0.f; + float min = 0; + float max = 0; + + for (unsigned i = 0; i < v1.nelements(); ++i) + { + min += std::min(v1[i], v2[i]); + max += std::max(v1[i], v2[i]); + } + + if (max) + res = min / max; + + res = 1 - res; + res = res * 4095; + + return (int) res; + } +} dist3; + + + + + + +int usage(const char* bin) +{ + std::cout << "Usage: " << bin << " input.dump min dist_type output.dump" << std::endl; + return 1; +} + +int main(int argc, char* argv[]) +{ + if (argc != 5) + return usage(argv[0]); + + float min = atof(argv[2]); + unsigned dist_type = atoi(argv[3]); + + + // Initialization. + image3d<float> input; + io::dump::load(input, argv[1]); + typedef image2d<util::array<float> > I; + I ima_arr; + initialize(ima_arr, slice(input, 0)); + for (unsigned int i = 0; i < input.nslices(); ++i) + { + image2d<float> tmp_slice = duplicate(slice(input, i)); + mln_piter_(image2d<float>) p(tmp_slice.domain()); + for_all(p) + ima_arr(p).append(tmp_slice(p) - min); // We set the minimum value to 0. + } + + + // Edges image creation. + typedef image_if<I, world::inter_pixel::is_pixel> Ix; + Ix imax = world::inter_pixel::immerse(ima_arr); + + // Edges distance computation. + image_if<image2d<int_u12>, world::inter_pixel::is_separator> edges; + if (dist_type == 1) + edges = world::inter_pixel::compute(imax, dist); + else if (dist_type == 2) + edges = world::inter_pixel::compute(imax, dist2); + else if (dist_type == 3) + edges = world::inter_pixel::compute(imax, dist3); + + io_save_edges_int_u12(edges, 0, "dist.pgm"); + + io::dump::save(edges.unmorph_(), argv[4]); + + return 0; +}