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

Last change on this file since 852 was 270, checked in by ansari, 26 years ago

Portage sur Linux (instances template) Reza 27/04/99
File size: 10.2 KB
Line 
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-- */
64DynCCD::DynCCD(int typ, float min, float max, float g,
65 float ron, float rf, float rfs)
66{
67if ( (typ >= kConstantNoise) && (typ <= kSqrtADUNoise) ) TypNoise = typ;
68else TypNoise = kConstantNoise;
69MinADU = min; MaxADU = max;
70Gain = g; RONoise = ron;
71RefFond = rf; RefSFond = rfs;
72}
73
74/* --Methode-- */
75void DynCCD::Set(int typ, float min, float max, float g,
76 float ron, float rf, float rfs)
77{
78if ( (typ >= kConstantNoise) && (typ <= kSqrtADUNoise) ) TypNoise = typ;
79MinADU = min; MaxADU = max;
80Gain = g; RONoise = ron;
81RefFond = rf; RefSFond = rfs;
82}
83
84/* --Methode-- */
85void DynCCD::Print()
86{
87printf("DynCCD: Type= %d Min/MaxADU= %g %g Gain/RoN= %g %g\n",
88 TypNoise, MinADU, MaxADU, Gain, RONoise);
89if (TypNoise == 2)
90 printf("... RefFond= %g RefSFond= %g \n", RefFond, RefSFond);
91return;
92}
93
94/* --Methode-- */
95float 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{
110float h,s,ronsq;
111float fond;
112
113if ( (pixel > MaxADU) || (pixel < MinADU) ) return(0.);
114if ( TypNoise == kConstantNoise) return(1.);
115if ( TypNoise == kSqrtADUNoise ) return(fsqrt(fabsf(pixel)));
116
117if ( TypNoise == kSigFondNoise)
118 { fond = RefFond;
119 ronsq = RefSFond * RefSFond; }
120else
121 { fond = 0;
122 ronsq = RONoise/Gain; ronsq *= ronsq; }
123
124h = (pixel>fond) ? (float)(pixel-fond) : 0.;
125s = ronsq+h/Gain;
126s = fsqrt(s);
127return(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 */
162template <class T>
163Image<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{
170float h,s,ronsq;
171float fond, min,max;
172int i, npix;
173float minois, offnois;
174
175if (pim == NULL) return(NULL);
176
177const T * pix = pim->ImagePtr();
178npix = pim->XSize()*pim->YSize();
179Image<T> * nois = new Image<T>(pim->XSize(), pim->YSize(), false);
180nois->SetOrg(pim->XOrg(), pim->YOrg());
181T * pno = nois->ImagePtr();
182
183min = dynccd->MinADU; max = dynccd->MaxADU;
184
185
186switch (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
252return(nois);
253}
254
255
256/* Nouvelle-Fonction */
257template <class T>
258void 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 */
270RzImage * 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{
276RzImage * nois = NULL;
277
278switch (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
314return(nois);
315}
316
317
318// ******** INSTANCES
319
320#if defined(__xlC) || defined(__aCC__)
321void instancetempdynccd(int n)
322{
323/* Cette fonction sert uniquement a forcer le compilo a instancier les
324 classes/fonctions template */
325
326ImageU2 iu2(n,n), *piu2;
327ImageI2 ii2(n,n), *pii2;
328ImageI4 ii4(n,n), *pii4;
329ImageR4 ir4(n,n), *pir4;
330DynCCD dyn(1,0.,65000., 4., 20.);
331
332piu2 = NoiseImage(&iu2, &dyn);
333pii2 = NoiseImage(&ii2, &dyn);
334pii4 = NoiseImage(&ii4, &dyn);
335pir4 = NoiseImage(&ir4, &dyn);
336
337ImgAddNoise(*piu2, dyn);
338ImgAddNoise(*pii2, dyn);
339ImgAddNoise(*pii4, dyn);
340ImgAddNoise(*pir4, dyn);
341
342return;
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__) )
363template Image<uint_2> * NoiseImage<uint_2>(Image<uint_2> const * , DynCCD const *);
364template Image< int_2> * NoiseImage< int_2>(Image< int_2> const * , DynCCD const *);
365template Image< int_4> * NoiseImage< int_4>(Image< int_4> const * , DynCCD const *);
366template Image< r_4> * NoiseImage< r_4>(Image< r_4> const * , DynCCD const *);
367
368template void ImgAddNoise<uint_2>(Image<uint_2>&, DynCCD const&);
369template void ImgAddNoise< int_2>(Image< int_2>&, DynCCD const&);
370template void ImgAddNoise< int_4>(Image< int_4>&, DynCCD const&);
371template void ImgAddNoise< r_4>(Image< r_4>&, DynCCD const&);
372
373#endif
374
375#if defined(__GNU_TEMPLATES__)
376template Image<uint_2> * NoiseImage(Image<uint_2> const *, DynCCD const *);
377template Image< int_2> * NoiseImage(Image< int_2> const *, DynCCD const *);
378template Image< int_4> * NoiseImage(Image< int_4> const *, DynCCD const *);
379template Image< r_4> * NoiseImage(Image< r_4> const *, DynCCD const *);
380
381template void ImgAddNoise(Image<uint_2>&, DynCCD const&);
382template void ImgAddNoise(Image< int_2>&, DynCCD const&);
383template void ImgAddNoise(Image< int_4>&, DynCCD const&);
384template void ImgAddNoise(Image< r_4>&, DynCCD const&);
385
386#endif
387
388
Note: See TracBrowser for help on using the repository browser.