Laboratoire RDP » Biophysics Team » lsm3d

echo 'compiling contour_init...'

$CCPP -o $bin_directory/contour_init contour_init.cpp $compiling_options

echo 'compiling lsm_contour_init...'

$CCPP -o $bin_directory/lsm_contour_init lsm_contour_init.cpp $compiling_options

echo 'compiling tiff2inr...'

$CCPP -o $bin_directory/tiff2inr tiff2inr.cpp $compiling_options -l:libtiff.so.5 

/*

Level-set Method to detect tissue contour (exterior shape)

- Sequential -

      Copyright 2016 ENS de Lyon

       File author(s):

           Typhaine Moreau, Annamaria Kiss <annamaria.kiss@ens-lyon.fr.fr>

       See accompanying file LICENSE.txt

To compile

 g++ -o contour_init contour_init.cpp -O2 -L/usr/X11R6/lib -lm -lpthread -lX11 -l:libtiff.so.5

 Need CImg.h and lsm_lib.h

To execute

 ./contour_init img t_up t_down smooth

 img : grayscale image of cells, .inr or .inr.gz

 t_up,t_down : linear threshold value (inr)

 smooth : amount of gaussian blur to apply to the image

*/

#include <iostream>

#include <math.h>

#include <sstream>

#include <vector>

#include <fstream>

#include "lsm_lib.h"

using namespace cimg_library;

using namespace std;

//------------------------------------------------------------------------------

//Main

//------------------------------------------------------------------------------

int main (int argc, char* argv[])

{

  clock_t begin=clock();

  if(argc!=5)

    {

      cout<<"!! wrong number of arguments"<<endl;

      cout<<"Usage : contour_init img t_up t_down smooth"<<endl;

      cout<<"Examples for parameter values:"<<endl;

      cout<<"------------------------------"<<endl;

      cout<<"img : grayscale image of cells, (.inr or .inr.gz)"<<endl;

      cout<<"Upper threshold : t_up = 20"<<endl;

      cout<<"Down threshold : t_down = 5"<<endl;

      cout<<"Gaussian filter : smooth = 1 (0, if image already filtered)"<<endl;

      return 0;

    }

  //ckeck filename and read image

  string filename=argv[1];

  float tailleVoxel[3] = {0};// resolution initialisation

  CImg<float> img=imread(filename, tailleVoxel); 

  //--------------------------------------------Parameters

  //linear threshold and smoothing

  int t_up=atoi(argv[2]);

  int t_down=atoi(argv[3]);

  float smooth=atof(argv[4]);

  //-------------------------------------------Names and directories

  //new name with arguments

  string ar2=argv[2];

  string ar3=argv[3];

  string ar4=argv[4];

  string insert="_initial"+ar2+"-"+ar3+"s"+ar4;

  filename.insert(filename.size()-4,insert);

  //-----------------------------------------Image Pre-processing

  //smooth image

  img.blur(smooth);

  //-------------------------------------------Initialization

  //initialize level-set

  CImg<unsigned char> segmented=threshold_linear_alongZ(img,t_up,t_down);

  int s=imsave(filename, tailleVoxel, segmented);

  return 0;

}

//Write an image image

//----------------------------------------------------------------------------

int imsave(string filename, float tailleVoxel[3], CImg<unsigned char> const & img)

{

	  img.save_inr(filename.c_str(),tailleVoxel);

	  string zip="gzip -f "+filename;

	  if(system(zip.c_str()));

  cout<<"Image saved in file :"<<filename<<".gz"<<endl;

  cout<<" - Voxel size : ("<<tailleVoxel[0]<<","<<tailleVoxel[1]<<","<<tailleVoxel[2]<<")"<<endl;

  cout<<" - Image size : ("<<img.width()<<","<<img.height()<<","<<img.depth()<<")"<<endl;

  return 1;

}

//Read an image image

//----------------------------------------------------------------------------

CImg<unsigned char> imread(string &filename, float (&tailleVoxel)[3])

{

	CImg<unsigned char> img_prev;

	cout<<"Image read : "<<filename<<endl;

    if(filename.compare(filename.size()-4,4,".inr")==0)

    {

      img_prev.get_load_inr(filename.c_str(),tailleVoxel); // reads resolution

    }

  else if(filename.compare(filename.size()-7,7,".inr.gz")==0)

    {

      string oldname = filename;

      filename.erase(filename.size()-3);

      string zip="gunzip -c "+oldname+" > "+filename;

      if(system(zip.c_str())); // decompress image file

      img_prev.load(filename.c_str()); //read image

      img_prev.get_load_inr(filename.c_str(),tailleVoxel); // read resolution

      zip="rm "+filename; 

      if(system(zip.c_str())); //removes decompressed image    

    }

  else

    {cout<<"!! wrong file extension : "<<filename<<endl;

	}

  cout<<" - Voxel size : ("<<tailleVoxel[0]<<","<<tailleVoxel[1]<<","<<tailleVoxel[2]<<")"<<endl;

  cout<<" - Image size : ("<<img_prev.width()<<","<<img_prev.height()<<","<<img_prev.depth()<<")"<<endl;

  return img_prev;

}

