Laboratoire RDP » Biophysics Team » lsm3d

/*     

       Copyright 2016 ENS de Lyon

       File author(s):

           Typhaine Moreau, Annamaria Kiss <annamaria.kiss@ens-lyon.fr.fr>

       See accompanying file LICENSE.txt

*/

#include <iostream>

#include <math.h>

#include <sstream>

#include <vector>

#include <fstream>

#include <algorithm>

#include "CImg.h"

using namespace cimg_library;

using namespace std;

//**************************************** GLOBAL IMAGE *************************************************************************

//Write an image image

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

int imsave(string filename, float tailleVoxel[3], CImg<unsigned char> const & img)

{

          img.save_inr(filename.c_str(),tailleVoxel);

          string zip="gzip -f "+filename;

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

  cout<<"Image saved in file :"<<filename<<".gz"<<endl;

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

  cout<<" - Image size : ("<<img.width()<<","<<img.height()<<","<<img.depth()<<")"<<endl;

  return 1;

}

//Read an image image

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

CImg<unsigned char> imread(string &filename, float (&tailleVoxel)[3])

{

        CImg<unsigned char> img_prev;

        cout<<"Image read : "<<filename<<endl;

    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;

        }

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

  cout<<" - Image size : ("<<img_prev.width()<<","<<img_prev.height()<<","<<img_prev.depth()<<")"<<endl;

  return img_prev;

}

//Invert the image

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

CImg<unsigned char> invert_image(CImg<unsigned char> const & img)

{

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

  CImg<unsigned char> inverted = unity - img;

  return inverted;

}

//Add black slices on each side of the image

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

CImg<float> add_side_slices(CImg<float> const & img, int side)

{

  int s=side*2;

  CImg<float> newImg(img._width+s,img._height+s,img._depth+s,1,1);

  cimg_forXYZ(img,x,y,z)

    {

      newImg(x+side,y+side,z+side)=img(x,y,z);   

    }

  return newImg;

}

//Remove slices

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

CImg<float> remove_side_slices(CImg<float> const & img, int side)

{

  int s=side*2;

  CImg<float> newImg(img._width-s,img._height-s,img._depth-s,1,0);

  cimg_forXYZ(newImg,x,y,z)

    {

      newImg(x,y,z)=img(x+side,y+side,z+side);   

    }

  return newImg;

}

//Threshold linear along Z   1:background  0:cells

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

CImg<unsigned char> threshold_linear_alongZ(CImg<float> const & img, int t1, int t2)

{

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

  for(int z=0; z<img._depth ; z++)

    {

      cimg_forXY(outside,x,y)

        {

          float thresh = t1 + (t2-t1)*(z*1.0/img._depth);

          if(img(x,y,z)<thresh)

            {outside(x,y,z)=1;}        

        }

    }

  outside.label();

  cimg_forXYZ(outside,x,y,z)

    {

      if(outside(x,y,z)>0) {outside(x,y,z)=0;}

      else {outside(x,y,z)=1;}

    }

  return outside;

}

//LSM contour init +c0 background, -c0 cells

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

CImg<float> lsm_contour_init(CImg<unsigned char> const & region,int c0)

{

  CImg<float> initLSM(region._width,region._height,region._depth,1,c0);

  initLSM=initLSM-2*c0*region;

  return initLSM;

}

//Gradient 3D with central finite differences

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

CImgList<float> gradient(CImg<float> const & img)

{

  //List of 3 images the same size as img

  CImgList<float> grad(3,img._width,img._height,img._depth,1,0);

  //Image loop with a 3*3*3 neighborhood

  CImg_3x3x3(I,float);

  cimg_for3x3x3(img,x,y,z,0,I,float)

    {

      grad(0,x,y,z) = (Incc - Ipcc)/2; //grad x = (img(x+1,y,z)-img(x-1,y,z))/2

      grad(1,x,y,z) = (Icnc - Icpc)/2; //grad y

      grad(2,x,y,z) = (Iccn - Iccp)/2; //grad z

    }

  return grad;

}

//Norm of the gradient

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

CImg<float> gradient_norm(CImg<float> const & img)

