#include "sopnamsp.h"
#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>

#include "simplesort.h"
#include "nbmath.h"
#include "histos.h"
#include "datatypes.h"
#include "imageop.h"
#include "nbtri.h"


//++
//  Module	Operations sur les images (C++)
//  Lib	Images++
//  include	imageop.h
//
//	Operations sur les images.
//--


// Outils prives.
template <class T>
float RzCmvSigmaCiel(T const * pix,int nx,int ny,int ndata,int bin,
		     int taille,float frac,int deg,
		     float low,float high,int deb);

#ifdef IMGRGCHECK
#define CHECKIMG2(_img_,_x_,_y_)                                  \
      if ((_x_ >= (_img_).siz_x) || (_y_ >= (_img_).siz_y) ||              \
          (_x_ < 0) || (_y_ < 0)) throw RangeCheckError("imageop.cc - CHECKIMG2") ;  \
      if (!(_img_).ImagePtr()) throw NullPtrError("imageop.cc - CHECKIMG2")
#else
#if defined(IMGVOIDCHECK)
#define CHECKIMG2(_img_,_x_,_y_) \
      if (!(_img_).ImagePtr()) THROW(NullPtrErr)
#else
#define CHECKIMG2(_img_,_x_,_y_)
#endif
#endif /* IMGRGCHECK */

// Les fonctions

/* Nouvelle-Fonction */

template <class T>
void CleanImagesOld(int nImg, Image<T>** imgTab, double nSig,
                 double lowbad, double highbad)

/*  Nettoie une serie d'image a partir de l'histo de la  */
/* serie de valeurs pour chaque pixel (Voir E. Aubourg)  */

{
  if (lowbad < MinRange((T)0))
    lowbad = MinRange((T)0);

  if (highbad > MaxRange((T)0))
    highbad = MaxRange((T)0);

  if (lowbad >= highbad)
    {lowbad = MinRange((T)0); highbad = MaxRange((T)0);}
  Histo h(lowbad, highbad, 500);

  for (int i=0; i<imgTab[0]->XSize(); i++)
    for (int j=0; j<imgTab[0]->YSize(); j++) {
      h.Zero();
      int k;
      for (k=0; k<nImg; k++)
        h.Add((*imgTab[k])(i,j));
      float x0, s0=0;
      x0 = h.CleanedMean(s0);
      float low = x0 - s0*nSig; if (low < lowbad) low = lowbad;
      float hi  = x0 + s0*nSig; // if (hi  > highbad) hi = highbad;
      for (k=0; k<nImg; k++)
        if ((*imgTab[k])(i,j) > hi || (*imgTab[k])(i,j) < low)
          (*imgTab[k])(i,j) = 0;
    }
}

/* Nouvelle-Fonction */

template <class T>
//++
void CleanImages(int nImg, Image<T>** imgTab, double nSig,
                 double lowbad, double highbad)
//
//	Nettoie une serie de `nImg' images `imgTab' a partir
//	d'une moyenne tronquee a `nSig' sigma
//	de la serie de valeurs pour chaque pixel.
//	Les valeurs des pixels sont prises entre `lowbad' et `highbad'.
//--
{
  if (lowbad < MinRange((T)0))
    lowbad = MinRange((T)0);

  if (highbad > MaxRange((T)0))
    highbad = MaxRange((T)0);

  if (lowbad >= highbad)
    {lowbad = MinRange((T)0); highbad = MaxRange((T)0);}
  
  double sx, sx2, sig, mean;
  int n;

  for (int i=0; i<imgTab[0]->XSize(); i++)
    for (int j=0; j<imgTab[0]->YSize(); j++) {
      sig = 1e20; mean = 0;
      for (int l = 0; l<3; l++) {
        double low = mean - nSig*sig;
        if (low < lowbad) low = lowbad;
        double high = mean + nSig*sig;
        if (high > highbad) high = highbad;
        sx = sx2 = n = 0;
        for (int k=0; k<nImg; k++) {
          double p = (*imgTab[k])(i,j);
          if (p > low && p < high) {
            sx += p;
            sx2 += p*p;
            n++;
          }
        }
        if (n>0) {
          mean = sx / n;
          sig  = sqrt( sx2/n - (sx/n)*(sx/n) );
        } // sinon on laisse tel quel, ca sera coupe...
      }
      double low = mean - nSig*sig;
      if (low < lowbad) low = lowbad;
      double high = mean + nSig*sig;
      if (high > highbad) high = highbad;
      for (int k=0; k<nImg; k++)
        if ((*imgTab[k])(i,j) > high || (*imgTab[k])(i,j) < low)
          (*imgTab[k])(i,j) = 0;
    }
}

/* Nouvelle-Fonction */

template <class T>
//++
float ImgFondCielHisto(Image<T> const& img, float binWidth, int degree, float frac,
                       float low, float high, int debug)
//
//	Calcul du fond de ciel par histogrammation. L'histogramme
//	va de `low' a `high' avec une largeur de bin `binWidth'.
//	La valeur la plus probalbe est fittee par un polynome de degre
//	`degree' en utilisant les bins entourant le maximum
//	a `frac' pourcent de sa hauteur.
//--
{
  // CHECKIMG2(img,0,0);

  if (degree < 2 || degree > 3)
    degree = 3;

  if (frac <= 0 || frac >= 1)
    frac = 0.5;

  if (low < MinRange((T)0)+2)
    low = MinRange((T)0)+2;

  if (high > MaxRange((T)0)-2)
    high = MaxRange((T)0)-2;

  if (low >= high)
    {low = MinRange((T)0); high = MaxRange((T)0);}

// verifier binwidth

//  if ( binWidth < 65536.0/img.XSize()*img.YSize() )
//    binWidth = 65536.0/img.XSize()*img.YSize();

  if (binWidth<0)
    binWidth = 1;

  if (debug > 1)
    cout << __FUNCTION__ << " binWidth =" << binWidth <<
      ", degree = " << degree <<
      ", frac = "   << frac   <<
      ", low = "    << low    <<
      ", high = "   << high   << endl;

  int nBins = int(((double)high-(double)low)/(double)binWidth + 1);
  Histo hist(low, high, nBins);

  int nPix = img.XSize()*img.YSize();
  for (int i=0; i<nPix; i++)
    hist.Add(img[i]);

  hist(0) = 0;

  float fd = hist.FitMax(degree, frac, debug);

  (float&) img.fond = fd;

  return fd;
}


/* Nouvelle-Fonction */

template <class T>
float ImgSigmaCiel(Image<T> const& img, int ndata, int binWidth,
		   int taille, float frac, int degree,
		   float low, float high, int debug)
{
  (float&) img.sigmaFond = RzCmvSigmaCiel(img.ImagePtr(), img.XSize(), img.YSize(),
				 ndata, binWidth, taille, frac,
				 degree, low, high, debug);
  return img.sigmaFond;
}


