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

#include "fmath.h"
#include "perandom.h"

#include "cimage.h"

#include "dynccd.h"

//++
// Class	DynCCD
// Lib		Images++ 
// include	dynccd.h
// 
//	Cette classe permet la specification des parametres 
//	definissant la dynamique du CCD, et doit etre utilise
//      pour le calcul des images de bruit.
//	- TypNoise = 1 :
//	Bruit = Sqrt( (RONoise/Gain)**2 + ValPix/Gain )      
//	- TypNoise = 2 :                                     
//	Bruit = Sqrt( RefSFond**2 + (ValPix-RefFond)/Gain )  
//	- TypNoise = 0 
//	Bruit = 1  (Constant pour tous les pixels) 
//	- TypNoise = 3         
//	Bruit = Sqrt(Abs(ValPix)) 
//
//	Les pixels hors dynamique sont marques (Bruit = 0)
//	( < MinADU ou > MaxADU) pour toutes valeurs de TypNoise.
//	Gain est exprime en electron par ADU, RONoise en electron,
//	Bruit, ValPix et RefFond en ADU.
//--

//++
// Links	Parents
// PPersist
//--
//++
// Links	Autres
// RzImage
// Image<T> 
//--
//++
// Titre	Methodes
//--
//++
// DynCCD(int TypNoise=0, float MinADU=-9.e19, float MaxADU=9.e19, float Gain=1., float RONoise=0., float RefFond=0., float RefSFond=0.);
//	Creation d'un objet DynCCD ("typ" determine la methode de calcul du bruit)
//	|Test verbatim 
//
// void Set(int TypNoise=0, float MinADU=-9.e19, float MaxADU=9.e19, float Gain=1., float RONoise=0., float RefFond=0., float RefSFond=0.);
//	Modification des parametres de la dynamique
// void Print()
// float Noise(float pixel) const
//	Renvoie la valeur du bruit pour "pixel" en ADU.
//--

/* ............................................................ */
/* :::::::::::::  methode de la classe  DynCCD :::::::::::::::: */
/* ............................................................ */

/* --Methode-- */
DynCCD::DynCCD(int typ, float min, float max, float g,
               float ron, float rf, float rfs)
{
if ( (typ >= kConstantNoise) && (typ <= kSqrtADUNoise) )  TypNoise = typ;
else  TypNoise = kConstantNoise;
MinADU = min;  MaxADU = max;
Gain = g;  RONoise = ron;
RefFond = rf;  RefSFond = rfs;
}

/* --Methode-- */
void DynCCD::Set(int typ, float min, float max, float g,
               float ron, float rf, float rfs)
{
if ( (typ >= kConstantNoise) && (typ <= kSqrtADUNoise) )  TypNoise = typ;
MinADU = min;  MaxADU = max;
Gain = g;  RONoise = ron;
RefFond = rf;  RefSFond = rfs;
}

/* --Methode-- */
void DynCCD::Print()
{
printf("DynCCD: Type= %d Min/MaxADU= %g %g Gain/RoN= %g %g\n",
       TypNoise, MinADU, MaxADU, Gain, RONoise);
if (TypNoise == 2)
  printf("... RefFond= %g  RefSFond= %g \n", RefFond, RefSFond);
return;
}

/* --Methode-- */
float DynCCD::Noise(float pixel) const

/* Cette fonction calcule la valeur du bruit pour pixel     */
/*  Si TypNoise = 1 :                                       */
/*     Bruit = Sqrt( (RONoise/Gain)**2 + ValPix/Gain )      */
/*  Si TypNoise = 2 :                                       */
/*     Bruit = Sqrt( RefSFond**2 + (ValPix-RefFond)/Gain )  */
/*  Si TypNoise = 0                                         */
/*     Bruit = 1  (Constant pour tous les pixels)           */
/*  Si TypNoise = 3                                         */
/*     Bruit = Sqrt(Abs(PixelADU))                          */
/*  Les pixels hors dynamique sont marques (Bruit = 0)      */
/*  ( < MinADU ou > MaxADU) pour tout valeur de TypNoise    */