{

  CImgList<float> grad = gradient(img);

  CImg<float> res(grad[0]._width,grad[0]._height,grad[0]._depth,1,0);

  float f=0;

  cimg_forXYZ(res,x,y,z)

    {

      f=grad(0,x,y,z)*grad(0,x,y,z) + grad(1,x,y,z)*grad(1,x,y,z) + grad(2,x,y,z)*grad(2,x,y,z);

      res(x,y,z)=sqrt(f);

    }

  return res;

}

//Edge indicator

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

CImg<float> edge_indicator(CImgList<float> const & grad)

{

  CImg<float> res(grad[0]._width,grad[0]._height,grad[0]._depth,1,0);

  float f=0;

  cimg_forXYZ(res,x,y,z)

    {

      f=grad(0,x,y,z)*grad(0,x,y,z) + grad(1,x,y,z)*grad(1,x,y,z) + grad(2,x,y,z)*grad(2,x,y,z);

      res(x,y,z)=1.0/(1.0+f);

      //res(x,y,z)=exp(-f);

    }

  return res;

}

//Edge indicator 1 (gradient)

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

CImg<float> edge_indicator1(CImg<float> const & img)

{

  CImgList<float> grad;

  grad = gradient(img);

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

  float f=0;

  cimg_forXYZ(res,x,y,z)

    {

      f=grad(0,x,y,z)*grad(0,x,y,z) + grad(1,x,y,z)*grad(1,x,y,z) + grad(2,x,y,z)*grad(2,x,y,z);

      res(x,y,z)=1.0/(1.0+f);

    }

  return res;

}

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

{

  CImgList<float> grad;

  grad = gradient(img);

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

  float f=0;

  cimg_forXYZ(res,x,y,z)

    {

      f=grad(0,x,y,z)*grad(0,x,y,z) + grad(1,x,y,z)*grad(1,x,y,z) + grad(2,x,y,z)*grad(2,x,y,z);

      res(x,y,z)=1.0/(1.0+f/gamma/gamma);

    }

  return res;

}

//Edge indicator 2 (hessienne, computed without smoothing)

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

CImg<float> edge_indicator2(CImg<float> const & img)

{

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

  float f=0;

  cimg_forXYZ(res,x,y,z)

    {

    int xp=x-1;

        if (x==0){xp=x;}

        int xn=x+1;

        if (x==img._width-1){xn=x;}

        int yp=y-1;

        if (y==0){yp=y;}

        int yn=y+1;

        if (y==img._height-1){yn=y;}

        int zp=z-1;

        if (z==0){zp=z;}

        int zn=z+1;

        if (z==img._depth-1){zn=z;}

        CImg<float> hess(3,3,1,1,0);

        hess(0,0) = img(xp,y,z) + img(xn,y,z) - 2*img(x,y,z);          // Hxx

        hess(1,1) = img(x,yp,z) + img(x,yn,z) - 2*img(x,y,z);          // Hyy

        hess(2,2) = img(x,y,zp) + img(x,y,zn) - 2*img(x,y,z);          // Hzz

        hess(1,0) = hess(0,1) = (img(xp,yp,z) + img(xn,yn,z) - img(xp,yn,z) - img(xn,yp,z))/4; // Hxy = Hyx

        hess(2,0) = hess(0,2) = (img(xp,y,zp) + img(xn,y,zn) - img(xp,y,zn) - img(xn,y,zp))/4; // Hxz = Hzx

        hess(2,1) = hess(1,2) = (img(x,yp,zp) + img(x,yn,zn) - img(x,yp,zn) - img(x,yn,zp))/4; // Hyz = Hzy

        //val : 3 eigen values as an array, decreasing order

        //vec : 3 eigen vectors as a 3*3 image

        //have to be double or create NaN

        CImg<double> val,vec;

        hess.symmetric_eigen(val,vec);

        if (val[2]>0) val[2]=0;

        val[2]=-val[2];

        //res(x,y,z)=1.0/(1.0 +val[2]*(50/3.));  //h=g test4

    //res(x,y,z)=1.0/(1.0 +val[2]*val[2]*(50/3.)*(50/3.) );

    //res(x,y,z)=1.0/(1.0 +exp(-val[2]));

    //res(x,y,z)=exp(-val[2]);

    //res(x,y,z)=exp(-val[2]*val[2]);

    res(x,y,z)=1.0/(1.0 +val[2]*val[2]);

    }

  return res;

}

//Edge indicator 2 (hessienne, computed with smoothing)

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

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

