Laboratoire RDP » Biophysics Team » lsm3d

/*

Level-set Method to detect cell's contours

Parallel standard version

To compile :

 g++ -o lsm_cells lsm_cells.cpp -O2 -L/usr/X11R6/lib -lm -lpthread -lX11 -fopenmp

 Need CImg.h and lsm_lib.h

To execute :

 ./lsm_cells img img_wat img_contour erosion

 image in .inr or .inr.gz, save in .inr.gz

 img : grayscale image of cells (unsigned char)

 img_wat : watershed 16bits image (unsigned short) with background label 1 and cells labels >1

 img_contour : binary image with background=1, all cells=0 (unsigned char)

 erosion : amount of erosion for each cell in watershed image (int)

*/

#include <iostream>

#include <math.h>

#include <sstream>

#include <vector>

#include <fstream>

#include <omp.h>

#include "CImg.h"

#include "lsm_lib.h"

using namespace cimg_library;

using namespace std;

//Return list of cell label

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

vector<int> index(const CImg<unsigned short> & wat, int nbcells)

{

  vector<int> list;

  vector<bool> test(nbcells,false);

  cimg_forXYZ(wat,x,y,z)

    {

      int ind=wat(x,y,z);

      if((test[ind]==false)and(ind!=1)and(ind!=0))

        {

          list.push_back(ind);

          test[ind]=true;

        }

    }

  return list;

}

//Return crop psi image and it's min coordinates for one cell 

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

CImg<float> box_cell(const CImg<unsigned short> & wat, int marge,int erosion, int c0,int indice, int &xmin, int &ymin, int &zmin)

{

  //search box (min and max coordinates with a marge)

  xmin=wat._width;

  ymin=wat._height;

  zmin=wat._depth;

  int xmax=0;

  int ymax=0;

  int zmax=0;

  cimg_forXYZ(wat,x,y,z)

    {

      if(wat(x,y,z)==indice)

        {

          if(x>xmax){xmax=x;}else if(x<xmin){xmin=x;}

          if(y>ymax){ymax=y;}else if(y<ymin){ymin=y;}

          if(z>zmax){zmax=z;}else if(z<zmin){zmin=z;}

        }

    }

  //marge and border conditions

  xmin=xmin-marge;

  if(xmin<0){xmin=0;}

  ymin=ymin-marge;

  if(ymin<0){ymin=0;}

  zmin=zmin-marge;

  if(zmin<0){zmin=0;}

  xmax=xmax+marge;

  if(xmax>=wat._width){xmax=wat._width-1;}

  ymax=ymax+marge;

  if(ymax>=wat._height){ymax=wat._height-1;}

  zmax=zmax+marge;

  if(zmax>=wat._depth){zmax=wat._depth-1;}

  //crop wat image to the size of box, make binary

  CImg<unsigned short>binary=wat.get_crop(xmin,ymin,zmin,0,xmax,ymax,zmax,0);

  cimg_forXYZ(binary,x,y,z)

    {

      if(binary(x,y,z)==indice){binary(x,y,z)=1;}

      else {binary(x,y,z)=0;}

    }

  //erode binary but not completely (vol stay >0)

  int vol=binary.sum();

  int nb_ero=0;

  while((vol>0)and(nb_ero<abs(erosion)))

  {

    CImg<unsigned char> binary_erode=binary.get_erode(3,3,3); 

        if (erosion<0) 

        {

                for(int i=0;i<3;i++)

                {

                binary=binary.get_dilate(2,2,2);

                }

                binary_erode=binary;

        }

      vol=binary_erode.sum();

      if(vol>0)

        {

          binary=binary_erode;

          nb_ero+=1;

        }

    }

  //initalize psi

  CImg<float>psi=binary;

  cimg_forXYZ(psi,x,y,z)

    {

      if(binary(x,y,z)==0){psi(x,y,z)=c0;}

      else {psi(x,y,z)=-c0;}

    }

  return psi;

}

//LSM segment : edge detection

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

CImg<float> lsm_segment2(CImg<float> psi,CImg<float> const & g,CImgList<float> const & gg,CImg<float> const & h, int lam, float mu,float alf, float beta, float epsilon, int dt, vector<int> min_list)

