/*-----------------------------------------------------------------------------

    File        : greycstoration.cpp

    Description : GREYCstoration - A tool to denoise, inpaint and resize images.

    The GREYCstoration algorithm is an implementation of diffusion tensor-directed
    diffusion PDE's for image regularization and interpolation, published in

    "Fast Anisotropic Smoothing of Multi-Valued Images 
    using Curvature-Preserving PDE's"
    (D. Tschumperlé)
    Research Report "Cahiers du GREYC 05/01", February 2005),
    GREYC, CNRS UMR 6072
    (see also http://www.greyc.ensicaen.fr/~dtschump/greycstoration)
    
    "Vector-Valued Image Regularization with PDE's : A Common Framework
    for Different Applications"
    (D. Tschumperlé, R. Deriche).
    IEEE Transactions on Pattern Analysis and Machine Intelligence,
    Vol 27, No 4, pp 506-517, April 2005.

    Copyright  : David Tschumperle - http://www.greyc.ensicaen.fr/~dtschump/

  This software is governed by the CeCILL  license under French law and
  abiding by the rules of distribution of free software.  You can  use, 
  modify and/ or redistribute the software under the terms of the CeCILL
  license as circulated by CEA, CNRS and INRIA at the following URL
  "http://www.cecill.info". 
  
  As a counterpart to the access to the source code and  rights to copy,
  modify and redistribute granted by the license, users are provided only
  with a limited warranty  and the software's author,  the holder of the
  economic rights,  and the successive licensors  have only  limited
  liability. 
  
  In this respect, the user's attention is drawn to the risks associated
  with loading,  using,  modifying and/or developing or reproducing the
  software by the user in light of its specific status of free software,
  that may mean  that it is complicated to manipulate,  and  that  also
  therefore means  that it is reserved for developers  and  experienced
  professionals having in-depth computer knowledge. Users are therefore
  encouraged to load and test the software's suitability as regards their
  requirements in conditions enabling the security of their systems and/or 
  data to be ensured and,  more generally, to use and operate it in the 
  same conditions as regards security. 
  
  The fact that you are presently reading this means that you have had
  knowledge of the CeCILL license and that you accept its terms.

  ------------------------------------------------------------------------------*/

#define greycstoration_version 2.4
#include "../CImg.h"
// The line below is necessary for compilation with non-standart C++ compilers.
#if (( defined(_MSC_VER) && _MSC_VER<=1200 ) || defined(__DMC__))
#define std
#endif

using namespace cimg_library;
typedef unsigned char uchar;
unsigned int iter = 0, nb_iter = 0;
CImgDisplay *disp = 0;

// get_geom() : read geometry from a string ('320x256' or '200%x200%' are valid geometries).
//------------
void get_geom(const char *geom, int &geom_w, int &geom_h) {
  char tmp[16];
  std::sscanf(geom,"%d%7[^0-9]%d%7[^0-9]",&geom_w,tmp,&geom_h,tmp+1);
  if (tmp[0]=='%') geom_w=-geom_w;
  if (tmp[1]=='%') geom_h=-geom_h;
}

