//  Classe image generique      
//                  E.Aubourg , E. Lesquoy  01-03/95 
//                  Modifs R. Ansari   04/95

// LAL (Orsay) / IN2P3-CNRS  DAPNIA/SPP (Saclay) / CEA

//  $Id: rzimage.cc,v 1.3 1999-04-27 15:03:02 ansari Exp $

#include "machdefs.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "rzimage.h"
#include "dvlist.h"
#include "generalfit.h"


/* .........  Methodes de la classe  RzVect    .........   */
//++
// Class 	RzVect<T>
// Lib 		Images++
// include	rzvect.h
//
//	Cette classe permet de construire des vecteurs de types T 
//	"(short, int, float, double, ...)". Cette classe est utilisee
//	par la classe *RzImage*.
//--

//++
// Titre	Constructeurs
//--
//++
// RzVect<T>(int n)
//	Constructeur avec allocation d'espace pour "n" elements
//	"RzVect<uint_1>, RzVect<uint_2>, RzVect<int_2>, RzVect<int_4>"
//	" RzVect<r_4>, RzVect<r_8>"
// RzVect<T>(int n, T* data)
//	Constructeur avec specification du tableau de stockage des elements
//	*Attention* "data" doit etre alloue a la bonne taille ("n")
// ~RzVect<T>()
//	Destructeur - Libere l'espace de stockage des elements si l'objet a
//	ete construit par le premier constructeur
//	
//--

/* --Methode-- */
template <class T>
RzVect<T>::RzVect(int n)
: p(new T[n]), nElt(n), nImgRef(1), shareable(1), ownTable(1)
{END_CONSTRUCTOR}

/* --Methode-- */
template <class T>
RzVect<T>::RzVect(int n, T* data)
: p(data), nElt(n), nImgRef(0), shareable(0), ownTable(0)
{END_CONSTRUCTOR}


/* --Methode-- */
template <class T>
RzVect<T>::~RzVect()
{
  if (ownTable) delete[] p;
  #ifdef DEBUG
  p = (T*) -1;
  #endif
}

//++
// Titre	Acces aux elements
//--
//++
// T& operator[](int i)
//	Acces a l'element "i" (lecture/ecriture)
// int  IElem(int k)
//      Valeur converti en entier "(int)" de l'element "k" (lecture)
// float  FElem(int k)
//      Valeur converti en reel "(float)" de l'element "k" (lecture)
// double  DElem(int k)
//      Valeur converti en reel "(double)" de l'element "k" (lecture)
//--

/* ..............   Classe  RzImage    ...............   */

//++
// Class 	RzImage
// Lib 		Images++
// include	rzimage.h
//
//	Cette classe permet de gerer les donnees correspondantes a une image.
//	Une image est un tableau rectangulaire de valeurs de pixels.
//	*RzImage* offre la possibilite de creer des images de valeurs de
//	pixels entier ou flottant "(unsigned char uint_1, short int_2," 
//	"unsigned short uint_2, int int_4, float r_4, double r_8)"
//	les objets de la classe RzImage sont faiblement types. La classe
//      *Image<T>* qui derive de RzImage permet la creation d'objets
//	types. Chaque image peut contenir un objet *DVList* permettant 
//	a l'utilisateur d'associer un ensemble de valeurs a l'image. 
//--
//++
// Links	Parents
// PPersist
//--
//++
// Links	Descendants
// Image<T>
// FitsImage<T>
//--
//++
// Links	Voir aussi
// DVList
//--

