[Olena-patches] r3576: Update Fabien' sandbox

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-03-25 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; Update Fabien' sandbox. * fabien/binarization/Makefile: Update. * fabien/binarization/test.cc: Update. * fabien/igr/matlab.cc: Update. * fabien/magick/Makefile: Update. * fabien/magick/save.cc: New test file for IM support. --- binarization/Makefile | 13 +-- binarization/test.cc | 123 +++++++++++++++++++++----------- igr/matlab.cc | 187 +++++++++++++++++++++++++++++++++++++++++++++----- magick/Makefile | 14 +++ magick/save.cc | 23 ++++++ 5 files changed, 292 insertions(+), 68 deletions(-) Index: trunk/milena/sandbox/fabien/igr/matlab.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/matlab.cc (revision 3575) +++ trunk/milena/sandbox/fabien/igr/matlab.cc (revision 3576) @@ -14,8 +14,7 @@ #include <mln/accu/sum.hh> #include <mln/accu/image/all.hh> #include <mln/arith/minus.hh> -#include <mln/arith/plus.hh> -#include <mln/arith/times.hh> +#include <mln/arith/all.hh> #include <mln/data/fill.hh> #include <mln/data/paste.hh> #include <mln/level/compute.hh> @@ -95,10 +94,8 @@ 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; @@ -106,40 +103,43 @@ image2d<bool> masque; 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 + data::fill((masque | pw::value(datasli) > pw::cst(seuil)).rw(), true); + // si on a choisi une région avec roi + int nargin = 0; + image2d<bool> roi; // FIXME: init this ROI, should be a domain if (nargin>2) - data::fill(masque | pw::value(roi2) == false, false); + data::fill((masque | pw::value(roi) == false).rw(), false); // On applique le masque sur image et imasoustraite for (int k = 0; k < arr_ima.nelements(); ++k) { - data::fill(arr_ima[k] | pw::value(masque) == false, 0.0); - data::fill(arr_sous[k] | pw::value(masque) == false, 0.0); + data::fill((arr_ima[k] | pw::value(masque) == false).rw(), 0.0); + data::fill((arr_sous[k] | pw::value(masque) == false).rw(), 0.0); } // On regarde si le seuillage et le masquage sont OK // FIXME: imagesc(abs(image(:,:,2))); - // On cherche les maxi (intensité c et temps T) des courbes signal(temps) + // On cherche les maxi (intensité ima_c et temps ima_t) des courbes signal(temps) // On masque toute la série d'images // Essai de lissage à parir de fin de montée vaculaire ini2 pour ameliorer // les seuillages a 10 50 90 NON EVALUE RIGOUREUSEMENT - for (int k = init2 - 1; k < dim3 - 1; ++k) - arr_ima[k] = 0.5 * arr_ima[k] + 0.25 * arr_ima[k - 1] + 0.25 * arr_ima[k + 1]; + util::array<image2d<double> > arr_smooth; + arr_smooth.append(arr_ima[0] * 1.0); + for (int k = ini2 - 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); for_all(p) - { if (masque(p) == false) { - for (int k = 0; k < arr_ima.nelements(); ++k) - arr_ima[k](p) = 0.0; + for (int k = 0; k < arr_smooth.nelements(); ++k) + arr_smooth[k](p) = 0.0; ima_c(p) = 0.0; ima_t(p) = 0.0; ++kk; } - } // FIXME: [c,T]=max(image,[ ],3); @@ -157,5 +157,160 @@ // Conversion du temps du pic en secondes à partir fin période de base tmax = acqui * (ima_t - ini); + // calcul des temps 10 et 90 et de la pente correspondante + for (int i = 0; i < dim1; ++i) + for (int j = 0; j < dim2; ++j) + { + [C10(i,j),idx_10(i,j)]=max(image(i,j,:)>=0.1*c(i,j)); + [C90(i,j),idx_90(i,j)]=max(image(i,j,:)>=0.9*c(i,j)); + [C50(i,j),idx_50(i,j)]=max(image(i,j,:)>=0.5*c(i,j)); + if (idx_90(i,j)-idx_10(i,j)==0) + masque2(i,j)=0; + else + masque2(i,j)=1; + } + + // calcul MTT + + // QUAND le signal redescent significativement après le maximum + // elimination des voxels où le signal est encore 90// du maximum en fin + // d'acquisition + int kk2 = 0; + for (int i = 0; i < dim1; ++i) + for (int j = 0; j < dim2; ++j) + if (image(i,j,dim3)>=0.9*c(i,j)) + { + masque3(i,j)=0; + kk2=kk2+1; + } + else + masque3(i,j)=1; + + // pour les voxels non masqués on inverse t et on cherche le max 50// a partir de la fin + for (int i = 0; i < dim1; ++i) + for (int j = 0; j < dim2; ++j) + for (int k = 0; k < dim3; ++k) + imageinv(i,j,k)=masque3(i,j)*image(i,j,dim3+1-k); + + for (int i = 0; i < dim1; ++i) + for (int j = 0; j < dim2; ++j) + if (masque3(i,j)==0) + { + idx_51(i,j)=0; + C51(i,j)=0; + } + else + [C51(i,j),idx_51(i,j)]=max(imageinv(i,j,:)>=0.8*c(i,j)); + + std::cout << "kk2 = " << kk2 << std::endl; + + // Calcul MTT defini comme largeur a mi hauteur de la courbe signal(temps) + // soit t50 montant et t50 descendant + for (int i = 0; i < dim1; ++i) + for (int j = 0; j < dim2; ++j) + if (masque(i,j)*masque2(i,j)*masque3(i,j)==0) + MTT(i,j)=0; + else + MTT(i,j)=acqui*(dim3-idx_51(i,j)-idx_50(i,j)); + + // calcul de la pente 10/90 (avec delta signal=80// du maximum c(i,j) et + // delta temps =t90 -t10 + for (int i = 0; i < dim1; ++i) + for (int j = 0; j < dim2; ++j) + if (masque(i,j)*masque2(i,j)==0) + pente(i,j)=0; + else + pente(i,j)=(0.8*c(i,j))/(acqui*(idx_90(i,j)-idx_10(i,j))); + + // calcul pente normalisée rehaussement normalisé et AUC normalisée + for (int i = 0; i < dim1; ++i) + for (int j = 0; j < dim2; ++j) + if (masque(i,j)==0) + { + pentenorm(i,j)=0; + reh(i,j)=0; + AUCnorm(i,j)=0; + } + else + { + pentenorm(i,j)=pente(i,j)/imageini(i,j); + reh(i,j)=c(i,j)/imageini(i,j); + AUCnorm(i,j)=AUC(i,j)/imageini(i,j); + } + + /*figure (3) + imagesc(AUCnorm); + title('image aire normalisée sous la courbe'); + // l'aire sous la courbe est négative là où il y a eu du bougé (bord tube) + + + figure (4) + imagesc(pentenorm); + title('pente normalisée10-90'); + + figure(5) + imagesc(reh); + title('image rehaussement relatif'); + caxis([0 5]); + // au dela de 5 il faut réviser la dose d'agent de contraste car surdosage + + figure(6) + imagesc(MTT); + title('image MTT'); + + figure(9) + imagesc(idx_51); + title('idx_51'); + + figure(10) + imagesc(C51); + + figure(7) + imagesc(tmax); + title('image Tmaximum'); + + figure(8) + subplot(2,2,1); + imagesc(pentenorm); + title('pente normalisée10-90'); + subplot(2,2,2) + imagesc(reh); + title('image rehaussement relatif'); + caxis([0 3]); + subplot(2,2,3) + imagesc(tmax); + title('image Tmaximum');*/ + + // COURBE S(t)sur pixel choisi graphiquement + // Permet de voir la courbe de rehaussement d'un voxel pointé + // SI LE VOXEL EST AU BORD DU TUBE TEMOIN: UN DETECTEUR DE MOUVEMENTS + // coordonx=1:size(image,3); + x = size(image, 1) / 2; + y = size(image, 2) / 2; + while (x < size(image, 1) && y < size(image, 2)) + { + // la saisie du point de coordonnées x,y se fait sur la figure ouverte en + // dernier + // pour sortir: cliquer dans le gris... + // pour affichage interactif coordonnées pixel, mais en fait pas utile + h=figure(8) + pixval; + // saisie coordonnées du pixel pour lequel on veut afficher la courbe de + // décroissance et le fit + [x,y]=ginput(1); + // on ajuste a l'entier supérieur pour trouver les "vraies" coordonnées + // en plus il faut inverser entre xy et ij... + coordx=ceil(y); + coordy=ceil(x); + // test sortie du programme + if ((coordx > size(image, 1)) | (coordy > size(image, 2))) + break; + fprintf(1,'\n// d\t// d\t// f\t// f\n',coordx,coordy,reh(coordx,coordy),image(coordx,coordy,1)); + titi = (abs(image(coordx, coordy, :))); + tata = squeeze(titi); + subplot(2,2,4) + plot(tata); + } + return 0; } Index: trunk/milena/sandbox/fabien/binarization/test.cc =================================================================== --- trunk/milena/sandbox/fabien/binarization/test.cc (revision 3575) +++ trunk/milena/sandbox/fabien/binarization/test.cc (revision 3576) @@ -1,30 +1,37 @@ #include <iostream> #include <fstream> +#include <mln/accu/median_h.hh> + #include <mln/algebra/vec.hh> +#include <mln/convert/from_to.hh> + #include <mln/core/image/image1d.hh> #include <mln/core/image/image2d.hh> #include <mln/core/alias/neighb1d.hh> #include <mln/draw/line.hh> -#include <mln/io/magick/load.hh> -#include <mln/io/pgm/save.hh> -#include <mln/io/pbm/save.hh> +#include <mln/io/magick/all.hh> +#include <mln/level/convert.hh> #include <mln/level/stretch.hh> +#include <mln/literal/colors.hh> + #include <mln/make/box2d.hh> #include <mln/morpho/closing/volume.hh> #include <mln/morpho/watershed/flooding.hh> +#include <mln/util/array.hh> #include <mln/util/couple.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/world/binary_2d/projected_histo.hh> #include <mln/world/binary_2d/subsample.hh> @@ -34,75 +41,109 @@ using value::int_u8; using value::int_u12; using value::label_16; +using value::rgb8; + + +int draw_lines(image2d<bool>& ima, int col_min, int col_max) +{ + box2d box = make::box2d(0, col_min, ima.nrows() - 1, col_max); + util::couple<image1d<float>, image1d<float> > histos = world::binary_2d::projected_histo(ima | box); + image1d<float> row_histo = histos.first(); + + /*int i = col_histo.ninds(); + while (i > 0 && col_histo.at_(i) < 60.0) + { + col_histo.at_(i) = 100.0; + --i; + }*/ + + row_histo = morpho::closing::volume(row_histo, c2(), 41); + + // Watershed + label_16 nbasins; + image1d<label_16> row_labels = morpho::watershed::flooding(row_histo, c2(), nbasins); + + int median = 0; + util::array<int> inter_lines; + int begin = 0; + int end = 0; + + // Draw lines + for (unsigned int i = 0; i < row_labels.ninds(); ++i) + if (row_labels.at_(i) == 0u) + { + if (end == 0) + { + begin = i; + end = i; + } + else if (end == begin) + { + end = i; + inter_lines.append(end - begin); + begin = end; + } + } + + image1d<int> ima_lines; + convert::from_to(inter_lines, ima_lines); + accu::median_h<int_u12> accu_med; + median = level::compute(accu_med, ima_lines); + + // Gnuplot files creation + /*std::ofstream fout_row("row.plot"); + for (unsigned int i = 0; i < row_histo.ninds(); ++i) + fout_row << i << " " << row_histo.at_(i) << std::endl; + + std::ofstream fout_col("col.plot"); + for (unsigned int i = 0; i < col_histo.ninds(); ++i) + fout_col << i << " " << col_histo.at_(i) << std::endl;*/ + return median; +} int main(int argc, char* argv[]) { if (argc != 3) { - std::cout << "Usage: " << argv[0] << "input.ext output.pgm" << std::endl; + std::cout << "Usage: " << argv[0] << "input output" << std::endl; return 1; } image2d<bool> ima; io::magick::load(ima, argv[1]); + image2d<rgb8> ima_color = level::convert(rgb8(), ima); //image2d<int_u8> ima_sampled = world::binary_2d::subsample(ima, 3); //io::pgm::save(ima_sampled, argv[2]); - box2d left_box = make::box2d(0, 0, ima.nrows() - 1, (ima.ncols() / 3)); - util::couple<image1d<float>, image1d<float> > histos_left = world::binary_2d::projected_histo(ima | left_box); - - box2d middle_box = make::box2d(0, (ima.ncols() / 3), ima.nrows() - 1, (2 * ima.ncols() / 3)); - util::couple<image1d<float>, image1d<float> > histos_middle = world::binary_2d::projected_histo(ima | middle_box); - - box2d right_box = make::box2d(0, (2 * ima.ncols() / 3), ima.nrows() - 1, ima.ncols() - 1); - util::couple<image1d<float>, image1d<float> > histos_right = world::binary_2d::projected_histo(ima | right_box); - - image1d<float> row_histo = histos_right.first(); - image1d<float> col_histo = histos_right.second(); - - int i = col_histo.ninds(); - while (i > 0 && col_histo.at_(i) < 60.0) - { - col_histo.at_(i) = 100.0; - --i; - } + int median_left = draw_lines(ima, 0, (ima.ncols() / 3)); + int median_middle = draw_lines(ima, (ima.ncols() / 3), (2 * ima.ncols() / 3)); + int median_right = draw_lines(ima, (2 * ima.ncols() / 3), ima.ncols() - 1); - // Closing - row_histo = morpho::closing::volume(row_histo, c2(), 31); + int median = (median_left + median_middle + median_right) / 3; - // Watershed - label_16 nbasins; - image1d<label_16> row_labels = morpho::watershed::flooding(row_histo, c2(), nbasins); + std::cout << "median = " << median << std::endl; // Draw lines - for (unsigned int i = 0; i < row_labels.ninds(); ++i) - { - if (row_labels.at_(i) == 0u) - { algebra::vec<2, unsigned int> vmin; algebra::vec<2, unsigned int> vmax; + for (unsigned int i = 0; i < ima.nrows(); ++i) + { + if (i % median == 0) + { vmin[0] = i; vmin[1] = 0; vmax[0] = i; vmax[1] = ima.ncols() - 1; mln_site_(image2d<bool>) pbeg(vmin); mln_site_(image2d<bool>) pend(vmax); - draw::line(ima, pbeg, pend, 0); + draw::line(ima_color, pbeg, pend, literal::red); } } - io::pbm::save(ima, argv[2]); - // Gnuplot files creation - std::ofstream fout_row("row.plot"); - for (unsigned int i = 0; i < row_histo.ninds(); ++i) - fout_row << i << " " << row_histo.at_(i) << std::endl; - - std::ofstream fout_col("col.plot"); - for (unsigned int i = 0; i < col_histo.ninds(); ++i) - fout_col << i << " " << col_histo.at_(i) << std::endl; + io::magick::save(ima_color, argv[2]); return 0; } Index: trunk/milena/sandbox/fabien/binarization/Makefile =================================================================== --- trunk/milena/sandbox/fabien/binarization/Makefile (revision 3575) +++ trunk/milena/sandbox/fabien/binarization/Makefile (revision 3576) @@ -1,10 +1,5 @@ -GCCFLAGS = -DNDEBUG -O1 -LLVMFLAGS = -DNDEBUG -O4 +CXX = llvm-g++ +CXXFLAGS = -DNDEBUG -O4 -all: gcc llvm - -gcc: test.cc - g++ -I../../../ `Magick++-config --cppflags --cxxflags --ldflags --libs` ${GCCFLAGS} test.cc -o test_gcc - -llvm: test.cc - llvm-g++ -I../../../ `Magick++-config --cppflags --cxxflags --ldflags --libs` ${LLVMFLAGS} test.cc -o test_llvm +all: test.cc + ${CXX} -I../../../ `Magick++-config --cppflags --cxxflags --ldflags --libs` ${CXXFLAGS} test.cc -o test Index: trunk/milena/sandbox/fabien/magick/save.cc =================================================================== --- trunk/milena/sandbox/fabien/magick/save.cc (revision 0) +++ trunk/milena/sandbox/fabien/magick/save.cc (revision 3576) @@ -0,0 +1,23 @@ +#include <mln/core/image/image2d.hh> + +#include <mln/io/magick/all.hh> +#include <mln/io/pbm/all.hh> +#include <mln/io/pgm/all.hh> +#include <mln/io/ppm/all.hh> + +#include <mln/value/int_u8.hh> +#include <mln/value/rgb8.hh> + +using namespace mln; + +int main(int argc, char* argv[]) +{ + if (argc != 2) + return 1; + + image2d<value::rgb8> ima; + io::magick::load(ima, argv[1]); + io::magick::save(ima, "out_magick.png"); + io::ppm::save(ima, "out_pnm.ppm"); + return 0; +} Index: trunk/milena/sandbox/fabien/magick/Makefile =================================================================== --- trunk/milena/sandbox/fabien/magick/Makefile (revision 3575) +++ trunk/milena/sandbox/fabien/magick/Makefile (revision 3576) @@ -1,2 +1,12 @@ -all: magick.cc - g++ -I../../../ `Magick++-config --cppflags --cxxflags --ldflags --libs` magick.cc -o mln_magick +CXX = llvm-g++ +CXXFLAGS = -DNDEBUG -O4 +LIBS = `Magick++-config --cppflags --cxxflags --ldflags --libs` +INC = -I../../../ + +all: magick + +magick: magick.cc + ${CXX} ${INC} ${LIBS} $^ -o mln_magick + +save: save.cc + ${CXX} ${INC} ${LIBS} $^ -o magick_save