[Olena-patches] r3585: Update IGR related code

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-03-27 Fabien Freling &lt;fabien.freling(a)lrde.epita.fr&gt; Update IGR related code. * fabien/igr/Makefile: Update. * fabien/igr/matlab.cc: Update. * fabien/igr/time_max.cc: Create tmax image. --- TODO | 5 igr/Makefile | 46 +++----- igr/matlab.cc | 296 ++++++++++++++++++++++++++++++++++---------------------- igr/time_max.cc | 113 +++++++++++++++++++++ 4 files changed, 318 insertions(+), 142 deletions(-) Index: trunk/milena/sandbox/fabien/igr/matlab.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/matlab.cc (revision 3584) +++ trunk/milena/sandbox/fabien/igr/matlab.cc (revision 3585) @@ -36,6 +36,8 @@ mln_precondition(first >=0 && first < arr_ima.nelements()); mln_precondition(last >=0 && last < arr_ima.nelements()); + initialize(result, arr_ima[first]); + accu::image::init(result); for (int i = first; i < last; ++i) accu::image::take(result, arr_ima[i]); @@ -58,10 +60,10 @@ 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 + const int acqui = 3; // durée acqui et temps interimages; en secondes + const int fseuil = 4; // ratio entre signaux trop faibles, éliminés, et bruit de fond mesuré sur l image + const int ini = 9; // nombre images ligne de base + const 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); @@ -86,10 +88,19 @@ data::fill(prodsignal1, 0.0); data::paste(datasli | pw::value(roi_noise) == pw::cst(true), prodsignal1); - // moyenne + ///////////// + // // + // Moyenne // + // // + ///////////// accu::mean<float> accu_mean; float moysignal1 = level::compute(accu_mean, prodsignal1 | pw::value(roi_noise) == pw::cst(true)); - // ecart type + + //////////////// + // // + // Ecart type // + // // + //////////////// int som1 = 0; int kk = 0; mln_piter_(image2d<float>) p(datasli.domain()); @@ -99,11 +110,16 @@ float ectys = std::sqrt(som1 / (nbrpix1 - 1)); float seuil = fseuil * ectys; - // Calcul du masque + ////////////////////// + // // + // Calcul du masque // + // // + ////////////////////// image2d<bool> masque; initialize(masque, datasli); data::fill(masque, false); 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 @@ -123,6 +139,12 @@ // On cherche les maxi (intensité ima_c et temps ima_t) des courbes signal(temps) // On masque toute la série d'images + ///////////// + // // + // Lissage // + // // + ///////////// + // Essai de lissage à parir de fin de montée vaculaire ini2 pour ameliorer // les seuillages a 10 50 90 NON EVALUE RIGOUREUSEMENT util::array<image2d<double> > arr_smooth; @@ -131,17 +153,33 @@ 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 // + // // + /////////////////////////////////// + image2d<float> ima_c; + initialize(ima_c, arr_smooth[0]); + data::fill(ima_c, 0.0); + image2d<unsigned> ima_t; + initialize(ima_t, ima_c); for_all(p) if (masque(p) == false) { 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; + ima_t(p) = 0; ++kk; } - // FIXME: [c,T]=max(image,[ ],3); + for_all(p) + for (unsigned k = 0; k < dim3; ++k) + if (ima_c(p) < arr_smooth[k](p)) + { + ima_c(p) = arr_smooth[k](p); + ima_t(p) = k; + } // Ou 'ima_c' est la valeur du max et 'ima_t' son index, le long de la dimension 3 // kk est le nombre de points masques @@ -149,151 +187,179 @@ std::cout << "kk = " << kk << std::endl; image2d<accu::sum<float> > accu_sum; - accu::image::init(result); - for (int i = first; i < last; ++i) - accu::image::take(result, arr_ima[i]); - ima_auc accu::image::to_result(result); + accu::image::init(accu_sum); + for (int k = 0; k < dim3; ++k) + accu::image::take(accu_sum, arr_smooth[k]); + image2d<float> ima_auc = accu::image::to_result(accu_sum); // Conversion du temps du pic en secondes à partir fin période de base - tmax = acqui * (ima_t - ini); + image2d<int> 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) + image2d<float> ima_c10; + initialize(ima_c10, arr_smooth[0]); + data::fill(ima_c10, 0.0); + image2d<unsigned> ima_idx10; + initialize(ima_idx10, ima_c10); + bool set_10 = false; + + image2d<float> ima_c90; + initialize(ima_c90, arr_smooth[0]); + data::fill(ima_c90, 0.0); + image2d<unsigned> ima_idx90; + initialize(ima_idx90, ima_c90); + bool set_90 = false; + + image2d<float> ima_c50; + initialize(ima_c50, arr_smooth[0]); + data::fill(ima_c50, 0.0); + image2d<unsigned> ima_idx50; + initialize(ima_idx50, ima_c50); + bool set_50 = false; + + image2d<bool> masque2; + initialize(masque2, arr_smooth[0]); + + // FIXME: This could be done with a bkw_piter. + for_all(p) { - [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; + for (unsigned k = 0; k < dim3; ++k) + { + if (!set_10 && arr_smooth[k](p) >= 0.1 * ima_c(p)) + { + ima_c10(p) = arr_smooth[k](p); + ima_idx10(p) = k; + set_10 = true; + } + if (!set_90 && arr_smooth[k](p) >= 0.9 * ima_c(p)) + { + ima_c90(p) = arr_smooth[k](p); + ima_idx90(p) = k; + set_90 = true; + } + if (!set_50 && arr_smooth[k](p) >= 0.5 * ima_c(p)) + { + ima_c50(p) = arr_smooth[k](p); + ima_idx50(p) = k; + set_50 = true; + } + } + masque2(p) = (ima_idx90(p) != ima_idx10(p)); } - // calcul MTT + //////////////// + // // + // 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 + image2d<bool> masque3; + initialize(masque3, masque2); 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)) + for_all(p) + if (arr_smooth[dim3 - 1](p) >= 0.9 * ima_c(p)) { - masque3(i,j)=0; - kk2=kk2+1; + masque3(p) = false; + ++kk2; } else - masque3(i,j)=1; + masque3(p) = true; // 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) + util::array<image2d<float> > arr_inv; for (int k = 0; k < dim3; ++k) - imageinv(i,j,k)=masque3(i,j)*image(i,j,dim3+1-k); + { + arr_inv.append(arr_ima[dim3 - 1 - k]); + data::fill((arr_inv[k] | pw::value(masque3) == false).rw(), 0.0); + } - for (int i = 0; i < dim1; ++i) - for (int j = 0; j < dim2; ++j) - if (masque3(i,j)==0) + image2d<float> ima_c51; + initialize(ima_c51, arr_smooth[0]); + data::fill(ima_c51, 0.0); + image2d<unsigned> ima_idx51; + initialize(ima_idx51, ima_c51); + bool set_51 = false; + for_all(p) + { + if (masque3(p)) { - idx_51(i,j)=0; - C51(i,j)=0; + for (unsigned k = 0; k < dim3; ++k) + if (!set_50 && arr_smooth[k](p) >= 0.8 * ima_c(p)) + { + ima_c51(p) = arr_inv[k](p); + ima_idx51(p) = k; + set_51 = true; + } } else - [C51(i,j),idx_51(i,j)]=max(imageinv(i,j,:)>=0.8*c(i,j)); + { + ima_c51(p) = 0; + ima_idx51(p) = 0; + } + } 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)); + image2d<int> ima_mtt; + initialize(ima_mtt, masque); + data::fill(ima_mtt, 0); + for_all(p) + if (masque(p) && masque2(p) && masque3(p)) + ima_mtt(p) = acqui * (dim3 - ima_idx51(p) - ima_idx50(p)); // 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))); + image2d<float> ima_pente; + initialize(ima_pente, masque); + data::fill(ima_pente, 0); + for_all(p) + if (masque(p) && masque2(p)) + ima_pente(p) = (0.8 * ima_c(p)) / (acqui * (ima_idx90(p) - ima_idx10(p))); - // 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 + //////////////////////////////////////////////// + // // + // Normalisation de la pente, du rehaussement // + // et de l'AUC // + // // + //////////////////////////////////////////////// + image2d<float> ima_pentenorm; + initialize(ima_pentenorm, masque); + data::fill(ima_pentenorm, 0); + + image2d<float> ima_reh; + initialize(ima_reh, masque); + data::fill(ima_reh, 0); + + image2d<float> ima_aucnorm; + initialize(ima_aucnorm, masque); + data::fill(ima_aucnorm, 0); + + for_all(p) + if (masque(p)) { - 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');*/ + ima_pentenorm(p) = ima_pente(p) / imageini(p); + ima_reh(p) = ima_c(p) / imageini(p); + ima_aucnorm(p) = ima_auc(p) / imageini(p); + } // 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)) + float x = masque.nrows() / 2; + float y = masque.ncols() / 2; + /*while (x < masque.nrows() && y < masque.ncols()) { // 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) + h=figure(8); pixval; // saisie coordonnées du pixel pour lequel on veut afficher la courbe de // décroissance et le fit @@ -305,12 +371,12 @@ // 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)); + 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) + subplot(2,2,4); plot(tata); - } + }*/ return 0; } Index: trunk/milena/sandbox/fabien/igr/time_max.cc =================================================================== --- trunk/milena/sandbox/fabien/igr/time_max.cc (revision 0) +++ trunk/milena/sandbox/fabien/igr/time_max.cc (revision 3585) @@ -0,0 +1,113 @@ +#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/dicom/load.hh> +#include <mln/io/magick/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; + + +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()); + + initialize(result, arr_ima[first]); + + 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 (unsigned int i = 0; i < input.nslices(); ++i) + arr_ima.append(duplicate(cast_image<float>(slice(input, i)))); + + // Calcul signal initialmoyen : 8 images de 2 à 9 moyennées=ligne de base + image2d<float> imageini = mean_slices(arr_ima, 1, 8); + + 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 // + // // + /////////////////////////////////// + image2d<float> ima_c; + initialize(ima_c, arr_smooth[0]); + data::fill(ima_c, 0.0); + image2d<unsigned> ima_t; + initialize(ima_t, ima_c); + data::fill(ima_c, 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; + } + + image2d<int_u8> result_c = level::stretch(int_u8(), ima_c); + image2d<int_u8> result_t = level::stretch(int_u8(), ima_t); + io::magick::save(result_c, "result_max.png"); + io::magick::save(result_t, "result_time.png"); + + return 0; +} Index: trunk/milena/sandbox/fabien/igr/Makefile =================================================================== --- trunk/milena/sandbox/fabien/igr/Makefile (revision 3584) +++ trunk/milena/sandbox/fabien/igr/Makefile (revision 3585) @@ -1,38 +1,30 @@ CXX = llvm-g++ -DCM_SRC=/Users/HiSoKa/Downloads/gdcm-2.0.10 -DCM_BIN=/Users/HiSoKa/Downloads/gdcmbin -DICOM_INC = -I${DCM_SRC}/Source/Common/ \ - -I${DCM_BIN}/Source/Common/ \ - -I${DCM_SRC}/Source/DataDictionary/ \ - -I${DCM_SRC}/Source/MediaStorageAndFileFormat/ \ - -I${DCM_SRC}/Source/DataStructureAndEncodingDefinition/ -DICOM_LIBS = -L${DCM_BIN}/bin \ - -lgdcmCommon -lgdcmDICT -lgdcmDSED -lgdcmIOD -lgdcmMSFF -lgdcmexpat -lgdcmjpeg12 -lgdcmjpeg16 -lgdcmjpeg8 -lgdcmopenjpeg -lgdcmuuid -lgdcmzlib \ +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 -CXXFLAGS = -DNDEBUG -O1 +DICOM = ${DICOM_INC} ${DICOM_LIBS} -all: 2d 3d wsd2d wsd3d grad +MAGICK = `Magick++-config --cppflags --ldflags --libs` -2d: seg_vol_irm.hh seg2d.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} seg2d.cc -o seg2d +seg2d: seg_vol_irm.hh seg2d.cc + ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} seg2d.cc -o seg2d -3d: seg_vol_irm.hh seg3d.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} seg3d.cc -o seg3d +seg3d: seg_vol_irm.hh seg3d.cc + ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} seg3d.cc -o seg3d wsd2d: watershed.hh watershed2d.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} $^ -o wsd2d + ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} $^ -o wsd2d wsd3d: watershed.hh watershed3d.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} $^ -o wsd3d - -wsd3dg: watershed.hh watershed3d.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} -g $^ -o wsd3dg + ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} $^ -o wsd3d nbasins: nbasins_finder.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} $^ -o nbasins_finder + ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} $^ -o nbasins_finder grad: grad.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} $^ -o grad + ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} $^ -o grad clo_vol: clo_vol.cc ${CXX} -I../../../ ${CXXFLAGS} $^ -o clo_vol @@ -41,16 +33,20 @@ ${CXX} -I../../../ ${CXXFLAGS} $^ -o wst graph: graph.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} $^ -o graph + ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} $^ -o graph med: med.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} $^ -o med + ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} $^ -o med thres: thres.cc ${CXX} -I../../../ ${CXXFLAGS} $^ -o thres matlab: matlab.cc - ${CXX} -I../../../ ${DICOM_INC} ${DICOM_LIBS} ${CXXFLAGS} -lm $^ -o matlab + ${CXX} -I../../../ ${DICOM} ${CXXFLAGS} -lm $^ -o matlab + +time_max: time_max.cc + ${CXX} -I../../../ ${DICOM} ${MAGICK} ${CXXFLAGS} -lm $^ -o time_max 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 Index: trunk/milena/sandbox/fabien/TODO =================================================================== --- trunk/milena/sandbox/fabien/TODO (revision 3584) +++ trunk/milena/sandbox/fabien/TODO (revision 3585) @@ -36,10 +36,11 @@ [X] 3D [X] Print nb bg regions // nb fg objets [ ] Profile for performance -[ ] ImageMagick support +[X] ImageMagick support [ ] Integrate external libraries (GDCM, IM) [ ] Send result images to lrde account -[ ] Translate Matlab code +[X] Translate Matlab code [X] Subsample binary images [X] Fast projected histogram [ ] Triple histogram +[ ] Create plot for each kind of point with each method