//++
// Titre	Constructeurs
//--
//++
// RzImage()
//	Creation d'une image vide, non typee.
// RzImage(PBaseDataTypes dType, int sizx, int sizy, -
//         int imgId = 0, char const * imgName = 0)
//	Creation d'une image de taille "sizx * sizy" pixels , de type "dType".
//      (dType = kuint_1, kuint_2, kint_2, kint_4, kr_4, kr_8).
//      "imgId" et "imgName" permettent d'indiquer une identification optionnelle.
// RzImage(const RzImage&, int sharePixels=0)
//	Constructeur par copie. Permet de partager le tableau des pixels (RzVect)
//	si "sharePixels" non nul (<> 0)
// RzImage(char *flnm)
//	Les RzImage etant des objets persistants (*PPersist*), ils peuvent etre
//	instancies a partir d'un fichier (nom = "flnm") ecrit par la methode Write()
// ~RzImage()
//	Destructeur. 
//  operator= (RzImage const &)
//	Copie les valeurs de pixels d'une image a l'autre.
//	Si la 1ere image a son propre tableau de pixel, il faut que 
//	les deux images aient le meme type de pixels et la meme taille.
//	Sinon, une exception ("typeMismatchErr" ou "sizeMismatchErr") est 
//	generee.  
//--

/* --Methode-- */
RzImage::RzImage(PBaseDataTypes dType)
: dataType(dType)
{
  voidP = NULL;
  isFits = 0;
  siz_x = siz_y = 0;
  vect.u2 = 0;
  mInfo = NULL;
  SetNameId(0, "");
  SetOrg(0, 0);
  SetPxSize();
  SetAtt();
  END_CONSTRUCTOR
}

/* --Methode-- */
RzImage::RzImage(PBaseDataTypes dType, int sizx, int sizy, int imgId, char const* imgName)
: dataType(dType)
{
  isFits = 0;
  siz_x = siz_y = 0;
  vect.u2 = 0;
  voidP = NULL;
  mInfo = NULL;
  Allocate(dType, sizx, sizy);
  SetNameId(imgId, imgName);
  SetOrg(0, 0);
  SetPxSize();
  SetAtt();
  END_CONSTRUCTOR
}

/* --Methode-- */
RzImage::RzImage(const RzImage& src, int sharePixels)
: dataType(kpbdt_unknown)
{
  isFits = 0;
  siz_x = siz_y = 0;
  vect.u2 = 0;
  voidP = NULL;
  mInfo = NULL;
  if (sharePixels) {
    ImgVectP vs = src.vect;
    Allocate(src.PixelType(), src.XSize(), src.YSize(), &vs); 
  }
  else { 
    Allocate(src.PixelType(), src.XSize(), src.YSize(), NULL);
    memcpy(voidP, src.voidP, siz_x*siz_y*DataSize(dataType)); 
  }

  SetNameId(src.id, src.name);
  SetOrg(src.XOrg(), src.YOrg());
  SetPxSize(src.XPxSize(), src.YPxSize());    
  SetAtt(src.nbNul, src.nbSat,
       src.minPix, src.maxPix, src.moyPix, src.sigPix,
       src.fond, src.sigmaFond);
  END_CONSTRUCTOR
}

/* --Methode-- */
RzImage::RzImage(char *flnm)
: dataType(kpbdt_unknown)
{
  isFits = 0;
  vect.u2 = 0;
  voidP = NULL;
  mInfo = NULL;
  PInPersist s(flnm);
  Read(s);
  END_CONSTRUCTOR
}

/* --Methode-- */
RzImage::~RzImage()
{
if (mInfo) delete mInfo;
DeAllocPixels();
}

/* --Methode-- */
RzImage& RzImage::operator = (const RzImage& src)
{

if (PixelType() >= 0)
  {
  if (src.PixelType() != PixelType())
    THROW(typeMismatchErr);
  if (siz_x != src.siz_x || siz_y != src.siz_y) 
    THROW(sizeMismatchErr);
  } 
else  Allocate(src.PixelType(), src.siz_x, src.siz_y);
  
memcpy(voidP, src.voidP, siz_x*siz_y*DataSize(dataType));
SetNameId(src.id, src.name);
SetOrg(src.XOrg(), src.YOrg());
SetPxSize(src.XPxSize(), src.YPxSize());    
SetAtt(src.nbNul, src.nbSat,
       src.minPix, src.maxPix, src.moyPix, src.sigPix,
       src.fond, src.sigmaFond);
 
return(*this); 
}