{

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

  // make a minimal blur before taking derivatives

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

  float f=0;

  cimg_forXYZ(res,x,y,z)

    {

    int xp=x-1;

        if (x==0){xp=x;}

        int xn=x+1;

        if (x==img._width-1){xn=x;}

        int yp=y-1;

        if (y==0){yp=y;}

        int yn=y+1;

        if (y==img._height-1){yn=y;}

        int zp=z-1;

        if (z==0){zp=z;}

        int zn=z+1;

        if (z==img._depth-1){zn=z;}

        CImg<float> hess(3,3,1,1,0);

        hess(0,0) = imgBlur(xp,y,z) + imgBlur(xn,y,z) - 2*imgBlur(x,y,z);          // Hxx

        hess(1,1) = imgBlur(x,yp,z) + imgBlur(x,yn,z) - 2*imgBlur(x,y,z);          // Hyy

        hess(2,2) = imgBlur(x,y,zp) + imgBlur(x,y,zn) - 2*imgBlur(x,y,z);          // Hzz

        hess(1,0) = hess(0,1) = (imgBlur(xp,yp,z) + imgBlur(xn,yn,z) - imgBlur(xp,yn,z) - imgBlur(xn,yp,z))/4; // Hxy = Hyx

        hess(2,0) = hess(0,2) = (imgBlur(xp,y,zp) + imgBlur(xn,y,zn) - imgBlur(xp,y,zn) - imgBlur(xn,y,zp))/4; // Hxz = Hzx

        hess(2,1) = hess(1,2) = (imgBlur(x,yp,zp) + imgBlur(x,yn,zn) - imgBlur(x,yp,zn) - imgBlur(x,yn,zp))/4; // Hyz = Hzy

        //val : 3 eigen values as an array, decreasing order

        //vec : 3 eigen vectors as a 3*3 image

        //have to be double or create NaN

        CImg<double> val,vec;

        hess.symmetric_eigen(val,vec);

        if (val[2]>0) val[2]=0;

        val[2]=-val[2];

    res(x,y,z)=1.0/(1.0 +val[2]*val[2]/gamma/gamma);

    }

  return res;

}

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

{

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

  float f=0;

  cimg_forXYZ(res,x,y,z)

    {

    int xp=x-1;

        if (x==0){xp=x;}

        int xn=x+1;

        if (x==img._width-1){xn=x;}

        int yp=y-1;

        if (y==0){yp=y;}

        int yn=y+1;

        if (y==img._height-1){yn=y;}

        int zp=z-1;

        if (z==0){zp=z;}

        int zn=z+1;

        if (z==img._depth-1){zn=z;}

        CImg<float> hess(3,3,1,1,0);

        hess(0,0) = img(xp,y,z) + img(xn,y,z) - 2*img(x,y,z);          // Hxx

        hess(1,1) = img(x,yp,z) + img(x,yn,z) - 2*img(x,y,z);          // Hyy

        hess(2,2) = img(x,y,zp) + img(x,y,zn) - 2*img(x,y,z);          // Hzz

        hess(1,0) = hess(0,1) = (img(xp,yp,z) + img(xn,yn,z) - img(xp,yn,z) - img(xn,yp,z))/4; // Hxy = Hyx

        hess(2,0) = hess(0,2) = (img(xp,y,zp) + img(xn,y,zn) - img(xp,y,zn) - img(xn,y,zp))/4; // Hxz = Hzx

        hess(2,1) = hess(1,2) = (img(x,yp,zp) + img(x,yn,zn) - img(x,yp,zn) - img(x,yn,zp))/4; // Hyz = Hzy

        //val : 3 eigen values as an array, decreasing order

        //vec : 3 eigen vectors as a 3*3 image

        //have to be double or create NaN

        CImg<double> val,vec;

        hess.symmetric_eigen(val,vec);

        if (val[2]>0) val[2]=0;

        //val[2]=-val[2];

    res(x,y,z)=1.0/(1.0 +val[2]*val[2]/gamma/gamma);

    }

  return res;

}

//Edge indicator 3 (image based)

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

CImg<float> edge_indicator3(CImg<float> const & img)

{

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

  cimg_forXYZ(res,x,y,z)

    {

      res(x,y,z)=1.0/(1.0+img(x,y,z)*img(x,y,z));

    }

  return res;

}

CImg<float> edge_indicator3(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)=1.0/(1.0+(img(x,y,z)/gamma)*(img(x,y,z)/gamma));

    }

  return res;

}

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;

}