/* Nouvelle-Fonction */

template <class T>
//++
float RzCmvSigmaCiel(T const* pix,int nx,int ny,int ndata,int bin,
		     int taille,float frac,int deg,
		     float low,float high,int deb)
//
//	Calcul du fond de ciel par histogrammation
//
//| INPUT:
//|   - pix   : tableau pix[nx][ny]
//|   - nx,ny : dimensions du tableau pix
//|   - ndata : nombre de donnees a utiliser
//|   - bin   : largeur du bin en ADU
//|   - taille: 1/2 dimension de l'histogramme en ADU
//|   - frac  : fraction du maximum au dessus de laquelle on fite
//|   - deg   : degre fit polynome sous la gaussienne (0 ou 2)
//|   - low   : valeur basse limite en ADU
//|   - high  : valeur haute limite en ADU
//|   - deb   : flag de debug
//--
//++
//| OUTPUT:
//|   - retourne : le sigma du fond de ciel si OK
//|                -1 si probleme definition
//|                -2 si probleme allocation
//|                -3 si probleme remplissage histogramme
//|                -4 si probleme recherche maximum
//|                -5 si probleme calcul fraction du maximum
//|                -6 si probleme nombre de bin
//|                -7 si probleme nombre de fit
//| 
//--
//                // EA seul changement par rapport a la version C :
//                // pix est un T*, en template. -> C++ necessaire.
{
FLOATDOUBLE X,Y,EY;
double (*Fun) (double ,double *,double *);
int_4 npar;
double par[6], epar[6],minpar[6],maxpar[6],stepar[6],stochi2;
int_4 i,j,i1,j1,j2,n,ip,nbin,inc,ibin0,max,imax,lmax,b1,b2;
float *histo,*xbin,*ehisto,fwhm,rc;
double vp,pxv;

/* set reasonable parameters */
Fun = Gauss1DPolF ;
if ( deg<=1 ) deg=0;  else deg=2;
low = ( low < 0 ) ? 0 : low;
high = ( high > 65530 ) ? 65530 : high;
if ( low >= high ) { low=0; high = 65530;}
bin = ( bin <= 0 ) ? 1 : bin;
nbin = int(2.*taille /bin + 1);
if ( nbin < 20 ) nbin = 20;
if ( nbin >= 65530 ) nbin = 65530;
ibin0 = int(-nbin*bin/2.);
inc = int(nx*ny / ndata); inc = int(sqrt( (double) inc ));
inc = (inc<1) ? 1 : inc;
if(deb>1)
  { printf("Sigma_Ciel: pix[%d,%d] ndata=%d bin=%d taille=%d"
          ,nx,ny,ndata,bin,taille);
  printf(" frac=%f low=%g high=%g nbin=%d ibin0=%d inc=%d\n"
        ,frac,low,high,nbin,ibin0,inc); }
if ( nx < 5 || ny < 5 ) return (-1.);

/* allocate place for histogramme */
if ( (histo  = (float *) calloc(nbin, sizeof(float)) ) == NULL )
       return(-2.);
if ( (xbin   = (float *) calloc(nbin, sizeof(float)) ) == NULL )
       {free(histo); return(-2.);}
if ( (ehisto = (float *) calloc(nbin, sizeof(float)) ) == NULL )
       {free(histo); free(xbin); return(-2.);}

/* fill histogramme */
ndata=0;
for (j=1;j<ny-1;j=j+inc) { for (i=1;i<nx-1;i=i+inc) {
  ip = (int)(*(pix+j*nx+i));
  if ( ip < low || ip > high ) continue;
  vp = (double) ip;
  pxv = 0.; n=0;
  for(j1=j-1;j1<=j+1;j1++) { for(i1=i-1;i1<=i+1;i1++) {
    if ( i1 == i && j1 == j ) continue;
    ip = (int)(*(pix+j1*nx+i1));
    if ( ip < low || ip > high ) continue;
    pxv += (double) ip;
    n++;
  } }
  if ( n != 8 ) continue;
  vp = vp - pxv/n;
  ip = int((vp-ibin0)/bin);
  if( ip<0 || ip >= nbin ) continue;
  histo[ip]++;
  ndata++;
} }
for(i=0;i<nbin;i++)
  { xbin[i] = (float) ibin0 + ( (float) i + 0.5 ) * (float) bin;
  ehisto[i] = (histo[i]<1.) ? -1. : sqrt((double) histo[i]);}

if(deb>1) printf("     nombre de data corrects = %d\n",ndata);
if( ndata < 1 ) {rc= -3.; goto END;}

/* find maximum */
max = -1; imax = -1;
for(i=nbin-1;i>=0;i--) if(histo[i]>max) { max=int(histo[i]); imax=i;}
if(deb>1) printf("     maxi histo = %d (bin=%d)\n",max,imax);
if ( imax == -1 || imax == 0 || imax == nbin-1 ) {rc= -4.; goto END;}
if(deb>2) {
  for(i=imax-10;i<=imax+10;i++) if(i>=0 && i<nbin) printf(" %5f",xbin[i]);
  printf("\n");
  for(i=imax-10;i<=imax+10;i++) if(i>=0 && i<nbin) printf(" %5d",(int) histo[i]);
  printf("\n"); }

/* find limites a half maximum */
lmax = (int) ( (float) max * 0.5 );
for(b1=imax;b1>0;b1--)    if(histo[b1]<=lmax) break;
for(b2=imax;b2<nbin;b2++)  if(histo[b2]<=lmax) break;
fwhm = (xbin[b2]-xbin[b1]);
if(deb>1) printf("     fwhm = %f (%d,%d)\n",fwhm,b1,b2);

/* find limites a frac*max */
lmax = (int) ( (float) max * frac ) + 1;
if( lmax < 1 )
   {if(deb>1) printf("     lmax=%d\n",lmax); rc= -5.; goto END;}
for(b1=imax;b1>0;b1--)    if(histo[b1]<=lmax) break;
for(b2=imax;b2<nbin;b2++)  if(histo[b2]<=lmax) break;

n=0;
for(i=b1;i<=b2;i++) if( histo[i] > 0 ) n++;
if(deb>1) printf("     limites bin %d %d (v=%d) bin non nuls=%d\n"
                ,b1,b2,lmax,n);
/* pas assez de bin, on en ajoute artificiellement */
/* 8 est une borne large pour fiter gauss + parab (3+3) */
if( n < 8 ) {
  j1 = b1; j2 = b2;
  for(i=1;i<nbin;i++) {
    if ( b1-i >= 0 ) j1=b1-i;
    if ( b2+i < nbin ) j2=b2+i;
    if ( histo[j1] > 0 ) n++;
    if ( histo[j2] > 0 ) n++;
    if( n >= 6 ) break;
  }
  b1 = j1; b2 = j2;
  if(deb>1) printf("     Re-update des limites bin %d %d (v=%d) bin non nuls=%d\n"
                  ,b1,b2,lmax,n);
}
if( n < 4 ) {rc= -6.; goto END;}
ndata = n;

/* fit dans les limites */
npar = 3+deg+1; Set_FitFunDPol(deg);
if(ndata<8) {npar=4; Set_FitFunDPol(0);}
if(ndata<5) {npar=3; Set_FitFunDPol(-1);}
X.fx  = xbin;
Y.fx  = histo;
EY.fx = ehisto;
stochi2 = 0.0001;
n = nbin;
par[0] = max;
par[1] = histo[imax-1]*xbin[imax-1]+histo[imax]*xbin[imax]+histo[imax+1]*xbin[imax+1];
par[1] = par[1]/(histo[imax-1]+histo[imax]+histo[imax+1]);
par[2] = (fwhm<1.) ? 1. : fwhm/2.36;
par[3] = par[4] = par[5] = 0.;
stepar[0] = par[0]/20.; stepar[1] = 1.; stepar[2] = par[2]/5.;
stepar[3] = 1.; stepar[5] = 0.1;
stepar[4] =  0.;
for(i=0;i<6;i++) {minpar[i] = 1.; maxpar[i] = -1.;}
minpar[2] = 0.1; maxpar[2] = 65530.;
if(deb>1) printf("     choix de fit a %d parametres b1=%d,b2=%d,ndata=%d)\n"
                ,npar,b1,b2,ndata);

i = (deb>=3) ? deb-2 : 0 ;
vp = FitFun(Fun,FITFUN_FLOAT,X,Y,EY,&n,b1,b2,par
           ,epar,stepar,minpar,maxpar,npar,stochi2,50,i);

if( vp<=0. && vp != -4. ) {rc= -7.; goto END;}
rc = par[2];
if(deb>1) printf("     sigma=%f (h=%f centrage=%f) xi2=%f (%f pour %d)\n"
                ,rc,par[0],par[1],vp,vp/(n-npar),npar);

END:
free(histo); free(xbin); free(ehisto);
return(rc);
}

