Laboratoire RDP » Biophysics Team » lsm3d

//Evolution AK2 contour INTERRACT

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

CImg<float> evolution_AK2_contour_interact(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, CImg<unsigned char> const & free, 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;

  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(free(xb,yb,zb)==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;

	    }	

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

	}

    }

  return u;

}