{

  int timestep_max=2000;

  int evolution_min=1;

  //int evolution_max=-1;

  int it=0;

  int it_stop=0;

  bool contour_evolve=true;

  int backsegm=0;

  CImg<float> psi_old=psi;

  cimg_forXYZ(psi,x,y,z)

    {

      if(psi(x,y,z)>=0){backsegm+=1;}

    }

  while((it<timestep_max)and(contour_evolve==true)and(backsegm>30))

    {

      psi_old=psi;

      //evolution

      psi=evolution_AK2_contour_cells(psi,g,gg,h,lam,mu,alf,beta,epsilon,dt,min_list);

      //new segmentation

      int new_backsegm=0;

      cimg_forXYZ(psi,x,y,z)

        {

          if(psi(x,y,z)>=0){new_backsegm+=1;}

        }

      //Stop criteria

      // * if the cell would disappear, we stop

      if(new_backsegm==0)

                        {psi=psi_old;

                     contour_evolve=false;}

      int bg_evolution=abs(new_backsegm-backsegm);

      // * if the evolution is less then evolution_min 3 consecutive times

      //if((it>10) and (bg_evolution<=evolution_min) and (bg_evolution>=evolution_max))

      if((it>10) and (bg_evolution<=evolution_min) )

        {

          it_stop+=1;

          if(it_stop>3)

            {contour_evolve=false;

                cout<<"stopit "<<endl;}

        }

      else

        {

          it_stop=0;

        }

      it+=1;

      backsegm=new_backsegm;

      //cout<<"backsegm= "<<backsegm<<endl;     

    }

//  cimg_forXYZ(psi,x,y,z)

//    {

//      if(psi(x,y,z)>=0){psi(x,y,z)=4;}

//      else{psi(x,y,z)=-4;}

//    }

  return psi;

}

//Reconstruct full segmented image from cell's images

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

CImg<unsigned short> reconstruct(const CImg<unsigned char> & background,const vector<CImg<float> > & psi_list, const vector< vector<int> > & min_list, int nbcells,const vector<int> & list)

{

  CImg<unsigned short> res=background;

  for(int i=0;i<nbcells;i++)

    {

      cimg_forXYZ(psi_list[i],xb,yb,zb)

        {

          int x=xb+min_list[i][0];

          int y=yb+min_list[i][1];

          int z=zb+min_list[i][2];

          if(psi_list[i](xb,yb,zb)>=0)

            {res(x,y,z)=list[i];}

        }

    }

  return res;

}

//Reconstruct segmented image from cell's images and treat overlap

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

CImg<unsigned short> reconstruct_overlap(CImg<unsigned char> const &  background,const vector<CImg<float> > & psi_list, const vector< vector<int> > & min_list, int nbcells, CImg<unsigned char> &free, const vector<int> & list)

{

  CImg<unsigned short> res=background;

  free=background;

  for(int i=0;i<nbcells;i++)

    {

      cimg_forXYZ(psi_list[i],xb,yb,zb)

        {

          int x=xb+min_list[i][0];

          int y=yb+min_list[i][1];

          int z=zb+min_list[i][2];

          if(psi_list[i](xb,yb,zb)>=0)

            {

              res(x,y,z)=list[i];

              free(x,y,z)+=1;

            }

        }

    }

  cimg_forXYZ(free,x,y,z)

    {

      if(free(x,y,z)>1)

        {

          if(background(x,y,z)==1)

            {

              free(x,y,z)=1;

              res(x,y,z)=1;

            }

          else

            {

              free(x,y,z)=0;

              res(x,y,z)=0;

            }

        }

    }

  return res;

}

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

//******************************************  MAIN **********************************************

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

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