/* --Methode-- */
void RzImage::Allocate(PBaseDataTypes dType, int sizx, int sizy, ImgVectP* pvpsh)
// ATTENTION : Si pvpsh non nul, on partage les pixels 
{
  DeAllocPixels(); // EA 090998 -- memory leak...
  siz_x = sizx;
  siz_y = sizy;
  if ((int)dataType<=0)  dataType = dType; // Sinon ERREUR. Exception?
  switch (dataType) {
  case kuint_1:
    vect.u1 = NULL;
    if (pvpsh)
      if (pvpsh->u1)  
        { ASSERT(pvpsh->u1->Shareable());
        vect.u1 = pvpsh->u1;  vect.u1->NImgRef()++; }
    if (!vect.u1)  vect.u1 = new RzVect<uint_1> (siz_x * siz_y);
    voidP   = vect.u1->p;
    IVal_p = &RzImage::IVal_u1;
    FVal_p = &RzImage::FVal_u1;
    DVal_p = &RzImage::DVal_u1;
    break;
  case kuint_2:
    vect.u2 = NULL;
    if (pvpsh)
      if (pvpsh->u2)  
        { ASSERT(pvpsh->u2->Shareable());
        vect.u2 = pvpsh->u2;  vect.u2->NImgRef()++; }
    if (!vect.u2)  vect.u2 = new RzVect<uint_2> (siz_x * siz_y);
    voidP   = vect.u2->p;
    IVal_p = &RzImage::IVal_u2;
    FVal_p = &RzImage::FVal_u2;
    DVal_p = &RzImage::DVal_u2;
    break;
  case kint_2:
    vect.i2 = NULL;
    if (pvpsh)
      if (pvpsh->i2)  
        { ASSERT(pvpsh->i2->Shareable());
        vect.i2 = pvpsh->i2;  vect.i2->NImgRef()++; }
    if (!vect.i2)  vect.i2 = new RzVect<int_2> (siz_x * siz_y);
    voidP   = vect.i2->p;
    IVal_p = &RzImage::IVal_i2;
    FVal_p = &RzImage::FVal_i2;
    DVal_p = &RzImage::DVal_i2;
    break;
  case kint_4:
    vect.i4 = NULL;
    if (pvpsh)
      if (pvpsh->i4)  
        { ASSERT(pvpsh->i4->Shareable());
        vect.i4 = pvpsh->i4;  vect.i4->NImgRef()++; }
    if (!vect.i4)  vect.i4 = new RzVect<int_4> (siz_x * siz_y);
    voidP   = vect.i4->p;
    IVal_p = &RzImage::IVal_i4;
    FVal_p = &RzImage::FVal_i4;
    DVal_p = &RzImage::DVal_i4;
    break;
  case kr_4:
    vect.r4 = NULL;
    if (pvpsh)
      if (pvpsh->r4)  
        { ASSERT(pvpsh->r4->Shareable());
        vect.r4 = pvpsh->r4;  vect.r4->NImgRef()++; }
    if (!vect.r4)  vect.r4 = new RzVect<r_4> (siz_x * siz_y);
    voidP   = vect.r4->p;
    IVal_p = &RzImage::IVal_r4;
    FVal_p = &RzImage::FVal_r4;
    DVal_p = &RzImage::DVal_r4;
    break;
  case kr_8:
    vect.r8 = NULL;
    if (pvpsh)
      if (pvpsh->r8)  
        { ASSERT(pvpsh->r8->Shareable());
        vect.r8 = pvpsh->r8;  vect.r8->NImgRef()++; }
    if (!vect.r8)  vect.r8 = new RzVect<r_8> (siz_x * siz_y);
    voidP   = vect.r8->p;
    IVal_p = &RzImage::IVal_r8;
    FVal_p = &RzImage::FVal_r8;
    DVal_p = &RzImage::DVal_r8;
    break;
  default:
    printf("RzImage::Allocate Error - Unsupported DataType= %d \n", (int)dataType);
    break;    
  } 
}

