[Olena-patches] r3695: Various update for IGR demo and new fixed segmentation

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-04-21 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; Various update for IGR demo and new fixed segmentation. * fabien/igr/Makefile: Add 'seg_fixed' target. * fabien/igr/check.sh: Update. * fabien/igr/fun_labels.cc: Update. * fabien/igr/fun_labels.sh: Update. * fabien/igr/med.cc: Update. * fabien/igr/norm.cc: Update. * fabien/igr/seg_fixed.cc: New. * fabien/igr/space_smooth/Makefile: Update. * fabien/igr/space_smooth/linear.cc: Update. * fabien/igr/thres.cc: Update. * fabien/igr/time_max_norm.cc: New tool for tmax map. * fabien/igr/time_smooth/Makefile: Update. * fabien/igr/time_smooth/linear.cc: Update. --- TODO | 18 +++++- igr/Makefile | 12 +++- igr/check.sh | 11 ++- igr/fun_labels.cc | 36 ++++++------ igr/fun_labels.sh | 36 +++++++++--- igr/med.cc | 15 ++--- igr/norm.cc | 14 +--- igr/seg_fixed.cc | 135 +++++++++++++++++++++++++++++++++++++++++++++ igr/space_smooth/Makefile | 2 igr/space_smooth/linear.cc | 21 ++++--- igr/thres.cc | 2 igr/time_max_norm.cc | 116 ++++++++++++++++++++++++++++++++++++++ igr/time_smooth/Makefile | 2 igr/time_smooth/linear.cc | 42 ++++++++++---- 14 files changed, 388 insertions(+), 74 deletions(-) Index: trunk/milena/sandbox/fabien/igr/seg_fixed.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/seg_fixed.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/seg_fixed.cc (revision 3695) @@ -0,0 +1,135 @@ +#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/routine/duplicate.hh> + +#include <mln/io/dump/all.hh> + +#include <mln/value/int_u12.hh> + +#include <mln/accu/sum.hh> +#include <mln/data/fill.hh> +#include <mln/convert/from_to.hh> +#include <mln/level/compute.hh> +#include <mln/util/array.hh> + + +using namespace mln; +using value::int_u12; + + + +inline +point2d +get_edge_location(point2d p, point2d n) +{ + point2d res(p); + + res.to_vec()[0] *= 2; + res.to_vec()[1] *= 2; + + if (n.to_vec()[0] > p.to_vec()[0]) + ++(res.to_vec()[0]); + if (n.to_vec()[0] < p.to_vec()[0]) + --(res.to_vec()[0]); + + if (n.to_vec()[1] > p.to_vec()[1]) + ++(res.to_vec()[1]); + if (n.to_vec()[1] < p.to_vec()[1]) + --(res.to_vec()[1]); + + return res; +} + + +template <typename I> +inline +float +compute_dist(image2d<util::array<I> > ima_arr, + image2d<float> ima_sum, + point2d p, point2d n) +{ + float res = 0.f; + + for (unsigned i = 0; i < ima_arr(p).nelements(); ++i) + res += std::min(ima_arr(p)[i], ima_arr(n)[i]); + res /= std::max(ima_sum(p), ima_sum(n)); + + return res; +} + +template <typename I> +inline +void +compute_sum_arrays(image2d<float> ima_sum, image2d<util::array<I> > ima_arr) +{ + mln_piter(image2d<float>) p(ima_sum.domain()); + for_all(p) + { + image1d<I> tmp_ima; + convert::from_to(ima_arr(p), tmp_ima); + accu::sum<I> accu_sum; + ima_sum(p) = level::compute(accu_sum, tmp_ima); + } +} + +template <typename I> +image2d<float> +dist_on_edges(image2d<util::array<I> > ima_arr) +{ + image2d<float> edges(2 * ima_arr.nrows(), 2 * ima_arr.ncols()); + data::fill(edges, -1.f); + image2d<float> ima_sum; + initialize(ima_sum, ima_arr); + compute_sum_arrays(ima_sum, ima_arr); + + mln_piter(image2d<util::array<I> >) p(ima_arr.domain()); + mln_niter(neighb2d) n(c4(), p); + for_all(p) + for_all(n) + { + point2d location = get_edge_location(p, n); + if (edges(location) != -1) + edges(location) = compute_dist(ima_arr, ima_sum, p, n); + } + + return edges; +} + + + + + +int usage(const char* bin) +{ + std::cout << "Usage: " << bin << " input.dump" << std::endl; + return 1; +} + +int main(int argc, char* argv[]) +{ + typedef int_u12 I; + + if (argc != 2) + return usage(argv[0]); + + image3d<I> input; + io::dump::save(input, argv[1]); + image2d<util::array<I> > ima_arr; + initialize(ima_arr, slice(input, 0)); + for (unsigned int i = 0; i < input.nslices(); ++i) + { + image2d<I> tmp_slice = duplicate(slice(input, i)); + mln_piter_(image2d<I>) p(tmp_slice.domain()); + for_all(p) + ima_arr(p).append(tmp_slice(p)); + } + + image2d<float> edges = dist_on_edges(ima_arr); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/time_max_norm.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/time_max_norm.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/time_max_norm.cc (revision 3695) @@ -0,0 +1,116 @@ +#include <sstream> + +#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_u8.hh> +#include <mln/value/int_u12.hh> + +#include <mln/io/dump/load.hh> +#include <mln/io/pgm/save.hh> + +#include <mln/accu/count.hh> +#include <mln/accu/mean.hh> +#include <mln/accu/sum.hh> +#include <mln/accu/image/all.hh> + +#include <mln/arith/minus.hh> +#include <mln/arith/all.hh> + +#include <mln/data/fill.hh> +#include <mln/data/paste.hh> + +#include <mln/level/compute.hh> +#include <mln/level/stretch.hh> + +#include <mln/math/abs.hh> +#include <mln/pw/all.hh> +#include <mln/trait/concrete.hh> +#include <mln/util/array.hh> + +#include <cmath> + +using namespace mln; +using value::int_u8; + + + +// FIXME: Protoype is wrong. Should be Image<I>. +template <typename I> +inline +void +save_tmax_map(util::array<image2d<I> >& arr_ima, unsigned dim3, const char* desc) +{ + image2d<float> ima_c; + initialize(ima_c, arr_ima[0]); + data::fill(ima_c, 0.0); + image2d<unsigned> ima_t; + initialize(ima_t, ima_c); + data::fill(ima_t, 0); + + mln_piter(image2d<float>) p(ima_c.domain()); + for_all(p) + for (unsigned k = 0; k < dim3; ++k) + if (ima_c(p) < arr_ima[k](p)) + { + ima_c(p) = arr_ima[k](p); + ima_t(p) = k; + } + + std::ostringstream smax; + smax << desc << "_max.pgm"; + std::ostringstream stime; + stime << desc << "_time.pgm"; + + image2d<int_u8> result_c = level::stretch(int_u8(), ima_c); + image2d<int_u8> result_t = level::stretch(int_u8(), ima_t); + io::pgm::save(result_c, smax.str()); + io::pgm::save(result_t, stime.str()); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 2) + { + std::cout << "Usage: " << argv[0] << " input" << std::endl; + return 1; + } + + image3d<float> input; + io::dump::load(input, argv[1]); + util::array<image2d<float> > arr_ima; + for (unsigned int i = 0; i < input.nslices(); ++i) + arr_ima.append(duplicate(cast_image<float>(slice(input, i)))); + + unsigned dim3 = arr_ima.nelements(); + + + ///////////// + // // + // Lissage // + // // + ///////////// + + util::array<image2d<double> > arr_smooth; + arr_smooth.append(arr_ima[0] * 1.0); + for (unsigned k = 1; k < dim3 - 1; ++k) + arr_smooth.append(arr_ima[k] * 0.5 + arr_ima[k - 1] * 0.25 + arr_ima[k + 1] * 0.25); + arr_smooth.append(arr_ima[dim3 - 1] * 1.0); + + + /////////////////////////////////// + // // + // Calcul image max et temps max // + // // + /////////////////////////////////// + + save_tmax_map(arr_ima, dim3, "tmax"); + save_tmax_map(arr_smooth, dim3, "tmax_smooth"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/check.sh =================================================================== --- trunk/milena/sandbox/fabien/igr/check.sh (revision 3694) +++ trunk/milena/sandbox/fabien/igr/check.sh (revision 3695) @@ -2,19 +2,20 @@ process_file () { - echo "Processing $3..." + echo "Processing $1..." | cowsay dist_max=10 input=$1 dim=$2 - ./grad $input $dim + ./crop $input 0 50 90 149 230 170 crop.dump + ./grad crop.dump $dim for lambda_closing in 50 100 500 1000 5000 10000 50000; do echo " for lambda_closing = ${lambda_closing}"; ./clo_vol grad.dump $dim ${lambda_closing} nbasins=`./wst clo_vol.dump $dim` echo " nbasins = $nbasins" - ../bin/dumpl32_to_colorize wst.dump $dim $nbasins results/colorize_${3}_${lambda_closing}.ppm - median=`./med wst.dump $dim $input $nbasins` + ../bin/dumpl16_to_colorize wst.dump $dim $nbasins results/colorize_${3}_${lambda_closing}.ppm + median=`./med wst.dump $dim crop.dump $nbasins` echo " median = $median" threshold=$(($median / 2)) ../bin/dumpi12_to_pgm med.dump $dim results/median_${3}_${lambda_closing}.pgm @@ -38,6 +39,8 @@ # rm *.dump } +make crop grad clo_vol wst med thres + #process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0049.dcm" 2 "49" process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0052.dcm" 3 "52" #process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0055.dcm" 2 "55" Index: trunk/milena/sandbox/fabien/igr/fun_labels.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/fun_labels.cc (revision 3694) +++ trunk/milena/sandbox/fabien/igr/fun_labels.cc (revision 3695) @@ -71,7 +71,7 @@ box3d slice = make::box3d(sli, min_row, min_col, sli, max_row, max_col); mln_VAR(slice_label, vol_label | slice); - accu::mean<int_u12> accu_mean; + accu::mean<mln_value(I)> accu_mean; float mean = level::compute(accu_mean, slice_label); arr.append(mean); } @@ -148,14 +148,15 @@ void plot_point(Image<I>& ima, Image<L>& ima_labels, point2d point, const char* desc) { util::array<float> arr; - label_16 prev_lbl; + mln_value(L) prev_lbl; int start = 0; int count = 1; + int nslices = geom::nslis(ima); - for (unsigned sli = 0; sli < geom::nslis(ima); ++sli) + for (unsigned sli = 0; sli < nslices; ++sli) { - image2d<int_u12> ima_slice = duplicate(slice(ima, sli)); - image2d<label_16> lbl_slice = duplicate(slice(ima_labels, sli)); + image2d<mln_value(I)> ima_slice = duplicate(slice(ima, sli)); + image2d<mln_value(L)> lbl_slice = duplicate(slice(ima_labels, sli)); if (sli == 0) prev_lbl = lbl_slice(point); if (lbl_slice(point) != prev_lbl) @@ -166,13 +167,13 @@ start = sli + 1; } // Taking the median value of the region. - accu::median_h<int_u12> accu_med; - int_u12 median = level::compute(accu_med, ima_slice | pw::value(lbl_slice) == pw::cst(lbl_slice(point))); - arr.append(median); + accu::mean<mln_value(I)> accu_mean; + mln_value(I) mean = level::compute(accu_mean, ima_slice | pw::value(lbl_slice) == pw::cst(lbl_slice(point))); + arr.append(mean); prev_lbl = lbl_slice(point); // Saving a image of the selected label in the current slice for debug. - data::fill((ima_slice | pw::value(lbl_slice) == pw::cst(prev_lbl)).rw(), 1750); + /*data::fill((ima_slice | pw::value(lbl_slice) == pw::cst(prev_lbl)).rw(), 1750); std::ostringstream str_ima; str_ima << "debug_" << desc << "_"; if (sli < 100) @@ -180,7 +181,7 @@ if (sli < 10) str_ima << "0"; str_ima << sli << ".pgm"; - io::pgm::save(level::stretch(int_u8(), ima_slice), str_ima.str()); + io::pgm::save(level::stretch(int_u8(), ima_slice), str_ima.str());*/ } if (!arr.is_empty()) @@ -198,6 +199,7 @@ int main(int argc, char *argv[]) { typedef label_16 L; + typedef float I; if (argc != 4) { @@ -208,8 +210,8 @@ L nlabels = atoi(argv[3]); - point2d p_tumeur(156, 114); - point2d p_air(34, 94); + point2d p_tumeur(163, 114); + point2d p_air(54, 94); point2d p_poumon(122, 115); image3d<L> ima_labels; io::dump::load(ima_labels, argv[1]); @@ -220,15 +222,15 @@ //plot_all_labels(ima, ima_labels, nlabels); - L l = 0; + /*L l = 0; for (unsigned i = 0; i < nlabels; ++i, l = i) - plot_label(ima, ima_labels, l); + plot_label(ima, ima_labels, l);*/ io::dump::save(ima_labels, "labels.dump"); - //plot_point(ima, ima_labels, p_tumeur, "tumeur"); - //plot_point(ima, ima_labels, p_air, "air"); - //plot_point(ima, ima_labels, p_poumon, "poumon"); + plot_point(ima, ima_labels, p_tumeur, "tumeur"); + plot_point(ima, ima_labels, p_air, "air"); + plot_point(ima, ima_labels, p_poumon, "poumon"); /*util::set<L> lbl_set; for (int sli = 0; sli < geom::nslis(ima_labels); ++sli) Index: trunk/milena/sandbox/fabien/igr/norm.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/norm.cc (revision 3694) +++ trunk/milena/sandbox/fabien/igr/norm.cc (revision 3695) @@ -18,6 +18,8 @@ #include <mln/data/paste.hh> // DEBUG +#include <mln/core/var.hh> +#include <mln/estim/min_max.hh> /*#include <mln/accu/min.hh> #include <mln/accu/max.hh> #include <mln/level/compute.hh>*/ @@ -60,18 +62,10 @@ image2d<float> ima_ini = mean_slices(ima_f, first, last); image3d<float> ima_result; - initialize(ima_result, ima_f); + initialize(ima_result, ima); for (int i = min_sli; i <= max_sli; ++i) - data::paste((slice(ima_f, i) - ima_ini), slice(ima_result, i).rw()); - - // DEBUG - /*accu::min<float> accu_min; - accu::max<float> accu_max; - float min = level::compute(accu_min, ima_ini); - float max = level::compute(accu_max, ima_ini); - std::cout << "min = " << min << std::endl; - std::cout << "max = " << max << std::endl;*/ + data::paste((slice(ima_f, i) - ima_ini) / ima_ini, slice(ima_result, i).rw()); return ima_result; } Index: trunk/milena/sandbox/fabien/igr/thres.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/thres.cc (revision 3694) +++ trunk/milena/sandbox/fabien/igr/thres.cc (revision 3695) @@ -34,7 +34,7 @@ using value::int_u12; using value::label_16; using value::label_32; - typedef label_32 L; + typedef label_16 L; if (argc != 4) { Index: trunk/milena/sandbox/fabien/igr/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/Makefile (revision 3694) +++ trunk/milena/sandbox/fabien/igr/Makefile (revision 3695) @@ -28,7 +28,7 @@ ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} $^ -o graph med: med.cc - ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} $^ -o med + ${CXX} -I../../../ ${CXXFLAGS} $^ -o med thres: thres.cc ${CXX} -I../../../ ${CXXFLAGS} $^ -o thres @@ -39,6 +39,9 @@ time_max: time_max.cc ${CXX} -I../../../ ${DICOM} ${MAGICK} ${CXXFLAGS} -lm $^ -o time_max +time_max_norm: time_max_norm.cc + ${CXX} -I../../../ ${CXXFLAGS} -lm $^ -o time_max_norm + first_slice_dicom: first_slice_dicom.cc ${CXX} -I../../../ ${DICOM} ${MAGICK} ${CXXFLAGS} -lm $^ -o first_slice_dicom @@ -57,6 +60,13 @@ norm: norm.cc ${CXX} -I../../../ ${CXXFLAGS} $^ -o norm +##################### + +seg_fixed: seg_fixed.cc + ${CXX} -I../../../ ${CXXFLAGS} $^ -o seg_fixed + +##################### + 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/fun_labels.sh =================================================================== --- trunk/milena/sandbox/fabien/igr/fun_labels.sh (revision 3694) +++ trunk/milena/sandbox/fabien/igr/fun_labels.sh (revision 3695) @@ -1,5 +1,9 @@ #!/bin/zsh +######### +# Plots # +######### + reconstruct_plot () { rm -f ${1}_${2}.plot @@ -10,11 +14,16 @@ rename_label_plots () { - for i in label_*.plot; do + for i in *_label_*.plot; do mv $i ${1}_$i done } + +############# +# Animation # +############# + create_anim () { echo " Creating ${1} animation..." @@ -43,15 +52,25 @@ rm debug_${1}_*.png } + + + + + + + + + process_file () { - echo "Processing $3..." + echo "Processing $1..." | cowsay input=$1 dim=$2 ./crop $input 0 50 90 149 230 170 crop.dump + ./norm crop.dump norm.dump ./grad crop.dump $dim - for lambda_closing in 10000; do + for lambda_closing in 1000; do echo " for lambda_closing = ${lambda_closing}"; ./clo_vol grad.dump $dim ${lambda_closing} nbasins=`./wst clo_vol.dump $dim` @@ -61,24 +80,25 @@ #../bin/dumpi12_to_png mean_slices.dump $dim mean_slices_${3}_${lambda_closing}.png #../bin/dumpi12_to_pgm mean_slices.dump $dim mean_slices_${3}_${lambda_closing}.pgm - ./norm crop.dump norm.dump ./fun_labels wst.dump norm.dump $nbasins - rename_label_plots ${lambda_closing} #./all_labels2gif.sh crop.dump labels.dump $nbasins ${lambda_closing} #mv *.gif results/ #mv *.plot results/plots/ -#reconstruct_plot tumeur ${lambda_closing} + reconstruct_plot tumeur ${lambda_closing} + reconstruct_plot air ${lambda_closing} + reconstruct_plot poumon ${lambda_closing} + rename_label_plots ${lambda_closing} #create_anim tumeur ${lambda_closing} -#reconstruct_plot air ${lambda_closing} #create_anim air ${lambda_closing} -#reconstruct_plot poumon ${lambda_closing} #create_anim poumon ${lambda_closing} done } +make crop grad clo_vol fun_labels norm label2gif + # 3D (2D + t) images only process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0052.dcm" 3 52 Index: trunk/milena/sandbox/fabien/igr/space_smooth/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/space_smooth/Makefile (revision 3694) +++ trunk/milena/sandbox/fabien/igr/space_smooth/Makefile (revision 3695) @@ -2,7 +2,7 @@ linear: linear.cc - ${CXX} -I../../../../ ${DICOM} ${CXXFLAGS} $^ -o linear + ${CXX} -I../../../../ ${CXXFLAGS} $^ -o linear median: median.cc ${CXX} -I../../../../ ${DICOM} ${CXXFLAGS} $^ -o median Index: trunk/milena/sandbox/fabien/igr/space_smooth/linear.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/space_smooth/linear.cc (revision 3694) +++ trunk/milena/sandbox/fabien/igr/space_smooth/linear.cc (revision 3695) @@ -7,7 +7,7 @@ #include <mln/value/int_u12.hh> -#include <mln/io/dicom/load.hh> +#include <mln/io/dump/load.hh> #include <mln/io/plot/save.hh> #include <mln/accu/median_h.hh> @@ -29,17 +29,19 @@ return 1; } + typedef float I; + /////////////////// // // // 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) + image3d<I> input; + io::dump::load(input, argv[1]); + util::array<image2d<I> > arr_ima; + for (unsigned int i = 0; i < input.nslices(); ++i) { - image2d<int_u12> tmp_slice = duplicate(slice(input, i)); + image2d<I> tmp_slice = duplicate(slice(input, i)); arr_ima.append(tmp_slice); } @@ -61,14 +63,15 @@ // Outputs // // // ///////////// - image2d<util::array<float> > ima_result(input.nrows(), input.ncols()); + image2d<util::array<float> > ima_result; + initialize(ima_result, slice(input, 0)); mln_piter_(image2d<util::array<float> >) p(ima_linear[0].domain()); for_all(p) - for (int i = 0; i < ima_linear.nelements(); ++i) + for (unsigned int i = 0; i < ima_linear.nelements(); ++i) ima_result(p).append(ima_linear[i](p)); io::plot::save(ima_result(point2d(156, 114)), "linear_tumeur.plot"); - io::plot::save(ima_result(point2d(34, 94)), "linear_air.plot"); + io::plot::save(ima_result(point2d(54, 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/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/time_smooth/Makefile (revision 3694) +++ trunk/milena/sandbox/fabien/igr/time_smooth/Makefile (revision 3695) @@ -2,7 +2,7 @@ linear: linear.cc - ${CXX} -I../../../../ ${DICOM} ${CXXFLAGS} $^ -o linear + ${CXX} -I../../../../ ${CXXFLAGS} $^ -o linear median: median.cc ${CXX} -I../../../../ ${DICOM} ${CXXFLAGS} $^ -o median Index: trunk/milena/sandbox/fabien/igr/time_smooth/linear.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/time_smooth/linear.cc (revision 3694) +++ trunk/milena/sandbox/fabien/igr/time_smooth/linear.cc (revision 3695) @@ -7,11 +7,13 @@ #include <mln/value/int_u12.hh> -#include <mln/io/dicom/load.hh> +#include <mln/io/dump/all.hh> #include <mln/io/plot/save.hh> #include <mln/accu/median_h.hh> #include <mln/util/array.hh> +#include <mln/linear/convolve.hh> +#include <mln/make/w_window3d.hh> using namespace mln; @@ -26,33 +28,48 @@ return 1; } + typedef float I; + /////////////////// // // // 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) + image3d<I> input; + io::dump::load(input, argv[1]); + image2d<util::array<I> > ima_arr; + initialize(ima_arr, duplicate(slice(input, 0))); + for (unsigned 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()); + image2d<I> tmp_slice = duplicate(slice(input, i)); + mln_piter_(image2d<I>) p(tmp_slice.domain()); for_all(p) { ima_arr(p).append(tmp_slice(p)); } } + /*image3d<I> output; + initialize(output, input);*/ //////////////////////// // // // Linear convolution // // // //////////////////////// - image2d<util::array<int_u12> > ima_linear; + /*float ws[][] = {{0, 0, 0, + 0, 0.25, 0, + 0, 0, 0}, + {0, 0, 0, + 0, 0.5, 0, + 0, 0, 0}, + {0, 0, 0, + 0, 0.25, 0, + 0, 0, 0}}; + + output = linear::convolve(input, make::w_window3d(ws));*/ + image2d<util::array<I> > ima_linear; initialize(ima_linear, ima_arr); - mln_piter_(image2d<int_u12>) p(ima_linear.domain()); - accu::median_h<int_u12> accu_med; + mln_piter_(image2d<I>) p(ima_linear.domain()); for_all(p) { ima_linear(p).append(ima_arr(p)[0]); @@ -66,10 +83,11 @@ // Outputs // // // ///////////// + //io::dump::save(output, "time_linear.dump"); io::plot::save(ima_arr(point2d(156, 114)), "ref_tumeur.plot"); io::plot::save(ima_linear(point2d(156, 114)), "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(54, 94)), "ref_air.plot"); + io::plot::save(ima_linear(point2d(54, 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"); Index: trunk/milena/sandbox/fabien/igr/med.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/med.cc (revision 3694) +++ trunk/milena/sandbox/fabien/igr/med.cc (revision 3695) @@ -9,7 +9,6 @@ #include <mln/core/alias/window3d.hh> #include <mln/io/dump/all.hh> -#include <mln/io/dicom/load.hh> #include <mln/value/int_u8.hh> #include <mln/value/int_u12.hh> @@ -45,7 +44,7 @@ using value::int_u12; using value::label_16; using value::label_32; - typedef label_32 L; + typedef label_16 L; if (argc != 5) { @@ -70,12 +69,12 @@ { image2d<L> labels; io::dump::load(labels, argv[1]); - image2d<int_u12> dcm; - io::dicom::load(dcm, argv[3]); + image2d<int_u12> input; + io::dump::load(input, argv[3]); mln_VAR(wst_dilate, morpho::elementary::dilation(extend(labels | (pw::value(labels) == 0u), labels), c8())); data::fill((labels | (pw::value(labels) == 0u)).rw(), wst_dilate); - mln_VAR(wst_mean, labeling::mean_values(dcm, labels, nbasins)); + mln_VAR(wst_mean, labeling::mean_values(input, labels, nbasins)); histo::array<int_u12> histogram = histo::compute(wst_mean); int j = 0; @@ -106,12 +105,12 @@ { image3d<L> labels; io::dump::load(labels, argv[1]); - image3d<int_u12> dcm; - io::dicom::load(dcm, argv[3]); + image3d<int_u12> input; + io::dump::load(input, argv[3]); mln_VAR(wst_dilate, morpho::elementary::dilation(extend(labels | (pw::value(labels) == 0u), labels), c26())); data::fill((labels | (pw::value(labels) == 0u)).rw(), wst_dilate); - mln_VAR(wst_mean, labeling::mean_values(dcm, labels, nbasins)); + mln_VAR(wst_mean, labeling::mean_values(input, labels, nbasins)); histo::array<int_u12> histogram = histo::compute(wst_mean); int j = 0; Index: trunk/milena/sandbox/fabien/TODO =================================================================== --- trunk/milena/sandbox/fabien/TODO (revision 3694) +++ trunk/milena/sandbox/fabien/TODO (revision 3695) @@ -47,6 +47,20 @@ [ ] Fix mean slices values [X] Compatibility Windows [X] Fix spatial smooth -[ ] Crop volume 52 -[ ] Plot all labels (3, 4 labels for tumeur) +[X] Crop volume 52 +[X] Plot all labels (3, 4 labels for tumeur) [ ] Prepare data for IGR + [X] crop->norm->tmax + [X] crop->norm->space_smooth->tmax + + (comparo normalisé) + 1)brut + 2)time_smooth + 3)space_smooth + 4)current_label_smooth (lambda faible, pas de ligne watershed) + + [X] plot label + [X] gif animés labels + + [X] fausses couleurs watershed + [X] images moyennes