r3545: Add a new file translating IGR Matlab code

19 Mar 2009

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

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

	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;
 }

    

Fabien Freling

tags

participants (1)