/* --Methode-- */
void RzImage::SetPixels(PBaseDataTypes dType, int sizx, int sizy, void* data)
{
  if (voidP != data)  DeAllocPixels();
  siz_x = sizx;
  siz_y = sizy;
  //  if (dataType<0) (int&) dataType = (int)dType; // Sinon ERREUR. Exception?
  if (dataType<0)  dataType = dType; // Sinon ERREUR. Exception?
  switch (dataType) {
  case kuint_1:
    vect.u1 = new RzVect<uint_1> (siz_x * siz_y, (uint_1*)data);
    voidP   = vect.u1->p;
    IVal_p = &RzImage::IVal_u1;
    FVal_p = &RzImage::FVal_u1;
    DVal_p = &RzImage::DVal_u1;
    break;
  case kuint_2:
    vect.u2 = new RzVect<uint_2> (siz_x * siz_y, (uint_2*)data);
    voidP   = vect.u2->p;
    IVal_p = &RzImage::IVal_u2;
    FVal_p = &RzImage::FVal_u2;
    DVal_p = &RzImage::DVal_u2;
    break;
  case kint_2:
    vect.i2 = new RzVect<int_2> (siz_x * siz_y, (int_2*)data);
    voidP   = vect.i2->p;
    IVal_p = &RzImage::IVal_i2;
    FVal_p = &RzImage::FVal_i2;
    DVal_p = &RzImage::DVal_i2;
    break;
  case kint_4:
    vect.i4 = new RzVect<int_4> (siz_x * siz_y, (int_4*)data);
    voidP   = vect.i4->p;
    IVal_p = &RzImage::IVal_i4;
    FVal_p = &RzImage::FVal_i4;
    DVal_p = &RzImage::DVal_i4;
    break;
  case kr_4:
    vect.r4 = new RzVect<r_4> (siz_x * siz_y, (r_4*)data);
    voidP   = vect.r4->p;
    IVal_p = &RzImage::IVal_r4;
    FVal_p = &RzImage::FVal_r4;
    DVal_p = &RzImage::DVal_r4;
    break;
  case kr_8:
    vect.r8 = new RzVect<r_8> (siz_x * siz_y, (r_8*)data);
    voidP   = vect.r8->p;
    IVal_p = &RzImage::IVal_r8;
    FVal_p = &RzImage::FVal_r8;
    DVal_p = &RzImage::DVal_r8;
    break;
  default:
    printf("RzImage::SetPixels Error - Unsupported DataType= %d \n", (int)dataType);
    break;    
  } 
}

/* --Methode-- */
void RzImage::DeAllocPixels()
{
if (!voidP)  return;
  switch (dataType) {
  case kuint_1:
    if (! (--vect.u1->NImgRef()))
      delete vect.u1;
    vect.u1 = NULL;
    break;
  case kuint_2:
    if (! (--vect.u2->NImgRef()))
      delete vect.u2;
    vect.u2 = NULL;
    break;
  case kint_2:
    if (! (--vect.i2->NImgRef()))
      delete vect.i2;
    vect.i2 = NULL;
    break;
  case kint_4:
    if (! (--vect.i4->NImgRef()))
      delete vect.i4;
    vect.i4 = NULL;
    break;
  case kr_4:
    if (! (--vect.r4->NImgRef()))
      delete vect.r4;
    vect.r4 = NULL;
    break;
  case kr_8:
    if (! (--vect.r8->NImgRef()))
      delete vect.r8;
    vect.r8 = NULL;
    break;
  default:
    printf("RzImage::DeAllocPixels Error - Unsupported DataType= %d \n", (int)dataType);
    break;    
  }

voidP = NULL;
siz_x = siz_y = 0;
vect.u2 = NULL;
}

