[Olena-patches] r3594: Add time smoothing in IGR and drafts of space smoothing

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-03-30 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; Add time smoothing in IGR and drafts of space smoothing. * fabien/igr/Makefile.rules: New file to centralize rules. * fabien/igr/img/slice_7.pgm: Remove. * fabien/igr/img: Remove. * fabien/igr/space_smooth/Makefile: New. * fabien/igr/space_smooth/linear.cc: Implements linear space smoothing. * fabien/igr/space_smooth/median.cc: Draft. * fabien/igr/space_smooth/morpho.cc: Draft. * fabien/igr/time_smooth/Makefile: New. * fabien/igr/time_smooth/linear.cc: Implements linear time smoothing. * fabien/igr/time_smooth/median.cc: Implements median time smoothing. * fabien/igr/time_smooth/morpho.cc: Implements morpho time smoothing. --- Makefile.rules | 9 ++++ space_smooth/Makefile | 15 +++++++ space_smooth/linear.cc | 75 ++++++++++++++++++++++++++++++++++++ space_smooth/median.cc | 102 +++++++++++++++++++++++++++++++++++++++++++++++++ space_smooth/morpho.cc | 102 +++++++++++++++++++++++++++++++++++++++++++++++++ time_smooth/Makefile | 15 +++++++ time_smooth/linear.cc | 77 ++++++++++++++++++++++++++++++++++++ time_smooth/median.cc | 102 +++++++++++++++++++++++++++++++++++++++++++++++++ time_smooth/morpho.cc | 102 +++++++++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 599 insertions(+) Index: trunk/milena/sandbox/fabien/igr/space_smooth/median.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/space_smooth/median.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/space_smooth/median.cc (revision 3594) @@ -0,0 +1,102 @@ +#include <mln/core/image/image2d.hh> +#include <mln/core/image/image3d.hh> +#include <mln/core/image/image_if.hh> +#include <mln/core/image/cast_image.hh> +#include <mln/core/image/slice_image.hh> +#include <mln/core/routine/duplicate.hh> + +#include <mln/value/int_u12.hh> + +#include <mln/io/dicom/load.hh> +#include <mln/io/plot/save.hh> + +#include <mln/accu/median_h.hh> +#include <mln/util/array.hh> + + +using namespace mln; +using value::int_u12; + + +int main(int argc, char* argv[]) +{ + if (argc != 2) + { + std::cout << "Usage: " << argv[0] << "input" << std::endl; + return 1; + } + + /////////////////// + // // + // Image loading // + // // + /////////////////// + image3d<int_u12> input; + io::dicom::load(input, argv[1]); + image2d<util::array<int_u12> > ima_arr(input.nrows(), input.ncols()); + 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)); + } + } + + //////////// + // // + // 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; + 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()); + } + ima_median3(p).append(ima_arr(p)[ima_arr(p).nelements() - 1]); + + // 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()); + } + ima_median5(p).append(ima_arr(p)[ima_arr(p).nelements() - 2]); + ima_median5(p).append(ima_arr(p)[ima_arr(p).nelements() - 1]); + } + + ///////////// + // // + // 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"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/space_smooth/morpho.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/space_smooth/morpho.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/space_smooth/morpho.cc (revision 3594) @@ -0,0 +1,102 @@ +#include <mln/core/image/image1d.hh> +#include <mln/core/alias/neighb1d.hh> +#include <mln/core/image/image2d.hh> +#include <mln/core/image/image3d.hh> +#include <mln/core/image/image_if.hh> +#include <mln/core/image/cast_image.hh> +#include <mln/core/image/slice_image.hh> +#include <mln/core/routine/duplicate.hh> + +#include <mln/value/int_u12.hh> + +#include <mln/io/dicom/load.hh> +#include <mln/io/plot/save.hh> + +#include <mln/accu/median_h.hh> +#include <mln/convert/from_to.hh> +#include <mln/morpho/closing/area.hh> +#include <mln/morpho/opening/area.hh> +#include <mln/util/array.hh> + +/*#include <mln/draw/line.hh> +#include <mln/io/magick/save.hh> +#include <mln/level/convert.hh> +#include <mln/level/stretch.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> +#include <mln/literal/all.hh> +using mln::value::rgb8; +using mln::value::int_u8;*/ + + +using namespace mln; +using value::int_u12; + + +int main(int argc, char* argv[]) +{ + if (argc != 2) + { + std::cout << "Usage: " << argv[0] << "input" << std::endl; + return 1; + } + + /////////////////// + // // + // Image loading // + // // + /////////////////// + image3d<int_u12> input; + io::dicom::load(input, argv[1]); + image2d<util::array<int_u12> > ima_arr(input.nrows(), input.ncols()); + 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)); + } + } + + //////////// + // // + // Morpho // + // // + //////////// + image2d<image1d<int_u12> > ima_morpho; + initialize(ima_morpho, ima_arr); + mln_piter_(image2d<int_u12>) p(ima_morpho.domain()); + accu::median_h<int_u12> accu_med; + for_all(p) + { + image1d<int_u12> tmp_ima; + convert::from_to(ima_arr(p), tmp_ima); + tmp_ima = morpho::closing::area(tmp_ima, c2(), 3); + tmp_ima = morpho::opening::area(tmp_ima, c2(), 3); + ima_morpho(p) = tmp_ima; + } + + ///////////// + // // + // Outputs // + // // + ///////////// + /*image2d<rgb8> ima_color = level::convert(rgb8(), level::stretch(int_u8(), slice(input, 0))); + algebra::vec<2, unsigned int> vmin; + algebra::vec<2, unsigned int> vmax; + vmin[0] = 160; + vmin[1] = 120; + vmax[0] = 122; + vmax[1] = 115; + mln_site_(image2d<bool>) pbeg(vmin); + mln_site_(image2d<bool>) pend(vmax); + draw::line(ima_color, pbeg, pend, literal::red); + io::magick::save(ima_color, "test.png");*/ + + io::plot::save(ima_morpho(point2d(160, 120)), "morpho_tumeur.plot"); + io::plot::save(ima_morpho(point2d(34, 94)), "morpho_air.plot"); + io::plot::save(ima_morpho(point2d(122, 115)), "morpho_poumon.plot"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/space_smooth/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/space_smooth/Makefile (revision 0) +++ trunk/milena/sandbox/fabien/igr/space_smooth/Makefile (revision 3594) @@ -0,0 +1,15 @@ +include ../Makefile.rules + + +linear: linear.cc + ${CXX} -I../../../../ ${DICOM} ${CXXFLAGS} $^ -o linear + +median: median.cc + ${CXX} -I../../../../ ${DICOM} ${CXXFLAGS} $^ -o median + +morpho: morpho.cc + ${CXX} -I../../../../ ${DICOM} ${MAGICK} ${CXXFLAGS} $^ -o morpho + +clean: + rm -rf *.dump *.p?m *.plot *.log *.csv *.dSYM + rm linear median morpho Index: trunk/milena/sandbox/fabien/igr/space_smooth/linear.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/space_smooth/linear.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/space_smooth/linear.cc (revision 3594) @@ -0,0 +1,75 @@ +#include <mln/core/image/image2d.hh> +#include <mln/core/image/image3d.hh> +#include <mln/core/image/image_if.hh> +#include <mln/core/image/cast_image.hh> +#include <mln/core/image/slice_image.hh> +#include <mln/core/routine/duplicate.hh> + +#include <mln/value/int_u12.hh> + +#include <mln/io/dicom/load.hh> +#include <mln/io/plot/save.hh> + +#include <mln/accu/median_h.hh> +#include <mln/linear/convolve.hh> +#include <mln/make/w_window2d.hh> +#include <mln/util/array.hh> + + +using namespace mln; +using value::int_u12; + + + +int main(int argc, char* argv[]) +{ + if (argc != 2) + { + std::cout << "Usage: " << argv[0] << "input" << std::endl; + return 1; + } + + /////////////////// + // // + // Image loading // + // // + /////////////////// + image3d<int_u12> input; + io::dicom::load(input, argv[1]); + util::array<image2d<int_u12> > arr_ima; + for (int i = 0; i < input.nslices(); ++i) + { + image2d<int_u12> tmp_slice = duplicate(slice(input, i)); + arr_ima.append(tmp_slice); + } + + //////////////////////// + // // + // Linear convolution // + // // + //////////////////////// + float ws[] = { 0, 1/8, 0, + 1/8, 1/2, 1/8, + 0, 1/8, 0 }; + + util::array<image2d<float> > ima_linear; + for (unsigned i = 1; i < arr_ima.nelements(); ++i) + ima_linear.append(linear::convolve(arr_ima[i], make::w_window2d(ws))); + + ///////////// + // // + // Outputs // + // // + ///////////// + image2d<util::array<float> > ima_result(input.nrows(), input.ncols()); + mln_piter_(image2d<util::array<float> >) p(ima_linear[0].domain()); + for_all(p) + for (int i = 0; i < ima_linear.nelements(); ++i) + ima_result(p).append(ima_linear[i](p)); + + io::plot::save(ima_result(point2d(160, 120)), "linear_tumeur.plot"); + io::plot::save(ima_result(point2d(34, 94)), "linear_air.plot"); + io::plot::save(ima_result(point2d(122, 115)), "linear_poumon.plot"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/time_smooth/median.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/time_smooth/median.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/time_smooth/median.cc (revision 3594) @@ -0,0 +1,102 @@ +#include <mln/core/image/image2d.hh> +#include <mln/core/image/image3d.hh> +#include <mln/core/image/image_if.hh> +#include <mln/core/image/cast_image.hh> +#include <mln/core/image/slice_image.hh> +#include <mln/core/routine/duplicate.hh> + +#include <mln/value/int_u12.hh> + +#include <mln/io/dicom/load.hh> +#include <mln/io/plot/save.hh> + +#include <mln/accu/median_h.hh> +#include <mln/util/array.hh> + + +using namespace mln; +using value::int_u12; + + +int main(int argc, char* argv[]) +{ + if (argc != 2) + { + std::cout << "Usage: " << argv[0] << "input" << std::endl; + return 1; + } + + /////////////////// + // // + // Image loading // + // // + /////////////////// + image3d<int_u12> input; + io::dicom::load(input, argv[1]); + image2d<util::array<int_u12> > ima_arr(input.nrows(), input.ncols()); + 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)); + } + } + + //////////// + // // + // 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; + 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()); + } + ima_median3(p).append(ima_arr(p)[ima_arr(p).nelements() - 1]); + + // 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()); + } + ima_median5(p).append(ima_arr(p)[ima_arr(p).nelements() - 2]); + ima_median5(p).append(ima_arr(p)[ima_arr(p).nelements() - 1]); + } + + ///////////// + // // + // 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"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/time_smooth/morpho.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/time_smooth/morpho.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/time_smooth/morpho.cc (revision 3594) @@ -0,0 +1,102 @@ +#include <mln/core/image/image1d.hh> +#include <mln/core/alias/neighb1d.hh> +#include <mln/core/image/image2d.hh> +#include <mln/core/image/image3d.hh> +#include <mln/core/image/image_if.hh> +#include <mln/core/image/cast_image.hh> +#include <mln/core/image/slice_image.hh> +#include <mln/core/routine/duplicate.hh> + +#include <mln/value/int_u12.hh> + +#include <mln/io/dicom/load.hh> +#include <mln/io/plot/save.hh> + +#include <mln/accu/median_h.hh> +#include <mln/convert/from_to.hh> +#include <mln/morpho/closing/area.hh> +#include <mln/morpho/opening/area.hh> +#include <mln/util/array.hh> + +/*#include <mln/draw/line.hh> +#include <mln/io/magick/save.hh> +#include <mln/level/convert.hh> +#include <mln/level/stretch.hh> +#include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> +#include <mln/literal/all.hh> +using mln::value::rgb8; +using mln::value::int_u8;*/ + + +using namespace mln; +using value::int_u12; + + +int main(int argc, char* argv[]) +{ + if (argc != 2) + { + std::cout << "Usage: " << argv[0] << "input" << std::endl; + return 1; + } + + /////////////////// + // // + // Image loading // + // // + /////////////////// + image3d<int_u12> input; + io::dicom::load(input, argv[1]); + image2d<util::array<int_u12> > ima_arr(input.nrows(), input.ncols()); + 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)); + } + } + + //////////// + // // + // Morpho // + // // + //////////// + image2d<image1d<int_u12> > ima_morpho; + initialize(ima_morpho, ima_arr); + mln_piter_(image2d<int_u12>) p(ima_morpho.domain()); + accu::median_h<int_u12> accu_med; + for_all(p) + { + image1d<int_u12> tmp_ima; + convert::from_to(ima_arr(p), tmp_ima); + tmp_ima = morpho::closing::area(tmp_ima, c2(), 3); + tmp_ima = morpho::opening::area(tmp_ima, c2(), 3); + ima_morpho(p) = tmp_ima; + } + + ///////////// + // // + // Outputs // + // // + ///////////// + /*image2d<rgb8> ima_color = level::convert(rgb8(), level::stretch(int_u8(), slice(input, 0))); + algebra::vec<2, unsigned int> vmin; + algebra::vec<2, unsigned int> vmax; + vmin[0] = 160; + vmin[1] = 120; + vmax[0] = 122; + vmax[1] = 115; + mln_site_(image2d<bool>) pbeg(vmin); + mln_site_(image2d<bool>) pend(vmax); + draw::line(ima_color, pbeg, pend, literal::red); + io::magick::save(ima_color, "test.png");*/ + + io::plot::save(ima_morpho(point2d(160, 120)), "morpho_tumeur.plot"); + io::plot::save(ima_morpho(point2d(34, 94)), "morpho_air.plot"); + io::plot::save(ima_morpho(point2d(122, 115)), "morpho_poumon.plot"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/time_smooth/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/time_smooth/Makefile (revision 0) +++ trunk/milena/sandbox/fabien/igr/time_smooth/Makefile (revision 3594) @@ -0,0 +1,15 @@ +include ../Makefile.rules + + +linear: linear.cc + ${CXX} -I../../../../ ${DICOM} ${CXXFLAGS} $^ -o linear + +median: median.cc + ${CXX} -I../../../../ ${DICOM} ${CXXFLAGS} $^ -o median + +morpho: morpho.cc + ${CXX} -I../../../../ ${DICOM} ${MAGICK} ${CXXFLAGS} $^ -o morpho + +clean: + rm -rf *.dump *.p?m *.plot *.log *.csv *.dSYM + rm linear median morpho Index: trunk/milena/sandbox/fabien/igr/time_smooth/linear.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/time_smooth/linear.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/time_smooth/linear.cc (revision 3594) @@ -0,0 +1,77 @@ +#include <mln/core/image/image2d.hh> +#include <mln/core/image/image3d.hh> +#include <mln/core/image/image_if.hh> +#include <mln/core/image/cast_image.hh> +#include <mln/core/image/slice_image.hh> +#include <mln/core/routine/duplicate.hh> + +#include <mln/value/int_u12.hh> + +#include <mln/io/dicom/load.hh> +#include <mln/io/plot/save.hh> + +#include <mln/accu/median_h.hh> +#include <mln/util/array.hh> + + +using namespace mln; +using value::int_u12; + + +int main(int argc, char* argv[]) +{ + if (argc != 2) + { + std::cout << "Usage: " << argv[0] << "input" << std::endl; + return 1; + } + + /////////////////// + // // + // Image loading // + // // + /////////////////// + image3d<int_u12> input; + io::dicom::load(input, argv[1]); + image2d<util::array<int_u12> > ima_arr(input.nrows(), input.ncols()); + 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)); + } + } + + //////////////////////// + // // + // Linear convolution // + // // + //////////////////////// + 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; + for_all(p) + { + ima_linear(p).append(ima_arr(p)[0]); + for (unsigned i = 1; i < ima_arr(p).nelements() - 1; ++i) + ima_linear(p).append(0.25 * ima_arr(p)[i - 1] + 0.5 * ima_arr(p)[i] + 0.25 * ima_arr(p)[i + 1]); + ima_linear(p).append(ima_arr(p)[ima_arr(p).nelements() - 1]); + } + + ///////////// + // // + // Outputs // + // // + ///////////// + io::plot::save(ima_arr(point2d(160, 120)), "ref_tumeur.plot"); + io::plot::save(ima_linear(point2d(160, 120)), "linear_tumeur.plot"); + io::plot::save(ima_arr(point2d(34, 94)), "ref_air.plot"); + io::plot::save(ima_linear(point2d(34, 94)), "linear_air.plot"); + io::plot::save(ima_arr(point2d(122, 115)), "ref_poumon.plot"); + io::plot::save(ima_linear(point2d(122, 115)), "linear_poumon.plot"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/Makefile.rules =================================================================== --- trunk/milena/sandbox/fabien/igr/Makefile.rules (revision 0) +++ trunk/milena/sandbox/fabien/igr/Makefile.rules (revision 3594) @@ -0,0 +1,9 @@ +CXX = llvm-g++ +CXXFLAGS = -DNDEBUG -O4 + +DICOM_INC = -I/usr/local/include/gdcm-2.0 +DICOM_LIBS = -lgdcmCommon -lgdcmDICT -lgdcmDSED -lgdcmIOD -lgdcmMSFF -lgdcmexpat -lgdcmjpeg12 -lgdcmjpeg16 -lgdcmjpeg8 -lgdcmopenjpeg -lgdcmuuid -lgdcmzlib \ + -framework CoreFoundation +DICOM = ${DICOM_INC} ${DICOM_LIBS} + +MAGICK = `Magick++-config --cppflags --ldflags --libs`