// Do one iteration of a Curvature-preserving PDE's.
//--------------------------------------------------
template<typename T,typename t> 
inline CImg<T>& blur_anisotropicPDE(CImg<T>& img,const CImg<t>& G,
				    const float amplitude=10.0f, const float dl=0.8f,const float da=30.0f,
				    const float gauss_prec=2.0f, const unsigned int interpolation=0,
				    const bool fast_approx=true) {

  static unsigned long last_time = cimg::time();
  // Check arguments and init variables
  if (img.is_empty()) throw CImgArgumentException("GREYCstoration : Specified input image is empty.");
  if (amplitude==0) return img;

  CImg<float> dest(img.width,img.height,1,img.dim,0), tmp(img.dim), W(img.width,img.height,1,3);
  const float sqrt2amplitude = (float)std::sqrt(2*amplitude);
  const int dx1 = img.dimx()-1, dy1 = img.dimy()-1;
  int counter=0, countermax = (int)(img.width*img.height*(360/da)), N = 0;

  for (float theta=(360%(int)da)/2.0f; theta<360; (theta+=da),N++) {
    const float
      thetar = (float)(theta*cimg::PI/180),
      vx = (float)(std::cos(thetar)),
      vy = (float)(std::sin(thetar));
    
    const t *pa = G.ptr(0,0,0,0), *pb = G.ptr(0,0,0,1), *pc = G.ptr(0,0,0,2);
    t *pd0 = W.ptr(0,0,0,0), *pd1 = W.ptr(0,0,0,1), *pd2 = W.ptr(0,0,0,2);
    cimg_mapXY(G,xg,yg) {
      const t
	a = *(pa++), b = *(pb++), c = *(pc++), 
	u = a*vx + b*vy, v = b*vx + c*vy,
	n = (t)std::sqrt(1e-5+u*u+v*v),
	dln = dl/n;
      *(pd0++) = u*dln;
      *(pd1++) = v*dln;
      *(pd2++) = n;
    }
    
    cimg_mapXY(img,x,y) {
      
      if (!((counter++)%200000)) {
	const unsigned int pc = (unsigned int)(100.0*(counter+iter*countermax)/(nb_iter*countermax));
	std::fprintf(stderr,"\r- Processing : %u %% (iteration %u)              ",pc,iter+1);
	if (disp) {
	  disp->set_title("GREYCstoration (%u%%)",pc,iter+1);
	  if (disp->is_closed || disp->key==cimg::keyQ || disp->key==cimg::keyESC) return img;
	  const unsigned long curr_time = cimg::time();
	  if (disp && curr_time-last_time>=2500) {
	    last_time = curr_time;
	    disp->display(dest);
	  }	  
	}
      }

      tmp.fill(0);
      const t cu = W(x,y,0,0), cv = W(x,y,0,1), n = W(x,y,0,2);
      const float
	fsigma = (float)(n*sqrt2amplitude),
	length = gauss_prec*fsigma,
	fsigma2 = 2*fsigma*fsigma;
      float l, pu=cu, pv=cv, X=(float)x, Y=(float)y;
      
      if (fast_approx) {
	unsigned int S=0;
	switch (interpolation) {
	case 0: // Nearest-neighbor interpolation
	  for (l=0; l<length; l+=dl) {
	    const float 
	      Xn = X<0?0:(X>=dx1?dx1:X),
	      Yn = Y<0?0:(Y>=dy1?dy1:Y);
	    const int xi = (int)(Xn+0.5f), yi = (int)(Yn+0.5f);
	    t u = W(xi,yi,0,0), v = W(xi,yi,0,1);
	    if ((pu*u+pv*v)<0) { u=-u; v=-v; }
	    cimg_mapV(img,k) tmp[k]+=(t)img(xi,yi,0,k);
	    X+=(pu=u); Y+=(pv=v); ++S;
	  } break;
	case 1: // Linear interpolation
	  for (l=0; l<length; l+=dl) {
	    t u = (t)(W.linear_pix2d(X,Y,0,0)), v = (t)(W.linear_pix2d(X,Y,0,1));
	    if ((pu*u+pv*v)<0) { u=-u; v=-v; }
	    cimg_mapV(img,k) tmp[k]+=(t)img.linear_pix2d(X,Y,0,k);
	    X+=(pu=u); Y+=(pv=v); ++S;
	  } break;
	default: // 2nd-order Runge-kutta interpolation
	  for (l=0; l<length; l+=dl) {
	    t u0 = (t)(0.5f*dl*W.linear_pix2d(X,Y,0,0)), v0 = (t)(0.5f*dl*W.linear_pix2d(X,Y,0,1));
	    if ((pu*u0+pv*v0)<0) { u0=-u0; v0=-v0; }
	    t u = (t)(W.linear_pix2d(X+u0,Y+v0,0,0)), v = (t)(W.linear_pix2d(X+u0,Y+v0,0,1));
	    if ((pu*u+pv*v)<0) { u=-u; v=-v; }
	    cimg_mapV(img,k) tmp[k]+=(t)img.linear_pix2d(X,Y,0,k);
	    X+=(pu=u); Y+=(pv=v); ++S;
	  } break;
	}
	if (S>0) cimg_mapV(dest,k) dest(x,y,0,k)+=tmp[k]/S;
	else cimg_mapV(dest,k) dest(x,y,0,k)+=(t)(img(x,y,0,k));
      }
      else {
	float S=0;
	switch (interpolation) {
	case 0: // Nearest-neighbor interpolation
	  for (l=0; l<length; l+=dl) {
	    const float 
	      coef = (float)std::exp(-l*l/fsigma2),
	      Xn = X<0?0:(X>=dx1?dx1:X),
	      Yn = Y<0?0:(Y>=dy1?dy1:Y);
	    const int xi = (int)(Xn+0.5f), yi = (int)(Yn+0.5f);
	    t u = W(xi,yi,0,0), v = W(xi,yi,0,1);
	    if ((pu*u+pv*v)<0) { u=-u; v=-v; }
	    cimg_mapV(img,k) tmp[k]+=(t)(coef*img(xi,yi,0,k));
	    X+=(pu=u); Y+=(pv=v); S+=coef;
	  } break;
	case 1: // Linear interpolation
	  for (l=0; l<length; l+=dl) {
	    const float coef = (float)std::exp(-l*l/fsigma2);
	    t u = (t)(W.linear_pix2d(X,Y,0,0)), v = (t)(W.linear_pix2d(X,Y,0,1));
	    if ((pu*u+pv*v)<0) { u=-u; v=-v; }
	    cimg_mapV(img,k) tmp[k]+=(t)(coef*img.linear_pix2d(X,Y,0,k));
	    X+=(pu=u); Y+=(pv=v); S+=coef;
	  } break;
	default: // 2nd-order Runge-kutta interpolation
	  for (l=0; l<length; l+=dl) {
	    const float coef = (float)std::exp(-l*l/fsigma2);
	    t u0 = (t)(0.5f*dl*W.linear_pix2d(X,Y,0,0)), v0 = (t)(0.5f*dl*W.linear_pix2d(X,Y,0,1));
	    if ((pu*u0+pv*v0)<0) { u0=-u0; v0=-v0; }
	    t u = (t)(W.linear_pix2d(X+u0,Y+v0,0,0)), v = (t)(W.linear_pix2d(X+u0,Y+v0,0,1));
	    if ((pu*u+pv*v)<0) { u=-u; v=-v; }
	    cimg_mapV(img,k) tmp[k]+=(t)(coef*img.linear_pix2d(X,Y,0,k));
	    X+=(pu=u); Y+=(pv=v); S+=coef;
	  } break;
	}
	if (S>0) cimg_mapV(dest,k) dest(x,y,0,k)+=tmp[k]/S;
	else cimg_mapV(dest,k) dest(x,y,0,k)+=(t)(img(x,y,0,k));
      }
    }    
  }
  const float *ptrs = dest.data+dest.size(); cimg_map(img,ptrd,T) *ptrd = (T)(*(--ptrs)/N);
  return img;
}