/* ................................................................................. */
/*                  Filtrage ( Convolution ) d' images                               */
/* ................................................................................. */

/* Nouvelle-Fonction */
template <class T>
//++
Image<T> * FilterImage(Image<T> const * pim, Image<T> * pom, ImageR4 *matx)
//
//	Filtrage ( Convolution ) d'une image `pim' par une matrice carree
//	representee par l'ImageR4 `matx'.
//	Le resultat se retrouve dans `pom'.
//	Si `pom' = NULL, elle est cree.
//	`pim' et `pom' doivent avoir la meme taille.
//--
{
 
float vpix, *mxe;
bool fgcr;
int i,j,k,l,m,n;
int ndim, nd2;

fgcr = false;

if ((pim == NULL) || (matx == NULL))
  { printf(" FilterImage_Erreur : pim ou matx = NULL ! \n");
  return(NULL);  }
if (matx->XSize() != matx->YSize())
  { printf("FilterImage_Erreur: Filtrage par matrice carre uniquement (%d,%d)\n",
            matx->XSize(), matx->YSize());  return(NULL); }
ndim = matx->XSize();
if ( (ndim < 3) || ((ndim%2) == 0) || (ndim > 25))
  { printf("FilterImage_Erreur: Pb dimension matrice (%d) (Impair, entre 3 et 25)\n", ndim);
  return(NULL); }
if (pom == NULL)
  {  pom = new Image<T>(pim->XSize(), pim->YSize(), false);   fgcr = true; }

if ( (pim->XSize() != pom->XSize()) || (pim->YSize() != pom->YSize()) )
  { printf("FilterImage_Erreur : Taille pim et pom differentes ! \n") ;
  if (fgcr) delete pom;   return(NULL) ;}

nd2 = ndim/2;

/*  Traitement de la partie centrale */
for(j=nd2; j<(pom->YSize()-nd2); j++)
  {
  for(i=nd2; i<(pom->XSize()-nd2); i++)
    {
    vpix = 0.;   mxe = matx->ImagePtr();
    for(l=-nd2; l<=nd2; l++)
      for(k=-nd2; k<=nd2; k++)
         { vpix += (float)((*pim)((i+k), (j+l))) * (*mxe);
         mxe++; }
    (*pom)(i, j) = (T)vpix;
    }
  }

/*  Traitement bord en Y  en haut */
for(j=0; j<nd2; j++)
  {
  for(i=nd2; i<(pom->XSize()-nd2); i++)
    {
    vpix = 0.;  mxe = matx->ImagePtr();
    for(l=-nd2; l<=nd2; l++)

      for(k=-nd2; k<=nd2; k++)
        {
        m = i+k;  n = j+l;
        if (n < 0)  n = 0;
        vpix += (float)((*pim)(m, n)) * (*mxe);
        mxe++;
        }
    (*pom)(i, j) = (T)vpix;
    }
  }

/* Traitement bord en Y en bas  */
for(j=(pom->YSize()-nd2); j<pom->YSize(); j++)
  {
  for(i=nd2; i<(pom->XSize()-nd2); i++)
    {
    vpix = 0.;  mxe = matx->ImagePtr();
    for(l=-nd2; l<=nd2; l++)
      for(k=-nd2; k<=nd2; k++)
        {
        m = i+k;  n = j+l;
        if (n >= pim->YSize())  n = (pim->YSize()-1);
        vpix += (float)((*pim)(m, n)) * (*mxe);
        mxe++;
        }
    (*pom)(i, j) = (T)vpix;
    }
  }

/*  traitement bords en X et coins */
for(j=0; j<pom->YSize(); j++)
  {
  for(i=0; i<nd2; i++)
    {
    vpix = 0.;  mxe = matx->ImagePtr();
    for(l=-nd2; l<=nd2; l++)
      for(k=-nd2; k<=nd2; k++)
        {
        m = i+k;  n = j+l;
        if (m < 0)  m = 0;
        if (n < 0)  n = 0;
        if (n >= pim->YSize())  n = (pim->YSize()-1);
        vpix += (float)((*pim)(m, n)) * (*mxe);
        mxe++;
        }
    (*pom)(i, j) = (T)vpix;
    }
  for(i=(pom->XSize()-nd2); i<pom->XSize(); i++)
    {
    vpix = 0.;  mxe = matx->ImagePtr();
    for(l=-nd2; l<=nd2; l++)
      for(k=-nd2; k<=nd2; k++)
        {
        m = i+k;  n = j+l;
        if (m >= pim->XSize())  m = (pim->XSize()-1);
        if (n < 0)  n = 0;
        if (n >= pim->YSize())  n = (pim->YSize()-1);
        vpix += (float)((*pim)(m, n)) * (*mxe);
        mxe++;
        }
    (*pom)(i, j) = (T)vpix;
    }
  }


return(pom);

}