//++
// Titre	Informations globales
//--
//++
// int_4  XSize(), YSize()
//	Taille de l'image en X, Y
// int_4  XOrg(),  YOrg()
//	Origines (X,Y) de l'image. Cette origine permet en particulier 
//	de conserver l'information de decalage, quand l'image est 
//	extraite d'une image plus grande. Les valeurs de l'origine
//	peuvent etre modifiees par la methode "SetOrg()"
// void SetOrg(int_4 orgx, int_4 orgy)
//	Modifie l'origine de l'image
// float  XPxSize()  YPxSize()
//	Taille des pixels (X,Y) de l'image. Peuvent etre modifiees
//	par la methode "SetPxSize()"
// void SetPxSize()
// PBaseDataTypes   PixelType()
//	Type des pixels de l'image "(kint_2, kint_4, kr_4, ...)"
// char * PixelNomType()
//	Nom du Type des pixels de l'image. Chaine de caracteres = kint_2, kint_4, ...
// char * Nom()
//	Nom de l'image. "Ident()" retourne l' "id" associe a l'image
// void SetNameId(int imgid = 0, char const * nom = NULL)
//	Permet de modifier l'identificateur et le nom d'une image
// void SetAtt(int nbnul = -1, int nbsat = -1, -
//             float minpix = 0, float maxpix = -1., -
//             float moypix = 0, float sigpix = -1., -
//             float fnd = -1., float sigfnd = -1.)
//	Des donnees membres du type reels (float) 
//	"fond, sigmaFond, minPix, moyPix,sigPix" et entier (int) "nbNul, nbSat"
//	sont associees a chaque image. Leurs valeurs peuvent etre modifiees
//	a l'aide de la methode "SetAtt()" ou la methode "CheckDyn()" 
// int   CheckDyn(double min=-9.e19, double max=9.e19)
//	Cette methode retourne le nombre de pixels hors dynamique "( <min ou >max )"
//	Met a jour les donnees membres "min/maxPix nbNul/Sat".
//	Calcule aussi la moyenne et le sigma des pixels -> "moyPix, sigPix"
//--

//++
// Titre	I/O
//--
//++
// void  Write(string const& flnm)
// void  Save(char *flnm)
//	Ecrit l'image en format PPF (PPersist) dans le fichier "flnm"
// void  Read(string const& flnm)
//	Lit le fichier "flnm" (au format PPF)
//--

/* --Methode-- */
DVList&  RzImage::Info()
{
if (mInfo == NULL)  mInfo = new DVList;
return(*mInfo);
}

/* --Methode-- */
void RzImage::SetOrg(int_4 orgx, int_4 orgy)
{
org_x = orgx;   org_y = orgy;
return;
}

/* --Methode-- */
void RzImage::SetPxSize(float szx, float szy)
{
if (szx <= 0.)  szx = 1.;
if (szy <= 0.)  szy = 1.;
pxsz_x = szx;   pxsz_y = szy;
return;
}

/* --Methode-- */
void RzImage::SetAtt(int nbnul, int nbsat,
              float minpix, float maxpix,
              float moypix, float sigpix,
              float fnd, float sigfnd)
{
nbNul = nbnul;  nbSat = nbsat;
minPix = minpix;   maxPix = maxpix;
moyPix = moypix; sigPix = sigpix;
fond = fnd;  sigmaFond = sigfnd;
return;
}


/* --Methode-- */
void RzImage::SetNameId(int imgid, char const * nom)
{
id = imgid;
if (nom)  { strncpy(name, nom, 32);  name[31] = '\0'; }
else name[0] = '\0';
return;
}

/* --Methode-- */
void RzImage::Save(char *flnm)
{
POutPersist s(flnm);
Write(s);
}

/* --Methode-- */
void RzImage::WriteSelf(POutPersist& s) const

/*  Ecriture de la structure image ds un fichier de maniere */ 
/* brute - (Dump memoire)                                   */
/* Pour utilisation par PPersist::Write()                   */

