PTdecode/CImg-1.3.0/examples/greycstoration.cpp

changeset 5
1204ebf9340d
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/PTdecode/CImg-1.3.0/examples/greycstoration.cpp	Mon Aug 03 14:09:20 2009 +0100
     1.3 @@ -0,0 +1,576 @@
     1.4 +/*
     1.5 + #
     1.6 + #  File        : greycstoration.cpp
     1.7 + #                ( C++ source file )
     1.8 + #
     1.9 + #  Description : GREYCstoration - A tool to denoise, inpaint and resize images.
    1.10 + #                This file is a part of the CImg Library project.
    1.11 + #                ( http://cimg.sourceforge.net )
    1.12 + #                See also the GREYCstoration web page
    1.13 + #                ( http://www.greyc.ensicaen.fr/~dtschump/greycstoration )
    1.14 + #
    1.15 + #    The GREYCstoration algorithm is an implementation of tensor-directed and
    1.16 + #    patch-based diffusion PDE's for image regularization and interpolation,
    1.17 + #    published in
    1.18 + #
    1.19 + #    "Defining Some Variational Methods on the Space of Patches :
    1.20 + #     Application to Multi-Valued Image Denoising and Registration"
    1.21 + #    (D. Tschumperle, L. Brun)
    1.22 + #    Rapport de recherche : Les cahiers du GREYC No 08-01, Mars 2008.
    1.23 + #
    1.24 + #    "Fast Anisotropic Smoothing of Multi-Valued Images
    1.25 + #    using Curvature-Preserving PDE's"
    1.26 + #    (D. Tschumperle)
    1.27 + #    International Journal of Computer Vision, May 2006.
    1.28 + #
    1.29 + #    "Vector-Valued Image Regularization with PDE's : A Common Framework
    1.30 + #    for Different Applications"
    1.31 + #    (D. Tschumperle, R. Deriche).
    1.32 + #    IEEE Transactions on Pattern Analysis and Machine Intelligence,
    1.33 + #    Vol 27, No 4, pp 506-517, April 2005.
    1.34 + #
    1.35 + #  Copyright   : David Tschumperle
    1.36 + #                ( http://www.greyc.ensicaen.fr/~dtschump/ )
    1.37 + #
    1.38 + #  License     : CeCILL v2.0
    1.39 + #                ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )
    1.40 + #
    1.41 + #  This software is governed by the CeCILL  license under French law and
    1.42 + #  abiding by the rules of distribution of free software.  You can  use,
    1.43 + #  modify and/ or redistribute the software under the terms of the CeCILL
    1.44 + #  license as circulated by CEA, CNRS and INRIA at the following URL
    1.45 + #  "http://www.cecill.info".
    1.46 + #
    1.47 + #  As a counterpart to the access to the source code and  rights to copy,
    1.48 + #  modify and redistribute granted by the license, users are provided only
    1.49 + #  with a limited warranty  and the software's author,  the holder of the
    1.50 + #  economic rights,  and the successive licensors  have only  limited
    1.51 + #  liability.
    1.52 + #
    1.53 + #  In this respect, the user's attention is drawn to the risks associated
    1.54 + #  with loading,  using,  modifying and/or developing or reproducing the
    1.55 + #  software by the user in light of its specific status of free software,
    1.56 + #  that may mean  that it is complicated to manipulate,  and  that  also
    1.57 + #  therefore means  that it is reserved for developers  and  experienced
    1.58 + #  professionals having in-depth computer knowledge. Users are therefore
    1.59 + #  encouraged to load and test the software's suitability as regards their
    1.60 + #  requirements in conditions enabling the security of their systems and/or
    1.61 + #  data to be ensured and,  more generally, to use and operate it in the
    1.62 + #  same conditions as regards security.
    1.63 + #
    1.64 + #  The fact that you are presently reading this means that you have had
    1.65 + #  knowledge of the CeCILL license and that you accept its terms.
    1.66 + #
    1.67 +*/
    1.68 +
    1.69 +#define cimg_plugin "plugins/greycstoration.h"
    1.70 +#ifndef cimg_debug
    1.71 +#if defined(_MSC_VER) || defined(WIN32)  || defined(_WIN32) || defined(__WIN32__) \
    1.72 + || defined(WIN64)    || defined(_WIN64) || defined(__WIN64__)
    1.73 +#define cimg_debug 2
    1.74 +#else
    1.75 +#define cimg_debug 1
    1.76 +#endif
    1.77 +#endif
    1.78 +#include "CImg.h"
    1.79 +#if cimg_OS!=2
    1.80 +#include <pthread.h>
    1.81 +#endif
    1.82 +#define gprintf if (verbose) std::fprintf
    1.83 +using namespace cimg_library;
    1.84 +
    1.85 +// The lines below are necessary when using a non-standard compiler as visualcpp6.
    1.86 +#ifdef cimg_use_visualcpp6
    1.87 +#define std
    1.88 +#endif
    1.89 +#ifdef min
    1.90 +#undef min
    1.91 +#undef max
    1.92 +#endif
    1.93 +
    1.94 +//-----------
    1.95 +// get_geom() : read geometry from a string (for instance '320x256' or '200%x200%').
    1.96 +//-----------
    1.97 +void get_geom(const char *geom, int &geom_w, int &geom_h) {
    1.98 +  char tmp[16];
    1.99 +  std::sscanf(geom,"%d%7[^0-9]%d%7[^0-9]",&geom_w,tmp,&geom_h,tmp+1);
   1.100 +  if (tmp[0]=='%') geom_w=-geom_w;
   1.101 +  if (tmp[1]=='%') geom_h=-geom_h;
   1.102 +}
   1.103 +
   1.104 +//--------------------------
   1.105 +// GREYCstoration main code
   1.106 +//--------------------------
   1.107 +template<typename T> void greycstoration(int argc, char **argv, T pixel_type) {
   1.108 +  pixel_type = (T)0;
   1.109 +  cimg_usage(" Open Source Algorithms for Image Denoising and Interpolation.");
   1.110 +  cimg_help("-------------------------------------------------------------------------\n"
   1.111 +            " GREYCstoration v3.0, by David Tschumperle.                              \n"
   1.112 +            " ------------------------------------------------------------------------\n"
   1.113 +            " This program allows to denoise, inpaint and resize 2D color images.     \n"
   1.114 +            " It is the result of the research work done in the IMAGE group of the    \n"
   1.115 +            " GREYC Lab (CNRS, UMR 6072) (http://www.greyc.ensicaen.fr/EquipeImage/)  \n"
   1.116 +            " by David Tschumperle (http://www.greyc.ensicaen.fr/~dtschump/). This    \n"
   1.117 +            " program has been primarily released to help people processing image data\n"
   1.118 +            " and to allow comparisons between regularization algorithms. This is an  \n"
   1.119 +            " open source software, distributed within the CImg Library package       \n"
   1.120 +            " (http://cimg.sourceforge.net), and submitted to the CeCILL License. If  \n"
   1.121 +            " you are interested to distribute it in a closed-source product, please  \n"
   1.122 +            " read the licence file carefully. If you are using 'GREYCstored' images  \n"
   1.123 +            " in your own publications, be kind to reference the GREYCstoration web   \n"
   1.124 +            " site or the related scientific papers. More informations available at : \n"
   1.125 +            "            ** http://cimg.sourceforge.net/greycstoration/ **            \n"
   1.126 +            "-------------------------------------------------------------------------\n");
   1.127 +
   1.128 +  // Read Global parameters
   1.129 +  //------------------------
   1.130 +  cimg_help("  + Global parameters\n    -----------------------");
   1.131 +  const char        *restore_mode    = cimg_option("-restore",(char*)0,"Restore image specified after '-restore'");
   1.132 +  const char        *inpaint_mode    = cimg_option("-inpaint",(char*)0,"Inpaint image specified after '-inpaint'");
   1.133 +  const char        *resize_mode     = cimg_option("-resize",(char*)0,"Resize image specified after '-resize'");
   1.134 +  const char        *clean_mode      = cimg_option("-clean",(char*)0,"Clean image specified after '-clean'");
   1.135 +  const char        *reference_image = cimg_option("-ref",(char*)0,"Reference image to compare with");
   1.136 +  const char         nb_bits         = cimg_option("-bits",8,"Define input value type (8='uchar', 16='ushort', 32='float')");
   1.137 +  const unsigned int value_factor    = cimg_option("-fact",0,"Define input value factor (0=auto)");
   1.138 +  const float        noise_g         = cimg_option("-ng",0.0f,"Add Gaussian noise before denoising");
   1.139 +  const float        noise_u         = cimg_option("-nu",0.0f,"Add Uniform noise before denoising");
   1.140 +  const float        noise_s         = cimg_option("-ns",0.0f,"Add Salt&Pepper noise before denoising");
   1.141 +  const unsigned int color_base      = cimg_option("-cbase",0,"Processed color base (0=RGB, 1=YCbCr)");
   1.142 +  const char        *channel_range   = cimg_option("-crange",(char*)0,"Processed channels (ex: '0-2')");
   1.143 +  const unsigned int saving_step     = cimg_option("-save",0,"Iteration saving step (>=0)");
   1.144 +  const bool         visu            = cimg_option("-visu",cimg_display?true:false,"Enable/Disable visualization windows (0 or 1)");
   1.145 +  const char        *file_o          = cimg_option("-o",(char*)0,"Filename of output image");
   1.146 +  const bool         append_result   = cimg_option("-append",false,"Append images in output file");
   1.147 +  const bool         verbose         = cimg_option("-verbose",true,"Verbose mode");
   1.148 +  const unsigned int jpg_quality     = cimg_option("-quality",100,"Output compression quality (if JPEG format)");
   1.149 +  unsigned int       nb_iterations   = cimg_option("-iter",(restore_mode||clean_mode)?1:(inpaint_mode?1000:3),
   1.150 +                                                   "Number of iterations (>0)");
   1.151 +  const float        sdt             = cimg_option("-sharp",(restore_mode||clean_mode)?0.0f:(inpaint_mode?0.0f:10.0f),
   1.152 +                                                   "Sharpening strength (activate sharpening filter, >=0)");
   1.153 +  const float        sp              = cimg_option("-se",(restore_mode||clean_mode)?0.5f:(inpaint_mode?0.5f:3.0f),
   1.154 +                                                   "Sharpening edge threshold (>=0)");
   1.155 +  const unsigned int tile_size       = cimg_option("-tile",512,"Activate tiled mode and set tile size (>=0)");
   1.156 +  const unsigned int tile_border     = cimg_option("-btile",4,"Size of tile borders (>=0)");
   1.157 +  const unsigned int nb_threads      = cimg_option("-threads",1,"Number of threads used (*experimental*, tile mode only, >0)");
   1.158 +  const bool         fast_approx     = cimg_option("-fast",true,"Try faster algorithm (true or false)");
   1.159 +
   1.160 +  // Declare specific algorithm parameters
   1.161 +  //--------------------------------------
   1.162 +  float amplitude = 0, sharpness = 0, anisotropy = 0, alpha = 0, sigma = 0, gauss_prec = 0, dl = 0, da = 0, sigma_s = 0, sigma_p = 0;
   1.163 +  unsigned int interpolation = 0, patch_size = 0, lookup_size = 0;
   1.164 +
   1.165 +  if (argc==1 ||
   1.166 +      (!restore_mode && !inpaint_mode && !resize_mode && !clean_mode) ||
   1.167 +      (restore_mode && inpaint_mode) || (restore_mode && resize_mode) || (restore_mode && clean_mode) ||
   1.168 +      (inpaint_mode && resize_mode) || (inpaint_mode && clean_mode)) {
   1.169 +    std::fprintf(stderr,"\n%s : You must specify (only) one of the '-restore', '-inpaint', '-resize' or '-clean' options.\n"
   1.170 +                 "(try option '-h', '-h -restore','-h -inpaint', '-h -resize' or '-h -clean' to get options relative to specific actions\n\n",
   1.171 +                 cimg::basename(argv[0]));
   1.172 +    std::exit(0);
   1.173 +  }
   1.174 +
   1.175 +  // Init variables
   1.176 +  //----------------
   1.177 +  CImg<T> img0, img, imgr;
   1.178 +  CImg<unsigned char> mask;
   1.179 +  CImgDisplay disp;
   1.180 +
   1.181 +  // Read specific parameters for image restoration
   1.182 +  //------------------------------------------------
   1.183 +  if (restore_mode) {
   1.184 +    cimg_help("\n  + Restoration mode parameters\n    ---------------------------");
   1.185 +    amplitude      = cimg_option("-dt",40.0f,"Regularization strength per iteration (>=0)");
   1.186 +    sharpness      = cimg_option("-p",0.9f,"Contour preservation (>=0)");
   1.187 +    anisotropy     = cimg_option("-a",0.15f,"Smoothing anisotropy (0<=a<=1)");
   1.188 +    alpha          = cimg_option("-alpha",0.6f,"Noise scale (>=0)");
   1.189 +    sigma          = cimg_option("-sigma",1.1f,"Geometry regularity (>=0)");
   1.190 +    gauss_prec     = cimg_option("-prec",2.0f,"Computation precision (>0)");
   1.191 +    dl             = cimg_option("-dl",0.8f,"Spatial integration step (0<=dl<=1)");
   1.192 +    da             = cimg_option("-da",30.0f,"Angular integration step (0<=da<=90)");
   1.193 +    interpolation  = cimg_option("-interp",0,"Interpolation type (0=Nearest-neighbor, 1=Linear, 2=Runge-Kutta)");
   1.194 +
   1.195 +    gprintf(stderr,"- Image Restoration mode\n");
   1.196 +    if (!cimg::strcmp("-restore",restore_mode)) {
   1.197 +      std::fprintf(stderr,"%s : You must specify a valid input image filename after the '-restore' flag.\n\n",cimg::basename(argv[0]));
   1.198 +      std::exit(0);
   1.199 +    }
   1.200 +    gprintf(stderr,"- Load input image '%s'...",cimg::basename(restore_mode));
   1.201 +    img.load(restore_mode);
   1.202 +    gprintf(stderr,"\r- Input image : '%s' (size %dx%d, value range [%g,%g])\n",
   1.203 +            cimg::basename(restore_mode),img.dimx(),img.dimy(),(double)img.min(),(double)img.max());
   1.204 +    if (noise_g || noise_u || noise_s) {
   1.205 +      img0 = img;
   1.206 +      img.noise(noise_g,0).noise(noise_u,1).noise(noise_s,2);
   1.207 +      gprintf(stderr,"\r- Noisy image : value range [%g,%g], PSNR Noisy / Original : %g\n",
   1.208 +              (double)img.min(),(double)img.max(),img.PSNR(img0));
   1.209 +    }
   1.210 +  }
   1.211 +
   1.212 +  // Specific parameters for image inpainting
   1.213 +  //-----------------------------------------
   1.214 +  if (inpaint_mode) {
   1.215 +    cimg_help("\n  + Inpainting mode parameters\n    --------------------------");
   1.216 +    const char *file_m        = cimg_option("-m",(char*)0,"Input inpainting mask");
   1.217 +    const unsigned int dilate = cimg_option("-dilate",0,"Inpainting mask dilatation (>=0)");
   1.218 +    const unsigned int init   = cimg_option("-init",4,"Inpainting init (0=black, 1=white, 2=noise, 3=unchanged, 4=smart)");
   1.219 +    amplitude                 = cimg_option("-dt",20.0f,"Regularization strength per iteration (>=0)");
   1.220 +    sharpness                 = cimg_option("-p",0.3f,"Contour preservation (>=0)");
   1.221 +    anisotropy                = cimg_option("-a",1.0f,"Smoothing anisotropy (0<=a<=1)");
   1.222 +    alpha                     = cimg_option("-alpha",0.8f,"Noise scale (>=0)");
   1.223 +    sigma                     = cimg_option("-sigma",2.0f,"Geometry regularity (>=0)");
   1.224 +    gauss_prec                = cimg_option("-prec",2.0f,"Computation precision (>0)");
   1.225 +    dl                        = cimg_option("-dl",0.8f,"Spatial integration step (0<=dl<=1)");
   1.226 +    da                        = cimg_option("-da",30.0f,"Angular integration step (0<=da<=90)");
   1.227 +    interpolation             = cimg_option("-interp",0,"Interpolation type (0=Nearest-neighbor, 1=Linear, 2=Runge-Kutta)");
   1.228 +
   1.229 +    gprintf(stderr,"- Image Inpainting mode\n");
   1.230 +    if (!cimg::strcmp("-inpaint",inpaint_mode)) {
   1.231 +      std::fprintf(stderr,"%s : You must specify a valid input image filename after the '-inpaint' flag.\n\n",
   1.232 +                   cimg::basename(argv[0]));
   1.233 +      std::exit(0);
   1.234 +    }
   1.235 +    gprintf(stderr,"- Load input image '%s'...",cimg::basename(inpaint_mode));
   1.236 +    img.load(inpaint_mode);
   1.237 +    gprintf(stderr,"\r- Input image : '%s' (size %dx%d, value range [%g,%g])\n",
   1.238 +            cimg::basename(inpaint_mode),img.dimx(),img.dimy(),(double)img.min(),(double)img.max());
   1.239 +    if (noise_g || noise_u || noise_s) {
   1.240 +      img0 = img;
   1.241 +      img.noise(noise_g,0).noise(noise_u,1).noise(noise_s,2);
   1.242 +      gprintf(stderr,"\r- Noisy image : value range [%g,%g], PSNR Noisy / Original : %g\n",
   1.243 +              (double)img.min(),(double)img.max(),img.PSNR(img0));
   1.244 +    }
   1.245 +    if (!file_m) {
   1.246 +      std::fprintf(stderr,"%s : You need to specify a valid inpainting mask filename after the '-m' flag.\n\n",
   1.247 +                   cimg::basename(argv[0]));
   1.248 +      std::exit(0);
   1.249 +    }
   1.250 +    if (cimg::strncasecmp("block",file_m,5)) {
   1.251 +      gprintf(stderr,"- Load inpainting mask '%s'...",cimg::basename(file_m));
   1.252 +      mask.load(file_m);
   1.253 +      gprintf(stderr,"\r- Inpainting mask : '%s' (size %dx%d)\n",
   1.254 +              cimg::basename(file_m),mask.dimx(),mask.dimy());
   1.255 +    }
   1.256 +    else {
   1.257 +      unsigned int l = 16; std::sscanf(file_m,"block%u",&l);
   1.258 +      mask.assign(img.dimx()/l,img.dimy()/l);
   1.259 +      cimg_forXY(mask,x,y) mask(x,y) = (x+y)%2;
   1.260 +      img0 = img;
   1.261 +    }
   1.262 +    mask.resize(img.dimx(),img.dimy(),1,1);
   1.263 +    if (dilate) mask.dilate(dilate);
   1.264 +    switch (init) {
   1.265 +    case 0 : { cimg_forXYV(img,x,y,k) if (mask(x,y)) img(x,y,k) = 0; } break;
   1.266 +    case 1 : { cimg_forXYV(img,x,y,k) if (mask(x,y)) img(x,y,k) = 255; } break;
   1.267 +    case 2 : { cimg_forXYV(img,x,y,k) if (mask(x,y)) img(x,y,k) = (T)(255*cimg::rand()); } break;
   1.268 +    case 3 : break;
   1.269 +    default: {
   1.270 +      typedef unsigned char uchar;
   1.271 +      CImg<unsigned char> tmask(mask), ntmask(tmask);
   1.272 +      CImg_3x3(M,uchar); Mpp = Mnp = Mpn = Mnn = 0;
   1.273 +      CImg_3x3(I,T); Ipp = Inp = Icc = Ipn = Inn = 0;
   1.274 +      while (ntmask.max()>0) {
   1.275 +        cimg_for3x3(tmask,x,y,0,0,M) if (Mcc && (!Mpc || !Mnc || !Mcp || !Mcn)) {
   1.276 +          const float
   1.277 +            ccp = Mcp?0.0f:1.0f, cpc = Mpc?0.0f:1.0f,
   1.278 +            cnc = Mnc?0.0f:1.0f, ccn = Mcn?0.0f:1.0f,
   1.279 +            csum = ccp + cpc + cnc + ccn;
   1.280 +          cimg_forV(img,k) {
   1.281 +            cimg_get3x3(img,x,y,0,k,I);
   1.282 +            img(x,y,k) = (T)((ccp*Icp + cpc*Ipc + cnc*Inc + ccn*Icn)/csum);
   1.283 +          }
   1.284 +          ntmask(x,y) = 0;
   1.285 +        }
   1.286 +        tmask = ntmask;
   1.287 +      }
   1.288 +    } break;
   1.289 +    }
   1.290 +  }
   1.291 +
   1.292 +  // Specific parameters for image resizing
   1.293 +  //----------------------------------------
   1.294 +  if (resize_mode) {
   1.295 +    cimg_help("\n  + Resizing mode parameters\n    ------------------------");
   1.296 +    const char *geom0         = cimg_option("-g0",(char*)0,"Input image geometry");
   1.297 +    const char *geom          = cimg_option("-g",(char*)0,"Output image geometry");
   1.298 +    const bool anchor         = cimg_option("-anchor",true,"Anchor original pixels (keep their original values)");
   1.299 +    const unsigned int init   = cimg_option("-init",3,"Initial estimate (1=block, 3=linear, 4=Moving average, 5=bicubic)");
   1.300 +    amplitude                 = cimg_option("-dt",20.0f,"Regularization strength per iteration (>=0)");
   1.301 +    sharpness                 = cimg_option("-p",0.2f,"Contour preservation (>=0)");
   1.302 +    anisotropy                = cimg_option("-a",0.9f,"Smoothing anisotropy (0<=a<=1)");
   1.303 +    alpha                     = cimg_option("-alpha",0.1f,"Noise scale (>=0)");
   1.304 +    sigma                     = cimg_option("-sigma",1.5f,"Geometry regularity (>=0)");
   1.305 +    gauss_prec                = cimg_option("-prec",2.0f,"Computation precision (>0)");
   1.306 +    dl                        = cimg_option("-dl",0.8f,"Spatial integration step (0<=dl<=1)");
   1.307 +    da                        = cimg_option("-da",30.0f,"Angular integration step (0<=da<=90)");
   1.308 +    interpolation             = cimg_option("-interp",0,"Interpolation type (0=Nearest-neighbor, 1=Linear, 2=Runge-Kutta)");
   1.309 +
   1.310 +    gprintf(stderr,"- Image Resizing mode\n");
   1.311 +    if (!geom && !geom0) {
   1.312 +      std::fprintf(stderr,"%s : You need to specify an output geometry or an input geometry (option -g or -g0).\n\n",
   1.313 +                   cimg::basename(argv[0]));
   1.314 +      std::exit(0);
   1.315 +    }
   1.316 +    if (!cimg::strcmp("-resize",resize_mode)) {
   1.317 +      std::fprintf(stderr,"%s : You must specify a valid input image filename after the '-resize' flag.\n\n",
   1.318 +                   cimg::basename(argv[0]));
   1.319 +      std::exit(0);
   1.320 +    }
   1.321 +    gprintf(stderr,"- Load input image '%s'...",cimg::basename(resize_mode));
   1.322 +    img.load(resize_mode);
   1.323 +    gprintf(stderr,"\r- Input image : '%s' (size %dx%d, value range [%g,%g])\n",
   1.324 +            cimg::basename(resize_mode),img.dimx(),img.dimy(),(double)img.min(),(double)img.max());
   1.325 +    if (noise_g || noise_u || noise_s) {
   1.326 +      img0 = img;
   1.327 +      img.noise(noise_g,0).noise(noise_u,1).noise(noise_s,2);
   1.328 +      gprintf(stderr,"\r- Noisy image : value range [%g,%g], PSNR Noisy / Original : %g\n",
   1.329 +              (double)img.min(),(double)img.max(),img.PSNR(img0));
   1.330 +    }
   1.331 +    int w, h;
   1.332 +    if (geom0) {
   1.333 +      int w0, h0;
   1.334 +      get_geom(geom0,w0,h0);
   1.335 +      w0 = w0>0?w0:-w0*img.dimx()/100;
   1.336 +      h0 = h0>0?h0:-h0*img.dimy()/100;
   1.337 +      gprintf(stderr,"- Reducing geometry to %dx%d using %s interpolation.\n",w0,h0,
   1.338 +              init==1?"bloc":(init==3?"linear":(init==5?"bicubic":"moving average")));
   1.339 +      img0.assign(img);
   1.340 +      w = img.dimx();
   1.341 +      h = img.dimy();
   1.342 +      img.resize(w0,h0,-100,-100,5);
   1.343 +    }
   1.344 +    if (geom) {
   1.345 +      get_geom(geom,w,h);
   1.346 +      w = w>0?w:-w*img.dimx()/100;
   1.347 +      h = h>0?h:-h*img.dimy()/100;
   1.348 +    }
   1.349 +    mask.assign(img.dimx(),img.dimy(),1,1,255);
   1.350 +    if (!anchor) mask.resize(w,h,1,1,1); else mask = !mask.resize(w,h,1,1,4);
   1.351 +    img.resize(w,h,1,-100,init);
   1.352 +    if (img0) { gprintf(stderr,"\r- PSNR Original / Thumbnail : %g\n",img.PSNR(img0)); }
   1.353 +  }
   1.354 +
   1.355 +  // Specific parameters for image cleaning
   1.356 +  //----------------------------------------
   1.357 +  if (clean_mode) {
   1.358 +    cimg_help("\n  + Cleaning mode parameters\n    ------------------------");
   1.359 +    patch_size  = cimg_option("-p",4,"Patch size (>0)");
   1.360 +    sigma_s     = cimg_option("-ss",15.0f,"Spatial sigma (>0)");
   1.361 +    sigma_p     = cimg_option("-sp",10.0f,"Patch sigma (>0)");
   1.362 +    lookup_size = cimg_option("-r",7,"Size of the lookup region (>0)");
   1.363 +
   1.364 +    gprintf(stderr,"- Image Cleaning mode\n");
   1.365 +    if (!cimg::strcmp("-clean",clean_mode)) {
   1.366 +      std::fprintf(stderr,"%s : You must specify a valid input image filename after the '-clean' flag.\n\n",
   1.367 +                   cimg::basename(argv[0]));
   1.368 +      std::exit(0);
   1.369 +    }
   1.370 +    gprintf(stderr,"- Load input image '%s'...",cimg::basename(clean_mode));
   1.371 +    img.load(clean_mode);
   1.372 +    gprintf(stderr,"\r- Input image : '%s' (size %dx%d, value range [%g,%g])\n",
   1.373 +            cimg::basename(clean_mode),img.dimx(),img.dimy(),(double)img.min(),(double)img.max());
   1.374 +    if (noise_g || noise_u || noise_s) {
   1.375 +      img0 = img;
   1.376 +      img.noise(noise_g,0).noise(noise_u,1).noise(noise_s,2);
   1.377 +      gprintf(stderr,"\r- Noisy image : value range [%g,%g], PSNR Noisy / Original : %g\n",
   1.378 +              (double)img.min(),(double)img.max(),img.PSNR(img0));
   1.379 +    }
   1.380 +  }
   1.381 +
   1.382 +  // Load reference image if any specified
   1.383 +  //--------------------------------------
   1.384 +  if (reference_image) {
   1.385 +    gprintf(stderr,"- Load reference image '%s'...",cimg::basename(reference_image));
   1.386 +    imgr.load(reference_image);
   1.387 +    gprintf(stderr,"\r- Reference image : '%s' (size %dx%d, value range [%g,%g])",
   1.388 +            cimg::basename(reference_image),imgr.dimx(),imgr.dimy(),(double)imgr.min(),(double)imgr.max());
   1.389 +    if (img0) { imgr.resize(img0); gprintf(stderr,", PSNR Reference / Original : %g dB\n",imgr.PSNR(img0)); }
   1.390 +    else { imgr.resize(img); gprintf(stderr,"\n"); }
   1.391 +  }
   1.392 +
   1.393 +  // Init images and display
   1.394 +  //-------------------------
   1.395 +  CImg<T> dest(img);
   1.396 +  unsigned int crange_beg = 0, crange_end = dest.dimv()-1U;
   1.397 +  if (color_base) {
   1.398 +    switch (nb_bits) {
   1.399 +    case 8: dest.RGBtoYCbCr(); break;
   1.400 +    case 16: (dest/=256).RGBtoYCbCr(); break;
   1.401 +    default: std::fprintf(stderr,"\n%s : YCbCr color base is not authorized for 32bits float-valued images.\n\n",
   1.402 +                          cimg::basename(argv[0])); std::exit(0);
   1.403 +    }
   1.404 +  }
   1.405 +  if (channel_range) std::sscanf(channel_range,"%u%*c%u",&crange_beg,&crange_end);
   1.406 +  gprintf(stderr,"- Color base : %s, Channels range : %u-%u\n",color_base?"YCbCr":"RGB",crange_beg,crange_end);
   1.407 +  if (!visu && !append_result) img.assign();
   1.408 +  if (visu) {
   1.409 +    const int sx = 2*CImgDisplay::screen_dimx()/3, sy = 2*CImgDisplay::screen_dimy()/3;
   1.410 +    int nwidth = dest.dimx(), nheight = dest.dimy();
   1.411 +    if (nwidth>sx) { nheight = nheight*sx/nwidth; nwidth = sx; }
   1.412 +    if (nheight>sy) { nwidth = nwidth*sy/nheight; nheight = sy; }
   1.413 +    disp.assign(nwidth,nheight,"GREYCstoration");
   1.414 +    if (color_base) {
   1.415 +      if (nb_bits==16) (dest.get_YCbCrtoRGB()*=256).display(disp);
   1.416 +      else dest.get_YCbCrtoRGB().display(disp);
   1.417 +    }
   1.418 +    else dest.display(disp);
   1.419 +  }
   1.420 +  const float gfact = (value_factor>0)?value_factor:((sizeof(T)==2)?1.0f/256:1.0f);
   1.421 +
   1.422 +  //---------------------------------
   1.423 +  // Begin GREYCstoration iterations
   1.424 +  //---------------------------------
   1.425 +  bool stop_all = false;
   1.426 +  for (unsigned int iter=0; iter<nb_iterations && !stop_all; iter++) {
   1.427 +    bool stop_iteration = false;
   1.428 +
   1.429 +    // Run one iteration of the GREYCstoration filter
   1.430 +    //------------------------------------------------
   1.431 +    CImg<T> dest_range = dest.get_shared_channels(crange_beg,crange_end);
   1.432 +    if (restore_mode)
   1.433 +      dest_range.greycstoration_run(amplitude,sharpness,anisotropy,alpha,sigma,gfact,dl,da,gauss_prec,
   1.434 +                                    interpolation,fast_approx,tile_size,tile_border,nb_threads);
   1.435 +    if (clean_mode)
   1.436 +      dest_range.greycstoration_patch_run(patch_size,sigma_s,sigma_p,lookup_size,
   1.437 +                                          tile_size,tile_border,nb_threads,fast_approx);
   1.438 +    if (inpaint_mode || resize_mode)
   1.439 +      dest_range.greycstoration_run(mask,amplitude,sharpness,anisotropy,alpha,sigma,gfact,dl,da,gauss_prec,
   1.440 +                                    interpolation,fast_approx,tile_size,tile_border,nb_threads);
   1.441 +    do {
   1.442 +      const unsigned int progress = (unsigned int)dest_range.greycstoration_progress();
   1.443 +      gprintf(stderr,"\r- Processing : Iteration %u/%u (%u%%)\t\t",1+iter,nb_iterations,progress);
   1.444 +      if (disp) {
   1.445 +        if (disp.is_resized) disp.resize();
   1.446 +        disp.set_title("Processing : Iteration %u/%u (%u%%)",1+iter,nb_iterations,progress);
   1.447 +        if (disp.is_closed || disp.is_keyQ || disp.is_keyESC) {
   1.448 +          dest_range.greycstoration_stop();
   1.449 +          stop_iteration = true;
   1.450 +          iter = nb_iterations-1;
   1.451 +        }
   1.452 +      }
   1.453 +      cimg::wait(200);
   1.454 +    } while (dest_range.greycstoration_is_running());
   1.455 +    if (!stop_iteration && sdt>0) dest_range.sharpen(sdt,sp,alpha/3,sigma/3);
   1.456 +    dest_range.cut(cimg::type<T>::min(),cimg::type<T>::max());
   1.457 +
   1.458 +    // Prepare for next iteration
   1.459 +    //---------------------------
   1.460 +    CImg<T> tmp_rgb = color_base?(nb_bits==16?dest.get_YCbCrtoRGB()*=256:dest.get_YCbCrtoRGB()):CImg<T>(),
   1.461 +      &dest_rgb = color_base?tmp_rgb:dest;
   1.462 +    if (disp && visu) dest_rgb.display(disp);
   1.463 +    if (file_o && saving_step && !(iter%saving_step)) dest_rgb.save(file_o,iter);
   1.464 +
   1.465 +    // Display result and allows user interaction if needed.
   1.466 +    //-------------------------------------------------------
   1.467 +    if (iter==nb_iterations-1) {
   1.468 +      gprintf(stderr,"\r- Processing : Done !                          \n");
   1.469 +      if (img0) { gprintf(stderr,"- PSNR Restored / Original : %g dB\n",dest_rgb.PSNR(img0)); }
   1.470 +      if (disp) {
   1.471 +        static bool first_time = true;
   1.472 +        if (first_time) {
   1.473 +          first_time = false;
   1.474 +          gprintf(stderr,
   1.475 +                  "- GREYCstoration interface :\n"
   1.476 +                  " > You can now zoom to a particular rectangular region,\n"
   1.477 +                  "   or press one of the following key on the display window :\n"
   1.478 +                  "   SPACE : Swap views.\n"
   1.479 +                  "   S     : Save a snapshot of the current image.\n"
   1.480 +                  "   I     : Run another iteration.\n"
   1.481 +                  "   Q     : Quit GREYCstoration.\n");
   1.482 +        }
   1.483 +
   1.484 +        CImgList<T> visu;
   1.485 +        visu.insert(img0,~0,true).insert(img,~0,true).insert(dest_rgb,~0,true).insert(imgr,~0U,true);
   1.486 +        const char *titles[4] = { "original", "noisy", "restored", "reference"};
   1.487 +        unsigned int visupos = 2;
   1.488 +        CImgDisplay dispz;
   1.489 +        CImg<T> zoom;
   1.490 +        int snb = 0;
   1.491 +        bool stop_interact = false;
   1.492 +        while (!stop_interact) {
   1.493 +          disp.show().set_title("GREYCstoration (%s)",titles[visupos]);
   1.494 +          const CImg<int> s = visu(visupos).get_select(disp,2);
   1.495 +          if (disp.is_closed) stop_interact = true;
   1.496 +          switch (disp.key) {
   1.497 +          case cimg::keySPACE: do { visupos = (visupos+1)%visu.size; } while (!visu(visupos)); break;
   1.498 +          case cimg::keyBACKSPACE: do { visupos = (visupos-1+visu.size)%visu.size; } while (!visu(visupos)); break;
   1.499 +          case cimg::keyQ: stop_interact = stop_all = true; break;
   1.500 +          case cimg::keyI:
   1.501 +            stop_interact = true;
   1.502 +            gprintf(stderr,"- Perform iteration %u...\n",++nb_iterations);
   1.503 +            dest_rgb.display(disp);
   1.504 +            break;
   1.505 +          case cimg::keyS:
   1.506 +            if (!snb) {
   1.507 +              if (!append_result) dest_rgb.save(file_o?file_o:"GREYCstoration.bmp");
   1.508 +              else CImgList<T>(img,dest_rgb).get_append('x').save(file_o?file_o:"GREYCstoration.bmp");
   1.509 +            }
   1.510 +            if (zoom) zoom.save(file_o?file_o:"GREYCstoration.bmp",snb);
   1.511 +            gprintf(stderr,"- Snapshot %u : '%s' saved\n",snb++,file_o?file_o:"GREYCstoration.bmp");
   1.512 +            break;
   1.513 +          }
   1.514 +          disp.key = 0;
   1.515 +          if (disp.is_resized) disp.resize().display(visu(visupos));
   1.516 +          if (dispz && dispz.is_resized) dispz.resize().display(zoom);
   1.517 +          if (dispz && dispz.is_closed) dispz.assign();
   1.518 +
   1.519 +          if (s[0]>=0 && s[1]>=0 && s[3]>=0 && s[4]>=0) {
   1.520 +            const int x0 = s[0], y0 = s[1], x1 = s[3], y1 = s[4];
   1.521 +            if (cimg::abs(x0-x1)>4 && cimg::abs(y0-y1)>4) {
   1.522 +              CImgList<T> tmp(img.get_crop(x0,y0,x1,y1), dest_rgb.get_crop(x0,y0,x1,y1));
   1.523 +              if (img0) tmp.insert(img0.get_crop(x0,y0,x1,y1),0);
   1.524 +              if (imgr) tmp.insert(imgr.get_crop(x0,y0,x1,y1));
   1.525 +              zoom = tmp.get_append('x','c');
   1.526 +              if (!dispz) {
   1.527 +                const int sx = 5*CImgDisplay::screen_dimx()/6, sy = 5*CImgDisplay::screen_dimy()/6;
   1.528 +                int nwidth = zoom.dimx(), nheight = zoom.dimy();
   1.529 +                if (nwidth>nheight) { nheight = nheight*sx/nwidth; nwidth = sx; }
   1.530 +                else { nwidth = nwidth*sy/nheight; nheight = sy; }
   1.531 +                dispz.assign(zoom.get_resize(nwidth,nheight));
   1.532 +                dispz.set_title("GREYCstoration (zoom) : - %s %s %s %s",
   1.533 +                                img0?"original -":"",
   1.534 +                                img?"noisy -":"",
   1.535 +                                dest?"restored -":"",
   1.536 +                                imgr?"reference -":"");
   1.537 +              } else dispz.resize(dispz.dimx(),dispz.dimx()*zoom.dimy()/zoom.dimx(),false);
   1.538 +              dispz.display(zoom).show();
   1.539 +            }
   1.540 +          }
   1.541 +        }
   1.542 +      }
   1.543 +    }
   1.544 +  }
   1.545 +
   1.546 +  // Save result and exit
   1.547 +  //----------------------
   1.548 +  if (file_o) {
   1.549 +    CImg<T> tmp_rgb = color_base?(nb_bits==16?dest.get_YCbCrtoRGB()*=256:dest.get_YCbCrtoRGB()):CImg<T>(),
   1.550 +      &dest_rgb = color_base?tmp_rgb:dest;
   1.551 +    if (jpg_quality) {
   1.552 +      gprintf(stderr,"\n- Saving output image '%s' (JPEG quality = %u%%)\n",file_o,jpg_quality);
   1.553 +      if (!append_result) dest_rgb.save_jpeg(file_o,jpg_quality);
   1.554 +      else CImgList<T>(img,dest_rgb).get_append('x').save_jpeg(file_o,jpg_quality);
   1.555 +    } else {
   1.556 +      gprintf(stderr,"\n- Saving output image '%s'\n",file_o);
   1.557 +      if (!append_result) dest_rgb.save(file_o);
   1.558 +      else CImgList<T>(img,dest_rgb).get_append('x').save(file_o);
   1.559 +    }
   1.560 +  }
   1.561 +  gprintf(stderr,"\n- Quit\n\n");
   1.562 +}
   1.563 +
   1.564 +
   1.565 +/*-----------------
   1.566 +  Main procedure
   1.567 +  ----------------*/
   1.568 +int main(int argc,char **argv) {
   1.569 +  const unsigned int color_base = cimg_option("-cbase",0,0);
   1.570 +  switch (cimg_option("-bits",8,0)) {
   1.571 +  case 32: { float pixel_type = 0; greycstoration(argc,argv,pixel_type); } break;
   1.572 +  case 16: {
   1.573 +    if (!color_base) { float pixel_type = 0; greycstoration(argc,argv,pixel_type); }
   1.574 +    else { unsigned short pixel_type = 0; greycstoration(argc,argv,pixel_type); }
   1.575 +  } break;
   1.576 +  default: { unsigned char pixel_type = 0; greycstoration(argc,argv,pixel_type); } break;
   1.577 +  }
   1.578 +  return 0;
   1.579 +}