Laboratoire RDP » Biophysics Team » lsm3d

echo 'compiling edge_indicator...'

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

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

/*

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;

}

CImg<float> edge_indicator3exp(CImg<float> const & img, float gamma)

{

  CImg<float> res(img._width,img._height,img._depth,1,0);

  cimg_forXYZ(res,x,y,z)

    {

		res(x,y,z)=exp(-(img(x,y,z)/gamma)*(img(x,y,z)/gamma));

    }

  return res;

}