{
float h,s,ronsq;
float fond;

if ( (pixel > MaxADU) || (pixel < MinADU) )   return(0.);
if ( TypNoise == kConstantNoise)  return(1.);
if ( TypNoise == kSqrtADUNoise )  return(fsqrt(fabsf(pixel)));

if ( TypNoise == kSigFondNoise)
  { fond = RefFond;
  ronsq = RefSFond * RefSFond; }
else
  { fond = 0;
  ronsq = RONoise/Gain;  ronsq *= ronsq; }

h = (pixel>fond) ? (float)(pixel-fond) : 0.;
s = ronsq+h/Gain;
s = fsqrt(s);
return(s);
}

/* --------------------------------------------------------------  */
/*      Quelques fonctions pour manipuler des images de bruit      */
/* --------------------------------------------------------------  */

//++
// Module	Images de bruit 
// Lib		Images++ 
// include	dynccd.h
//
//	Ces fonctions permettent le calcul d'image de bruit a partir d'une 
//	image (RzImage ou Image<T>) et d'un objet DynCCD 
//--
//++
// Links 	Voir classes
// DynCCD
// RzImage
// Image<T>
//--
//++
// Titre	Les fonctions 
//--

//++
// Function	RzImage * NoiseImage(RzImage const *pim, DynCCD const * dynccd)
//	Construit et renvoie l'image du bruit pour l'image "*pim" (RzImage)
// Function	Image<T> * NoiseImage(Image<T> const * pim, DynCCD const * dynccd)
//      Meme fonctionalite pour une image typee (ImageU2, ImageR4, ...)
// Function	ImgAddNoise(Image<T>&, DynCCD const&)
//	Calcule l'image du bruit et le rajoute a l'image originale
//--

/* Nouvelle-Fonction */
template <class T>
Image<T> * NoiseImage(Image<T> const * pim, DynCCD const * dynccd)

/* Creation et Calcul d'une image de bruit a partir de l'image */
/* pim et de dynccd. Voir la methode DynCCD::Noise() pour la   */
/* description du calcul                                       */

{
float h,s,ronsq;
float fond, min,max;
int i, npix;
float minois, offnois;

if (pim == NULL)  return(NULL);

const T * pix = pim->ImagePtr();
npix = pim->XSize()*pim->YSize();
Image<T> * nois = new Image<T>(pim->XSize(), pim->YSize(), false);
nois->SetOrg(pim->XOrg(), pim->YOrg());
T * pno = nois->ImagePtr();

min = dynccd->MinADU;   max = dynccd->MaxADU;


switch (dynccd->TypNoise)
  {
  case kConstantNoise :
    for(i=0; i<npix; i++)
      {
      if ( (*pix <= max) && (*pix >= min) )  *pno = 1;
      else *pno = 0;
      pix++;   pno++;
      }
    break;

  case kSqrtADUNoise :
    for(i=0; i<npix; i++)
      {
      if ( (*pix <= max) && (*pix >= min) )  *pno = 1;
      else *pno = (T) fsqrt(fabsf((float)(*pix)));    
      pix++;   pno++;
      }
    break;

  case kPhotonNoise :
  case kSigFondNoise :
    if ( dynccd->TypNoise == kSigFondNoise)
      { fond = dynccd->RefFond;
      ronsq = dynccd->RefSFond * dynccd->RefSFond; }
    else
      { fond = 0;
      ronsq = dynccd->RONoise/dynccd->Gain;  ronsq *= ronsq; }

// Calcul de minois / offnois pour obtenir un bruit correct malgre 
// les conversions (float) -> (entier)

    switch(pim->PixelType())
      {
      case kuint_2:
      case kint_2:
      case kint_4:
        minois = 1.001;
        offnois = 0.5;
        break;
      case kr_4:
      case kr_8:
        minois = 1.e-9;
        offnois = 0.;
        break;
      default:
        minois = 1.e-9;
        offnois = 0.;
        break;
      }

    for(i=0; i<npix; i++)
      {
      if ( (*pix <= max) && (*pix >= min) ) 
        {
        h = (*pix>fond) ? (float)(*pix)-fond : 0.;
        s = ronsq+h/dynccd->Gain;
        s = fsqrt(s)+offnois;
        *pno = (s > minois) ? (T)s : (T)minois;
        } 
      else *pno = 0;
      pix++;   pno++;
      }
    break;
  }

return(nois);
}