{

  double begin=omp_get_wtime();

  if(argc<5)

    {

    cout<<"!! wrong number of arguments ("<<argc<<")"<<endl;

    cout<<"Usage : lsm_cells img img_wat img_contour erosion [a b smooth lsm_type]"<<endl;

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

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

        cout<<"img_wat : image of seeds, (.inr or .inr.gz)"<<endl;

        cout<<"img_contour : mask, where cells do not evolve, (.inr or .inr.gz)"<<endl;

        cout<<"              if 'None', then cells can evolve on the whole image"<<endl;

        cout<<"erosion : amount of erosion of seeds for initialisation (uint8) --> -2, 0, 2"<<endl;

        cout<<"              if 0, then no erosion or dilation"<<endl;

        cout<<"              if negative, then a dilation is performed"<<endl;

        cout<<"a : area term (float) --> 0 or 0.5 or 1 (the default is 0.5)"<<endl;

        cout<<"              if negative, the object retracts"<<endl;

        cout<<"              if positive, the object inflates"<<endl;

        cout<<"b : curvature term (float) --> 0 or 1 (the default is 0)"<<endl;

        cout<<"gamma : scale parameter (float>0) --> 0.5 or 1 (the default is 1)"<<endl;

        cout<<"smooth : gaussian blur to apply to the image (int) --> 0 or 1 (the default is 0)"<<endl;

    cout<<"lsm_type : image, gradient or hessien based evolution --> 'i', 'g' or 'h' (the default is g)"<<endl;  

      return 0;

    }

  //----------------------------------------------read images and check the names

  //Original image

  CImg<char> description;

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

  //float tailleVoxel[3] = {0.195177,0.195177,0.195177};

  bool gzipped=false;

  string filename_img=argv[1];

  CImg<unsigned char> img_prev;

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

    {

      img_prev.load(filename_img.c_str());

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

    }

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

    {

          gzipped = true;

      string oldname = filename_img;

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

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

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

      img_prev.load(filename_img.c_str());

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

      zip="rm "+filename_img;

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

    }

  else

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

      return 0;}

  CImg<float> img=img_prev;

  img_prev.assign();

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

  cout<<"size : "<<img.size()<<endl;

  cout<<"size of original image : "<< img._width<<' '<< img._height <<' '<< img._depth<<' '<< endl;

  //Watershed

  string filename=argv[2];

  CImg<unsigned short> wat;

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

    {

      wat.load(filename.c_str());

    }

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

    {

      wat.load_gzip_external(filename.c_str());

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

    }

  else

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

      return 0;}

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

cout<<"size : "<<wat.size()<<endl;

  cout<<"size of original image : "<< wat._width<<' '<< wat._height <<' '<< wat._depth<<' '<< endl;

  //Background

  string filename_bg=argv[3];

  CImg<unsigned char> background(img._width, img._height, img._depth,1,0);

   if(filename_bg.compare("None")==0)

    {

       cout<<"!! no background entered !!"<<endl;

    }

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

    {

      background.load(filename_bg.c_str());

    }

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

    {

      background.load_gzip_external(filename_bg.c_str());

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

    }

  else

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

      return 0;}

  cout<<"background image : "<<filename_bg<<endl;

  if((wat.size()!=img.size())or(background.size()!=img.size()))

    {

      cout<<"!! images are not the same size"<<endl;

      return 0;

    }

  else cout<<"size : "<<img.size()<<endl;

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

  //model parameters

  int c0=-4;

  int erosion=atoi(argv[4]);

  int marge=10; //10

  int lam=10;

  float alf=0.5;//1

  float beta=0;

  float gamma=1;

  float smooth=0;

  string lsm_type="g";

  string ar=argv[4];

  string insert="_cellLSM-d"+ar;

  if(argc>5)

    {

    alf=atof(argv[5]);

    ar=argv[5];

    insert=insert+"-a"+ar;

    }

  if(argc>6)

    {

    beta=atof(argv[6]);

    ar=argv[6];

    insert+="-b"+ar;

    }

  if(argc>7)

    {

    gamma=atof(argv[7]);

    ar=argv[7];

    insert+="-g"+ar;

        }

  if(argc>8)

    {

    smooth=atof(argv[8]);

    ar=argv[8];

    insert+="-s"+ar;

        }

  if(argc>9)

    {

        lsm_type=argv[9];

        ar=argv[9];

        insert+="-"+ar;

        }

  //numerical parameters

  float epsilon=0.5;

  int dt=1;   //it was 50, changed into 1 the 2015/10/16

  float mu=0.1/dt;

  int timestep_max=10000; // it was 1000, changed into 10000 the 2015/10/16

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

  //new name with arguments

  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;

  int syst=system(mkdir.c_str()); 

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

  syst=system(mkdir.c_str());

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

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

  filename=outputdir+"/"+filename;

  string wat_eroded_name=filename;

  wat_eroded_name.insert(filename.size()-4,"_eroded");

  string edge_detection_name=filename;

  edge_detection_name.insert(filename.size()-4,"_evoEdge");

  string edge_evolve_name=filename;

  edge_evolve_name.insert(filename.size()-4,"_final");

  //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<<"watershed : "<<argv[2]<<endl;

  file<<"background : "<<argv[3]<<endl;

  file<<"_________________________________"<<endl;

  file<<"Parameters"<<endl;

  file<<"lsm_type : "<<lsm_type<<endl;

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

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

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

  file<<"gamma :"<<gamma<<endl;

  //file<<"alphabis : "<<alfabis<<endl;

  //file<<"betabis : "<<betabis<<endl;

  file<<"erosion : "<<erosion<<endl;

  file<<"marge : "<<marge<<endl;

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

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

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

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

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

  //define cut

  int cut=img._depth/2;

  cut = 82;

  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;

  if (smooth>0)

  {

  img.blur(smooth);

  }

  //-------------------------------------Initialization with erosion

  //compute fixed terms

  CImg<float> g;

  if(lsm_type.compare("g")==0)

    {

  g=edge_indicator1(img, gamma);

    }

  else if (lsm_type.compare("h")==0)

   {

  g=edge_indicator2s(img, gamma);

   }

  else if (lsm_type.compare("i")==0)

   {

  g=edge_indicator3(img, gamma);

   }

  else

  cout<<"Wrong lsm type given :'"<<lsm_type<<"'('g' for gradient-based, 'h' for Hessien-based) "<<endl;

  CImgList<float> gg=gradient(g);

  CImg<float> h=g;

  img.assign();

  CImg<float> gout=g;

  gout.crop(0,0,cut,0,g._width,g._height,cut,0);

  gout.normalize(0,255).save((filename_cut+"_edge_indicator.png").c_str());

  gout.assign();

  //initialize psi for every cell

  int maxcells=wat.max()+1; //indice maximum  

  vector<int> list=index(wat,maxcells);

  int nbcells=list.size();

  cout<<"number of cells : "<<nbcells<<endl;

  file<<"number of cells : "<<nbcells<<endl;

  vector<CImg<float> > psi_list(nbcells);

  vector<vector<int> > min_list(nbcells,vector<int>(3,0));

  vector<vector<int> > max_list(nbcells,vector<int>(3,0));

