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

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

#include "cimage.h"

#include "dynccd.h"

/*! 
  \class SOPHYA::DynCCD
  \ingroup NTools
  
  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.

  \sa SOPHYA::Image
*/


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

/* --Methode-- */
//!  Creation d'un objet DynCCD ("typ" determine la methode de calcul du bruit)
DynCCD::DynCCD(int_4 typ, r_8 min, r_8 max, r_8 g,
               r_8 ron, r_8 rf, r_8 rfs)
{
if ( (typ >= kConstantNoise) && (typ <= kSqrtADUNoise) )  TypNoise = typ;
else  TypNoise = kConstantNoise;
MinADU = min;  MaxADU = max;
Gain = g;  RONoise = ron;
RefFond = rf;  RefSFond = rfs;
}

/* --Methode-- */
//!	Modification des parametres de la dynamique
void DynCCD::Set(int_4 typ, r_8 min, r_8 max, r_8 g,
               r_8 ron, r_8 rf, r_8 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;
}

//!	Renvoie la valeur du bruit pour "pixel" en ADU.
/*!
  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  
*/
/* --Methode-- */
r_8 DynCCD::Noise(r_8 pixel) const
{
r_8 h,s,ronsq;
r_8 fond;

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

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

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

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


/* Nouvelle-Fonction */
/*!
  Construit et renvoie l'image du bruit pour l'image "*pim" 
  
  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                                       
*/
template <class T>
Image<T> NoiseImage(Image<T> const & pim, DynCCD const & dynccd)

{
r_8 h,s,ronsq;
r_8 fond, min,max;
int_4 k, npix;
r_8 minois, offnois;

npix = pim.XSize()*pim.YSize();
Image<T> nois(pim.XSize(), pim.YSize());
nois.SetOrg(pim.XOrg(), pim.YOrg());

min = dynccd.MinADU;   max = dynccd.MaxADU;
for(k=0; k<npix; k++)      {
  if ( (pim[k] <= max) && (pim[k] >= min) )  nois[k] = (T)1;
  else nois[k] = 0;
}


return(nois);
}


/* Nouvelle-Fonction */
//!	Calcule l'image du bruit et le rajoute a l'image originale
template <class T>
void ImgAddNoise(Image<T>& img, DynCCD const& dyn)

{
  int_4 nPix = img.XSize() * img.YSize();
	
  for (int_4 i=0; i<nPix; i++)
    img[i] += (T) (dyn.Noise(img[i])*NorRand());
}





#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template NoiseImage<uint_2>
#pragma define_template NoiseImage<int_4>
#pragma define_template NoiseImage<r_4>

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

#if defined(ANSI_TEMPLATES) || defined(__GNU_TEMPLATES__)
#if defined( __IBMCPP__ )
namespace SOPHYA {
#endif
template Image<uint_2> NoiseImage<uint_2>(Image<uint_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_4>(Image< int_4>&, DynCCD const&);
template void ImgAddNoise<   r_4>(Image<   r_4>&, DynCCD const&);
#if defined( __IBMCPP__ )
}
#endif

#endif