/* Nouvelle-Fonction */
template <class T>
void ImgAddNoise(Image<T>& img, DynCCD const& dyn)
{
	T* p = img.ImagePtr();
	int nPix = img.XSize() * img.YSize();
	
	for (int i=0; i<nPix; i++, p++)
		*p += (T) (dyn.Noise(*p)*NorRand());
}



/* Nouvelle-Fonction */
RzImage * NoiseImage(RzImage const * pim, DynCCD const * dynccd)

/* Creation et Calcul d'une image de bruit a partir de l'image */
/* pim et de dynccd. Voir la methode DynCCD::Noise() pour la   */
/* description du calcul                                       */
{
RzImage * nois = NULL;

switch (pim->PixelType()) 
  {
  case kuint_2: 
    {
    ImageU2 pix((RzImage&)(*pim));
    nois = NoiseImage(&pix, dynccd);
    break;
    }
  case kint_2:
    {
    ImageI2 pix((RzImage&)(*pim));
    nois = NoiseImage(&pix, dynccd);
    break;
    }
  case kint_4:
    {
    ImageI4 pix((RzImage&)(*pim));
    nois = NoiseImage(&pix, dynccd);
    break;
    }
  case kr_4:
    {
    ImageR4 pix((RzImage&)(*pim));
    nois = NoiseImage(&pix, dynccd);
    break;
    }  
  case kuint_4:
  case kr_8:
    cerr << "NoiseImage(RzImage, ...) Unsupported image type (kuint_4/kr_8) " 
         << int(pim->PixelType()) << "\n" ; 
    break;
  default:
    cerr << "NoiseImage(RzImage, ...) Unknown image type !! " << int(pim->PixelType()) << "\n" ;
    break;
  }   

return(nois);
}


// ******** INSTANCES

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

ImageU2  iu2(n,n), *piu2;
ImageI2  ii2(n,n), *pii2;
ImageI4  ii4(n,n), *pii4;
ImageR4  ir4(n,n), *pir4;
DynCCD dyn(1,0.,65000., 4., 20.);

piu2 = NoiseImage(&iu2, &dyn);
pii2 = NoiseImage(&ii2, &dyn);
pii4 = NoiseImage(&ii4, &dyn);
pir4 = NoiseImage(&ir4, &dyn);

ImgAddNoise(*piu2, dyn);
ImgAddNoise(*pii2, dyn);
ImgAddNoise(*pii4, dyn);
ImgAddNoise(*pir4, dyn);

return;
}
#endif


#ifdef __CXX_PRAGMA_TEMPLATES__

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

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

#endif


#if ( defined(ANSI_TEMPLATES) && !defined(__aCC__) )
template Image<uint_2> * NoiseImage<uint_2>(Image<uint_2> const * , DynCCD const *);
template Image< int_2> * NoiseImage< int_2>(Image< int_2> const * , DynCCD const *);
template Image< int_4> * NoiseImage< int_4>(Image< int_4> const * , DynCCD const *);
template Image<   r_4> * NoiseImage<   r_4>(Image<   r_4> const * , DynCCD const *);

template void ImgAddNoise<uint_2>(Image<uint_2>&, DynCCD const&);
template void ImgAddNoise< int_2>(Image< int_2>&, DynCCD const&);
template void ImgAddNoise< int_4>(Image< int_4>&, DynCCD const&);
template void ImgAddNoise<   r_4>(Image<   r_4>&, DynCCD const&);

#endif

#if defined(__GNU_TEMPLATES__)
template Image<uint_2> * NoiseImage(Image<uint_2> const *, DynCCD const *);
template Image< int_2> * NoiseImage(Image< int_2> const *, DynCCD const *);
template Image< int_4> * NoiseImage(Image< int_4> const *, DynCCD const *);
template Image<   r_4> * NoiseImage(Image<   r_4> const *, DynCCD const *);

template void ImgAddNoise(Image<uint_2>&, DynCCD const&);
template void ImgAddNoise(Image< int_2>&, DynCCD const&);
template void ImgAddNoise(Image< int_4>&, DynCCD const&);
template void ImgAddNoise(Image<   r_4>&, DynCCD const&);

#endif