// Do one iteration of Shock-Filter PDE
//--------------------------------------
template<typename T,typename t> 
inline CImg<T>& shock_filterPDE(CImg<T>& img, const CImg<t>& G, const float dt=10.0f, const float edge=1.0f) {
  CImg<> val,vec,U(img.dimx(),img.dimy(),1,3),veloc(img);
  CImg_3x3(I,float);
  const float nedge = 0.5f*edge;
  cimg_mapXY(img,x,y) {
    if (G(x,y,0)==G(x,y,1) && G(x,y,1)==G(x,y,2) && G(x,y,2)==0) 
      U(x,y,0)=U(x,y,1)=U(x,y,2)=0;
    else {
      G.get_tensor_at(x,y).symmetric_eigen(val,vec);
      U(x,y,0) = vec(0,0);
      U(x,y,1) = vec(0,1);
      U(x,y,2) = 1.0f-(float)std::pow(1.0f+val[0]+val[1],-nedge);
    }
  } 
  cimg_mapV(img,k) cimg_map3x3(img,x,y,0,k,I) {
    const float 
      a = U(x,y,0), b = U(x,y,1), norm = U(x,y,2),
      ixx = Inc+Ipc-2*Icc, iyy = Icn+Icp-2*Icc, ixy = 0.25f*(Inn+Ipp-Inp-Ipn),
      ixf = Inc-Icc, ixb = Icc-Ipc, iyf = Icn-Icc, iyb = Icc-Icp,
      itt = a*a*ixx + 2*a*b*ixy + b*b*iyy,
      it = a*cimg::minmod(ixf,ixb) + b*cimg::minmod(iyf,iyb);
    veloc(x,y,k) = -norm*cimg::sign(itt)*cimg::abs(it);
  }
  const CImgStats stats(veloc,false);
  const float vmax = (float)cimg::max(cimg::abs(stats.min),cimg::abs(stats.max));
  if (vmax>1e-5f) img+=veloc*dt/vmax;
  return img;
}

/*-----------------
  Main procedure
  ----------------*/
