Laboratoire RDP » Biophysics Team » lsm3d

/*

Anisotropic blur : removes noise while keeping cell walls

Parallel version

3D only

To compile :

 g++ -o edge_indicator edge_indicator.cpp -O2 -L/usr/X11R6/lib -lm -lpthread -lX11 -fopenmp -l:libtiff.so.4

 Needs CImg.h

To execute :

 ./edge_indicator img K sigma [filtertype]

 img : image in .inr or .inr.gz

 K : proportional to width of cell walls (0.1 - 0.3)

 sigma : amount of gaussian blur (float)

 filtertype (optional) : by default 1/(1+x/K) type edge function. If "exp" is marked, then exp(-x/K) edge function is used

 * authors : Annamaria Kiss (RDP)

 * ENS-Lyon

 * 

*/

#include <iostream>

#include <math.h>

#include <sstream>

#include <fstream>

#include "CImg.h"

#include <omp.h>

using namespace cimg_library;

using namespace std;

#include "lsm_lib.h"

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

//Main

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

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

{

 // double begin_main=omp_get_wtime();

  if(argc<4 || argc>5)

    {

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

      cout<<"how to execute : edge_indicator img K sigma gamma upDiter nbIter [type]"<<endl;

      cout<<" img : image in .inr or .inr.gz"<<endl;

          cout<<" K : proportional to width of cell walls (0.1 - 0.3)"<<endl;

          cout<<" sigma : amount of gaussian blur (float)"<<endl;

          cout<<" filtertype (optional) : by default 'hill' that is 1/(1+x/K) type edge function. If 'exp' is marked, then exp(-x/K) edge function is used"<<endl; 

      return 0;

    }

  string filter_type="hill";

  if (argc==5)

        {

        if (string(argv[4]).compare("exp")==0)

                {filter_type="exp";}

        }

  //Open image with name in argument

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

  string name=argv[1];

  CImg<> img1;

  CImg<char> description;

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

  bool gzipped = false;

  bool tiffile = false;

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

    {   

      img1.load(name.c_str());

      img1.get_load_inr(name.c_str(),tailleVoxel); // reads resolution

    }

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

    {

      gzipped = true;

      string oldname = name;

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

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

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

      img1.load(name.c_str()); //read image

      img1.get_load_inr(name.c_str(),tailleVoxel); // read resolution

      zip="rm "+name; 

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

    }

  else if(name.compare(name.size()-4,4,".tif")==0)

        {

          tiffile = true;

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

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

        }

  else

    {cout<<"!! wrong file extension (not an inr/inr.gz/tif)  : "<<name<<endl;

      return 0;}

  CImg<float> img=img1;

  img1.assign();

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

  cout<<"gzipped : "<<gzipped<<endl;

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

  if(img.depth()<=1)

    {

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

      return  0;

    }

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

  float dt=0.1;

  float  K= atof(argv[2]);

  float sigma=atof(argv[3]);

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

  //--------------------------------------------------------Name and directories

  string insert=string("_edgeind-K")+argv[2]+"s"+argv[3]+"-type-"+filter_type;

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

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

  if(test!=name.npos)

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

  //---------------------------------------------------------------------Processing

  //Evolution during nbIter

      //-------------------------------------------------------------------------3D

      cout<<"3D image"<<endl;

      cout <<"Computing Hessian"<<endl;

      int nthreads;

      CImg<float> fx;

      CImg<float> edgeind;

      // make a blur before taking derivatives

          CImg<float> imgBlur=img.get_blur(sigma);

      fx=edge_lambda3(imgBlur);

    if (filter_type.compare("hill")==0)

                {

                        edgeind=edge_indicator3(img,K);

                }

                else

                {

                        edgeind=edge_indicator3exp(img,K);

                } 

        float maxi=edgeind.max();

        float mini=edgeind.min();

        edgeind=maxi-edgeind;

        maxi=edgeind.max();

        mini=edgeind.min();

        edgeind=edgeind*255./(maxi-mini);

    //edgeind=(edgeind-maxi)*255./(mini-maxi);

          //Saving results

              string newname=name;

              CImg<unsigned char> imgbis=edgeind;

              //CImg<float> imgbis=edgeind;

              cout<<"Save image"<<endl;

              cout<<"image : "<<newname<<endl;

              if (tiffile)

                        {

                                imgbis.save_tiff(newname.c_str(),0,tailleVoxel,description.data());

                        }

                  else

                        {

                          imgbis.save_inr(newname.c_str(),tailleVoxel);

                          if (gzipped) 

                          {

                                string zip="gzip -f "+newname;

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

                                }

                  }

              imgbis.assign();

  cout << "Done!"<<endl;

  return 0;

}