[Olena-patches] r3576: Update Fabien' sandbox

25 Mar 2009

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox

ChangeLog:
2009-03-25  Fabien Freling  <fabien.freling@lrde.epita.fr>

	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

    