/* Nouvelle-Fonction */
//++
ImageR4* NoiseFiltImage(ImageR4& img, ImageR4& filtre, DynCCD& dynccd)
//	Calcule une image de bruit, a partir de "img" et de "dynccd"
//	prenant en compte le filtre de convolution "filtre".
//	L'image de bruit renvoyee contient les correlations de bruit entre
//	pixels, du au filtrage (convolution) de l'image "img" par le 
//	filtre "filtre".
//--
{
  ImageR4* ImgBrt = NoiseImage(&img,&dynccd);

  int i,j;  
  for (i=0; i<ImgBrt->XSize(); i++)
    for (j=0; j<ImgBrt->YSize(); j++)  (*ImgBrt)(i,j) =(*ImgBrt)(i,j) *(*ImgBrt)(i,j);

  ImageR4 Filter2(filtre.XSize(), filtre.YSize()) ;
  for (i=0; i<filtre.XSize(); i++)
    for (j=0; j<filtre.YSize(); j++) (Filter2)(i,j) = filtre(i,j) * filtre(i,j);  
 
  ImageR4* out2;
  out2 = new ImageR4(ImgBrt->XSize(), ImgBrt->YSize());
  FilterImage(ImgBrt, out2, &Filter2);   
  for (i=0; i<out2->XSize(); i++)
    for (j=0; j<out2->YSize(); j++)  (*out2)(i,j) =sqrt((*out2)(i,j));
  delete ImgBrt;
  return(out2);
}


//////////////////////////////////////////////////////////////////////////////
/* Nouvelle-Fonction    CMV 14/10/97 */
template <class T>
//++
Image<T> * FilterStat(Image<T> const * pim, Image<T> * pom
                     ,int N,int M,float lowbad,float highbad)
//
//	Filtrage median d'une image `(*pim)'.
//	Le resultat se retrouve dans `(*pom)'.
//	Si pom = NULL l'image de sortie est cree automatiquement.
//	pim et pom doivent avoir la meme taille.
//	Le masque de calcul du filtre median est un rectangle `N'x`M'
//	La dynamique de l'image est `lowbad,highbad'.
//	Si lowbad>=highbad dynamique infinie.
//--
{
if(pim == NULL)
  {printf(" FilterStat_Erreur : pim = NULL ! \n");
  return(NULL);}
if( N<=0 || M<=0 )
  {printf("FilterStat_Erreur: Filtrage par matrice (%d,%d)\n",N,M);
  return(NULL);}
bool fgcr = false;
if(pom == NULL)
  {pom = new Image<T>(pim->XSize(), pim->YSize(), false);   fgcr = true;}
if( (pim->XSize() != pom->XSize()) || (pim->YSize() != pom->YSize()) )
  { printf("FilterImage_Erreur : Taille pim et pom differentes ! \n");
  if(fgcr) delete pom;   return(NULL);}

T low  = (T) lowbad;
T high = (T) highbad;
if(high>low) fgcr = true; else fgcr = false;

N /= 2; M /= 2;
double * tab = new double[(2*N+1)*(2*M+1)];

for(int j=0;j<pom->YSize();j++) {
  for(int i=0;i<pom->XSize();i++) {
    int n=0;
    for(int jj=j-M;jj<=j+M;jj++) {
      if( jj<0 || jj >= pom->YSize() ) continue;
      for(int ii=i-N;ii<=i+N;ii++) {
        if( ii<0 || ii >= pom->XSize() ) continue;
          if( fgcr && ( (*pim)(ii,jj)<=low || high<=(*pim)(ii,jj) ) ) continue;
          tab[n] = (double) (*pim)(ii,jj);
          n++;
      }
    }
    if(n==0) (*pom)(i,j) = 0;
    else if(n==1) (*pom)(i,j) = (T) tab[0];
    else {
      HeapSort(n,tab);
      if(n%2==1) (*pom)(i,j) = (T) tab[(n-1)/2];
      else (*pom)(i,j) = (T) (tab[n/2]+tab[n/2-1])/2.;
    }
  }
}

delete [] tab;
return(pom);
}

//////////////////////////////////////////////////////////////////////////////
/* Nouvelle-Fonction    CMV 7/8/98 */
template <class T>
//++
Image<T> * CSplineImage(Image<T> const & f1,double pxsz,double x0,double y0,int natural)
//
//	Fonction retournant un pointeur sur une image qui contient
//	le Spline 3 a 2 dimensions de l'images f1 (NULL est retourne si probleme).
//	`pxsz' est la taille du pixel de l'image retournee en considerant
//	que l'image `f1' en entree a une taille de pixel egale a 1.
//	`x0,y0' est le decalage eventuel a appliquer de l'image origine a
//	l'image retournee, ce decalage etant mesure en unites de pixels
//	de l'image origine (1).
//	`natural' a la signification decrite dans la classe `CSpline2'.
//	Attention, la taille de l'image retournee sera `nx/pxsz' et `ny/pxsz'.
//	Pour eviter les problemes d'extrapolation sur les bords de l'image
//	d'origine, il leurs est assigne la valeur du centre du pixel auquel ils
//	appartiennent.
//--
{
//// open for result
int nx = f1.XSize();
int ny = f1.YSize();
float dum;
dum = nx/pxsz;
int dx = (int) (dum+0.5);
dum = ny/pxsz;
int dy = (int) (dum+0.5);
cout<<"Taille image en entree: "<<nx<<" "<<ny<<endl
    <<"Taille image en sortie: "<<dx<<" "<<dy<<" (pix size="<<pxsz<<")"<<endl
    <<"Decalage en x="<<x0<<" en y="<<y0<<endl;
if(nx<=1 || ny<=1) return NULL;
if(dx<=0 || dy<=0) return NULL;

Image<T> *f2 = NULL;
TRY {
  f2 = new Image<T>(dx,dy);
} CATCHALL {
  cout<<"Impossible d'ouvrir le fichier resultat"<<endl;
  return NULL;
} ENDTRY

//// Spline 3 2D ... nx+2 pour tenir compte des effets de bords !
double *x = new double[nx+2];
double *y = new double[ny+2];
double *z = new double[(nx+2)*(ny+2)];
{
int k;
// Centre pixel (0,0) en (0.5,0.5) et on ajoute les coordonnees
// des bords pour extrapolation sur tout l'intervalle.
// Les coordonnees vont de [0,nx] et [0,ny] (valeurs extremes comprises).
x[0] = 0.; x[nx+1] = nx;
y[0] = 0.; y[ny+1] = ny;
for(k=1;k<=nx;k++) x[k] = k-0.5;
for(k=1;k<=ny;k++) y[k] = k-0.5;
// Remplissage des valeurs avec effet de bord [0] et [nx+1],[ny+1].
k = 0;
int ii,jj;
for(int j=0;j<ny+2;j++) {
  jj = j-1;
  if(jj<0) jj=0; else if(jj>=ny) jj=ny-1;
  for(int i=0;i<nx+2;i++) {
    ii = i-1;
    if(ii<0) ii=0; else if(ii>=nx) ii=nx-1;
    z[k] = (double) f1(ii,jj);
    k++;
  }
}
}

CSpline2 spl(nx+2,x,ny+2,y,z,natural,true);
spl.ComputeCSpline();

*f2 = 0.;
{
// Attention au decalage possible !
// Les valeurs vont de [0.5*pxsz,(dx-0.5)*pxsz] idem pour y...
// a corriger eventuellement du decalage x0,y0.
double xmin = -pxsz/10000., xmax = nx+pxsz/10000.;
double ymin = -pxsz/10000., ymax = ny+pxsz/10000.;
for(int j=0;j<dy;j++) {
  double yy = (j+0.5)*pxsz - y0;
  if( yy<=ymin || yy>=ymax ) continue;
  for(int i=0;i<dx;i++) {
    double xx = (i+0.5)*pxsz - x0;
    if( xx<=xmin || xx>=xmax ) continue;
    (*f2)(i,j) = (T) spl.CSplineInt(xx,yy);
}}}

//// nettoyage
delete [] x;
delete [] y;
delete [] z;

return f2;
}
//////////////////////////////////////////////////////////////////////////////