//Edge image (based on lambda3 of the hessien)

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

CImg<float> edge_lambda3(CImg<float> const & img)

{

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

  float f=0;

  cimg_forXYZ(res,x,y,z)

    {

    int xp=x-1;

        if (x==0){xp=x;}

        int xn=x+1;

        if (x==img._width-1){xn=x;}

        int yp=y-1;

        if (y==0){yp=y;}

        int yn=y+1;

        if (y==img._height-1){yn=y;}

        int zp=z-1;

        if (z==0){zp=z;}

        int zn=z+1;

        if (z==img._depth-1){zn=z;}

        CImg<float> hess(3,3,1,1,0);

        hess(0,0) = img(xp,y,z) + img(xn,y,z) - 2*img(x,y,z);          // Hxx

        hess(1,1) = img(x,yp,z) + img(x,yn,z) - 2*img(x,y,z);          // Hyy

        hess(2,2) = img(x,y,zp) + img(x,y,zn) - 2*img(x,y,z);          // Hzz

        hess(1,0) = hess(0,1) = (img(xp,yp,z) + img(xn,yn,z) - img(xp,yn,z) - img(xn,yp,z))/4; // Hxy = Hyx

        hess(2,0) = hess(0,2) = (img(xp,y,zp) + img(xn,y,zn) - img(xp,y,zn) - img(xn,y,zp))/4; // Hxz = Hzx

        hess(2,1) = hess(1,2) = (img(x,yp,zp) + img(x,yn,zn) - img(x,yp,zn) - img(x,yn,zp))/4; // Hyz = Hzy

        //val : 3 eigen values as an array, decreasing order

        //vec : 3 eigen vectors as a 3*3 image

        //have to be double or create NaN

        CImg<double> val,vec;

        hess.symmetric_eigen(val,vec);

        if (val[2]>0) val[2]=0;

    res(x,y,z)=-val[2];

    }

  return res;

}

//Wall indicator

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

CImg<float> wall_indicator(CImg<float> const & img)

{

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

  cimg_forXYZ(res,x,y,z)

    {

      res(x,y,z)=1.0/(1.0+(img(x,y,z)*img(x,y,z)));

    }

  return res;

}

//Dirac function

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

CImg<float> Dirac(CImg<float> const & img, float sigma)

{

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

  float f=0;

  cimg_forXYZ(img,x,y,z)

    {

      if(abs(img(x,y,z))<=sigma)

        {

          f=(1.0/(2*sigma))*(1+cos(M_PI*img(x,y,z)/sigma));

        }

      else {f=0;}

      fImg(x,y,z)=f;

    }

  return fImg;

} 

//Normal vector

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

CImgList<float> normal_vector(CImg<float> const & img)

{

  CImgList<float> normal_vector(3,img._width,img._height,img._depth,1,0); 

  CImg_3x3x3(I,float);        

  cimg_for3x3x3(img,x,y,z,0,I,float) 

    {

      float nx = (Incc-Ipcc)/2.0;    

      float ny = (Icnc-Icpc)/2.0;   

      float nz = (Iccn-Iccp)/2.0;

      //border

      if((x==0)or(y==0)or(z==0))

        {nx=ny=nz=0;}

      if((x==img.width()-1)or(y==img.height()-1)or(z==img.depth()-1))

        {nx=ny=nz=0;}

      float norm=1;

      if((nx!=0)or(ny!=0)or(nz!=0))

        {norm = sqrt(nx*nx+ny*ny+nz*nz);}

      normal_vector(0,x,y,z)=nx/norm;

      normal_vector(1,x,y,z)=ny/norm;

      normal_vector(2,x,y,z)=nz/norm;

    }

  return normal_vector;

}

//Curvature

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

CImg<float> curvature2(CImgList<float> const & N)

