Laboratoire RDP » Biophysics Team » lsm3d

/*

Level-set Method to detect contour (exterior shape)

Sequential

Stop criteria parameters in command lines

To compile

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

 Need CImg.h and lsm_lib.h

To execute

 ./lsm_contour img t_up t_down a b smooth perUp perDown

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

 t_up,t_down : linear threshold value (inr)

 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!=9)

    {

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

      cout<<"Usage : lsm_contour img t_up t_down 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<<"Upper threshold : t_up = 20"<<endl;

      cout<<"Down threshold : t_down = 5"<<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;

    }

  //ckeck filename and read image

  string filename=argv[1];

  CImg<unsigned char> img_prev;

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

  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;

      return 0;}

  CImg<float> img=img_prev;

  img_prev.assign();

  cout<<"original image : "<<filename<<endl;

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

  //model parameters

  int lam=10;

  int alf=atoi(argv[4]);

  int beta=atoi(argv[5]);

  //numerical parameters

  float epsilon=1.5;

  int dt=100;

  float mu=0.1/dt;

  int timestep_max=2000;

  //linear threshold

  int t_up=atoi(argv[2]);

  int t_down=atoi(argv[3]);

  float smooth=atof(argv[6]);

  float perUp=atof(argv[7]);

  float perDown=atof(argv[8]);

  cout<<"Voxel size : ("<<tailleVoxel[0]<<","<<tailleVoxel[1]<<","<<tailleVoxel[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 ar6=argv[6];

  string insert="_LSMcont"+ar2+"-"+ar3+"a"+ar4+"b"+ar5+"s"+ar6;

  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<<"_________________________________"<<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<<"\nthreshold up : "<<t_up<<endl;

  file<<"threshold down : "<<t_down<<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);

  //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<unsigned char> segmented=threshold_linear_alongZ(img,t_up,t_down);

  string segmentedName=filename+".gz";

  segmentedName.insert(filename.size()-4,"_initial");

  segmented.save_gzip_external(segmentedName.c_str());

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

          segSave.save_inr(result_name.c_str(),tailleVoxel);

          segSave.assign();

          string zip="gzip -f "+result_name;

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

        }

      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;

}