/* Nouvelle-Fonction */
int FondSigmaCiel(RzImage *pim, float low, float high,int pvsz
              ,float nbsig1,float nbsig2,float nbsig3
              ,float binsg,float frac,int modu,int deb)
{
int rc;

  switch (pim->PixelType())
    {
    case kuint_2:
      {
      ImageU2 imt((RzImage&)(*pim));
      rc = imt.FondSigmaCiel(low, high, pvsz, nbsig1, nbsig2, nbsig3, binsg, frac, modu, deb);
      pim->fond = imt.fond;
      pim->sigmaFond = imt.sigmaFond;
      break;
      }
    case kint_2:
      {
      ImageI2 imt((RzImage&)(*pim));
      rc = imt.FondSigmaCiel(low, high, pvsz, nbsig1, nbsig2, nbsig3, binsg, frac, modu, deb);
      pim->fond = imt.fond;
      pim->sigmaFond = imt.sigmaFond;
      break;
      }
    case kint_4:
      {
      ImageI4 imt((RzImage&)(*pim));
      rc = imt.FondSigmaCiel(low, high, pvsz, nbsig1, nbsig2, nbsig3, binsg, frac, modu, deb);
      pim->fond = imt.fond;
      pim->sigmaFond = imt.sigmaFond;
      break;
      }
    case kr_4:
      {
      ImageR4 imt((RzImage&)(*pim));
      rc = imt.FondSigmaCiel(low, high, pvsz, nbsig1, nbsig2, nbsig3, binsg, frac, modu, deb);
      pim->fond = imt.fond;
      pim->sigmaFond = imt.sigmaFond;
      break;
      }
    default:
      {
      cout<<"FindSigmaCiel: type d'image non-prevu "<<(int) pim->PixelType()<<endl;
      exit(-1);
      }
    }

return(rc);
}

/* Nouvelle-Fonction */
int CalFondImage(RzImage *pim, double low, double high, int action)

/*  ATTENTION: Adaptation a partir d'une version C              */ 
/*             CalculFond() de imgstat.c                        */
 
/*  Cette fonction calcule le fond et le sigma du ciel a la CMV */
/*  CMV / E.Aubourg / Reza                                      */
/*  action Bit 0  CheckDyn()                                    */
/*         Bit 1  ImgFondCielHisto()                            */
/*         Bit 2  ImgSigmaCiel                                  */
/*         Bit 3  FondSig2Ciel()                                */
/*         Bit 1+2+3 (=14)   Fond_Ciel() + Sigma_Ciel() Si OK   */
/*                            FondSig2Ciel()  sinon             */
/*  action = 0  ou  >= 15  CheckDyn + Bit 1+2+3 (=14)           */
/*  Rc=0 OK  Erreur Sinon                                       */
/*  Sortie : Mise a jour de la struct RzImage                   */

{
int n,ok,rc,fga;
float fnd, sg;

int nbsig, vitesse, nbin;
float minaff, maxaff, minpix, maxpix;
float fracmax,fracnul,fracsat;
char *meth;

#define DBG  0

if (action >= 15)   action = 15;
if (action <= 0)  action = 15;
fga = action;
fnd = sg = -1.;
ok = 1;

if (action & 1)  
  { pim->CheckDyn(low, high); 
   action--;  }

n = pim->XSize()*pim->YSize();

if (action & 6)
  {
  int npu = (n/4 < 5000) ? 5000 : n/4;
  meth = (char *)"Fond/Sig";
  switch (pim->PixelType())
    {
    case kuint_2:
      {
      ImageU2 imt((RzImage&)(*pim));
      if (action & 2) 
        fnd = ImgFondCielHisto(imt, 1.0f, 2, 0.5f, (float)low, (float)high, DBG);
      if (action & 4)
        {  sg = ImgSigmaCiel(imt, npu, 1, 200, 0.1f, 2, 
			     (float)low, (float)high, DBG);
          if (sg < -1.e-3 )  ok = 0; }
      break;
      }
    case kint_2:
      {
      ImageI2 imt((RzImage&)(*pim));
      if (action & 2) 
        fnd = ImgFondCielHisto(imt, 1.0f, 2, 0.5f, (float)low, (float)high, DBG);
      if (action & 4)
        {  sg = ImgSigmaCiel(imt, npu, 1, 200, 0.1f, 2,
			     (float)low, (float)high, DBG);
          if (sg < -1.e-3 )  ok = 0; }
      break;
      }
    case kint_4:
      {
      ImageI4 imt((RzImage&)(*pim));
      if (action & 2) 
        fnd = ImgFondCielHisto(imt, 1.0f, 2, 0.5f, (float)low, (float)high, DBG);
      if (action & 4)
        {  sg = ImgSigmaCiel(imt, npu, 1, 200, 0.1f, 2, 
			     (float)low, (float)high, DBG);
          if (sg < -1.e-3 )  ok = 0; }
      break;
      }
    case kr_4:
      {
      ImageR4 imt((RzImage&)(*pim));
      if (action & 2) 
        fnd = ImgFondCielHisto(imt, 1.0f, 2, 0.5f, (float)low, (float)high, DBG);
      if (action & 4)
        {  sg = ImgSigmaCiel(imt, npu, 1, 200, 0.1f, 2, 
			     (float)low,(float) high, DBG);
          if (sg < -1.e-3 )  ok = 0; }
      break;
      }

    case kuint_4:
    case kr_8:
      cerr << "CalFondImage(RzImage *, ...) Unsupported image type (kuint_4/kr_8) " 
           << (int) pim->PixelType() << "\n" ; 
      break;
    default:
      cerr << "CalFondImage(RzImage *, ...) Unknown image type !! " << (int) pim->PixelType() << "\n" ;
    break;
    }

  }

if ( (action == 8) || (!ok && (action & 8)) )
  {
  fracmax = 0.85;   nbsig = 1;     vitesse = 3;
  minpix = low;   maxpix = high;  nbin = 10000;
  rc = RzFondSig2Ciel(pim, vitesse, fracmax, nbsig, nbin,  
                      &minpix, &maxpix, &fnd, &sg,
                      &minaff, &maxaff, &fracnul, &fracsat);
  if (rc != 0)   ok = 0;
  meth = (char *)"FondSig2";
  }

if (ok > 0)     /*  Mise a jour */ 
  {
  if (fnd > -1.)  pim->fond = fnd; 
  if (sg > -1. )  pim->sigmaFond = sg;     
  printf("CalculFond_Info: Meth=%s (%d) Fond/Sig= %5g %4g Nul/Sat=%d %d\n",
         meth, fga, fnd, sg, pim->nbNul, pim->nbSat);
  return(0);
  }
else
  {
  printf("CalculFond_Erreur: Action=%d  NbNul/Sat=%d %d\n",fga,
          pim->nbNul, pim->nbSat);
  return(1);
  }

}