/*

Level-set Method to detect tissue contour (exterior shape)

- Sequential -

      Copyright 2016 ENS de Lyon

       File author(s):

           Typhaine Moreau, Annamaria Kiss <annamaria.kiss@ens-lyon.fr.fr>

       See accompanying file LICENSE.txt

To compile

 g++ -o lsm_contour_init lsm_contour_init.cpp -O2 -L/usr/X11R6/lib -lm -lpthread -lX11 -l:libtiff.so.5

 Need CImg.h and lsm_lib.h

To execute

 ./lsm_contour_init img initial_contour a b smooth perUp perDown

 img : grayscale image of cells (.inr or .inr.gz)

 initial_contour : 8bit image, with values 0=tissue, 1=background (.inr or .inr.gz)

 a : area term (float) --> 0.5, 1

 b : curvature term (float)

 smooth : amount of gaussian blur to apply to the image

 perUp, perDown : the algorithm stops when 10 succesive iteration are between perUp and perDown (in % of background growth)

*/

#include <iostream>

#include <math.h>

#include <sstream>

#include <vector>

#include <fstream>

#include "lsm_lib.h"

using namespace cimg_library;

using namespace std;

//------------------------------------------------------------------------------

//Main

//------------------------------------------------------------------------------

int main (int argc, char* argv[])

