Line data Source code
1 : /** \file ccdDetector.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Provides a class to simulate a CCD.
4 : * \ingroup mxAO_files
5 : *
6 : */
7 :
8 : #ifndef ccdDetector_hpp
9 : #define ccdDetector_hpp
10 :
11 : // #include <Eigen/Dense>
12 :
13 : #include "../../improc/imageTransforms.hpp"
14 : #include "../../math/randomT.hpp"
15 :
16 : #include "wavefront.hpp"
17 :
18 : #include <fcntl.h>
19 : #include <iostream>
20 :
21 : #ifndef BREAD_CRUMB
22 : #define MXLIB_CCD_DETECTOR_LOCAL_BREAD_CRUMB
23 : #ifdef DEBUG
24 : #define BREAD_CRUMB std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
25 : #else
26 : #define BREAD_CRUMB
27 : #endif
28 : #endif
29 :
30 : namespace mx
31 : {
32 : namespace AO
33 : {
34 : namespace sim
35 : {
36 :
37 : /// A simulated CCD detector
38 : /** A simulated CCD detector, including an optional EMCCD noise model.
39 : *
40 : * \ingroup mxAOSim
41 : */
42 : template <typename _realT>
43 : class ccdDetector
44 : {
45 :
46 : public:
47 : typedef _realT realT;
48 :
49 : typedef wavefront<realT> wavefrontT;
50 :
51 : //typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imageT;
52 :
53 : typedef mx::math::randomT<realT, std::mt19937_64, std::normal_distribution<realT>> norm_distT;
54 :
55 : typedef mx::math::randomT<int, std::mt19937_64, std::poisson_distribution<int>> poisson_distT;
56 :
57 : typedef mx::math::randomT<realT, std::mt19937_64, std::gamma_distribution<realT>> gamma_distT;
58 :
59 : ccdDetector();
60 :
61 : protected:
62 : norm_distT m_normVar; ///< Gets normal-distributed variates
63 : poisson_distT m_poissonVar; ///< Gets Poisson distributed variates
64 : gamma_distT m_gammaVar; ///< Gets gamma distributed variates
65 :
66 : realT m_qe{ 1 }; ///< The quantum efficiency.
67 :
68 : realT m_darkCurrent{ 0 }; ///< The dark current, per pixel per second
69 : realT m_ron{ 0 }; ///< The readout noise, electrons per pixel per read.
70 :
71 : realT m_cic{ 0 }; ///< EMCCD clock induced charge, electrons per pixel per read.
72 : realT m_gain{ 1 }; ///< Electron multiplication gain. If >1, then EMCCD is modeled.
73 :
74 : realT m_expTime{ 1 }; ///< The exposure time, in seconds.
75 :
76 : int m_rows{ 0 }; ///< The detector size, in rows.
77 : int m_cols{ 0 }; ///< The detector size, in columns.
78 :
79 : bool m_noNoise{ false }; ///< If true no noise is added to the exposed image.
80 :
81 : // This is probably vestigial, commented out on 2018-07-23
82 : // Delete after testing with an aoSim compile.
83 : // long get_seed();
84 :
85 : public:
86 : /// Get the current value of qe.
87 : /**
88 : * \returns the current value of m_qe
89 : */
90 : realT qe();
91 :
92 : void qe( const realT &q );
93 :
94 : realT darkCurrent();
95 :
96 : void darkCurrent( const realT &dc );
97 :
98 : realT ron();
99 :
100 : void ron( const realT &r );
101 :
102 : realT cic();
103 :
104 : void cic( const realT &c );
105 :
106 : realT gain();
107 :
108 : void gain( const realT &g );
109 :
110 : realT expTime();
111 :
112 : void expTime( const realT &dt );
113 :
114 : int rows();
115 :
116 : void rows( const int &r );
117 :
118 : int cols();
119 :
120 : void cols( const int &c );
121 :
122 : void setSize( const int &r, const int &c );
123 :
124 : bool noNoise();
125 :
126 : void noNoise( bool nn );
127 :
128 : /// Rebin and add noise to the input image, placing the result in the output image.
129 : /** The output image must be the same size or smaller than the input image.
130 : * The output image is resized only if necessary.
131 : * The input image is multiplied by expTime, so its flux should be in photons/sec.
132 : * Noise is modeled as follows:
133 : * -# The number of input photons for pixel (i,j) is \f$ n_{ie} = (in(i,j) + darkCurrent)*expTime + cic \f$
134 : * -# The Poisson distribution with mean and variance \f$ n_{ie} \f$ is sampled.
135 : * -# If gain > 1, the gamma distribution is sampled
136 : * -# Read noise is applied
137 : * -# Result is rounded to nearest whole integer.
138 : */
139 : template<typename imageTout, typename imageTin>
140 : void exposeImage( imageTout &out, ///< [out] The output image, after all above steps applied.
141 : imageTin &in ///< [in] the input image, in photons/sec flux units.
142 : );
143 : };
144 :
145 : template <typename realT>
146 1 : ccdDetector<realT>::ccdDetector()
147 : {
148 1 : m_normVar.seed();
149 1 : m_poissonVar.seed();
150 1 : m_gammaVar.seed();
151 1 : }
152 :
153 : /*template<typename realT>
154 : long ccdDetector<realT>::get_seed()
155 : {
156 : long int seedval;
157 :
158 : int fd;
159 :
160 : fd = open("/dev/urandom", O_RDONLY);
161 :
162 : int rv =read(fd, &seedval, sizeof(long int));
163 :
164 : close(fd);
165 :
166 :
167 : if(rv < 0)
168 : {
169 : std::cerr << "read from /dev/urandom returned error\n";
170 :
171 : return 0;
172 : }
173 :
174 : return seedval;
175 : }
176 : */
177 :
178 : template <typename realT>
179 0 : realT ccdDetector<realT>::qe()
180 : {
181 0 : return m_qe;
182 : }
183 :
184 : template <typename realT>
185 0 : void ccdDetector<realT>::qe( const realT &q )
186 : {
187 0 : m_qe = q;
188 0 : }
189 :
190 : template <typename realT>
191 0 : realT ccdDetector<realT>::darkCurrent()
192 : {
193 0 : return m_darkCurrent;
194 : }
195 :
196 : template <typename realT>
197 0 : void ccdDetector<realT>::darkCurrent( const realT &dc )
198 : {
199 0 : m_darkCurrent = dc;
200 0 : }
201 :
202 : template <typename realT>
203 0 : realT ccdDetector<realT>::ron()
204 : {
205 0 : return m_ron;
206 : }
207 :
208 : template <typename realT>
209 0 : void ccdDetector<realT>::ron( const realT &r )
210 : {
211 0 : m_ron = r;
212 0 : }
213 :
214 : template <typename realT>
215 0 : realT ccdDetector<realT>::cic()
216 : {
217 0 : return m_cic;
218 : }
219 :
220 : template <typename realT>
221 0 : void ccdDetector<realT>::cic( const realT &c )
222 : {
223 0 : m_cic = c;
224 0 : }
225 :
226 : template <typename realT>
227 0 : realT ccdDetector<realT>::gain()
228 : {
229 0 : return m_gain;
230 : }
231 :
232 : template <typename realT>
233 0 : void ccdDetector<realT>::gain( const realT &g )
234 : {
235 0 : m_gain = g;
236 0 : }
237 :
238 : template <typename realT>
239 0 : realT ccdDetector<realT>::expTime()
240 : {
241 0 : return m_expTime;
242 : }
243 :
244 : template <typename realT>
245 1 : void ccdDetector<realT>::expTime( const realT &dt )
246 : {
247 1 : m_expTime = dt;
248 1 : }
249 :
250 : template <typename realT>
251 0 : int ccdDetector<realT>::rows()
252 : {
253 0 : return m_rows;
254 : }
255 :
256 : template <typename realT>
257 0 : void ccdDetector<realT>::rows( const int &r )
258 : {
259 0 : m_rows = r;
260 0 : }
261 :
262 : template <typename realT>
263 0 : int ccdDetector<realT>::cols()
264 : {
265 0 : return m_cols;
266 : }
267 :
268 : template <typename realT>
269 0 : void ccdDetector<realT>::cols( const int &c )
270 : {
271 0 : m_cols = c;
272 0 : }
273 :
274 : template <typename realT>
275 3 : void ccdDetector<realT>::setSize( const int &r, const int &c )
276 : {
277 3 : m_rows = r;
278 3 : m_cols = c;
279 3 : }
280 :
281 : template <typename realT>
282 0 : bool ccdDetector<realT>::noNoise()
283 : {
284 0 : return m_noNoise;
285 : }
286 :
287 : template <typename realT>
288 0 : void ccdDetector<realT>::noNoise( bool nn )
289 : {
290 0 : m_noNoise = nn;
291 0 : }
292 :
293 : template <typename realT>
294 : template<typename imageTout, typename imageTin>
295 : void ccdDetector<realT>::exposeImage( imageTout &out, imageTin &in )
296 : {
297 :
298 : using poisson_param_t = typename std::poisson_distribution<int>::param_type;
299 : using gamma_param_t = typename std::gamma_distribution<realT>::param_type;
300 :
301 : BREAD_CRUMB
302 :
303 : out.resize( m_rows, m_cols );
304 :
305 : BREAD_CRUMB
306 :
307 : improc::imageDownSample( out, in );
308 :
309 : BREAD_CRUMB
310 :
311 : if( m_noNoise )
312 : return;
313 :
314 : for( int j = 0; j < m_cols; ++j )
315 : {
316 : for( int i = 0; i < m_rows; ++i )
317 : {
318 : realT charge = ( out( i, j ) * m_qe + m_darkCurrent ) * m_expTime + m_cic;
319 :
320 : if( charge > 1000 )
321 : {
322 : out( i, j ) = charge + m_normVar * sqrt( charge );
323 : }
324 : else
325 : {
326 : m_poissonVar.distribution.param( poisson_param_t{ charge } );
327 : out( i, j ) = m_poissonVar;
328 : }
329 :
330 : if( m_gain > 1 )
331 : {
332 : m_gammaVar.distribution.param( gamma_param_t{ out( i, j ), m_gain } );
333 :
334 : out( i, j ) = m_gammaVar;
335 : }
336 :
337 : if( m_ron > 0 )
338 : {
339 : out( i, j ) += m_normVar * m_ron;
340 : }
341 :
342 : if( m_gain > 1 )
343 : {
344 : out( i, j ) /= m_gain;
345 : }
346 :
347 : // out(i,j) = round(out(i,j));
348 : }
349 : }
350 : }
351 :
352 : } // namespace sim
353 : } // namespace AO
354 : } // namespace mx
355 :
356 : #ifdef MXLIB_CCD_DETECTOR_LOCAL_BREAD_CRUMB
357 : #undef BREAD_CRUMB
358 : #undef MXLIB_CCD_DETECTOR_LOCAL_BREAD_CRUMB
359 : #endif
360 :
361 : #endif // ccdDetector_hpp
|
