Laboratoire RDP » Biophysics Team » lsm3d

/*

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

Sequential standard version

* background : 0; object : 255

To compile

 g++ -o lsm_detect_contour-mgx-tiffonly.exe lsm_detect_contour-mgx-tiffonly.cpp -O2 -L/usr/X11R6/lib -lm -lpthread -lX11 -l:libtiff.so.4

 Needs CImg.h and lsm_lib.h

To execute

 ./lsm_detect_contour-mgx.exe img t_up t_down b smooth

 img : grayscale image of cells, (.tif)

 t_up,t_down : linear threshold value (int)

 b : curvature term (float)

 smooth : amount of gaussian blur to apply to the image

 * implemented by Typhaine Moreau during her M1 internship (April-July 2014)

 * supervisors Annamaria Kiss (RDP) and Cerasela Calugaru (CBP)

 * ENS-Lyon

 * 

 * Reference : Li et al. Distance Regularized Level Set Evolution and Its Application to Image Segmentation, 

 * IEEE Transactions on Image Processing, Vol.19. No. 12, December 2010

*/

#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!=6)

    {

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

      cout<<"how to execute : ./lsm_detect_contour.exe img t_up t_down b smooth"<<endl;

      return 0;

    }

  //Open image with name in argument

  //image is convert to unsigned char (0 to 255) and then to float

  string name=argv[1];

  //CImg<unsigned char> img_prev;

  CImg<> img1;

  CImg<char> description;

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

          cout<<"file to read : "<<name<<endl;

          img1.load_tiff(name.c_str(),0,~0U,1,tailleVoxel,&description);

  CImg<float> img=img1;

  img1.assign();

  cout<<"image : "<<name<<endl;

  cout<<"image depth : "<<img.depth()<<endl;

  if(img.depth()<=1)

    {

      cout<<"This is not a 3D image"<<endl;

      return  0;

    }

  string filename = name;

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

  //model parameters

  int lam=10;

  int alf=0;

  int beta=atoi(argv[4]);

  //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[5]);

  //-------------------------------------------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"+ar2+"-"+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())); 

  mkdir="mkdir -p "+outputdir+"/evolution";

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

  string filename_cut=outputdir+"/evolution/"+filename;

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

  filename=outputdir+"/"+filename;

  string result_name=filename;

  //txt files 

  ofstream file;

  string txt_name=filename_cut+".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;

  ofstream bg_file;

  string bg_name=filename_cut+"_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);

  //define and save cut

  int cut=img._depth/2;

  cout <<"cut z= "<<cut<<endl;

  file <<"\ncut z= "<<cut<<endl;

  CImg<float> imgCut=img.get_crop(0,0,cut,0,img._width,img._height,cut,0);

  imgCut.save((filename_cut+".png").c_str());

  //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<0.002) and (bg100)>-0.002)

        {

          it_stop+=1;

          if(it_stop>9)

            {contour_evolves=false;}

        }

      else

        {

          it_stop=0;

        }

      //Graphics to follow evolution

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

        {

          ostringstream oss;

          oss << it;

          string iterstring="";

          if(it<10) {iterstring="00"+oss.str();}

          else if (it<100) {iterstring="0"+oss.str();}

          else {iterstring=oss.str();}

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

          cout<<" *** time step for graphics : "<<iterstring<<endl;

          //Save cut

          CImg<float>segCut=segmented.get_crop(0,0,cut,0,img._width,img._height,cut,0);

          string temp_name=filename_cut+"it"+iterstring+".png";

          segCut.normalize(0,255).save(temp_name.c_str());

        }

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

                {

                  //Save result

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

                  segSave = invert_image(segSave)*255;

                        segSave.save_tiff(result_name.c_str(),0,tailleVoxel,description.data());

                  segSave.assign();

                }

      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;

}