{
char buf[256];
unsigned long sz;
int i,j;
int_4 itab[5];

sprintf(buf,"\nRzImage(%s)[%d*%d]  Org: %d %d  PxSz: %g %g (Id=%d Nom=%s)\n",
        PixelNomType(),siz_x, siz_y, org_x, org_y, pxsz_x, pxsz_y, id, name);
i = strlen(buf); 
sprintf(buf+i,"NbSat,Nul=%d %d Min,Max= %g %g Moy,Sig= %g %g Fond,Sig= %g %g\n",
        nbSat, nbNul, minPix, maxPix, moyPix, sigPix, fond, sigmaFond);
i += strlen(buf+i);
for(j=i; j<255; j++)  buf[j] = ' ';  buf[255] = '\0';
s.PutBytes(buf, 256);

itab[0] = (int)dataType;
itab[1] = siz_x;  itab[2] = siz_y;
itab[3] = itab[4] = 0; 
if (mInfo)  itab[3] = 1;    // Has DVList as info
s.PutI4s(itab, 5);

s.PutR4(fond);
s.PutR4(sigmaFond);
s.PutI4(nbNul);
s.PutI4(nbSat);
s.PutR4s(&minPix, 4); // minPix -> sigPix

s.PutI4(org_x);
s.PutI4(org_y);
s.PutR4(pxsz_x);
s.PutR4(pxsz_y);
s.PutI4(isFits);
s.PutBytes(name, 32);
s.PutI4(id);

sz = XSize()*YSize();
switch (dataType) {
case kuint_1:
  s.PutBytes(voidP, sz);
  break;
case kuint_2:
  s.PutU2s((uint_2*)voidP, sz);
  break;
case kint_2:
  s.PutI2s((int_2*)voidP, sz);
  break;
case kint_4:
  s.PutI4s((int_4*)voidP, sz);
  break;
case kr_4:
  s.PutR4s((r_4*)voidP, sz);
  break;
case kr_8:
  s.PutR8s((r_8*)voidP, sz);
  break;
default:
  printf("RzImage::WriteSelf Error - Unsupported DataType= %d \n", (int)dataType);
  break;    
} 

if (mInfo) mInfo->Write(s);
return;
}


/* --Methode-- */
void RzImage::ReadSelf(PInPersist& s)
/*  Relecture de fichier cree par WriteSelf()    */
/*  Pour utilisation par PPersist::Read()        */
{
char buf[256];
unsigned long sz;
int_4 itab[5];
int_4 opt;

s.GetBytes(buf, 256);
s.GetI4s(itab, 5);
opt = (int_4)PixelType();
if ((opt > 0) && (opt != itab[0]))  { 
    cerr << "RzImage::ReadSelf() Image type (=" << opt 
         << ") mismatch (type in file=" << itab[0] << ")" << endl;
    THROW(typeMismatchErr); }
Allocate((PBaseDataTypes)itab[0], itab[1], itab[2]);

s.GetR4(fond);
s.GetR4(sigmaFond);
s.GetI4(nbNul);
s.GetI4(nbSat);
s.GetR4s(&minPix, 4); // minPix -> sigPix

s.GetI4(org_x);
s.GetI4(org_y);
s.GetR4(pxsz_x);
s.GetR4(pxsz_y);
s.GetI4(isFits);
s.GetBytes(name, 32);
s.GetI4(id);

sz = XSize()*YSize();
switch (dataType) {
case kuint_1:
  s.GetBytes(voidP, sz);
  break;
case kuint_2:
  s.GetU2s((uint_2*)voidP, sz);
  break;
case kint_2:
  s.GetI2s((int_2*)voidP, sz);
  break;
case kint_4:
  s.GetI4s((int_4*)voidP, sz);
  break;
case kr_4:
  s.GetR4s((r_4*)voidP, sz);
  break;
case kr_8:
  s.GetR8s((r_8*)voidP, sz);
  break;
default:
  printf("RzImage::ReadSelf Error - Unsupported DataType= %d \n", (int)dataType);
  break;    
} 

if (itab[3])  // has info
  { if (mInfo == NULL)  mInfo = new DVList;
  mInfo->Read(s); }

return;
}


//++
// Titre	Impression
//--
//++
// void  Print(ostream& os)
//	Imprime sur le flot "os" les informations globales concernant l'image
// void  Print()
//	Imprime (sur cout) les informations globales concernant l'image
// ostream&   operator << (ostream& s, RzImage const & img)
//      sortie sur flot "s" (Appel a "Print(s)").
//--