#pragma omp parallel shared(wat,marge,erosion,c0,psi_list,min_list,list)

  {

  int xmin=0;

  int ymin=0;

  int zmin=0;

#pragma omp for schedule(dynamic)

  for(int i=0;i<nbcells;i++) 

    {

      int ind=list[i];

      psi_list[i]=box_cell(wat,marge,erosion,c0,ind,xmin,ymin,zmin);

      min_list[i][0]=xmin;

      min_list[i][1]=ymin;

      min_list[i][2]=zmin;

      cout <<"cell "<<ind<<endl;

    }

  }

  wat.assign();

  //reconstruct image of eroded cell

  CImg<unsigned short> wat_eroded=reconstruct(background,psi_list,min_list,nbcells,list);

  wat_eroded.save_inr(wat_eroded_name.c_str(),tailleVoxel);

  string zip="gzip -f "+wat_eroded_name;

  syst=system(zip.c_str());

  wat_eroded.crop(0,0,cut,0,wat_eroded._width,wat_eroded._height,cut,0);

  wat_eroded.normalize(0,255).save((filename_cut+"_eroded.png").c_str());

  wat_eroded.assign();

  //Segmented inital = background segmentation

  CImg<unsigned char> segmented=background;

  double end1=omp_get_wtime();

  double time1=double(end1-begin);

  cout<<"initialization with erosion : "<<time1<<endl;

  file<<"\ninitialization with erosion : "<<time1<<endl;

  cout<<"Evolving cells..... "<<endl;

  //---------------------------------------------------------Edge detection

  //evolve each cell one by one, attract to maximal gradient

#pragma omp parallel shared(psi_list,min_list,g,gg,h,lam,mu,alf,beta,epsilon,dt,list)

  {

#pragma omp for schedule(dynamic)

  for(int i=0;i<nbcells;i++) 

    {

      psi_list[i]=lsm_segment2(psi_list[i],g,gg,h,lam,mu,alf,beta,epsilon,dt,min_list[i]);

      cout <<"cell evolution "<<list[i]<<endl;

    }

  }

  //reconstruct image of edge detection, overlap=not segmented

  CImg<unsigned char>free=background;

  CImg<unsigned short> edge=reconstruct_overlap(background,psi_list,min_list,nbcells,free,list);

  edge.save_inr(edge_detection_name.c_str(),tailleVoxel);

  zip="gzip -f "+edge_detection_name;

  syst=system(zip.c_str());

  edge.crop(0,0,cut,0,edge._width,edge._height,cut,0);

  edge.normalize(0,255).save((filename_cut+"_evoEdge.png").c_str());

  double end2=omp_get_wtime();

  double time2=double(end2-begin);

  cout<<"edge detection : "<<time2<<endl;

  file<<"edge detection : "<<time2<<endl;

  double end=omp_get_wtime();

  double time=double(end-begin);

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

  cout<<"initialization with erosion : "<<time1<<endl;

  cout<<"edge detection : "<<time2<<endl;

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

  file<<"initialization with erosion : "<<time1<<endl;

  file<<"edge detection : "<<time2<<endl;

  file.close();

  return 0;

}