| 1 | #include "machdefs.h"
|
|---|
| 2 | #include <stdlib.h>
|
|---|
| 3 | #include <stdio.h>
|
|---|
| 4 |
|
|---|
| 5 | #include "fmath.h"
|
|---|
| 6 | #include "perandom.h"
|
|---|
| 7 |
|
|---|
| 8 | #include "cimage.h"
|
|---|
| 9 |
|
|---|
| 10 | #include "dynccd.h"
|
|---|
| 11 |
|
|---|
| 12 | //++
|
|---|
| 13 | // Class DynCCD
|
|---|
| 14 | // Lib Images++
|
|---|
| 15 | // include dynccd.h
|
|---|
| 16 | //
|
|---|
| 17 | // Cette classe permet la specification des parametres
|
|---|
| 18 | // definissant la dynamique du CCD, et doit etre utilise
|
|---|
| 19 | // pour le calcul des images de bruit.
|
|---|
| 20 | // - TypNoise = 1 :
|
|---|
| 21 | // Bruit = Sqrt( (RONoise/Gain)**2 + ValPix/Gain )
|
|---|
| 22 | // - TypNoise = 2 :
|
|---|
| 23 | // Bruit = Sqrt( RefSFond**2 + (ValPix-RefFond)/Gain )
|
|---|
| 24 | // - TypNoise = 0
|
|---|
| 25 | // Bruit = 1 (Constant pour tous les pixels)
|
|---|
| 26 | // - TypNoise = 3
|
|---|
| 27 | // Bruit = Sqrt(Abs(ValPix))
|
|---|
| 28 | //
|
|---|
| 29 | // Les pixels hors dynamique sont marques (Bruit = 0)
|
|---|
| 30 | // ( < MinADU ou > MaxADU) pour toutes valeurs de TypNoise.
|
|---|
| 31 | // Gain est exprime en electron par ADU, RONoise en electron,
|
|---|
| 32 | // Bruit, ValPix et RefFond en ADU.
|
|---|
| 33 | //--
|
|---|
| 34 |
|
|---|
| 35 | //++
|
|---|
| 36 | // Links Parents
|
|---|
| 37 | // PPersist
|
|---|
| 38 | //--
|
|---|
| 39 | //++
|
|---|
| 40 | // Links Autres
|
|---|
| 41 | // RzImage
|
|---|
| 42 | // Image<T>
|
|---|
| 43 | //--
|
|---|
| 44 | //++
|
|---|
| 45 | // Titre Methodes
|
|---|
| 46 | //--
|
|---|
| 47 | //++
|
|---|
| 48 | // DynCCD(int TypNoise=0, float MinADU=-9.e19, float MaxADU=9.e19, float Gain=1., float RONoise=0., float RefFond=0., float RefSFond=0.);
|
|---|
| 49 | // Creation d'un objet DynCCD ("typ" determine la methode de calcul du bruit)
|
|---|
| 50 | // |Test verbatim
|
|---|
| 51 | //
|
|---|
| 52 | // 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.);
|
|---|
| 53 | // Modification des parametres de la dynamique
|
|---|
| 54 | // void Print()
|
|---|
| 55 | // float Noise(float pixel) const
|
|---|
| 56 | // Renvoie la valeur du bruit pour "pixel" en ADU.
|
|---|
| 57 | //--
|
|---|
| 58 |
|
|---|
| 59 | /* ............................................................ */
|
|---|
| 60 | /* ::::::::::::: methode de la classe DynCCD :::::::::::::::: */
|
|---|
| 61 | /* ............................................................ */
|
|---|
| 62 |
|
|---|
| 63 | /* --Methode-- */
|
|---|
| 64 | DynCCD::DynCCD(int typ, float min, float max, float g,
|
|---|
| 65 | float ron, float rf, float rfs)
|
|---|
| 66 | {
|
|---|
| 67 | if ( (typ >= kConstantNoise) && (typ <= kSqrtADUNoise) ) TypNoise = typ;
|
|---|
| 68 | else TypNoise = kConstantNoise;
|
|---|
| 69 | MinADU = min; MaxADU = max;
|
|---|
| 70 | Gain = g; RONoise = ron;
|
|---|
| 71 | RefFond = rf; RefSFond = rfs;
|
|---|
| 72 | }
|
|---|
| 73 |
|
|---|
| 74 | /* --Methode-- */
|
|---|
| 75 | void DynCCD::Set(int typ, float min, float max, float g,
|
|---|
| 76 | float ron, float rf, float rfs)
|
|---|
| 77 | {
|
|---|
| 78 | if ( (typ >= kConstantNoise) && (typ <= kSqrtADUNoise) ) TypNoise = typ;
|
|---|
| 79 | MinADU = min; MaxADU = max;
|
|---|
| 80 | Gain = g; RONoise = ron;
|
|---|
| 81 | RefFond = rf; RefSFond = rfs;
|
|---|
| 82 | }
|
|---|
| 83 |
|
|---|
| 84 | /* --Methode-- */
|
|---|
| 85 | void DynCCD::Print()
|
|---|
| 86 | {
|
|---|
| 87 | printf("DynCCD: Type= %d Min/MaxADU= %g %g Gain/RoN= %g %g\n",
|
|---|
| 88 | TypNoise, MinADU, MaxADU, Gain, RONoise);
|
|---|
| 89 | if (TypNoise == 2)
|
|---|
| 90 | printf("... RefFond= %g RefSFond= %g \n", RefFond, RefSFond);
|
|---|
| 91 | return;
|
|---|
| 92 | }
|
|---|
| 93 |
|
|---|
| 94 | /* --Methode-- */
|
|---|
| 95 | float DynCCD::Noise(float pixel) const
|
|---|
| 96 |
|
|---|
| 97 | /* Cette fonction calcule la valeur du bruit pour pixel */
|
|---|
| 98 | /* Si TypNoise = 1 : */
|
|---|
| 99 | /* Bruit = Sqrt( (RONoise/Gain)**2 + ValPix/Gain ) */
|
|---|
| 100 | /* Si TypNoise = 2 : */
|
|---|
| 101 | /* Bruit = Sqrt( RefSFond**2 + (ValPix-RefFond)/Gain ) */
|
|---|
| 102 | /* Si TypNoise = 0 */
|
|---|
| 103 | /* Bruit = 1 (Constant pour tous les pixels) */
|
|---|
| 104 | /* Si TypNoise = 3 */
|
|---|
| 105 | /* Bruit = Sqrt(Abs(PixelADU)) */
|
|---|
| 106 | /* Les pixels hors dynamique sont marques (Bruit = 0) */
|
|---|
| 107 | /* ( < MinADU ou > MaxADU) pour tout valeur de TypNoise */
|
|---|
| 108 |
|
|---|
| 109 | {
|
|---|
| 110 | float h,s,ronsq;
|
|---|
| 111 | float fond;
|
|---|
| 112 |
|
|---|
| 113 | if ( (pixel > MaxADU) || (pixel < MinADU) ) return(0.);
|
|---|
| 114 | if ( TypNoise == kConstantNoise) return(1.);
|
|---|
| 115 | if ( TypNoise == kSqrtADUNoise ) return(fsqrt(fabsf(pixel)));
|
|---|
| 116 |
|
|---|
| 117 | if ( TypNoise == kSigFondNoise)
|
|---|
| 118 | { fond = RefFond;
|
|---|
| 119 | ronsq = RefSFond * RefSFond; }
|
|---|
| 120 | else
|
|---|
| 121 | { fond = 0;
|
|---|
| 122 | ronsq = RONoise/Gain; ronsq *= ronsq; }
|
|---|
| 123 |
|
|---|
| 124 | h = (pixel>fond) ? (float)(pixel-fond) : 0.;
|
|---|
| 125 | s = ronsq+h/Gain;
|
|---|
| 126 | s = fsqrt(s);
|
|---|
| 127 | return(s);
|
|---|
| 128 | }
|
|---|
| 129 |
|
|---|
| 130 | /* -------------------------------------------------------------- */
|
|---|
| 131 | /* Quelques fonctions pour manipuler des images de bruit */
|
|---|
| 132 | /* -------------------------------------------------------------- */
|
|---|
| 133 |
|
|---|
| 134 | //++
|
|---|
| 135 | // Module Images de bruit
|
|---|
| 136 | // Lib Images++
|
|---|
| 137 | // include dynccd.h
|
|---|
| 138 | //
|
|---|
| 139 | // Ces fonctions permettent le calcul d'image de bruit a partir d'une
|
|---|
| 140 | // image (RzImage ou Image<T>) et d'un objet DynCCD
|
|---|
| 141 | //--
|
|---|
| 142 | //++
|
|---|
| 143 | // Links Voir classes
|
|---|
| 144 | // DynCCD
|
|---|
| 145 | // RzImage
|
|---|
| 146 | // Image<T>
|
|---|
| 147 | //--
|
|---|
| 148 | //++
|
|---|
| 149 | // Titre Les fonctions
|
|---|
| 150 | //--
|
|---|
| 151 |
|
|---|
| 152 | //++
|
|---|
| 153 | // Function RzImage * NoiseImage(RzImage const *pim, DynCCD const * dynccd)
|
|---|
| 154 | // Construit et renvoie l'image du bruit pour l'image "*pim" (RzImage)
|
|---|
| 155 | // Function Image<T> * NoiseImage(Image<T> const * pim, DynCCD const * dynccd)
|
|---|
| 156 | // Meme fonctionalite pour une image typee (ImageU2, ImageR4, ...)
|
|---|
| 157 | // Function ImgAddNoise(Image<T>&, DynCCD const&)
|
|---|
| 158 | // Calcule l'image du bruit et le rajoute a l'image originale
|
|---|
| 159 | //--
|
|---|
| 160 |
|
|---|
| 161 | /* Nouvelle-Fonction */
|
|---|
| 162 | template <class T>
|
|---|
| 163 | Image<T> * NoiseImage(Image<T> const * pim, DynCCD const * dynccd)
|
|---|
| 164 |
|
|---|
| 165 | /* Creation et Calcul d'une image de bruit a partir de l'image */
|
|---|
| 166 | /* pim et de dynccd. Voir la methode DynCCD::Noise() pour la */
|
|---|
| 167 | /* description du calcul */
|
|---|
| 168 |
|
|---|
| 169 | {
|
|---|
| 170 | float h,s,ronsq;
|
|---|
| 171 | float fond, min,max;
|
|---|
| 172 | int i, npix;
|
|---|
| 173 | float minois, offnois;
|
|---|
| 174 |
|
|---|
| 175 | if (pim == NULL) return(NULL);
|
|---|
| 176 |
|
|---|
| 177 | const T * pix = pim->ImagePtr();
|
|---|
| 178 | npix = pim->XSize()*pim->YSize();
|
|---|
| 179 | Image<T> * nois = new Image<T>(pim->XSize(), pim->YSize(), false);
|
|---|
| 180 | nois->SetOrg(pim->XOrg(), pim->YOrg());
|
|---|
| 181 | T * pno = nois->ImagePtr();
|
|---|
| 182 |
|
|---|
| 183 | min = dynccd->MinADU; max = dynccd->MaxADU;
|
|---|
| 184 |
|
|---|
| 185 |
|
|---|
| 186 | switch (dynccd->TypNoise)
|
|---|
| 187 | {
|
|---|
| 188 | case kConstantNoise :
|
|---|
| 189 | for(i=0; i<npix; i++)
|
|---|
| 190 | {
|
|---|
| 191 | if ( (*pix <= max) && (*pix >= min) ) *pno = 1;
|
|---|
| 192 | else *pno = 0;
|
|---|
| 193 | pix++; pno++;
|
|---|
| 194 | }
|
|---|
| 195 | break;
|
|---|
| 196 |
|
|---|
| 197 | case kSqrtADUNoise :
|
|---|
| 198 | for(i=0; i<npix; i++)
|
|---|
| 199 | {
|
|---|
| 200 | if ( (*pix <= max) && (*pix >= min) ) *pno = 1;
|
|---|
| 201 | else *pno = (T) fsqrt(fabsf((float)(*pix)));
|
|---|
| 202 | pix++; pno++;
|
|---|
| 203 | }
|
|---|
| 204 | break;
|
|---|
| 205 |
|
|---|
| 206 | case kPhotonNoise :
|
|---|
| 207 | case kSigFondNoise :
|
|---|
| 208 | if ( dynccd->TypNoise == kSigFondNoise)
|
|---|
| 209 | { fond = dynccd->RefFond;
|
|---|
| 210 | ronsq = dynccd->RefSFond * dynccd->RefSFond; }
|
|---|
| 211 | else
|
|---|
| 212 | { fond = 0;
|
|---|
| 213 | ronsq = dynccd->RONoise/dynccd->Gain; ronsq *= ronsq; }
|
|---|
| 214 |
|
|---|
| 215 | // Calcul de minois / offnois pour obtenir un bruit correct malgre
|
|---|
| 216 | // les conversions (float) -> (entier)
|
|---|
| 217 |
|
|---|
| 218 | switch(pim->PixelType())
|
|---|
| 219 | {
|
|---|
| 220 | case kuint_2:
|
|---|
| 221 | case kint_2:
|
|---|
| 222 | case kint_4:
|
|---|
| 223 | minois = 1.001;
|
|---|
| 224 | offnois = 0.5;
|
|---|
| 225 | break;
|
|---|
| 226 | case kr_4:
|
|---|
| 227 | case kr_8:
|
|---|
| 228 | minois = 1.e-9;
|
|---|
| 229 | offnois = 0.;
|
|---|
| 230 | break;
|
|---|
| 231 | default:
|
|---|
| 232 | minois = 1.e-9;
|
|---|
| 233 | offnois = 0.;
|
|---|
| 234 | break;
|
|---|
| 235 | }
|
|---|
| 236 |
|
|---|
| 237 | for(i=0; i<npix; i++)
|
|---|
| 238 | {
|
|---|
| 239 | if ( (*pix <= max) && (*pix >= min) )
|
|---|
| 240 | {
|
|---|
| 241 | h = (*pix>fond) ? (float)(*pix)-fond : 0.;
|
|---|
| 242 | s = ronsq+h/dynccd->Gain;
|
|---|
| 243 | s = fsqrt(s)+offnois;
|
|---|
| 244 | *pno = (s > minois) ? (T)s : (T)minois;
|
|---|
| 245 | }
|
|---|
| 246 | else *pno = 0;
|
|---|
| 247 | pix++; pno++;
|
|---|
| 248 | }
|
|---|
| 249 | break;
|
|---|
| 250 | }
|
|---|
| 251 |
|
|---|
| 252 | return(nois);
|
|---|
| 253 | }
|
|---|
| 254 |
|
|---|
| 255 |
|
|---|
| 256 | /* Nouvelle-Fonction */
|
|---|
| 257 | template <class T>
|
|---|
| 258 | void ImgAddNoise(Image<T>& img, DynCCD const& dyn)
|
|---|
| 259 | {
|
|---|
| 260 | T* p = img.ImagePtr();
|
|---|
| 261 | int nPix = img.XSize() * img.YSize();
|
|---|
| 262 |
|
|---|
| 263 | for (int i=0; i<nPix; i++, p++)
|
|---|
| 264 | *p += (T) (dyn.Noise(*p)*NorRand());
|
|---|
| 265 | }
|
|---|
| 266 |
|
|---|
| 267 |
|
|---|
| 268 |
|
|---|
| 269 | /* Nouvelle-Fonction */
|
|---|
| 270 | RzImage * NoiseImage(RzImage const * pim, DynCCD const * dynccd)
|
|---|
| 271 |
|
|---|
| 272 | /* Creation et Calcul d'une image de bruit a partir de l'image */
|
|---|
| 273 | /* pim et de dynccd. Voir la methode DynCCD::Noise() pour la */
|
|---|
| 274 | /* description du calcul */
|
|---|
| 275 | {
|
|---|
| 276 | RzImage * nois = NULL;
|
|---|
| 277 |
|
|---|
| 278 | switch (pim->PixelType())
|
|---|
| 279 | {
|
|---|
| 280 | case kuint_2:
|
|---|
| 281 | {
|
|---|
| 282 | ImageU2 pix((RzImage&)(*pim));
|
|---|
| 283 | nois = NoiseImage(&pix, dynccd);
|
|---|
| 284 | break;
|
|---|
| 285 | }
|
|---|
| 286 | case kint_2:
|
|---|
| 287 | {
|
|---|
| 288 | ImageI2 pix((RzImage&)(*pim));
|
|---|
| 289 | nois = NoiseImage(&pix, dynccd);
|
|---|
| 290 | break;
|
|---|
| 291 | }
|
|---|
| 292 | case kint_4:
|
|---|
| 293 | {
|
|---|
| 294 | ImageI4 pix((RzImage&)(*pim));
|
|---|
| 295 | nois = NoiseImage(&pix, dynccd);
|
|---|
| 296 | break;
|
|---|
| 297 | }
|
|---|
| 298 | case kr_4:
|
|---|
| 299 | {
|
|---|
| 300 | ImageR4 pix((RzImage&)(*pim));
|
|---|
| 301 | nois = NoiseImage(&pix, dynccd);
|
|---|
| 302 | break;
|
|---|
| 303 | }
|
|---|
| 304 | case kuint_4:
|
|---|
| 305 | case kr_8:
|
|---|
| 306 | cerr << "NoiseImage(RzImage, ...) Unsupported image type (kuint_4/kr_8) "
|
|---|
| 307 | << int(pim->PixelType()) << "\n" ;
|
|---|
| 308 | break;
|
|---|
| 309 | default:
|
|---|
| 310 | cerr << "NoiseImage(RzImage, ...) Unknown image type !! " << int(pim->PixelType()) << "\n" ;
|
|---|
| 311 | break;
|
|---|
| 312 | }
|
|---|
| 313 |
|
|---|
| 314 | return(nois);
|
|---|
| 315 | }
|
|---|
| 316 |
|
|---|
| 317 |
|
|---|
| 318 | // ******** INSTANCES
|
|---|
| 319 |
|
|---|
| 320 | #if defined(__xlC) || defined(__aCC__)
|
|---|
| 321 | void instancetempdynccd(int n)
|
|---|
| 322 | {
|
|---|
| 323 | /* Cette fonction sert uniquement a forcer le compilo a instancier les
|
|---|
| 324 | classes/fonctions template */
|
|---|
| 325 |
|
|---|
| 326 | ImageU2 iu2(n,n), *piu2;
|
|---|
| 327 | ImageI2 ii2(n,n), *pii2;
|
|---|
| 328 | ImageI4 ii4(n,n), *pii4;
|
|---|
| 329 | ImageR4 ir4(n,n), *pir4;
|
|---|
| 330 | DynCCD dyn(1,0.,65000., 4., 20.);
|
|---|
| 331 |
|
|---|
| 332 | piu2 = NoiseImage(&iu2, &dyn);
|
|---|
| 333 | pii2 = NoiseImage(&ii2, &dyn);
|
|---|
| 334 | pii4 = NoiseImage(&ii4, &dyn);
|
|---|
| 335 | pir4 = NoiseImage(&ir4, &dyn);
|
|---|
| 336 |
|
|---|
| 337 | ImgAddNoise(*piu2, dyn);
|
|---|
| 338 | ImgAddNoise(*pii2, dyn);
|
|---|
| 339 | ImgAddNoise(*pii4, dyn);
|
|---|
| 340 | ImgAddNoise(*pir4, dyn);
|
|---|
| 341 |
|
|---|
| 342 | return;
|
|---|
| 343 | }
|
|---|
| 344 | #endif
|
|---|
| 345 |
|
|---|
| 346 |
|
|---|
| 347 | #ifdef __CXX_PRAGMA_TEMPLATES__
|
|---|
| 348 |
|
|---|
| 349 | #pragma define_template NoiseImage<uint_2>
|
|---|
| 350 | #pragma define_template NoiseImage<int_2>
|
|---|
| 351 | #pragma define_template NoiseImage<int_4>
|
|---|
| 352 | #pragma define_template NoiseImage<r_4>
|
|---|
| 353 |
|
|---|
| 354 | #pragma define_template ImgAddNoise<uint_2>
|
|---|
| 355 | #pragma define_template ImgAddNoise<int_2>
|
|---|
| 356 | #pragma define_template ImgAddNoise<int_4>
|
|---|
| 357 | #pragma define_template ImgAddNoise<r_4>
|
|---|
| 358 |
|
|---|
| 359 | #endif
|
|---|
| 360 |
|
|---|
| 361 |
|
|---|
| 362 | #if ( defined(ANSI_TEMPLATES) && !defined(__aCC__) )
|
|---|
| 363 | template Image<uint_2> * NoiseImage<uint_2>(Image<uint_2> const * , DynCCD const *);
|
|---|
| 364 | template Image< int_2> * NoiseImage< int_2>(Image< int_2> const * , DynCCD const *);
|
|---|
| 365 | template Image< int_4> * NoiseImage< int_4>(Image< int_4> const * , DynCCD const *);
|
|---|
| 366 | template Image< r_4> * NoiseImage< r_4>(Image< r_4> const * , DynCCD const *);
|
|---|
| 367 |
|
|---|
| 368 | template void ImgAddNoise<uint_2>(Image<uint_2>&, DynCCD const&);
|
|---|
| 369 | template void ImgAddNoise< int_2>(Image< int_2>&, DynCCD const&);
|
|---|
| 370 | template void ImgAddNoise< int_4>(Image< int_4>&, DynCCD const&);
|
|---|
| 371 | template void ImgAddNoise< r_4>(Image< r_4>&, DynCCD const&);
|
|---|
| 372 |
|
|---|
| 373 | #endif
|
|---|
| 374 |
|
|---|
| 375 | #if defined(__GNU_TEMPLATES__)
|
|---|
| 376 | template Image<uint_2> * NoiseImage(Image<uint_2> const *, DynCCD const *);
|
|---|
| 377 | template Image< int_2> * NoiseImage(Image< int_2> const *, DynCCD const *);
|
|---|
| 378 | template Image< int_4> * NoiseImage(Image< int_4> const *, DynCCD const *);
|
|---|
| 379 | template Image< r_4> * NoiseImage(Image< r_4> const *, DynCCD const *);
|
|---|
| 380 |
|
|---|
| 381 | template void ImgAddNoise(Image<uint_2>&, DynCCD const&);
|
|---|
| 382 | template void ImgAddNoise(Image< int_2>&, DynCCD const&);
|
|---|
| 383 | template void ImgAddNoise(Image< int_4>&, DynCCD const&);
|
|---|
| 384 | template void ImgAddNoise(Image< r_4>&, DynCCD const&);
|
|---|
| 385 |
|
|---|
| 386 | #endif
|
|---|
| 387 |
|
|---|
| 388 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