/* --Methode-- */
void RzImage::Print(ostream& os) const
{
char buff[256];
os << " ---- RzImage::Print() ----- \n";
sprintf(buff, "RzImage(%s)[%d*%d]  Org: %d %d  PxSz= %g %g (Id=%d Nom=%s)\n",
       PixelNomType(),siz_x, siz_y, org_x, org_y, pxsz_x, pxsz_y, id, name);
os << buff;
sprintf(buff, "NbSat,Nul=%d %d Min,Max= %g %g Moy,Sig= %g %g Fond,Sig= %g %g\n",
       nbSat, nbNul, minPix, maxPix, moyPix, sigPix, fond, sigmaFond);
os << buff;
if (mInfo) mInfo->Print(os);
os << endl;
}

/* --Methode-- */
int RzImage::CheckDyn(double min, double max)

/*  Cette methode retourne le nombre de pixels hors dynamique    */
/* ds le pave ( <min ou >max )                                   */
/* Met a jour la structure image Min/MaxPix nbNul/Sat            */
/* Calcule aussi la moyenne et le sigma des pixels, met a jour   */
/* moypix et sigpix dans la structure                            */

{
float pixv=0;
double dv;

float minadu = PutInRange(min, pixv); 
float maxadu = PutInRange(max, pixv);
float minpix = PutInRange(maxadu+1, pixv);
float maxpix = PutInRange(minadu-1, pixv);

double m = 0.;
double s = 0.;
int n9 = 0;
int n0 = 0;
int n = XSize()*YSize();
for(int i=0; i<n; i++) 
  {
  pixv = FValue(i);
  if (pixv <= maxadu && pixv >= minadu) 
    {  
    if (pixv < minpix)  minpix = pixv;
    else if (pixv > maxpix)  maxpix = pixv; 
    dv = (double)pixv;
    m += dv;  s += (dv*dv);
    }
  else 
    {  
    if (pixv > maxadu)  n9++;
    if (pixv < minadu)  n0++; 
    }
  }

nbNul = n0;
nbSat = n9;
minPix = minpix;
maxPix = maxpix;

n -= (n0+n9);
if (n > 0)
  { dv = (double)n;  m /= dv;  s = s/dv - m*m;
  if (s>0.)  s = sqrt(s);  else s = 0.;
  moyPix = m;  sigPix = s; }
else { moyPix = 0.;  sigPix = -1.; }

return(n0+n9);
}

////////////////////////////////////////////////////////////////
//++
RzImage* RzImage::FitResidus(GeneralFit& gfit)
//
//	Retourne une RzImage R4 contenant les residus du fit ``gfit''.
//--
{
if(XSize()<=0 || YSize()<=0) return NULL;
GeneralFunction* f = gfit.GetFunction();
if(f==NULL) return NULL;
Vector par = gfit.GetParm();
RzImage* img = new RzImage(DataType((r_4)0.),XSize(),YSize());
img->SetOrg(XOrg(),YOrg());
img->SetPxSize(XPxSize(),YPxSize());
r_4* ptr = (r_4 *) img->VoidP();
for(int i=0;i<XSize();i++) {
  double x[2];
  x[0] = XOrg()+(i+0.5)*XPxSize();
  for(int j=0;j<YSize();j++) {
    x[1] = YOrg()+(j+0.5)*YPxSize();
    *(ptr+i+j*img->siz_x) = FValue(i,j) - (float) f->Value(x,par.Data());
  }
}
return img;
}

//++
RzImage* RzImage::FitFunction(GeneralFit& gfit)
//
//	Retourne RzImage R4 contenant la fonction du fit ``gfit''.
//--
{
if(XSize()<=0 || YSize()<=0) return NULL;
GeneralFunction* f = gfit.GetFunction();
if(f==NULL) return NULL;
Vector par = gfit.GetParm();
RzImage* img = new RzImage(DataType((r_4)0.),XSize(),YSize());
img->SetOrg(XOrg(),YOrg());
img->SetPxSize(XPxSize(),YPxSize());
r_4* ptr = (r_4 *) img->VoidP();
for(int i=0;i<XSize();i++) {
  double x[2];
  x[0] = XOrg()+(i+0.5)*XPxSize();
  for(int j=0;j<YSize();j++) {
    x[1] = YOrg()+(j+0.5)*YPxSize();
    *(ptr+i+j*img->siz_x) = (float) f->Value(x,par.Data());
  }
}
return img;
}