int main(int argc,char **argv) {
  cimg_usage("Fast Anisotropic Smoothing of Multi-valued Images with Curvature-Preserving PDE's");
  std::fprintf(stderr,
	       "---------------------------------------------------------------------------------------------\n"
	       " GREYCstoration v%g :\n"
	       " ------------------------------\n"
	       " The GREYCstoration algorithm is the result of a research work done in the\n"
	       " IMAGE group of the GREYC Lab (CNRS, UMR 6072)/France, by David Tschumperle\n"
	       " ( http://www.greyc.ensicaen.fr/EquipeImage/ and http://www.greyc.ensicaen.fr/~dtschump/ ).\n"
	       " It can be used to denoise, inpaint or resize 2D color images.\n"
	       " This algorithm is intended to help people for processing image data and for comparison\n"
	       " purposes. This is an open source software, distributed within the CImg Library package\n"
	       " (http://cimg.sourceforge.net), and submitted to the CeCiLL License (view file LICENSE.txt).\n"
	       " If you are interested to distribute this algorithm in a closed-source product,\n"
	       " you are invited to contact David Tschumperle (mail available on his web page).\n"
	       " Please also cite the paper below when using results of the GREYCstoration\n"
	       " algorithm in your own publications :\n\n"
	       "    \"Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE's.\"\n"
	       "    ( D. Tschumperle )\n"
	       "    International Journal of Computer Vision, IJCV(68), No. 1, June 2006, pp. 65-82.\n"
	       " You can also take a look at the official GREYCstoration web page :\n"
	       "    http://www.greyc.ensicaen.fr/~dtschump/greycstoration/\n"
	       "---------------------------------------------------------------------------------------------\n\n",
	       greycstoration_version);

  // Read command line parameters
  //------------------------------
  const bool help = cimg_option("-h",false,(char*)0) || cimg_option("--help",false,(char*)0) || cimg_option("-help",false,(char*)0);

  if (help) std::fprintf(stderr,
			 "  + Input/Output parameters\n"
			 "    -----------------------\n");
  const char *restore = cimg_option("-restore",(const char*)0,"Restore the image specified after '-restore'");
  const char *inpaint = cimg_option("-inpaint",(const char*)0,"Inpaint the image specified after '-inpaint'");
  const char *resize = cimg_option("-resize",(const char*)0, "Resize the image specified after '-resize'");
  const char *refer = cimg_option("-reference",(const char*)0,"Reference image to compare with");
  const unsigned int save = cimg_option("-save",0,"Iteration saving step (>=0)");
  const bool visu = cimg_option("-visu",true,"Enable/Disable visualization (0 or 1)");
  const char *file_o = cimg_option("-o",(const char*)0,"Output image");
  const bool append_result = cimg_option("-append",false,"Append original image to output file");    

  if (help) std::fprintf(stderr,"\n"
			 "  + Global algorithm parameters\n"
			 "    ---------------------------\n");
  const float dt = cimg_option("-dt",restore?50.0f:(inpaint?20.0f:20.0f),"Regularization strength for one iteration (>=0)");
  nb_iter = cimg_option("-iter",restore?1:(inpaint?1000:3),"Number of regularization iterations (>0)");
  const float p = cimg_option("-p",restore?0.9f:(inpaint?0.3f:0.2f),"Contour preservation for regularization (>=0)");
  const float anisotropy = cimg_option("-a",restore?0.2f:(inpaint?1.0f:0.9f),"Regularization anisotropy (0<=a<=1)");
  const float alpha = cimg_option("-alpha",restore?0.5f:(inpaint?0.8f:0.1f),"Noise scale(>=0)");
  const float sigma = cimg_option("-sigma",restore?1.0f:(inpaint?2.0f:1.5f),"Geometry regularity (>=0)");
  const bool fast_approx = cimg_option("-fast",true,"Use fast approximation for regularization (0 or 1)");
  const float gauss_prec = cimg_option("-prec",2.0f,"Precision of the gaussian function for regularization (>0)");
  const float dl = cimg_option("-dl",0.8f,"Spatial integration step for regularization (0<=dl<=1)");
  const float da = cimg_option("-da",30.0f,"Angular integration step for regulatization (0<=da<=90)");
  const unsigned int interp = cimg_option("-interp",0,"Interpolation type (0=Nearest-neighbor, 1=Linear, 2=Runge-Kutta)");
  const float sdt = cimg_option("-sdt",restore?0.0f:(inpaint?0.0f:10.0f),"Sharpening strength for one iteration (>=0)");
  const float sp = cimg_option("-sp",restore?0.5f:(inpaint?0.5f:3.0f),"Sharpening edge threshold (>=0)");


  if (argc==1 || (!restore && !inpaint && !resize)) {
    cimg::dialog(cimg::basename(argv[0]),
		 "You must specify one of the '-restore', '-inpaint', or '-resize' options.\n"
		 "(try option '-h', '-h -restore','-h -inpaint' or '-h -resize' to get options relative to specific actions\n");
    std::exit(0);
  }
  
  // Init variables
  //----------------
  CImg<float> img0,img,imgr;
  CImg<uchar> mask;
  float valmin=0, valmax=0;
  
  // Init parameters and variables for image restoration
  //-----------------------------------------------------
  if (restore) { 
    if (help) std::fprintf(stderr,"\n"
			   "  + Restoration mode parameters\n"
			   "    ---------------------------\n");
    const float ng = cimg_option("-ng",0.0f,"Add synthetic gaussian noise before denoising (>=0)");
    const float nu = cimg_option("-nu",0.0f,"Add synthetic uniform noise before denoising (>=0)");
    const float ns = cimg_option("-ns",0.0f,"Add synthetic Salt&Pepper noise before denoising (>=0)");
    std::fprintf(stderr,"- Image Restoration mode\n");
    if (!cimg::strcmp("-restore",restore)) {
      cimg::dialog(cimg::basename(argv[0]),"You must specify a valid input image filename after the '-restore' flag.");
      std::exit(0);
    }
    std::fprintf(stderr,"- Load input image '%s'...",cimg::basename(restore));
    img.load(restore);
    const CImgStats stats(img,false);
    valmin = (float)stats.min; valmax = (float)stats.max;
    std::fprintf(stderr,"\r- Input image : '%s' (size %ux%u, value range [%g,%g])\n",
		 cimg::basename(restore),img.width,img.height,stats.min,stats.max);   
    if (ng>0 || nu>0 || ns>0) {
      img0 = img;
      if (ng>0) img.noise(ng,0);
      if (nu>0) img.noise(nu,1);
      if (ns>0) img.noise(ns,2);
      img.cut(valmin,valmax);
      const CImgStats stats(img,false);
      std::fprintf(stderr,"\r- Noisy image : value range [%g,%g], PSNR Noisy / Original : %g\n",
		   stats.min,stats.max,img.PSNR(img0));      
    }
  }
  
  // Init parameters and variables for image inpainting
  //----------------------------------------------------
  if (inpaint) {
    if (help) std::fprintf(stderr,"\n"
			   "  + Inpainting mode parameters\n"
			   "    --------------------------\n");
    const char *file_m = cimg_option("-m",(const char*)0,"Input inpainting mask");
    const unsigned int dilate = cimg_option("-dilate",0,"Inpainting mask dilatation (>=0)");
    const unsigned int init = cimg_option("-init",4,"Inpainting init (0=black, 1=white, 2=noise, 3=unchanged, 4=interpol)");
    std::fprintf(stderr,"- Image Inpainting mode\n");
    if (!cimg::strcmp("-inpaint",inpaint)) {
      cimg::dialog(cimg::basename(argv[0]),"You must specify a valid input image filename after the '-inpaint' flag.");
      std::exit(0);
    }
    std::fprintf(stderr,"- Load input image '%s'...",cimg::basename(inpaint));
    img.load(inpaint);
    const CImgStats stats(img,false);
    valmin = (float)stats.min; valmax = (float)stats.max;
    std::fprintf(stderr,"\r- Input image : '%s' (size %ux%u, value range [%g,%g])\n",
		 cimg::basename(inpaint),img.width,img.height,stats.min,stats.max);
    if (!file_m) {
      cimg::dialog(cimg::basename(argv[0]),"You need to specify a valid inpainting mask filename after the '-m' flag.");
      std::exit(0);
    }
    if (cimg::strncasecmp("block",file_m,5)) {
      std::fprintf(stderr,"- Load inpainting mask '%s'...",cimg::basename(file_m));
      mask.load(file_m);
      std::fprintf(stderr,"\r- Inpainting mask : '%s' (size %ux%u)\n",cimg::basename(file_m),mask.width,mask.height);
    }
    else {
      int l=16; std::sscanf(file_m,"block%d",&l);
      mask.assign(img.dimx()/l,img.dimy()/l);
      cimg_mapXY(mask,x,y) mask(x,y)=(x+y)%2;
      img0 = img;
    }
    mask.resize(img.dimx(),img.dimy(),1,1);
    if (dilate) mask.dilate(dilate);
    switch (init) {
    case 0 : { cimg_mapXYV(img,x,y,k) if (mask(x,y)) img(x,y,k) = 0; } break;
    case 1 : { cimg_mapXYV(img,x,y,k) if (mask(x,y)) img(x,y,k) = 255; } break;
    case 2 : { cimg_mapXYV(img,x,y,k) if (mask(x,y)) img(x,y,k) = (float)(255*cimg::rand()); } break;
    case 3 : break;
    default: {
      CImg<uchar> tmask(mask),ntmask(tmask);
      CImg_3x3(M,uchar);
      CImg_3x3(I,float);
      while (CImgStats(ntmask,false).max>0) {
	cimg_map3x3(tmask,x,y,0,0,M) if (Mcc && (!Mpc || !Mnc || !Mcp || !Mcn)) {
	  const float ccp = Mcp?0.0f:1.0f, cpc = Mpc?0.0f:1.0f,
	    cnc = Mnc?0.0f:1.0f, ccn = Mcn?0.0f:1.0f, csum = ccp + cpc + cnc + ccn;
	  cimg_mapV(img,k) {
	    cimg_get3x3(img,x,y,0,k,I);
	    img(x,y,k) = (ccp*Icp + cpc*Ipc + cnc*Inc + ccn*Icn)/csum;
	  }
	  ntmask(x,y) = 0;
	}
	tmask = ntmask;
      }
    } break;    
    }
  }
  
  // Init parameters and variables for image resizing
  //--------------------------------------------------
  if (resize) {
    if (help) std::fprintf(stderr,"\n"
			   "  + Resizing mode parameters\n"
			   "    ------------------------\n");
    const char *geom0 = cimg_option("-g0",(const char*)0,"Input image geometry");
    const char *geom  = cimg_option("-g",(const char*)0,"Output image geometry");
    const bool anchor = cimg_option("-anchor",true,"Anchor original pixels");
    const unsigned int init = cimg_option("-init",5,"Initial estimate (1=block, 3=linear, 5=bicubic)");
    std::fprintf(stderr,"- Image Resizing mode\n");
    if (!geom && !geom0) {
      cimg::dialog(cimg::basename(argv[0]),"You need to specify an output geometry or an input geometry (option -g or -g0)");
      std::exit(0); 
    }
    if (!cimg::strcmp("-resize",resize)) {
      cimg::dialog(cimg::basename(argv[0]),"You must specify a valid input image filename after the '-resize' flag.");
      std::exit(0);
    }
    std::fprintf(stderr,"- Load input image '%s'...",cimg::basename(resize));
    img.load(resize);
    const CImgStats stats(img,false);
    valmin = (float)stats.min; valmax = (float)stats.max;
    std::fprintf(stderr,"\r- Input image : '%s' (size %ux%u, value range [%g,%g])\n",
		 cimg::basename(resize),img.width,img.height,stats.min,stats.max);
    int w,h;
    if (geom0) {
      int w0,h0;
      get_geom(geom0,w0,h0);
      w0 = w0>0?w0:-w0*img.dimx()/100;
      h0 = h0>0?h0:-h0*img.dimy()/100;
      std::fprintf(stderr,"- Reducing geometry to %dx%d using %s interpolation.\n",w0,h0,
		   init==1?"bloc":(init==3?"linear":"bicubic"));
      img0.assign(img);
      w = img.dimx();
      h = img.dimy();
      img.resize(w0,h0,-100,-100,5);
    }
    if (geom) {
      get_geom(geom,w,h);
      w = w>0?w:-w*img.dimx()/100;
      h = h>0?h:-h*img.dimy()/100;
    }    
    mask.assign(img.dimx(),img.dimy(),1,1,255);
    if (!anchor) mask.resize(w,h,1,1,1); else mask = ~mask.resize(w,h,1,1,4);
    img.resize(w,h,1,-100,init);
    if (!img0.is_empty()) std::fprintf(stderr,"\r- PSNR Original / Thumbnail : %g\n",img.PSNR(img0));
    img.cut(valmin,valmax);
  }

  // Load reference image if any
  //-----------------------------
  if (refer) {
    std::fprintf(stderr,"- Load reference image '%s'...",cimg::basename(refer));
    imgr.load(refer);
    CImgStats stats(imgr,false);
    std::fprintf(stderr,"\r- Reference image : '%s' (size %ux%u, value range [%g,%g])",
		 cimg::basename(refer),imgr.width,imgr.height,stats.min,stats.max);
    if (!img0.is_empty()) {
      imgr.resize(img0);
      std::fprintf(stderr,", PSNR Reference / Original : %g dB\n",imgr.PSNR(img0));
    } else {
      imgr.resize(img);
      std::fprintf(stderr,"\n");
    }
  }
  
  // Init images and display  
  //-------------------------
  CImg<> dest(img);
  if (!visu && !append_result) img.empty();  
  if (visu) {
    const int sx = 2*CImgDisplay::screen_dimx()/3, sy = 2*CImgDisplay::screen_dimy()/3;
    int nwidth=dest.dimx(), nheight=dest.dimy();
    if (nwidth>sx) { nheight = nheight*sx/nwidth; nwidth = sx; }
    if (nheight>sy) { nwidth = nwidth*sy/nheight; nheight = sy; }
    disp = new CImgDisplay(img.get_resize(nwidth,nheight),"GREYCstoration",1);
    dest.display(*disp);
  }

  //-------------------------------------
  // Begin regularization PDE iterations
  //-------------------------------------
  std::fprintf(stderr,"\r- Processing : initialization...    \t\t");
  if (disp) disp->set_title("GREYCstoration (initialization)");

  bool stopflag = false;
  for (iter=0; iter<nb_iter && !stopflag && (!disp || !disp->is_closed); iter++) {
    
    // Compute smoothed structure tensor field G
    //-------------------------------------------
    CImg_3x3(I,float);
    CImg<> G(dest.dimx(),dest.dimy(),1,3,0);
    const CImg<> blur = dest.get_blur(alpha).normalize(valmin,valmax);
    cimg_mapV(blur,k) cimg_map3x3(blur,x,y,0,k,I) {
      const float ixf = Inc-Icc, iyf = Icn-Icc, ixb = Icc-Ipc, iyb = Icc-Icp;
      G(x,y,0) += 0.5f*(ixf*ixf+ixb*ixb);
      G(x,y,1) += 0.25f*(ixf*iyf+ixf*iyb+ixb*iyf+ixb*iyb);
      G(x,y,2) += 0.5f*(iyf*iyf+iyb*iyb);
    } 
    if (sigma>0) G.blur(sigma);    
    if (mask.data) cimg_mapXY(G,x,y) if (!mask(x,y)) G(x,y,0) = G(x,y,1) = G(x,y,2) = 0;

    // Compute tensor field sqrt(T)
    //------------------------------
    if (dt>0) {
      CImg<> T(G.dimx(),G.dimy(),1,3), val, vec;
      const float 
	np1 = 0.5f*p,
	np2 = np1/(1e-7+1-anisotropy);
      float nmax = 0;  
      cimg_mapXY(G,x,y) {
	G.get_tensor_at(x,y).symmetric_eigen(val,vec);
	if (G(x,y,0)==G(x,y,1) && G(x,y,1)==G(x,y,2) && G(x,y,2)==0) T(x,y,0)=T(x,y,1)=T(x,y,2)=0;
	else {
	  const float l1 = val[1], l2 = val[0], u = vec(1,0), v = vec(1,1),      
	    n1 = (float)std::pow(1.0f+l1+l2,-np1),
	    n2 = (float)std::pow(1.0f+l1+l2,-np2);
	  T(x,y,0) = n1*u*u + n2*v*v;
	  T(x,y,1) = u*v*(n1-n2);
	  T(x,y,2) = n1*v*v + n2*u*u;
	  if (n1>nmax) nmax=n1;
	}
      }    
      std::fprintf(stderr,"\n- Maximum smoothing amplitude = %g\n",nmax*dt);
      blur_anisotropicPDE(dest,T,dt,dl,da,gauss_prec,interp,fast_approx);
      dest.cut(valmin,valmax);
    }
    if (sdt>0) shock_filterPDE(dest,G,sdt,sp);

    // Prepare for next iteration
    //---------------------------
    if (disp && disp->is_resized) disp->resize();
    if (disp && visu) dest.display(*disp);
    if (file_o && save && !(iter%save)) dest.save(file_o,iter);

    // Display result and allows user interaction if needed.
    //-------------------------------------------------------
    if (disp && (disp->key==cimg::keyQ || disp->key==cimg::keyESC)) nb_iter=iter+1;

    if (iter==nb_iter-1) {
      std::fprintf(stderr,"\r- Processing : Done !                          \n");
      if (!img0.is_empty()) std::fprintf(stderr,"- PSNR Restored / Original : %g dB\n",dest.PSNR(img0));
      if (disp) {
	static bool first_time = true;
	if (first_time) {
	  first_time = false;
	  std::fprintf(stderr,
		       "- User interface :\n"
		       " > You can now zoom to a particular rectangular region, or press one of the following key :\n"
		       "   SPACE : Swap between original/restored view.\n"
		       "   S     : Save a snapshot of the image.\n"
		       "   I     : Run another iteration.\n"
		       "   Q     : Quit.\n");
	}

	CImgl<> visu;
	visu.insert_shared(img0).insert_shared(img).insert_shared(dest).insert_shared(imgr);
	const char *titles[4] = { "original", "noisy", "restored", "reference"};
	unsigned int visupos = 2;
	CImgDisplay *dispz = 0;
	CImg<> zoom;
	int s[6], snb=0;
	bool quit = false;
	while (!quit && !disp->is_closed) {
	  disp->set_title("GREYCstoration (%s)",titles[visupos]);
	  visu(visupos).feature_selection(s,2,*disp);
	  switch (disp->key) {
	  case cimg::keySPACE: do { visupos = (visupos+1)%visu.size; } while (visu(visupos).is_empty()); break;
	  case cimg::keyBACKSPACE: do { visupos = (visupos-1+visu.size)%visu.size; } while (visu(visupos).is_empty()); break;
	  case cimg::keyQ: disp->is_closed = quit = stopflag = true; break;
	  case cimg::keyI: 
	    quit = true;
	    std::fprintf(stderr,"- Perform iteration %u...\n",++nb_iter);
	    dest.display(*disp);
	    break;
	  case cimg::keyS:
	    if (!snb) dest.save(file_o?file_o:"GREYCstoration.bmp");
	    if (!zoom.is_empty()) zoom.save(file_o?file_o:"GREYCstoration.bmp",snb);
	    std::fprintf(stderr,"- Snapshot %u : '%s' saved\n",snb++,file_o?file_o:"GREYCstoration.bmp");
	    break;
	  }
	  disp->key = 0;	  
	  if (disp->is_resized) disp->resize().display(visu(visupos));
	  if (dispz && dispz->is_resized) dispz->resize().display(zoom);
	  if (dispz && dispz->is_closed) { delete dispz; dispz = 0; }

	  if (s[0]>=0 && s[1]>=0 && s[3]>=0 && s[4]>=0) {
	    const int x0 = s[0], y0 = s[1], x1 = s[3], y1 = s[4];
	    if (cimg::abs(x0-x1)>4 && cimg::abs(y0-y1)>4) {
	      CImgl<> tmp(img.get_crop(x0,y0,x1,y1), dest.get_crop(x0,y0,x1,y1));
	      if (!img0.is_empty()) tmp.insert(img0.get_crop(x0,y0,x1,y1),0);
	      if (!imgr.is_empty()) tmp.insert(imgr.get_crop(x0,y0,x1,y1));
	      zoom = tmp.get_append('x','c');
	      if (!dispz) {
		const int sx = 5*CImgDisplay::screen_dimx()/6, sy = 5*CImgDisplay::screen_dimy()/6;
		int nwidth = zoom.dimx(), nheight = zoom.dimy();
		if (nwidth>nheight) { nheight = nheight*sx/nwidth; nwidth = sx; }
		else { nwidth = nwidth*sy/nheight; nheight = sy; }
		dispz = new CImgDisplay(zoom.get_resize(nwidth,nheight),0,1);
		dispz->set_title("GREYCstoration (zoom) : - %s %s %s %s",
				 img0.is_empty()?"":"original -",
				 img.is_empty()?"":"noisy -",
				 dest.is_empty()?"":"restored -",
				 imgr.is_empty()?"":"reference -");
	      } else dispz->resize(dispz->dimx(),dispz->dimx()*zoom.dimy()/zoom.dimx(),false);
	      dispz->display(zoom).show();
	    }
	  }
	}
	if (dispz) delete dispz;
      }
    }
  }
  
  // Save result and exit
  //----------------------
  dest.cut(valmin,valmax);
  if (file_o) {
    std::fprintf(stderr,"\n- Saving output image '%s'\n",file_o);
    if (!append_result) dest.save(file_o);
    else CImgl<>(img,dest).get_append('x').save(file_o);
  }
  if (disp) delete disp;
  std::fprintf(stderr,"\n- Quit\n\n");

  return 0;
}
