[Olena-patches] r3597: Update spatial smoothing

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-03-31 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; Update spatial smoothing. * fabien/igr/Makefile: Update. * fabien/igr/space_smooth/linear.cc: Fix spatial smoothing. * fabien/igr/space_smooth/median.cc: Implement spatial smoothing. * fabien/igr/space_smooth/morpho.cc: Draft. * fabien/igr/time_smooth/linear.cc: Small fix. --- Makefile | 5 +- space_smooth/linear.cc | 8 +-- space_smooth/median.cc | 111 +++++++++++++++++++++++++++++++------------------ space_smooth/morpho.cc | 12 +---- time_smooth/linear.cc | 2 5 files changed, 83 insertions(+), 55 deletions(-) Index: trunk/milena/sandbox/fabien/igr/space_smooth/median.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/space_smooth/median.cc (revision 3596) +++ trunk/milena/sandbox/fabien/igr/space_smooth/median.cc (revision 3597) @@ -11,6 +11,7 @@ #include <mln/io/plot/save.hh> #include <mln/accu/median_h.hh> +#include <mln/core/alias/window2d.hh> #include <mln/util/array.hh> @@ -33,15 +34,11 @@ /////////////////// image3d<int_u12> input; io::dicom::load(input, argv[1]); - image2d<util::array<int_u12> > ima_arr(input.nrows(), input.ncols()); + util::array<image2d<int_u12> > arr_ima; for (int i = 0; i < input.nslices(); ++i) { image2d<int_u12> tmp_slice = duplicate(slice(input, i)); - mln_piter_(image2d<int_u12>) p(tmp_slice.domain()); - for_all(p) - { - ima_arr(p).append(tmp_slice(p)); - } + arr_ima.append(tmp_slice); } //////////// @@ -49,41 +46,52 @@ // Median // // // //////////// - image2d<util::array<int_u12> > ima_median3; - image2d<util::array<int_u12> > ima_median5; - initialize(ima_median3, ima_arr); - initialize(ima_median5, ima_arr); - mln_piter_(image2d<int_u12>) p(ima_median3.domain()); - accu::median_h<int_u12> accu_med; + window2d win_c12p; + win_c12p + .insert(-1, -1) + .insert( 0, -1) + .insert(+1, -1) + .insert(-1, 0) + .insert( 0, 0) + .insert(+1, 0) + .insert(-1, +1) + .insert( 0, +1) + .insert(+1, +1) + .insert(+2, +2) + .insert(-2, -2) + .insert(+2, -2) + .insert(-2, +2); + window2d win_c4 = win_c4p(); + util::array<image2d<int_u12> > ima_median3; + util::array<image2d<int_u12> > ima_median5; + accu::median_h<int_u12> accu_med3; + accu::median_h<int_u12> accu_med5; + image2d<int_u12> tmp; + initialize(tmp, arr_ima[0]); + mln_piter_(image2d<int_u12>) p(arr_ima[0].domain()); + mln_qiter_(window2d) q3(win_c4, p); + mln_qiter_(window2d) q5(win_c12p, p); + for (unsigned i = 0; i < arr_ima.nelements(); ++i) + { for_all(p) { // Median 3 - ima_median3(p).append(ima_arr(p)[0]); - for (unsigned i = 1; i < ima_arr(p).nelements() - 1; ++i) - { - accu_med.init(); - accu_med.take(ima_arr(p)[i - 1]); - accu_med.take(ima_arr(p)[i]); - accu_med.take(ima_arr(p)[i + 1]); - ima_median3(p).append(accu_med.to_result()); + accu_med3.init(); + for_all(q3) if (arr_ima[i].domain().has(q3)) + accu_med3.take(arr_ima[i](q3)); + tmp(p) = accu_med3.to_result(); } - ima_median3(p).append(ima_arr(p)[ima_arr(p).nelements() - 1]); + ima_median3.append(duplicate(tmp)); + for_all(p) + { // Median 5 - ima_median5(p).append(ima_arr(p)[0]); - ima_median5(p).append(ima_arr(p)[1]); - for (unsigned i = 2; i < ima_arr(p).nelements() - 2; ++i) - { - accu_med.init(); - accu_med.take(ima_arr(p)[i - 2]); - accu_med.take(ima_arr(p)[i - 1]); - accu_med.take(ima_arr(p)[i]); - accu_med.take(ima_arr(p)[i + 1]); - accu_med.take(ima_arr(p)[i + 2]); - ima_median5(p).append(accu_med.to_result()); + accu_med5.init(); + for_all(q5) if (arr_ima[i].domain().has(q5)) + accu_med5.take(arr_ima[i](q5)); + tmp(p) = accu_med5.to_result(); } - ima_median5(p).append(ima_arr(p)[ima_arr(p).nelements() - 2]); - ima_median5(p).append(ima_arr(p)[ima_arr(p).nelements() - 1]); + ima_median5.append(duplicate(tmp)); } ///////////// @@ -91,12 +99,35 @@ // Outputs // // // ///////////// - io::plot::save(ima_median3(point2d(160, 120)), "median3_tumeur.plot"); - io::plot::save(ima_median3(point2d(34, 94)), "median3_air.plot"); - io::plot::save(ima_median3(point2d(122, 115)), "median3_poumon.plot"); - io::plot::save(ima_median5(point2d(160, 120)), "median5_tumeur.plot"); - io::plot::save(ima_median5(point2d(34, 94)), "median5_air.plot"); - io::plot::save(ima_median5(point2d(122, 115)), "median5_poumon.plot"); + util::array<int_u12> arr_tumeur3; + util::array<int_u12> arr_tumeur5; + for (unsigned i = 0; i < ima_median3.nelements(); ++i) + { + arr_tumeur3.append(ima_median3[i](point2d(160, 120))); + arr_tumeur5.append(ima_median5[i](point2d(160, 120))); + } + io::plot::save(arr_tumeur3, "median3_tumeur.plot"); + io::plot::save(arr_tumeur5, "median5_tumeur.plot"); + + util::array<int_u12> arr_air3; + util::array<int_u12> arr_air5; + for (unsigned i = 0; i < ima_median3.nelements(); ++i) + { + arr_air3.append(ima_median3[i](point2d(34, 94))); + arr_air5.append(ima_median5[i](point2d(34, 94))); + } + io::plot::save(arr_air3, "median3_air.plot"); + io::plot::save(arr_air5, "median5_air.plot"); + + util::array<int_u12> arr_poumon3; + util::array<int_u12> arr_poumon5; + for (unsigned i = 0; i < ima_median3.nelements(); ++i) + { + arr_poumon3.append(ima_median3[i](point2d(122, 115))); + arr_poumon5.append(ima_median5[i](point2d(122, 115))); + } + io::plot::save(arr_poumon3, "median3_poumon.plot"); + io::plot::save(arr_poumon5, "median5_poumon.plot"); return 0; } Index: trunk/milena/sandbox/fabien/igr/space_smooth/morpho.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/space_smooth/morpho.cc (revision 3596) +++ trunk/milena/sandbox/fabien/igr/space_smooth/morpho.cc (revision 3597) @@ -48,15 +48,11 @@ /////////////////// image3d<int_u12> input; io::dicom::load(input, argv[1]); - image2d<util::array<int_u12> > ima_arr(input.nrows(), input.ncols()); + util::array<image2d<int_u12> > arr_ima; for (int i = 0; i < input.nslices(); ++i) { image2d<int_u12> tmp_slice = duplicate(slice(input, i)); - mln_piter_(image2d<int_u12>) p(tmp_slice.domain()); - for_all(p) - { - ima_arr(p).append(tmp_slice(p)); - } + arr_ima.append(tmp_slice); } //////////// @@ -64,10 +60,8 @@ // Morpho // // // //////////// - image2d<image1d<int_u12> > ima_morpho; - initialize(ima_morpho, ima_arr); + util::array<image2d<int_u12> > ima_morpho; mln_piter_(image2d<int_u12>) p(ima_morpho.domain()); - accu::median_h<int_u12> accu_med; for_all(p) { image1d<int_u12> tmp_ima; Index: trunk/milena/sandbox/fabien/igr/space_smooth/linear.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/space_smooth/linear.cc (revision 3596) +++ trunk/milena/sandbox/fabien/igr/space_smooth/linear.cc (revision 3597) @@ -48,12 +48,12 @@ // Linear convolution // // // //////////////////////// - float ws[] = { 0, 1/8, 0, - 1/8, 1/2, 1/8, - 0, 1/8, 0 }; + float ws[] = { 0, 1/8.f, 0, + 1/8.f, 1/2.f, 1/8.f, + 0, 1/8.f, 0 }; util::array<image2d<float> > ima_linear; - for (unsigned i = 1; i < arr_ima.nelements(); ++i) + for (unsigned i = 0; i < arr_ima.nelements(); ++i) ima_linear.append(linear::convolve(arr_ima[i], make::w_window2d(ws))); ///////////// Index: trunk/milena/sandbox/fabien/igr/time_smooth/linear.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/time_smooth/linear.cc (revision 3596) +++ trunk/milena/sandbox/fabien/igr/time_smooth/linear.cc (revision 3597) @@ -52,7 +52,7 @@ image2d<util::array<int_u12> > ima_linear; initialize(ima_linear, ima_arr); mln_piter_(image2d<int_u12>) p(ima_linear.domain()); - accu::median<int_u12> accu_med; + accu::median_h<int_u12> accu_med; for_all(p) { ima_linear(p).append(ima_arr(p)[0]); Index: trunk/milena/sandbox/fabien/igr/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/Makefile (revision 3596) +++ trunk/milena/sandbox/fabien/igr/Makefile (revision 3597) @@ -47,6 +47,9 @@ time_max: time_max.cc ${CXX} -I../../../ ${DICOM} ${MAGICK} ${CXXFLAGS} -lm $^ -o time_max +first_slice_dicom: first_slice_dicom.cc + ${CXX} -I../../../ ${DICOM} ${MAGICK} ${CXXFLAGS} -lm $^ -o first_slice_dicom + 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 + rm seg2d seg3d wsd2d wsd3d nbasins_finder grad clo_vol wst graph med thres matlab time_max first_slice_dicom