/* Nouvelle-Fonction */

int RzFondSig2Ciel(RzImage * pim, int vitesse, float FracMax, int NbSigmas, int NBin,
	           float *MinPix, float *MaxPix, float *Fond, float *SigFond,
                   float *MinAff, float *MaxAff, float *FracNul, float *FracSature)

/*   ATTENTION : Cette fonction a ete mis en C++, pour traiter des   */
/*               RzImage a partir de la fonction FondSig2Ciel()      */
/*               de imgstat.c                                        */
/*                                        R. Ansari  04/95           */

/*    Cette fonction calcule les valeurs suivantes pour une image    */
/*    Arguments d'entree :                                           */
/*     - pim : Pointeur RzImage                                      */
/*     - vitesse : 1 pixel / vitesse pris en compte                  */
/*     - FracMax : Parametre de calcul (voir ci-dessous)             */
/*     - NbSigmas: MinAff =  Fond + NbSigmas*SigFond                 */
/*       **Note : Si NbSigmas < 0 , on prend :                       */
/*                 MinAff =  Fond + NbSigmas*SigFond                 */
/*              et MaxAff =  Fond - NbSigmas*SigFond                 */
/*         Dans ce cas la gamme d'affichage fait ressortir le fond.  */
/*     - NBin : Nb de bins de l'histo                                */
/*     - La valeur mini des pixels et la valeur de saturation        */
/*       (ValNul et ValSat) doivent etre fournies ds Min/MaxPix      */
/*              sorties :                                            */
/*  - MinPix,MaxPix = Valeurs Min et Max des pixels  (<ValSat)       */
/*  - Fond =  Pic de l'histo des pixels                              */
/*  - SigFond = sigma du fond calcule par la relation                */
/*     NbPixel(Fond+SigFond) =  Pic histo * 0.606 (=Exp(-0.5))       */
/*  - MinAff = Bruit + NbSigmas * Sigma bruit                        */
/*  - MaxAff est calcule tel que la relation suivante soit verifie   */
/*     NbPixel(MinAff<Pixel<MaxAff) =  Nb de pixel entre Min-MaxAff  */
/*      = FracMax * NbPixel(MinAff<Pixel<ValSat                      */
/*  - FracNul = Nb pixel < ValNul / Nb total pixel                   */
/*  - FracSat = Nb pixel satures (>ValSat) / Nb total pixel          */

/*     return code : 0 OK , 1 Erreur calloc                          */
/*                   2 : Pic histos non trouve                       */
/*                   3,4 : Probleme de Min/MaxAff avec ValSat        */
/*                   5 : Gamme dynamique minimale non respectee      */
/*                         (MaxAff - MinAff > MINDYNAMIC             */