{

  clock_t begin=clock();

  if(argc!=8)

    {

      cout<<"!! wrong number of arguments"<<endl;

      cout<<"Usage : lsm_contour_init img initial_contour a b smooth perUp perDown"<<endl;

      cout<<"Examples for parameter values:"<<endl;

      cout<<"------------------------------"<<endl;

      cout<<"img : grayscale image of cells, (.inr or .inr.gz)"<<endl;

      cout<<"initial_contour : 8bit image, with values 0=tissue, 1=background (.inr or .inr.gz)"<<endl;

      cout<<"Area term : a = 0 (0.5, 1)"<<endl;

      cout<<"Curvature term : b = 0 (1)"<<endl;

      cout<<"Gaussian filter : smooth = 1 (0, if image already filtered)"<<endl;

      cout<<"Stop criteria : the contour evolution is in [perDown,perUp] for 10 consecutive iterations"<<endl;

      cout<<"     perUp = 0.002, perDown = -0.002"<<endl;

      return 0;

    }

  //check filename and read image

  //-----------------------------

  string filename=argv[1];

  float tailleVoxel[3] = {0};// resolution initialisation  

  CImg<float> img=imread(filename,tailleVoxel);

 //check filename and read the initial contour image

 //-----------------------------

  string filename_init=argv[2];

  float tailleVoxel_init[3] = {0};// resolution initialisation

  CImg<float> img_init=imread(filename_init,tailleVoxel_init);

  //--------------------------------------------Parameters

  //model parameters

  int lam=10;

  int alf=atoi(argv[3]);

  int beta=atoi(argv[4]);

  //numerical parameters

  float epsilon=1.5;

  int dt=100;

  float mu=0.1/dt;

  int timestep_max=2000;

  //smoothing and stop criteria

  float smooth=atof(argv[5]);

  float perUp=atof(argv[6]);

  float perDown=atof(argv[7]);

  cout<<"Voxel size : ("<<tailleVoxel[0]<<","<<tailleVoxel[1]<<","<<tailleVoxel[2]<<")"<<endl;

  cout<<"Voxel size in contour : ("<<tailleVoxel_init[0]<<","<<tailleVoxel_init[1]<<","<<tailleVoxel_init[2]<<")"<<endl;

  //-------------------------------------------Names and directories

  //new name with arguments

  string ar2=argv[2];

  string ar3=argv[3];

  string ar4=argv[4];

  string ar5=argv[5];

  string insert="_LSMcont_a"+ar3+"b"+ar4+"s"+ar5;

  filename.insert(filename.size()-4,insert);

  //create directories and update names

  size_t test=filename.rfind("/");

  if(test!=filename.npos)

    {filename.erase(0,test+1);}

  string outputdir=filename;

  outputdir.erase(filename.size()-4);

  string mkdir="mkdir -p "+outputdir;

  if(system(mkdir.c_str())); 

  string filename_txt=outputdir+"/"+filename;

  filename_txt.erase(filename_txt.size()-4);

  filename=outputdir+"/"+filename;

  string result_name=filename;

  //txt files 

  ofstream file;

  string txt_name=filename_txt+".txt";

  file.open(txt_name.c_str());

  file<<argv[0]<<endl;

  time_t t;

  struct tm * timeinfo;

  time(&t);

  timeinfo=localtime(&t);

  file<<asctime(timeinfo);

  file<<"image : "<<argv[1]<<endl;

  file<<"initial contour : "<<argv[2]<<endl;

  file<<"_________________________________"<<endl;

  file<<"Parameters"<<endl;

  file<<"lambda : "<<lam<<endl;

  file<<"alpha : "<<alf<<endl;

  file<<"epsilon : "<<epsilon<<endl;

  file<<"dt : "<<dt<<endl;

  file<<"mu : "<<mu<<endl;

  file<<"timestep_max : "<<timestep_max<<endl;

  file<<"beta : "<<beta<<endl;

  file<<"perUp : "<<perUp<<endl;

  file<<"perDown : "<<perDown<<endl;

  ofstream bg_file;

  string bg_name=filename_txt+"_BGgrowth.txt";

  bg_file.open(bg_name.c_str());

  bg_file<<"it\tbg_growth"<<endl;

  //-----------------------------------------Image Pre-processing

  //add slices

  img=add_side_slices(img,3);

  img_init=add_side_slices(img_init,3);

  CImg<unsigned char> segmented = img_init;

  //smooth image

  file<<"smooth : "<<smooth<<endl;

  img.blur(smooth);

  //-------------------------------------------Initialization

  //compute fixed terms

  CImg<float> g=edge_indicator(gradient(img));

  CImgList<float> gg=gradient(g);

  //initialize level-set

  int c0=-4;

  CImg<float> psi=lsm_contour_init(segmented,c0);

  int it=0;

  int it_stop=0;

  bool contour_evolves=true;

  int nb_pix=img.width()*img.height()*img.depth();

  double prev_backsegm=segmented.sum();

  //-------------------------------------------Time iterations

  while( (it<timestep_max) and (contour_evolves==true) )

    {

      //LSM

      psi=evolution_AK2_contour(psi,g,gg,g,lam,mu,alf,beta,epsilon,dt);

      //Update segmentation

      double backsegm=0;

      cimg_forXYZ(segmented,x,y,z)

	{

	  if(psi(x,y,z)>0)

	    {

	      segmented(x,y,z)=1;

	      backsegm+=1;

	    }

	  else

	    {segmented(x,y,z)=0;}

	}	

      //Background evolution

      double bg_evolution=backsegm-prev_backsegm;

      double bg100=(bg_evolution*1.0/nb_pix)*100;

      prev_backsegm=backsegm;

      cout<<"----------------------------------- it : "<<it<<endl;

      cout<<"bg growth : "<<bg_evolution<<endl;

      cout<<"% of bg growth : "<<bg100<<endl;

      bg_file<<it<<"\t"<<bg_evolution<<endl;

      //Stop criteria 

      if((it>10) and (bg100<perUp) and (bg100>perDown))

	{

	  it_stop+=1;

	  if(it_stop>9)

	    {contour_evolves=false;}

	}

      else

	{

	  it_stop=0;

	}

    //Save result

      if((((it%50)==0)and(it!=0))or(contour_evolves==false)or(it==timestep_max-1))

	{	  

	  CImg<unsigned char>segSave=remove_side_slices(segmented,3);

	  int sv=imsave(result_name, tailleVoxel, segSave);

	}

      it+=1;

    }

  clock_t end=clock();

  double time=double(end-begin)/CLOCKS_PER_SEC;

  cout <<"elapsed time : "<<time<<" sec ( ~ "<<time/60<<" mn ~ "<<time/60/60<<" h)"<<endl;

  file <<"last iteration : "<<it-1<<endl;

  file <<"elapsed time : "<<time<<" sec ( ~ "<<time/60<<" mn ~ "<<time/60/60<<" h)"<<endl;

  file<<"width "<<img.width()<<endl;

  file<<"height "<<img.height()<<endl;

  file<<"depth "<<img.depth()<<endl;

  file<<"number of pixel "<<img._width*img._height*img._depth<<endl;

  file.close();

  return 0;

}