source: Sophya/trunk/SophyaLib/NTools/dynccd.cc@ 2639

#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        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_4TypNoise=0, r_8 MinADU=-9.e19, r_8 MaxADU=9.e19, r_8 Gain=1., r_8 RONoise=0., r_8 RefFond=0., r_8 RefSFond=0.);
//      Creation d'un objet DynCCD ("typ" determine la methode de calcul du bruit)
//      |Test verbatim 
//
// void Set(int_4TypNoise=0, r_8 MinADU=-9.e19, r_8 MaxADU=9.e19, r_8 Gain=1., r_8 RONoise=0., r_8 RefFond=0., r_8 RefSFond=0.);
//      Modification des parametres de la dynamique
// void Print()
// r_8 Noise(r_8 pixel) const
//      Renvoie la valeur du bruit pour "pixel" en ADU.
//--

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

/* --Methode-- */
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-- */
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;
}

/* --Methode-- */
r_8 DynCCD::Noise(r_8 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    */

{
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      */
/* --------------------------------------------------------------  */

//++
// 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                                       */

{
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 */
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__)
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&);
#endif