{

#define MINNBIN 100
#define MAXNBIN 32765

int nbin;

float minhis = 0., maxhis = 65530. , binwidth;
int xdim, ydim;
int_4 under,over,tothis;
int_4 *phis;
float ValNul,ValSat;
int_4 *q;

float val;
int i;
float min, max;
int_4 nbpix,nbpix2, imin, imax;
int bin,binmax;
int sigl,sigr;
float fracfond,fracsignal;

float MINDYNAMIC = 5;          /* Gamme dynamique minimale (MaxAff-MinAff)*/

#define FRACPIC		0.6065     /*  Valeur de exp(-0.5)   */

xdim = pim->XSize();
ydim = pim->YSize();

ValNul = *MinPix;
ValSat = *MaxPix;
if ((ValSat < ValNul+0.01))  ValSat = ValNul+0.01;  

nbin = NBin;
if (nbin < MINNBIN)  nbin = MINNBIN;
if (nbin > MAXNBIN)  nbin = MAXNBIN;

if (vitesse < 1)  vitesse = 1;
if (vitesse > xdim*ydim/100)  vitesse = xdim*ydim/100;


/*  Je modifie min et max histo et nbin en fonction de valnul et valsat */
/*  je garde binwidth constant     */

binwidth = (ValSat-ValNul) / nbin;
minhis = ValNul;
maxhis = minhis + nbin*binwidth;


(*MinAff) = minhis;
(*MaxAff) = maxhis;
*MinPix = 0.;
*MaxPix = -1.;
*Fond = -1.;
*SigFond = -1.;
*FracNul = 0.0;
*FracSature = 0.0;

/* allocation dynamique de memoire */
phis = new int_4[nbin];

/* mise a zero de l'histogramme */
for ( q=phis ; q<phis+nbin ; *q++ = 0 )  ;
under = over = tothis = 0;



nbpix2 = nbpix = 0;
min = ValSat+1.;
max = ValNul-1;

/*   Boucle sur les pixels   */
/*     Pour aller plus vite, on prend un pixel sur 7   */
for ( i = 0 ; i < xdim*ydim ; i += vitesse )
  {
  val = pim->IValue(i);

   /*   Calcul PixMax/Min et NbPix sature  */
  if (val > ValSat)  nbpix++;
  else {
    if(val < ValNul) nbpix2++;
    else {
	    if (val > max) max = val;
	    else if (val < min) min = val;
	}
  }

   /* remplissage de l'histogramme */
  bin = (int)((val - minhis) / binwidth) ;
  if ( bin < 0 )  under++;
  else 
    {
    if (bin < nbin) { (*(phis+bin))++; tothis++; }
    else   over++; 
    }
  }     /*  Fin de for : Boucle sur les pixels   */

(*MinPix) = min; 
(*MaxPix) = max;
/*       Ne pas oublier le facteur vitesse pour le calcul FracSatue  */
(*FracNul) = (float)(nbpix2 * vitesse) / (float)(xdim*ydim); 
(*FracSature) = (float)(nbpix * vitesse) / (float)(xdim*ydim); 

/*  On verifie qu'il y a eu une entree ds l'histo */
if (tothis <= 0)
  { delete[] phis;  return(2); }
  
/* recherche du maximum */
binmax = -1 ; imax = -1 ; q = phis;
for ( i = 0 ; i < nbin ;  i++, q++ )
  if ( *(q) > imax ) {imax = *(q) ; binmax = i ;}

if ( (binmax < 0) || (binmax >= nbin) )
  { delete[] phis;  return(2); }

*Fond = minhis + binmax * binwidth;   

/*   On va chercher le sigma du fond  */
imin = (int_4) ( (float)(imax) * FRACPIC ) ;

/*
printf("StatsDebug : min/maxhis  %g  %g  (nbin= %d) w=%g\n",
        minhis,maxhis,nbin,binwidth);
printf(".binmax= %d  imax(=pic) = %d  , imin= %d  \n ",binmax,imax,imin);   
*/

nbpix = 0;

sigr = 0;
q = phis+binmax;
while ( ((*q) > imin) && (sigr < (nbin-binmax)) ) 
  { sigr++; nbpix += (*q) ; q++; }

sigl = 0;
q = phis+binmax;
while ( ((*q) > imin) && (sigl < binmax) ) 
  { sigl++; nbpix += (*q) ; q--; }
nbpix -= (*(phis+binmax));

(*SigFond) = (float)(sigl+sigr)*binwidth/2.; 
fracfond = (float)(nbpix) / (float)(tothis);

(*MinAff) = minhis+(binmax-sigl)*binwidth;
(*MaxAff) = minhis+(binmax+sigr)*binwidth;


/*   Check sur NbSigmas  :  */
if (NbSigmas == 0) NbSigmas = 1;
if (NbSigmas < 0)    /*   On considere le cas ou NbSigmas est negatif   */
  {
  if ( (bin = binmax+NbSigmas*sigr) < 0 )  bin = 0;
  (*MinAff) = minhis + (float)bin*binwidth;
  if ( (bin = binmax-NbSigmas*sigr) >= nbin ) bin = nbin-1;
  (*MaxAff) = minhis + (float)bin*binwidth;
  goto sortie;
  }

if ( (bin = binmax+NbSigmas*sigl) >= nbin ) 
  { delete[] phis;  return(4); }

(*MinAff) = minhis + bin*binwidth;

/*   Calcul du nb de pixel entre MinAff et ValSat  */
nbpix = 0;
for( q = phis+bin ; q < phis+nbin ; q++)  nbpix += (*q);

fracsignal = (float)(nbpix) / (float)(tothis+over+under);
nbpix = (int_4)((float)nbpix * FracMax);
nbpix2 = 0; q = phis+bin; i = bin;
while( (nbpix2 < nbpix) && (i < nbin) )  
  { nbpix2 += (*q); i++; q++; }

(*MaxAff) = minhis + (float)i*binwidth;


sortie:        /*   On sort en liberant l'espace   */

/*
printf("  Stats..Debug : Fond,sigl,sigr,SigFond = %g,%d,%d,%g \n",
(*Fond),sigl,sigr,(*SigFond));
 printf("    Min/MaxPix= %g  %g , Min/MaxAff= %g  %g \n",(*MinPix),
(*MaxPix),(*MinAff),(*MaxAff));
printf("    FracFond,FracSignal %g  %g  FracNul,FracSat  %g  %g\n",
       fracfond,fracsignal,(*FracNul),(*FracSature));   
*/

delete[] phis;
if ( (*MaxAff) < ((*MinAff)+MINDYNAMIC) )  return(5);
return(0);

}


// ********** INSTANCES

#if defined(__xlC) || defined(__aCC__)
void instancetempimageop(int n)
{
/* Cette fonction sert uniquement a forcer le compilo a instancier les
  classes/fonctions template   */

int m = n-50;
ImageU2  iu2(n,n), *tiu2[5], xiu2(n,n);
ImageI2  ii2(n,n), *tii2[5], xii2(n,n);
ImageI4  ii4(n,n), *tii4[5], xii4(n,n);
ImageR4  ir4(n,n), *tir4[5], xir4(n,n);

ImageR4  matx(5,5);

RzCmvSigmaCiel(iu2.ImagePtr(), n, n, 100, 50, 50, (float)0.5, 2, (float)0., (float)65000., 0);
RzCmvSigmaCiel(ii2.ImagePtr(), n, n, 100, 50, 50, (float)0.5, 2, (float)0., (float)65000., 0);
RzCmvSigmaCiel(ii4.ImagePtr(), n, n, 100, 50, 50, (float)0.5, 2, (float)0., (float)65000., 0);
RzCmvSigmaCiel(ir4.ImagePtr(), n, n, 100, 50, 50, (float)0.5, 2, (float)0., (float)65000., 0);

ImgFondCielHisto(iu2);
ImgFondCielHisto(ii2);
ImgFondCielHisto(ii4);
ImgFondCielHisto(ir4);

ImgSigmaCiel(iu2, 100, 1, 100);
ImgSigmaCiel(ii2, 100, 1, 100);
ImgSigmaCiel(ii4, 100, 1, 100);
ImgSigmaCiel(ir4, 100, 1, 100);

CleanImages(5, tiu2); 
CleanImages(5, tii2); 
CleanImages(5, tii4); 
CleanImages(5, tir4); 


FilterImage(&iu2, &xiu2, &matx);
FilterImage(&ii2, &xii2, &matx);
FilterImage(&ii4, &xii4, &matx);
FilterImage(&ir4, &xir4, &matx);

FilterStat(&iu2, &xiu2,(int)3,(int)3,(float)1.,(float)-1.);
FilterStat(&ii2, &xii2,(int)3,(int)3,(float)1.,(float)-1.);
FilterStat(&ii4, &xii4,(int)3,(int)3,(float)1.,(float)-1.);
FilterStat(&ir4, &xir4,(int)3,(int)3,(float)1.,(float)-1.);

CSplineImage(iu2,(double)1,(double)0,(double)0,(int)0);
CSplineImage(ii2,(double)1,(double)0,(double)0,(int)0);
CSplineImage(ii4,(double)1,(double)0,(double)0,(int)0);
CSplineImage(ir4,(double)1,(double)0,(double)0,(int)0);

return;
}

#endif


#ifdef __CXX_PRAGMA_TEMPLATES__

#pragma define_template RzCmvSigmaCiel<uint_2>
#pragma define_template RzCmvSigmaCiel<int_2>
#pragma define_template RzCmvSigmaCiel<int_4>
#pragma define_template RzCmvSigmaCiel<r_4>

#pragma define_template ImgFondCielHisto<uint_2>
#pragma define_template ImgFondCielHisto<int_2>
#pragma define_template ImgFondCielHisto<int_4>
#pragma define_template ImgFondCielHisto<r_4>