{

  CImg<float> K(N[0]._width,N[0]._height,N[0]._depth,1,0);

  cimg_forXYZ(N[0],x,y,z)

    {

      float kx=0;

      float ky=0;

      float kz=0;

      if(x==0)

        {kx=N[0](x+1,y,z)-N[0](x,y,z);}

      else if(x==N[0]._width-1)

        {kx=N[0](x,y,z)-N[0](x-1,y,z);}

      else

        {kx=(N[0](x+1,y,z)-N[0](x-1,y,z))/2.0;}

     if(y==0)

        {ky=N[1](x,y+1,z)-N[1](x,y,z);}

      else if(y==N[1]._height-1)

        {ky=N[1](x,y,z)-N[1](x,y-1,z);}

      else

        {ky=(N[1](x,y+1,z)-N[1](x,y-1,z))/2.0;}

     if(z==0)

        {kz=N[2](x,y,z+1)-N[2](x,y,z);}

      else if(z==N[2]._depth-1)

        {kz=N[2](x,y,z)-N[2](x,y,z-1);}

      else

        {kz=(N[2](x,y,z+1)-N[2](x,y,z-1))/2.0;}

     K(x,y,z)=kx+ky+kz;

    }

  return K;

}

// Compute image laplacian

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

CImg<float> laplacian(CImg<float> const & img)

{

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

  CImg_3x3x3(I,float);

  cimg_for3x3x3(img,x,y,z,0,I,float)

    {

      laplacian(x,y,z) =Incc+Ipcc+Icnc+Icpc+Iccn+Iccp - 6*Iccc;

    }

  return laplacian;

}

//Evolution AK2 contour

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

CImg<float> evolution_AK2_contour(CImg<float> u, CImg<float> const & g, CImgList<float> const & gg, CImg<float> const & h, int lam, float mu, float alf, float beta, float epsilon, int dt)

{

  CImg<float> diracU=Dirac(u,epsilon); 

  CImgList<float> N=normal_vector(u);

  CImg<float> K=curvature2(N);

  CImg<float> L=laplacian(u);

  float term=0;

  cimg_forXYZ(u,x,y,z)

    {

      term=0;

      if (diracU(x,y,z)!=0)

        {

          //weighted length term

          if(lam!=0)

            {term=( gg(0,x,y,z)*N(0,x,y,z) + gg(1,x,y,z)*N(1,x,y,z) + gg(2,x,y,z)*N(2,x,y,z) + (g(x,y,z)*K(x,y,z))) *lam;}

          //length term

          if(beta!=0)

            {term=term + K(x,y,z)*beta;}

          //weighted area term

          if(alf!=0)

            {term=term + h(x,y,z)*alf;}

          //regularization

          term=term*diracU(x,y,z);

          //penalizing term

          term=term + (L(x,y,z)-K(x,y,z))*mu;

        }

      else

        {

          //penalizing term

          term=(L(x,y,z)-K(x,y,z))*mu;

        }

      u(x,y,z)=u(x,y,z) + term*dt;

    }

  return u;

}

//Evolution AK2 contour cells

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

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

{

  CImg<float> diracU=Dirac(u,epsilon); 

  CImgList<float> N=normal_vector(u);

  CImg<float> K=curvature2(N);

  CImg<float> L=laplacian(u);

  float term=0;

  int xmax=u._width;

  int ymax=u._height;

  int zmax=u._depth;

  int c0=-4;

  cimg_forXYZ(u,x,y,z)

    {

      int xb=x+min_list[0];

      int yb=y+min_list[1];

      int zb=z+min_list[2];

      term=0;

      if (diracU(x,y,z)!=0)

        {

          //weighted length term

          if(lam!=0)

            {term=( gg(0,xb,yb,zb)*N(0,x,y,z) + gg(1,xb,yb,zb)*N(1,x,y,z) + gg(2,xb,yb,zb)*N(2,x,y,z) + (g(xb,yb,zb)*K(x,y,z))) *lam;}

          //length term

          if(beta!=0)

            {term=term + K(x,y,z)*beta;}

          //weighted area term

          if(alf!=0)

            {term=term + h(xb,yb,zb)*alf;}            

          //regularization

          term=term*diracU(x,y,z);

          //penalizing term

          term=term + (L(x,y,z)-K(x,y,z))*mu;

        }

      else

        {

          //penalizing term

          term=(L(x,y,z)-K(x,y,z))*mu;

        }

      // computing the evolution

      u(x,y,z)=u(x,y,z) + term*dt;

      // box boundary regularisation

      if ((x==0) or (x==xmax) or (x==1) or (x==xmax-1) or(y==0) or (y==ymax)or(y==1) or (y==ymax-1) or (z==0) or (z==zmax)or (z==1) or (z==zmax-1))

      {u(x,y,z)=c0;}

    }

  return u;

}