[Olena-patches] r3545: Add a new file translating IGR Matlab code

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-03-19 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; Add a new file translating IGR Matlab code. * fabien/bin/dump2pbm.cc: Minor fix. * fabien/igr/Makefile: Update. * fabien/igr/check.sh: Update. * fabien/igr/matlab.cc: New file translating matlab code. * fabien/igr/med.cc: Update. * fabien/igr/thres.cc: Update. * fabien/magick/load.hh: Moved to mln/io/magick. * fabien/magick/magick.cc: Update. --- TODO | 12 +++-- bin/dump2pbm.cc | 28 +++++++++--- igr/Makefile | 3 + igr/check.sh | 10 ++-- igr/matlab.cc | 126 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ igr/med.cc | 47 ++++++++++++++++++++ igr/thres.cc | 7 ++- magick/magick.cc | 9 +++ 8 files changed, 223 insertions(+), 19 deletions(-) Index: trunk/milena/sandbox/fabien/igr/matlab.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/matlab.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/matlab.cc (revision 3545) @@ -0,0 +1,126 @@ +#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/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/data/fill.hh> +#include <mln/data/paste.hh> +#include <mln/level/compute.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; + + +inline +image2d<float> mean_slices(util::array<image2d<float> >& arr_ima, int first, int last) +{ + image2d<accu::mean<float> > result; + + mln_precondition(first >=0 && first < arr_ima.nelements()); + mln_precondition(last >=0 && last < arr_ima.nelements()); + + accu::image::init(result); + for (int i = first; i < last; ++i) + accu::image::take(result, arr_ima[i]); + + return accu::image::to_result(result); +} + + +int main(int argc, char* argv[]) +{ + if (argc != 2) + { + std::cout << "Usage: " << argv[0] << "input" << std::endl; + return 1; + } + + image3d<value::int_u12> input; + io::dicom::load(input, argv[1]); + util::array<image2d<float> > arr_ima; + for (int i = 0; i < input.nslices(); ++i) + arr_ima.append(duplicate(cast_image<float>(slice(input, i)))); + + int acqui = 3; // durée acqui et temps interimages; en secondes + int fseuil = 4; // ratio entre signaux trop faibles, éliminés, et bruit de fond mesuré sur l image + int ini = 9; // nombre images ligne de base + int ini2 = 20; // a la fin de la montée vasculaire, à la dixieme image post injection: lissage des images par trois + + // Calcul signal initialmoyen : 8 images de 2 à 9 moyennées=ligne de base + image2d<float> imageini = mean_slices(arr_ima, 1, 8); + + int dim1 = imageini.nrows(); + int dim2 = imageini.ncols(); + int dim3 = arr_ima.nelements(); + + // Calcul auc aire sous la courbe + util::array<image2d<float> > arr_sous; + for (int k = 0; k < arr_ima.nelements(); ++k) + arr_sous.append(arr_ima[k] - imageini); + + // Mesure bruit de fond pour seuiller + // calculé sur la première image + image2d<bool> roi_noise; // FIXME: init this ROI, should be a domain + accu::count<bool> accu_nbrpix1; + unsigned nbrpix1 = level::compute(accu_nbrpix1, roi_noise); + image2d<float> datasli = arr_ima[0]; + + image2d<float> prodsignal1; + data::fill(prodsignal1, 0.0); + data::paste(datasli | pw::value(roi_noise) == pw::cst(true), prodsignal1); + + // moyenne + accu::mean<float> accu_mean; + float moysignal1 = level::compute(accu_mean, prodsignal1 | pw::value(roi_noise) == pw::cst(true)); + // ecart type + int som1 = 0; + int kk = 0; + mln_piter_(image2d<float>) p(datasli.domain()); + for_all(p) + { + if (roi_noise(p)) + som1 += std::pow(math::abs(datasli(p) - moysignal1), 2); + } + float ectys = std::sqrt(som1 / (nbrpix1 - 1)); + float seuil = fseuil * ectys; + + // Calcul du masque + image2d<bool> masque; // FIXME: init this ROI, should be a domain + initialize(masque, datasli); + data::fill(masque, false) + data::fill(masque | pw::value(datasli) > pw::cst(seuil), true); + // si on a choisi une région avec roi2 + if (nargin>2) + masque=masque.*roi2; + + //on applique le masque sur image et imasoustraite + + for k=2:dim3 + for i=1:dim1 + for j=1:dim2 + if(masque(i,j)==0) + image(i,j,k)=0; + imasoustraite(i,j,k)=0; + end; + end; + end; + end; + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/check.sh =================================================================== --- trunk/milena/sandbox/fabien/igr/check.sh (revision 3544) +++ trunk/milena/sandbox/fabien/igr/check.sh (revision 3545) @@ -15,11 +15,11 @@ echo " nbasins = $nbasins" ../bin/dumpl32_to_colorize wst.dump $dim $nbasins results/colorize_${3}_${lambda_closing}.ppm median=`./med wst.dump $dim $input $nbasins` - echo " median = $median" +#echo " median = $median" threshold=$(($median / 2)) ../bin/dumpi12_to_pgm med.dump $dim results/median_${3}_${lambda_closing}.pgm ./thres med.dump $dim $threshold - mv result.pbm results/result_${3}_${lambda_closing}.pbm + ../bin/dump2pbm result.dump $dim results/result_${3}_${lambda_closing}.pbm #if [ $2 -eq 2 ]; then # if [ ${lambda_closing} -eq 100 ]; then @@ -39,8 +39,8 @@ } 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_0052.dcm" 3 "52" process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0055.dcm" 2 "55" process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0058.dcm" 2 "58" -#process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0061.dcm" 3 "61" -#process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0064.dcm" 3 "64" +process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0061.dcm" 3 "61" +process_file "/Users/HiSoKa/Work/IGR/souris18/irm/IM_0064.dcm" 3 "64" Index: trunk/milena/sandbox/fabien/igr/med.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/med.cc (revision 3544) +++ trunk/milena/sandbox/fabien/igr/med.cc (revision 3545) @@ -63,6 +63,8 @@ } unsigned median = 0; + int bg = 0; + int fg = 0; if (dim == 2) { @@ -88,14 +90,57 @@ accu::median_h<int_u12> accu_med; median = level::compute(accu_med, ima_histo | pw::value(ima_histo) != pw::cst(0)); + for (int i = 0; i < histogram.nvalues(); ++i) + { + if (histogram[i]) + { + if (histogram[i] > median / 2) + ++fg; + else ++bg; + } + } + io::dump::save(wst_mean, "med.dump"); } else { - // FIXME + image3d<L> labels; + io::dump::load(labels, argv[1]); + image3d<int_u12> dcm; + io::dicom::load(dcm, 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)); + + histo::array<int_u12> histogram = histo::compute(wst_mean); + int j = 0; + int k = 0; + for (int i = 0; i < histogram.nvalues(); ++i) + { + if (histogram[i]) + histogram[i] = i; + } + image1d<unsigned> ima_histo; + convert::from_to(histogram, ima_histo); + accu::median_h<int_u12> accu_med; + median = level::compute(accu_med, ima_histo | pw::value(ima_histo) != pw::cst(0)); + + for (int i = 0; i < histogram.nvalues(); ++i) + { + if (histogram[i]) + { + if (histogram[i] > median / 2) + ++fg; + else ++bg; + } + } + + io::dump::save(wst_mean, "med.dump"); } std::cout << median << std::endl; + std::cerr << " [ " << bg << " <" << median / 2 << "> " << fg << " ]" << std::endl; return 0; } Index: trunk/milena/sandbox/fabien/igr/thres.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/thres.cc (revision 3544) +++ trunk/milena/sandbox/fabien/igr/thres.cc (revision 3545) @@ -56,11 +56,14 @@ image2d<int_u12> input; io::dump::load(input, argv[1]); image2d<bool> bin_result = binarization::threshold(input, threshold); - io::pbm::save(bin_result, "result.pbm"); + io::dump::save(bin_result, "result.dump"); } else { - // FIXME + image3d<int_u12> input; + io::dump::load(input, argv[1]); + image3d<bool> bin_result = binarization::threshold(input, threshold); + io::dump::save(bin_result, "result.dump"); } return 0; Index: trunk/milena/sandbox/fabien/igr/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/Makefile (revision 3544) +++ trunk/milena/sandbox/fabien/igr/Makefile (revision 3545) @@ -48,5 +48,8 @@ thres: thres.cc g++ -I../../../ ${CXXFLAGS} $^ -o thres +matlab: matlab.cc + g++ -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} -lm $^ -o matlab + clean: rm -rf *.dump *.p?m *.plot *.log *.csv *.dSYM Index: trunk/milena/sandbox/fabien/TODO =================================================================== --- trunk/milena/sandbox/fabien/TODO (revision 3544) +++ trunk/milena/sandbox/fabien/TODO (revision 3545) @@ -30,7 +30,13 @@ [X] Test 3D workflow on 2D images [X] Create colorized colors for graph [X] Diff with % -[ ] Find biggest dark regions (threshold value or median accu - median / 2 - ) - [ ] Learn regions value - [ ] Threshold +[X] Find biggest dark regions (threshold value or median accu - median / 2 - ) + [X] Learn regions value + [X] Threshold + [X] 3D + [X] Print nb bg regions // nb fg objets [ ] Profile for performance +[ ] ImageMagick support +[ ] Integrate external libraries (GDCM, IM) +[ ] Send result images to lrde account +[ ] Translate Matlab code Index: trunk/milena/sandbox/fabien/bin/dump2pbm.cc =================================================================== --- trunk/milena/sandbox/fabien/bin/dump2pbm.cc (revision 3544) +++ trunk/milena/sandbox/fabien/bin/dump2pbm.cc (revision 3545) @@ -1,5 +1,6 @@ - #include <mln/core/image/image2d.hh> +#include <mln/core/image/image3d.hh> +#include <mln/core/image/slice_image.hh> #include <mln/make/image3d.hh> #include <mln/debug/slices_2d.hh> @@ -12,7 +13,7 @@ void usage(char* argv[]) { - std::cerr << "usage: " << argv[0] << " input.dump output.pbm" << std::endl; + std::cerr << "usage: " << argv[0] << " input.dump dim output.pbm" << std::endl; abort(); } @@ -23,13 +24,28 @@ using namespace mln; using value::int_u16; - if (argc != 3) + if (argc != 4) usage(argv); + unsigned dim = atoi(argv[2]); + if (dim != 2 && dim != 3) + { + std::cout << "dimensions invalid" << std::endl; + return 1; + } + + if (dim == 2) + { + image2d<bool> ima; + io::dump::load(ima, argv[1]); + io::pbm::save(ima, argv[3]); + } + else + { image3d<bool> vol; io::dump::load(vol, argv[1]); - int_u16 bg = 0; - image2d<bool> ima = debug::slices_2d(vol, 1.f, bg); - io::pbm::save(ima, argv[2]); + image2d<bool> ima = debug::slices_2d(vol, 1.f, false); + io::pbm::save(ima, argv[3]); + } } Index: trunk/milena/sandbox/fabien/magick/load.hh (deleted) =================================================================== Index: trunk/milena/sandbox/fabien/magick/magick.cc =================================================================== --- trunk/milena/sandbox/fabien/magick/magick.cc (revision 3544) +++ trunk/milena/sandbox/fabien/magick/magick.cc (revision 3545) @@ -1,11 +1,15 @@ #include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> #include <mln/value/rgb8.hh> -#include "load.hh" +#include <mln/io/magick/load.hh> +#include <mln/io/pgm/save.hh> +#include <mln/io/pbm/save.hh> int main(int argc, char* argv[]) { using namespace mln; + using value::int_u8; using value::rgb8; if (argc != 2) @@ -14,9 +18,10 @@ return 1; } - image2d<rgb8> lena; + image2d<bool> lena; io::magick::load(lena, argv[1]); + io::pbm::save(lena, "result.pbm"); return 0; }