//++
// Titre	Acces aux pixels
//	*Attention:* Les methodes IValue() FValue() et DValue() sont lentes.
//	Ne pas les utiliser dans les parties du code ou la vitesse de traitement
//	est critique.
//--
//++
// void * VoidP()
//	Renvoie un pointeur void sur le tableau des pixels. Doit etre converti
//	dans le bon type suivant le contenu de l'image.
// int IValue(int i, int j)
//	Renvoie la valeur du pixel "(i,j)" converti en int
// float FValue(int i, int j)
//	Renvoie la valeur du pixel "(i,j)" converti en float
// double DValue(int i, int j)
//	Renvoie la valeur du pixel "(i,j)" converti en double
// int IValue(int k)
//	Renvoie la valeur du pixel numero "k" converti en int
// float FValue(int k)
//	Renvoie la valeur du pixel numero "k" converti en float
// double DValue(int k)
//	envoie la valeur du pixel numero "k" converti en double
//
//	Les pixels sont rangees ligne apres ligne 
//|	{(1,1), (2,1), ... (m, 1), (1,2), ... , (m,2), .... (m,n)}
//--

/* --Methodes-- */
int RzImage::IVal_u1(int k) const          {return vect.u1->IElem(k);}
int RzImage::IVal_u2(int k) const          {return vect.u2->IElem(k);}
int RzImage::IVal_i2(int k) const          {return vect.i2->IElem(k);}
int RzImage::IVal_i4(int k) const          {return vect.i4->IElem(k);}
int RzImage::IVal_r4(int k) const          {return vect.r4->IElem(k);}
int RzImage::IVal_r8(int k) const          {return vect.r8->IElem(k);}

float RzImage::FVal_u1(int k) const        {return vect.u1->FElem(k);}
float RzImage::FVal_u2(int k) const        {return vect.u2->FElem(k);}
float RzImage::FVal_i2(int k) const        {return vect.i2->FElem(k);}
float RzImage::FVal_i4(int k) const        {return vect.i4->FElem(k);}
float RzImage::FVal_r4(int k) const        {return vect.r4->FElem(k);}
float RzImage::FVal_r8(int k) const        {return vect.r8->FElem(k);}

double RzImage::DVal_u1(int k) const       {return vect.u1->DElem(k);}
double RzImage::DVal_u2(int k) const       {return vect.u2->DElem(k);}
double RzImage::DVal_i2(int k) const       {return vect.i2->DElem(k);}
double RzImage::DVal_i4(int k) const       {return vect.i4->DElem(k);}
double RzImage::DVal_r4(int k) const       {return vect.r4->DElem(k);}
double RzImage::DVal_r8(int k) const       {return vect.r8->DElem(k);}



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

RzVect<uint_1> rzvu1(n);
RzVect<uint_2> rzvu2(n);
RzVect<int_2> rzvi2(n);
RzVect<int_4> rzvi4(n);
RzVect<r_4> rzvr4(n);
RzVect<r_8> rzvr8(n);
return;
}

#endif

#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template RzVect<uint_1>
#pragma define_template RzVect<uint_2>
#pragma define_template RzVect<int_2>
#pragma define_template RzVect<int_4>
#pragma define_template RzVect<r_4>
#pragma define_template RzVect<r_8>
#endif

#ifdef GNU_TEMPLATES
template class RzVect<uint_1>;
template class RzVect<uint_2>;
template class RzVect<int_2>;
template class RzVect<int_4>;
template class RzVect<r_4>;
template class RzVect<r_8>;
#endif

#if defined(ANSI_TEMPLATES) 
template class RzVect<uint_1>;
template class RzVect<uint_2>;
template class RzVect<int_2>;
template class RzVect<int_4>;
template class RzVect<r_4>;
template class RzVect<r_8>;
#endif