URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-05-25 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; Change processing chain, now filter before normalization. * fabien/igr/Makefile: Add filter target. * fabien/igr/filter.cc: New tool for filtering input. * fabien/igr/fixed_seg/main.cc: Update closing to closing::sum * fabien/igr/point_filtering/main.cc: Minor update. --- Makefile | 5 + filter.cc | 232 ++++++++++++++++++++++++++++++++++++++++++++++++ fixed_seg/main.cc | 146 +++++++----------------------- point_filtering/main.cc | 29 +++--- 4 files changed, 287 insertions(+), 125 deletions(-) Index: trunk/milena/sandbox/fabien/igr/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/Makefile (revision 3870) +++ trunk/milena/sandbox/fabien/igr/Makefile (revision 3871) @@ -67,6 +67,11 @@ ##################### +filter: filter.cc + ${CXX} -I../../../ ${CXXFLAGS} $^ -o filter + +##################### + clean: rm -rf *.dump *.p?m *.plot *.log *.csv *.dSYM rm seg2d seg3d wsd2d wsd3d nbasins_finder grad clo_vol wst graph med thres matlab time_max first_slice_dicom Index: trunk/milena/sandbox/fabien/igr/point_filtering/main.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/point_filtering/main.cc (revision 3870) +++ trunk/milena/sandbox/fabien/igr/point_filtering/main.cc (revision 3871) @@ -9,6 +9,7 @@ #include <mln/morpho/closing/structural.hh> #include <mln/morpho/opening/structural.hh> #include <mln/util/array.hh> +#include <mln/win/segment1d.hh> using namespace mln; @@ -19,6 +20,7 @@ if (argc != 2) { std::cout << "Usage: " << argv[0] << "point.plot" << std::endl; + return 1; } util::array<float> arr; @@ -30,16 +32,17 @@ io::plot::save(ima, "raw.plot"); // Morpho filtering. - window1d win_1; - win_1 - .insert(-2) - .insert(-1) - .insert(0) - .insert(1) - .insert(2); + win::segment1d seg3(3); + win::segment1d seg5(5); + win::segment1d seg7(7); + win::segment1d seg9(9); + win::segment1d seg11(11); + win::segment1d seg13(13); + win::segment1d seg15(15); + win::segment1d seg21(21); - image1d<float> opening_ima = morpho::opening::structural(ima, win_1); - image1d<float> closing_ima = morpho::closing::structural(ima, win_1); + image1d<float> opening_ima = morpho::opening::structural(ima, seg21); + image1d<float> closing_ima = morpho::closing::structural(ima, seg21); image1d<accu::mean<float> > result; @@ -52,8 +55,8 @@ io::plot::save(ima_morpho, "morpho.plot"); // Morpho (again). - opening_ima = morpho::opening::structural(ima_morpho, win_1); - closing_ima = morpho::closing::structural(ima_morpho, win_1); + opening_ima = morpho::opening::structural(ima_morpho, seg13); + closing_ima = morpho::closing::structural(ima_morpho, seg13); initialize(result, ima_morpho); @@ -64,8 +67,8 @@ io::plot::save(ima_morpho2, "morpho2.plot"); // Morpho (the return of the revenge). - opening_ima = morpho::opening::structural(ima_morpho2, win_1); - closing_ima = morpho::closing::structural(ima_morpho2, win_1); + opening_ima = morpho::opening::structural(ima_morpho2, seg5); + closing_ima = morpho::closing::structural(ima_morpho2, seg5); initialize(result, ima_morpho2); Index: trunk/milena/sandbox/fabien/igr/fixed_seg/main.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/fixed_seg/main.cc (revision 3870) +++ trunk/milena/sandbox/fabien/igr/fixed_seg/main.cc (revision 3871) @@ -19,7 +19,6 @@ #include <mln/value/int_u12.hh> #include <mln/value/label_16.hh> #include <mln/value/rgb8.hh> -#include <mln/value/float01_8.hh> #include <mln/accu/sum.hh> #include <mln/accu/mean.hh> @@ -41,13 +40,14 @@ #include <mln/math/diff_abs.hh> #include <mln/morpho/dilation.hh> #include <mln/morpho/erosion.hh> -#include <mln/morpho/closing/area.hh> +#include <mln/morpho/closing/sum.hh> #include <mln/morpho/closing/structural.hh> #include <mln/morpho/opening/structural.hh> #include <mln/morpho/elementary/gradient.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> @@ -65,7 +65,6 @@ using value::int_u8; using value::int_u12; using value::label_16; -using value::float01_8; const float saturation = 1.0; @@ -74,18 +73,6 @@ namespace mln { - struct int_u12_from_float : Function_v2v< int_u12_from_float > - { - typedef value::int_u12 result; - result operator()(float f) const - { - mln_precondition(f >= 0.f && f <= 1.f); - unsigned i = f / saturation * 4095; - return i > 4095 ? 4095 : i; - } - }; - - template <typename I> void io_save_edges_int_u12(const I& input, value::int_u8 bg, @@ -107,21 +94,16 @@ // Mean image. //------------ + template <typename V> inline image1d<float> -mean_image(image1d<V>& input) +mean_image(image1d<V>& input, unsigned seg_size) { - window1d win_1; - win_1 - .insert(-2) - .insert(-1) - .insert(0) - .insert(1) - .insert(2); + win::segment1d seg(seg_size); - image1d<V> closing_ima = morpho::closing::structural(input, win_1); - image1d<V> opening_ima = morpho::opening::structural(input, win_1); + image1d<V> closing_ima = morpho::closing::structural(input, seg); + image1d<V> opening_ima = morpho::opening::structural(input, seg); image1d<accu::mean<float> > result; @@ -166,19 +148,21 @@ res /= std::max(sum_v1, sum_v2); res = 1 - res; - res = res * 4095; + res = res * (4095 / saturation); + if (res > 4095) + return 4095; return (int) res; } } dist; -/*struct dist_morpho_t : Function_vv2v<dist_morpho_t> +struct dist_morpho_t : Function_vv2v<dist_morpho_t> { - typedef float result; + typedef int_u12 result; template <typename V> - float operator()(util::array<V> v1, util::array<V> v2) const + int_u12 operator()(util::array<V> v1, util::array<V> v2) const { float res = 0.f; @@ -186,12 +170,16 @@ image1d<V> tmp_ima; convert::from_to(v1, tmp_ima); - image1d<float> morpho_ima = mean_image(tmp_ima); + image1d<float> morpho_ima = mean_image(tmp_ima, 15); + morpho_ima = mean_image(morpho_ima, 11); + morpho_ima = mean_image(morpho_ima, 7); float sum_v1 = level::compute(accu_sum, morpho_ima); image1d<V> tmp_ima2; convert::from_to(v2, tmp_ima2); - image1d<float> morpho_ima2 = mean_image(tmp_ima2); + image1d<float> morpho_ima2 = mean_image(tmp_ima2, 15); + morpho_ima2 = mean_image(tmp_ima2, 11); + morpho_ima2 = mean_image(tmp_ima2, 7); float sum_v2 = level::compute(accu_sum, morpho_ima2); mln_piter(image1d<float>) p(morpho_ima.domain()); @@ -202,10 +190,14 @@ return 1; res /= std::max(sum_v1, sum_v2); + res = 1 - res; + res = res * (4095 / saturation); + if (res > 4095) + return 4095; - return 1 - res; + return res; } -} dist_morpho;*/ +} dist_morpho; @@ -252,18 +244,16 @@ // Edges distance computation. image_if<image2d<int_u12>, world::inter_pixel::is_separator> edges; - //if (!is_smooth) + if (!is_smooth) edges = world::inter_pixel::compute(imax, dist); - //else - // edges = world::inter_pixel::compute(imax, dist_morpho); - //io::dump::save(edges.unmorph_(), "edges.dump"); - //mln_VAR(d, level::transform(edges, int_u12_from_float())); + else + edges = world::inter_pixel::compute(imax, dist_morpho); io_save_edges_int_u12(edges, 0, "dist.pgm"); // Closing. - /*mln_VAR(d_clo, morpho::closing::area(d, world::inter_pixel::e2e(), lambda)); - //io_save_edges_int_u12(d_clo, 0, "d_clo.pgm"); + mln_VAR(d_clo, morpho::closing::sum(edges, world::inter_pixel::e2e(), lambda)); + io_save_edges_int_u12(d_clo, 0, "d_clo.pgm"); // Watershed. @@ -271,11 +261,11 @@ L nbasins; mln_VAR(wst, morpho::watershed::flooding(d_clo, world::inter_pixel::e2e(), nbasins)); - std::cout << "nbasins: " << nbasins << std::endl;*/ + std::cout << "nbasins: " << nbasins << std::endl; - /*mln_VAR(w_all, wst.unmorph_()); - io::dump::save(w_all, "watershed_edges.dump"); + 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); @@ -284,7 +274,7 @@ 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");*/ + io::pgm::save(labeling::wrap(int_u8(), w_all), "watershed.pgm"); @@ -326,74 +316,6 @@ - // Deviation. - /*util::array<accu::stat::deviation<float> > arr_dev; - for (unsigned i = 0; i < means.nelements(); ++i) - arr_dev.append(accu::stat::deviation<float> (means[i])); - util::array<float> deviations = labeling::compute(arr_dev, e, wst, nbasins); - - // Display. - { - typedef image_if<image2d<float>, world::inter_pixel::is_separator> Fsx; - Fsx ima_dev; - initialize(ima_dev, wst); - data::paste(wst, ima_dev); - for (unsigned i = 1; i < deviations.nelements(); ++i) - data::fill((ima_dev | pw::value(ima_dev) == pw::cst(i)).rw(), deviations[i]); - mln_VAR(display_dev, world::inter_pixel::display_edge(ima_dev.unmorph_(), 0.0, 3)); - io::pgm::save(level::stretch(int_u8(), display_dev), "05_dev.pgm"); - }*/ - - - - // Plots labels. - /*image2d<L> w_simple = world::inter_pixel::full2image(w_all); - plot_label(input, w_simple, 191u); - plot_label(input, w_simple, 171u); - plot_label(input, w_simple, 188u); - - plot_label(input, w_simple, 16u); - - plot_label(input, w_simple, 187u);*/ - - - - - - - - - - /*mln_VAR(clo, morpho::closing::volume(edges | world::inter_pixel::dim2::is_edge(), world::inter_pixel::e2e(), atoi(argv[2]))); - - // Debug. - //debug::println("clo", clo); - - // Display. - image2d<float> display_ima2 = world::inter_pixel::display_edge(clo.unmorph_(), 0.0, 3); - io::pgm::save(level::stretch(int_u8(), display_ima2), "edges2.pgm"); - - // Watershed. - typedef label_16 L; - L nbasins; - mln_VAR(wst, morpho::watershed::flooding(clo, world::inter_pixel::e2e(), nbasins)); - - // Debug. - //debug::println("wst", wst); - - // Extension. - image2d<L> w_all = wst.unmorph_(); - // edges -> pixel - mln_VAR(w_pixels, w_all | world::inter_pixel::dim2::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); - - - // Save labels map. - std::cout << "nbasins: " << nbasins << std::endl; - io::dump::save(wst.unmorph_(), "watershed_fixed.dump");*/ return 0; } Index: trunk/milena/sandbox/fabien/igr/filter.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/filter.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/filter.cc (revision 3871) @@ -0,0 +1,232 @@ +#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/sum.hh> +#include <mln/morpho/closing/structural.hh> +#include <mln/morpho/opening/structural.hh> +#include <mln/morpho/elementary/gradient.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; + + + + +// Mean image. +//------------ + +template <typename V> +inline +image1d<float> +mean_image(image1d<V>& input, unsigned seg_size) +{ + win::segment1d seg(seg_size); + + image1d<V> closing_ima = morpho::closing::structural(input, seg); + image1d<V> opening_ima = morpho::opening::structural(input, seg); + + image1d<accu::mean<float> > result; + + initialize(result, input); + + accu::image::init(result); + accu::image::take(result, closing_ima); + accu::image::take(result, opening_ima); + + return accu::image::to_result(result); +} + + + +// Distance function. +//------------------- + +struct dist_t : Function_vv2v<dist_t> +{ + typedef int_u12 result; + + template <typename V> + int_u12 operator()(util::array<V> v1, 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; + if (res > 4095) + return 4095; + + return (int) res; + } +} dist; + + +struct dist_morpho_t : Function_vv2v<dist_morpho_t> +{ + typedef int_u12 result; + + template <typename V> + int_u12 operator()(util::array<V> v1, util::array<V> v2) const + { + float res = 0.f; + + accu::sum<V> accu_sum; + + image1d<V> tmp_ima; + convert::from_to(v1, tmp_ima); + image1d<float> morpho_ima = mean_image(tmp_ima, 15); + morpho_ima = mean_image(morpho_ima, 11); + morpho_ima = mean_image(morpho_ima, 7); + float sum_v1 = level::compute(accu_sum, morpho_ima); + + image1d<V> tmp_ima2; + convert::from_to(v2, tmp_ima2); + image1d<float> morpho_ima2 = mean_image(tmp_ima2, 15); + morpho_ima2 = mean_image(tmp_ima2, 11); + morpho_ima2 = mean_image(tmp_ima2, 7); + float sum_v2 = level::compute(accu_sum, morpho_ima2); + + mln_piter(image1d<float>) p(morpho_ima.domain()); + for_all(p) + res += std::min(morpho_ima(p), morpho_ima2(p)); + + if (sum_v1 == 0 && sum_v2 == 0) + return 1; + + res /= std::max(sum_v1, sum_v2); + res = 1 - res; + res = res * 4095; + if (res > 4095) + return 4095; + + return res; + } +} dist_morpho; + + + + + + +int usage(const char* bin) +{ + std::cout << "Usage: " << bin << " input.dump output.dump" << std::endl; + return 1; +} + +int main(int argc, char* argv[]) +{ + if (argc != 3) + return usage(argv[0]); + + // Initialization. + typedef int_u12 input_type; + image3d<input_type> input; + io::dump::load(input, argv[1]); + int min_slice = input.bbox().pmin().sli(); + int max_slice = input.bbox().pmax().sli(); + typedef image2d<util::array<input_type> > I; + I ima_arr; + initialize(ima_arr, slice(input, 0)); + for (int i = min_slice; i <= max_slice; ++i) + { + image2d<input_type> tmp_slice = duplicate(slice(input, i)); + mln_piter_(image2d<input_type>) p(tmp_slice.domain()); + for_all(p) + ima_arr(p).append(tmp_slice(p)); + } + + + image3d<float> output; + initialize(output, input); + + mln_piter_(I) p(ima_arr.domain()); + for_all(p) + { + image1d<input_type> tmp_ima; + convert::from_to(ima_arr(p), tmp_ima); + image1d<float> morpho_ima = mean_image(tmp_ima, 15); + morpho_ima = mean_image(morpho_ima, 11); + morpho_ima = mean_image(morpho_ima, 7); + + for (int i = min_slice; i <= max_slice; ++i) + output.at_(i, p.row(), p.col()) = morpho_ima.at_(i - min_slice); + } + + + io::dump::save(output, argv[2]); + + return 0; +}