#pragma define_template ImgSigmaCiel<uint_2>
#pragma define_template ImgSigmaCiel<int_2>
#pragma define_template ImgSigmaCiel<int_4>
#pragma define_template ImgSigmaCiel<r_4>

#pragma define_template CleanImages<uint_2>
#pragma define_template CleanImages<int_2>
#pragma define_template CleanImages<int_4>
#pragma define_template CleanImages<r_4>


#pragma define_template FilterImage<uint_2>
#pragma define_template FilterImage<int_2>
#pragma define_template FilterImage<int_4>
#pragma define_template FilterImage<r_4>

#pragma define_template FilterStat<uint_2>
#pragma define_template FilterStat<int_2>
#pragma define_template FilterStat<int_4>
#pragma define_template FilterStat<r_4>

#pragma define_template CSplineImage<uint_2>
#pragma define_template CSplineImage<int_2>
#pragma define_template CSplineImage<int_4>
#pragma define_template CSplineImage<r_4>

#endif


#if ( defined(ANSI_TEMPLATES) && !defined(__aCC__) )
template float RzCmvSigmaCiel<uint_2 const>(uint_2 const*, int, int, int, int, int, float, int, float, float, int);
template float RzCmvSigmaCiel< int_2 const>( int_2 const*, int, int, int, int, int, float, int, float, float, int);
template float RzCmvSigmaCiel< int_4 const>( int_4 const*, int, int, int, int, int, float, int, float, float, int);
template float RzCmvSigmaCiel<   r_4 const>(   r_4 const*, int, int, int, int, int, float, int, float, float, int);
template float ImgFondCielHisto<uint_2>(Image<uint_2> const&, float, int, float, float , float , int);
template float ImgFondCielHisto< int_2>(Image< int_2> const&, float, int, float, float , float , int);
template float ImgFondCielHisto< int_4>(Image< int_4> const&, float, int, float, float , float , int);
template float ImgFondCielHisto<   r_4>(Image<   r_4> const&, float, int, float, float , float , int);
template float ImgSigmaCiel<uint_2>(Image<uint_2> const& img, int, int, int, float, int, float , float, int);
template float ImgSigmaCiel< int_2>(Image< int_2> const& img, int, int, int, float, int, float , float, int);
template float ImgSigmaCiel< int_4>(Image< int_4> const& img, int, int, int, float, int, float , float, int);
template float ImgSigmaCiel<   r_4>(Image<   r_4> const& img, int, int, int, float, int, float , float, int);
template void CleanImages<uint_2>(int, Image<uint_2>**, double, double, double);
template void CleanImages< int_2>(int, Image< int_2>**, double, double, double);
template void CleanImages< int_4>(int, Image< int_4>**, double, double, double);
template void CleanImages<   r_4>(int, Image<   r_4>**, double, double, double);


template Image<uint_2>* FilterImage<uint_2>(Image<uint_2> const *, Image<uint_2>*, ImageR4 *);
template Image< int_2>* FilterImage< int_2>(Image< int_2> const *, Image< int_2>*, ImageR4 *);
template Image< int_4>* FilterImage< int_4>(Image< int_4> const *, Image< int_4>*, ImageR4 *);
template Image<   r_4>* FilterImage<   r_4>(Image<   r_4> const *, Image<   r_4>*, ImageR4 *);

template Image<uint_2>* FilterStat<uint_2>(Image<uint_2> const *, Image<uint_2>*,int,int,float,float);
template Image< int_2>* FilterStat< int_2>(Image< int_2> const *, Image< int_2>*,int,int,float,float);
template Image< int_4>* FilterStat< int_4>(Image< int_4> const *, Image< int_4>*,int,int,float,float);
template Image<   r_4>* FilterStat<   r_4>(Image<   r_4> const *, Image<   r_4>*,int,int,float,float);

template Image<uint_2>* CSplineImage<uint_2>(Image<uint_2> const &,double,double,double,int);
template Image< int_2>* CSplineImage< int_2>(Image< int_2> const &,double,double,double,int);
template Image< int_4>* CSplineImage< int_4>(Image< int_4> const &,double,double,double,int);
template Image<   r_4>* CSplineImage<   r_4>(Image<   r_4> const &,double,double,double,int);

#endif


#ifdef GNU_TEMPLATES
template float RzCmvSigmaCiel(uint_2 const*, int, int, int, int, int, float, int, float, float, int);
template float RzCmvSigmaCiel( int_2 const*, int, int, int, int, int, float, int, float, float, int);
template float RzCmvSigmaCiel( int_4 const*, int, int, int, int, int, float, int, float, float, int);
template float RzCmvSigmaCiel(   r_4 const*, int, int, int, int, int, float, int, float, float, int);
template float ImgFondCielHisto(Image<uint_2> const&, float, int, float, float , float , int);
template float ImgFondCielHisto(Image< int_2> const&, float, int, float, float , float , int);
template float ImgFondCielHisto(Image< int_4> const&, float, int, float, float , float , int);
template float ImgFondCielHisto(Image<   r_4> const&, float, int, float, float , float , int);
template float ImgSigmaCiel(Image<uint_2> const& img, int, int, int, float, int, float , float, int);
template float ImgSigmaCiel(Image< int_2> const& img, int, int, int, float, int, float , float, int);
template float ImgSigmaCiel(Image< int_4> const& img, int, int, int, float, int, float , float, int);
template float ImgSigmaCiel(Image<   r_4> const& img, int, int, int, float, int, float , float, int);
template void CleanImages(int, Image<uint_2>**, double, double, double);
template void CleanImages(int, Image< int_2>**, double, double, double);
template void CleanImages(int, Image< int_4>**, double, double, double);
template void CleanImages(int, Image<   r_4>**, double, double, double);


template Image<uint_2>* FilterImage(Image<uint_2> const *, Image<uint_2>*, ImageR4 *);
template Image< int_2>* FilterImage(Image< int_2> const *, Image< int_2>*, ImageR4 *);
template Image< int_4>* FilterImage(Image< int_4> const *, Image< int_4>*, ImageR4 *);
template Image<   r_4>* FilterImage(Image<   r_4> const *, Image<   r_4>*, ImageR4 *);

template Image<uint_2>* FilterStat(Image<uint_2> const *, Image<uint_2>*,int,int,float,float);
template Image< int_2>* FilterStat(Image< int_2> const *, Image< int_2>*,int,int,float,float);
template Image< int_4>* FilterStat(Image< int_4> const *, Image< int_4>*,int,int,float,float);
template Image<   r_4>* FilterStat(Image<   r_4> const *, Image<   r_4>*,int,int,float,float);

template Image<uint_2>* CSplineImage(Image<uint_2> const &,double,double,double,int);
template Image< int_2>* CSplineImage(Image< int_2> const &,double,double,double,int);
template Image< int_4>* CSplineImage(Image< int_4> const &,double,double,double,int);
template Image<   r_4>* CSplineImage(Image<   r_4> const &,double,double,double,int);

#endif