Line data Source code
1 : /** \file directPhaseSensor.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Declaration and definition of a direct phase sensor.
4 : * \ingroup mxAO_sim_files
5 : *
6 : */
7 :
8 : #ifndef directPhaseSensor_hpp
9 : #define directPhaseSensor_hpp
10 :
11 : #include "../../math/randomT.hpp"
12 :
13 : #include "../../ioutils/fits/fitsFile.hpp"
14 : #include "../../improc/eigenCube.hpp"
15 :
16 : // #include <mx/fraunhoferImager.hpp>
17 : #include "../../sigproc/psdFilter.hpp"
18 :
19 : #include "wavefront.hpp"
20 :
21 : #ifdef DEBUG
22 : #define BREAD_CRUMB std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
23 : #else
24 : #define BREAD_CRUMB
25 : #endif
26 :
27 : namespace mx
28 : {
29 : namespace AO
30 : {
31 : namespace sim
32 : {
33 :
34 : template <typename _realT>
35 : struct wfsImageT
36 : {
37 : typedef _realT realT;
38 :
39 : realT iterNo;
40 :
41 : /// The wavefront sensor detector image type
42 : typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imageT;
43 :
44 : imageT image;
45 : };
46 :
47 : template <typename _realT, typename _detectorT>
48 : class directPhaseSensor
49 : {
50 : public:
51 : typedef _realT realT;
52 :
53 : typedef std::complex<_realT> complexT;
54 :
55 : /// The wavefront data type
56 : typedef wavefront<_realT> wavefrontT;
57 :
58 : /// The pupil type
59 : typedef improc::eigenImage<realT> pupilT;
60 :
61 : /// The wavefront complex field type
62 : typedef mx::wfp::imagingArray<std::complex<_realT>, mx::wfp::fftwAllocator<std::complex<_realT>>, 0> complexFieldT;
63 :
64 : typedef _detectorT detectorT;
65 :
66 : typedef mx::math::randomT<realT, std::mt19937_64, std::normal_distribution<realT>> norm_distT;
67 :
68 : /// Alias for a poisson random variate
69 : typedef mx::math::randomT<long, std::mt19937_64, std::poisson_distribution<long>> poisson_distT;
70 :
71 : protected:
72 : /* Standard WFS Interface: */
73 : int m_wfSz{ 0 }; ///< Size of the wavefront in pixels
74 :
75 : int m_detRows{ 0 }; ///< The number of rows of the WFS m_detector. After forming the image the WFS detector plane
76 : ///< is re-binned to this.
77 : int m_detCols{ 0 }; ///< The number of columns of the WFS m_detector. After forming the image the WFS detector
78 : ///< plane is re-binned to this.
79 :
80 : realT m_lambda{ 0 }; ///< Central wavelength, in meters
81 :
82 : int m_iTime{ 0 }; ///< Integration time in loop steps
83 :
84 : int m_roTime{ 1 }; ///< Readout time in loop steps
85 :
86 : realT m_simStep{ 0.001 }; ///< The simulation stepsize in seconds.
87 :
88 : /* Direct Phase Sensor specific: */
89 :
90 : realT m_npix{ 0 }; ///< The number of pixels assumed actually used for WFS in the S/N calc.
91 : realT m_Fbg{ 0.0 }; ///< The background flux in photons/sec/pixel.
92 :
93 : int m_iTime_counter{ 0 };
94 :
95 : int m_reading{ 0 };
96 :
97 : int m_roTime_counter{ 0 };
98 :
99 : bool m_firstRun{ true };
100 :
101 : std::vector<wavefrontT> m_wavefronts;
102 :
103 : int m_lastWavefront;
104 :
105 : /// The image formed by the WFS
106 : wfsImageT<realT> m_wfsImage;
107 :
108 : public:
109 : bool m_poissonNoise{ true };
110 :
111 : /// Default c'tor
112 : directPhaseSensor();
113 :
114 : /* The standard WFS interface: */
115 :
116 : detectorT m_detector;
117 :
118 : /// The image on the detector, resized from wfsImage
119 : wfsImageT<realT> m_detectorImage;
120 :
121 : /// Get the wavefront size in pixels
122 : /**
123 : * \returns the wavefront size in pixels
124 : */
125 : int wfSz();
126 :
127 : /// Set the wavefront size in pixels.
128 : void wfSz( int sz /**< [in] The new size */ );
129 :
130 : /// Get the detector rows in pixels
131 : /**
132 : * \returns m_detRows
133 : */
134 : int detRows();
135 :
136 : /// Set the detector rows in pixels while preserving the current column count.
137 : void detRows( int sz /**< [in] The new row count. */ );
138 :
139 : /// Get the detector Cols in pixels
140 : /**
141 : * \returns m_detCols
142 : */
143 : int detCols();
144 :
145 : /// Set the detector columns in pixels while preserving the current row count.
146 : void detCols( int sz /**< [in] The new column count. */ );
147 :
148 : /// Set both detector dimensions in pixels.
149 : void detSize( int nrows, /**< [in] The new row count. */
150 : int ncols /**< [in] The new column count. */ );
151 :
152 : /// Get the PyWFS central wavelength
153 : /**
154 : * \returns the central wavelength in meters
155 : */
156 : realT lambda();
157 :
158 : /// Set the PyWFS central wavelength
159 : /**
160 : * \param d is the new central wavelength in meters
161 : */
162 : void lambda( realT l );
163 :
164 : /// Get the PyWFS integration time, in time steps
165 : int iTime();
166 :
167 : /// Set the PyWFS integration time, in time steps
168 : void iTime( int it );
169 :
170 : /// Get the PyWFS detector readout time, in time steps
171 : int roTime();
172 :
173 : /// Set the PyWFS detector readout time, in time steps
174 : void roTime( int rt );
175 :
176 : /// Get the simulation step-size, in seconds.
177 : realT simStep();
178 :
179 : /// Set the simulation step-size, in seconds.
180 : void simStep( realT st );
181 :
182 : /// Get the number of pixels used for noise calculations
183 : /**
184 : * \returns the number of pixels used by the wfs, m_npix.
185 : */
186 : realT npix();
187 :
188 : /// Set the number of pixels used fo rnoise calculations
189 : /**
190 : */
191 : void npix( realT np /**< [in] the number of pixels to be used by the WFS*/ );
192 :
193 : /// Get the background rate
194 : /**
195 : * \returns the background rate in photons/pixel/sec, the current value of m_Fbg.
196 : */
197 : realT Fbg();
198 :
199 : /// Set the background rate
200 : /**
201 : */
202 : void Fbg( realT bg /**< [in] the background rate in photons/pixel/sec.*/ );
203 :
204 : template <typename AOSysT>
205 : void linkSystem( AOSysT &AOSys );
206 :
207 : /// Record a wavefront without sensing
208 : /** This merely inserts the wavefront in the circular buffer but does no other processing
209 : * Retruns 0 on success
210 : */
211 : int recordWavefront( wavefrontT &pupilPlane );
212 :
213 : /// Sense the wavefront aberrations
214 : /** Returns true if a new wavefront measurement is ready.
215 : * Retruns false if still integrating.
216 : */
217 : bool senseWavefront( wavefrontT &pupilPlane );
218 :
219 : /// Sense the wavefront aberrations in a calibration mode
220 : /** Allows for faster calibrations.
221 : */
222 : bool senseWavefrontCal( wavefrontT &pupilPlane );
223 :
224 : void doSenseWavefront();
225 :
226 : /** \name Spatial Filtering
227 : * @{
228 : */
229 : protected:
230 : bool m_applyFilter{ false }; ///< Flag to control whether or not the spatial filter is applied.
231 :
232 : int m_filterWidth{
233 : 0 }; ///< The half-width of the filter, in Fourier domain pixels corresponding to 1/D spatial frequency units.
234 :
235 : sigproc::psdFilter<realT, 2> m_filter; ///< The spatial filter class object.
236 :
237 : public:
238 : /// Set the flag controlling whether or not the spatial filter is applied
239 : void applyFilter( bool af /**< [in] true to apply the filter, false to not apply*/ );
240 :
241 : /// Get the current value of the flag controlling whether or not the spatial filter is applied
242 : /**
243 : * \returns the current value of m_applyFilter
244 : */
245 : bool applyFilter();
246 :
247 : /// Set the spatial filter half-width
248 : /** This sets up the filter and loads it into m_filter. The filter is a hard-edged square of width 2*width+1 pixels
249 : * in the Fourier plane.
250 : *
251 : * m_filterWidth must still be set to true after calling this function.
252 : */
253 : void filterWidth( int width /**< [in] the desired half-width of the filter*/ );
254 :
255 : /// Get the current value of the filter half-width
256 : /**
257 : * \returns the current value of m_filterWidth
258 : */
259 : int filterWidth();
260 :
261 : /// Provides access to the filter for use by other parts of the simulation.
262 : /**
263 : * \returns a const reference to m_filter.
264 : */
265 : const sigproc::psdFilter<realT, 2> &filter();
266 :
267 : ///@}
268 :
269 : /** \name Photon Noise
270 : * To turn on the addition of photon noise to the WFS image, you must set beta_p to be greater than 0,
271 : * and must supply the pupil image.
272 : *
273 : * @{
274 : */
275 : protected:
276 : norm_distT m_normVar; ///< Gets normal-distributed variates
277 : poisson_distT m_poissonVar;
278 :
279 : /// The photon noise senstivity parameter.
280 : /** If 0 (default) no noise is applied. A value of 1 represents the ideal interferomter.
281 : * Note that this is constant vs. spatial frequency. See Guyon (2005) \cite guyon_2005 for more information about
282 : * \f$ \beta_p \f$.
283 : */
284 : realT m_beta_p{ 0 };
285 :
286 : /// The pupil is needed to properly normalize poisson noise.
287 : /** If null, then no noise will be applied.
288 : */
289 : pupilT *m_pupil{ nullptr };
290 :
291 : /// Array used internally to calculated noise, global to avoid re-allocations.
292 : typename wfsImageT<realT>::imageT m_noiseIm;
293 :
294 : public:
295 : /// Set the new value of the photon noise sensitivity parameter, beta_p.
296 : void beta_p( realT bp /**< [in] the new value of beta_p */ );
297 :
298 : /// Get the current value of the photon noise sensitivity parameter.
299 : /**
300 : * \returns the current value of m_beta_p.
301 : */
302 : realT beta_p();
303 :
304 : /// Set the pupil
305 : /** Sets the m_pupil pointer to the value provided.
306 : */
307 : void pupil( pupilT *pupil /**< [in] pointer to the pupil image */ );
308 :
309 : /// Set the pupil
310 : /** Sets the m_pupil pointer to the address of the image provided.
311 : */
312 : void pupil( pupilT &pupil /**< [in] the pupil image */ );
313 :
314 : /// Get the current pupil as a pointer
315 : /**
316 : * \returns the current value of m_pupil
317 : */
318 : pupilT *pupil();
319 :
320 : ///@}
321 : };
322 :
323 : template <typename _realT, typename _detectorT>
324 1 : directPhaseSensor<_realT, _detectorT>::directPhaseSensor()
325 : {
326 1 : iTime( 1 );
327 :
328 1 : m_normVar.seed();
329 1 : m_poissonVar.seed();
330 1 : }
331 :
332 : template <typename _realT, typename _detectorT>
333 0 : int directPhaseSensor<_realT, _detectorT>::wfSz()
334 : {
335 0 : return m_wfSz;
336 : }
337 :
338 : template <typename _realT, typename _detectorT>
339 0 : void directPhaseSensor<_realT, _detectorT>::wfSz( int sz )
340 : {
341 0 : if( m_wfSz == sz )
342 0 : return;
343 :
344 0 : m_wfSz = sz;
345 : }
346 :
347 : template <typename _realT, typename _detectorT>
348 3 : int directPhaseSensor<_realT, _detectorT>::detRows()
349 : {
350 3 : return m_detRows;
351 : }
352 :
353 : template <typename _realT, typename _detectorT>
354 1 : void directPhaseSensor<_realT, _detectorT>::detRows( int sz )
355 : {
356 1 : detSize( sz, m_detCols );
357 1 : }
358 :
359 : template <typename _realT, typename _detectorT>
360 3 : int directPhaseSensor<_realT, _detectorT>::detCols()
361 : {
362 3 : return m_detCols;
363 : }
364 :
365 : template <typename _realT, typename _detectorT>
366 1 : void directPhaseSensor<_realT, _detectorT>::detCols( int sz )
367 : {
368 1 : detSize( m_detRows, sz );
369 1 : }
370 :
371 : template <typename _realT, typename _detectorT>
372 3 : void directPhaseSensor<_realT, _detectorT>::detSize( int nrows, int ncols )
373 : {
374 3 : if( m_detRows == nrows && m_detCols == ncols )
375 0 : return;
376 :
377 3 : m_detRows = nrows;
378 3 : m_detCols = ncols;
379 :
380 3 : m_detector.setSize( m_detRows, m_detCols );
381 3 : m_detectorImage.image.resize( m_detRows, m_detCols );
382 : }
383 :
384 : template <typename _realT, typename _detectorT>
385 0 : _realT directPhaseSensor<_realT, _detectorT>::lambda()
386 : {
387 0 : return m_lambda;
388 : }
389 :
390 : template <typename _realT, typename _detectorT>
391 0 : void directPhaseSensor<_realT, _detectorT>::lambda( _realT l )
392 : {
393 0 : m_lambda = l;
394 0 : }
395 :
396 : template <typename _realT, typename _detectorT>
397 : template <typename AOSysT>
398 : void directPhaseSensor<_realT, _detectorT>::linkSystem( AOSysT &AOSys )
399 : {
400 : static_cast<void>( AOSys );
401 : }
402 :
403 : template <typename _realT, typename _detectorT>
404 0 : int directPhaseSensor<_realT, _detectorT>::iTime()
405 : {
406 0 : return m_iTime;
407 : }
408 :
409 : template <typename _realT, typename _detectorT>
410 1 : void directPhaseSensor<_realT, _detectorT>::iTime( int it )
411 : {
412 1 : if( it < 1 )
413 : {
414 0 : std::cerr << "iTime must be >= 1. Correcting\n";
415 0 : it = 1;
416 : }
417 :
418 1 : m_iTime = it;
419 :
420 1 : m_wavefronts.resize( m_iTime + 2 + 100 );
421 1 : m_lastWavefront = -1;
422 :
423 1 : m_detector.expTime( m_simStep * m_iTime );
424 1 : }
425 :
426 : template <typename _realT, typename _detectorT>
427 0 : int directPhaseSensor<_realT, _detectorT>::roTime()
428 : {
429 0 : return m_roTime;
430 : }
431 :
432 : template <typename _realT, typename _detectorT>
433 0 : void directPhaseSensor<_realT, _detectorT>::roTime( int rt )
434 : {
435 0 : if( rt < 1 )
436 : {
437 0 : std::cerr << "roTime must be >= 1. Correcting\n";
438 0 : rt = 1;
439 : }
440 :
441 0 : m_roTime = rt;
442 0 : }
443 :
444 : template <typename _realT, typename _detectorT>
445 0 : _realT directPhaseSensor<_realT, _detectorT>::simStep()
446 : {
447 0 : return m_simStep;
448 : }
449 :
450 : template <typename _realT, typename _detectorT>
451 0 : void directPhaseSensor<_realT, _detectorT>::simStep( _realT st )
452 : {
453 :
454 0 : m_simStep = st;
455 :
456 0 : m_detector.expTime( m_simStep * m_iTime );
457 0 : }
458 :
459 : template <typename realT, typename _detectorT>
460 0 : realT directPhaseSensor<realT, _detectorT>::npix()
461 : {
462 0 : return m_npix;
463 : }
464 :
465 : template <typename realT, typename _detectorT>
466 0 : void directPhaseSensor<realT, _detectorT>::npix( realT np )
467 : {
468 0 : m_npix = np;
469 0 : }
470 :
471 : template <typename realT, typename _detectorT>
472 0 : realT directPhaseSensor<realT, _detectorT>::Fbg()
473 : {
474 0 : return m_Fbg;
475 : }
476 :
477 : template <typename realT, typename _detectorT>
478 0 : void directPhaseSensor<realT, _detectorT>::Fbg( realT bg )
479 : {
480 0 : m_Fbg = bg;
481 0 : }
482 :
483 : template <typename _realT, typename _detectorT>
484 0 : int directPhaseSensor<_realT, _detectorT>::recordWavefront( wavefrontT &pupilPlane )
485 : {
486 :
487 0 : ++m_lastWavefront;
488 0 : if( (size_t)m_lastWavefront >= m_wavefronts.size() )
489 0 : m_lastWavefront = 0;
490 :
491 0 : wavefrontT pPlane = pupilPlane;
492 :
493 0 : m_wavefronts[m_lastWavefront].amplitude = pPlane.amplitude;
494 0 : m_wavefronts[m_lastWavefront].phase = pPlane.phase;
495 0 : m_wavefronts[m_lastWavefront].iterNo = pPlane.iterNo;
496 :
497 0 : return 0;
498 0 : }
499 :
500 : template <typename _realT, typename _detectorT>
501 0 : bool directPhaseSensor<_realT, _detectorT>::senseWavefront( wavefrontT &pupilPlane )
502 : {
503 : using poisson_param_t = typename std::poisson_distribution<long>::param_type;
504 :
505 0 : ++m_lastWavefront;
506 0 : if( (size_t)m_lastWavefront >= m_wavefronts.size() )
507 0 : m_lastWavefront = 0;
508 :
509 0 : wavefrontT pPlane = pupilPlane;
510 :
511 0 : m_wavefronts[m_lastWavefront].amplitude = pPlane.amplitude;
512 0 : m_wavefronts[m_lastWavefront].phase = pPlane.phase;
513 0 : m_wavefronts[m_lastWavefront].iterNo = pPlane.iterNo;
514 :
515 : // std::cerr << m_lastWavefront << " " << pPlane.iterNo << "\n";
516 : // Always skip the first one for averaging to center of iTime.
517 0 : if( m_firstRun )
518 : {
519 0 : m_firstRun = false;
520 0 : return false;
521 : }
522 :
523 0 : ++m_iTime_counter;
524 :
525 0 : bool rv = false;
526 :
527 : if( 0 ) // m_reading)
528 : {
529 : ++m_roTime_counter;
530 :
531 : if( m_roTime_counter >= m_roTime )
532 : {
533 : m_detectorImage.image = m_wfsImage.image.block(
534 : 0.5 * ( m_wfsImage.image.rows() - 1 ) - 0.5 * ( m_detectorImage.image.rows() - 1 ),
535 : 0.5 * ( m_wfsImage.image.cols() - 1 ) - 0.5 * ( m_detectorImage.image.cols() - 1 ),
536 : m_detectorImage.image.rows(),
537 : m_detectorImage.image.cols() );
538 :
539 : m_detectorImage.iterNo = m_wfsImage.iterNo;
540 :
541 : m_roTime_counter = 0;
542 : m_reading = 0;
543 : rv = true;
544 : }
545 : }
546 :
547 0 : if( m_iTime_counter >= m_iTime )
548 : {
549 0 : doSenseWavefront();
550 :
551 0 : m_iTime_counter = 0;
552 :
553 0 : m_reading = 1;
554 0 : m_roTime_counter = 0;
555 :
556 : // Just do the read
557 0 : m_detectorImage.image =
558 0 : m_wfsImage.image.block( 0.5 * ( m_wfsImage.image.rows() - 1 ) - 0.5 * ( m_detectorImage.image.rows() - 1 ),
559 0 : 0.5 * ( m_wfsImage.image.cols() - 1 ) - 0.5 * ( m_detectorImage.image.cols() - 1 ),
560 : m_detectorImage.image.rows(),
561 : m_detectorImage.image.cols() );
562 :
563 : //*** Spatial Filter:
564 0 : if( m_applyFilter )
565 : {
566 : // std::cerr << "Filtering . . . \n";
567 0 : m_filter.filter( m_detectorImage.image );
568 0 : if( m_pupil != nullptr )
569 0 : m_detectorImage.image *= *m_pupil;
570 : }
571 :
572 : //*** Adding Noise:
573 0 : if( m_beta_p > 0 && m_pupil != nullptr )
574 : {
575 0 : m_noiseIm.resize( m_detectorImage.image.rows(), m_detectorImage.image.cols() );
576 :
577 : // Calculate intensity at each pixel
578 0 : for( int c = 0; c < m_noiseIm.cols(); ++c )
579 : {
580 0 : for( int r = 0; r < m_noiseIm.rows(); ++r )
581 : {
582 0 : m_noiseIm( r, c ) =
583 0 : pow( pupilPlane.amplitude( r, c ), 2 ) * ( *m_pupil )( r, c ) * m_detector.expTime();
584 : }
585 : }
586 : // std::cerr << "Total Phots: " << m_noiseIm.sum() << "\n";
587 :
588 0 : realT sqrtFbg = sqrt( m_Fbg * m_detector.expTime() );
589 :
590 : // Add noise
591 0 : for( int c = 0; c < m_noiseIm.cols(); ++c )
592 : {
593 0 : for( int r = 0; r < m_noiseIm.rows(); ++r )
594 : {
595 0 : if( m_noiseIm( r, c ) == 0 ) // assuming branch prediction makes this worth it
596 : {
597 0 : continue;
598 : }
599 : else
600 : {
601 0 : realT err = sqrt( m_noiseIm( r, c ) ) * m_normVar;
602 0 : if( sqrtFbg > 0 )
603 0 : err += sqrtFbg * m_normVar;
604 0 : if( m_detector.ron() > 0 )
605 0 : err += m_detector.ron() * m_normVar;
606 :
607 0 : m_detectorImage.image( r, c ) +=
608 0 : m_beta_p * err / m_noiseIm( r, c ); // Fractional noise in photons is noise in radians
609 : }
610 : }
611 : }
612 :
613 : } // if(m_beta_p > 0 && m_pupil != nullptr)
614 :
615 0 : m_detectorImage.iterNo = m_wfsImage.iterNo;
616 :
617 0 : m_roTime_counter = 0;
618 0 : m_reading = 0;
619 0 : rv = true;
620 : }
621 :
622 : // std::cerr << "DPWFS: " << m_iTime_counter << " " << m_reading << " " << m_roTime_counter << " " << rv << "\n";
623 :
624 0 : return rv;
625 0 : }
626 :
627 : template <typename _realT, typename _detectorT>
628 0 : bool directPhaseSensor<_realT, _detectorT>::senseWavefrontCal( wavefrontT &pupilPlane )
629 : {
630 :
631 0 : m_lastWavefront = 1;
632 :
633 0 : m_wavefronts[0].amplitude = pupilPlane.amplitude;
634 0 : m_wavefronts[0].phase = pupilPlane.phase;
635 :
636 0 : m_wavefronts[1].amplitude = pupilPlane.amplitude;
637 0 : m_wavefronts[1].phase = pupilPlane.phase;
638 :
639 0 : doSenseWavefront();
640 :
641 0 : m_detectorImage.image =
642 0 : m_wfsImage.image.block( 0.5 * ( m_wfsImage.image.rows() - 1 ) - 0.5 * ( m_detectorImage.image.rows() - 1 ),
643 0 : 0.5 * ( m_wfsImage.image.cols() - 1 ) - 0.5 * ( m_detectorImage.image.cols() - 1 ),
644 : m_detectorImage.image.rows(),
645 : m_detectorImage.image.cols() );
646 :
647 0 : if( m_applyFilter )
648 : {
649 0 : m_filter.filter( m_detectorImage.image );
650 : }
651 :
652 0 : return true;
653 : }
654 :
655 : template <typename _realT, typename _detectorT>
656 0 : void directPhaseSensor<_realT, _detectorT>::doSenseWavefront()
657 : {
658 0 : wavefrontT pupilPlane;
659 :
660 : BREAD_CRUMB;
661 :
662 : /* Here make average wavefront for now */
663 0 : int _firstWavefront = m_lastWavefront - m_iTime;
664 0 : if( _firstWavefront < 0 )
665 0 : _firstWavefront += m_wavefronts.size();
666 :
667 : // std::cerr << m_lastWavefront << " " << m_iTime_counter << " " << _firstWavefront << "\n";
668 :
669 0 : pupilPlane.amplitude = 0.5 * m_wavefronts[_firstWavefront].amplitude;
670 0 : pupilPlane.phase = 0.5 * m_wavefronts[_firstWavefront].phase;
671 :
672 0 : pupilPlane.iterNo = 0.5 * m_wavefronts[_firstWavefront].iterNo;
673 :
674 : // std::cerr << "DPS Averaging: " << m_wavefronts[_firstWavefront].iterNo << " " ;
675 : BREAD_CRUMB;
676 :
677 0 : if( m_wavefronts[_firstWavefront].iterNo < m_iTime )
678 : {
679 0 : m_wfsImage.image = pupilPlane.phase;
680 0 : m_wfsImage.iterNo = pupilPlane.iterNo;
681 0 : return;
682 : }
683 :
684 0 : for( int i = 0; i < m_iTime - 1; ++i )
685 : {
686 0 : ++_firstWavefront;
687 0 : if( (size_t)_firstWavefront >= m_wavefronts.size() )
688 0 : _firstWavefront = 0;
689 :
690 : // std::cerr << m_lastWavefront << " " << m_iTime_counter << " " << _firstWavefront << "\n";
691 :
692 0 : pupilPlane.amplitude += m_wavefronts[_firstWavefront].amplitude;
693 0 : pupilPlane.phase += m_wavefronts[_firstWavefront].phase;
694 0 : pupilPlane.iterNo += m_wavefronts[_firstWavefront].iterNo;
695 : }
696 :
697 0 : ++_firstWavefront;
698 0 : if( (size_t)_firstWavefront >= m_wavefronts.size() )
699 0 : _firstWavefront = 0;
700 :
701 : // std::cerr << m_lastWavefront << " " << m_iTime_counter << " " << _firstWavefront << "\n";
702 :
703 0 : pupilPlane.amplitude += 0.5 * m_wavefronts[_firstWavefront].amplitude;
704 0 : pupilPlane.phase += 0.5 * m_wavefronts[_firstWavefront].phase;
705 0 : pupilPlane.iterNo += 0.5 * m_wavefronts[_firstWavefront].iterNo;
706 :
707 0 : pupilPlane.amplitude /= ( m_iTime );
708 0 : pupilPlane.phase /= ( m_iTime );
709 0 : pupilPlane.iterNo /= ( m_iTime );
710 :
711 0 : m_wfsImage.image = pupilPlane.phase;
712 0 : m_wfsImage.iterNo = pupilPlane.iterNo;
713 0 : }
714 :
715 : template <typename _realT, typename _detectorT>
716 0 : void directPhaseSensor<_realT, _detectorT>::applyFilter( bool af )
717 : {
718 0 : m_applyFilter = af;
719 0 : }
720 :
721 : template <typename _realT, typename _detectorT>
722 0 : bool directPhaseSensor<_realT, _detectorT>::applyFilter()
723 : {
724 0 : return m_applyFilter;
725 : }
726 :
727 : template <typename _realT, typename _detectorT>
728 0 : void directPhaseSensor<_realT, _detectorT>::filterWidth( int width )
729 : {
730 0 : int nr = m_detectorImage.image.rows();
731 0 : int nc = m_detectorImage.image.cols();
732 :
733 0 : typename wfsImageT<_realT>::imageT filterMask;
734 :
735 0 : filterMask.resize( nr, nc );
736 0 : filterMask.setZero();
737 :
738 0 : filterMask.block( 0, 0, width + 1, width + 1 ).setConstant( 1.0 );
739 0 : filterMask.block( 0, nc - width, width + 1, width ).setConstant( 1.0 );
740 0 : filterMask.block( nr - width, 0, width, width + 1 ).setConstant( 1.0 );
741 0 : filterMask.block( nr - width, nc - width, width, width ).setConstant( 1.0 );
742 :
743 0 : m_filter.psdSqrt( filterMask, static_cast<realT>( 1 ) / nr, static_cast<realT>( 1 ) / nc );
744 :
745 0 : m_filterWidth = width;
746 0 : }
747 :
748 : template <typename _realT, typename _detectorT>
749 0 : int directPhaseSensor<_realT, _detectorT>::filterWidth()
750 : {
751 0 : return m_filterWidth;
752 : }
753 :
754 : template <typename realT, typename detectorT>
755 0 : const sigproc::psdFilter<realT, 2> &directPhaseSensor<realT, detectorT>::filter()
756 : {
757 0 : return m_filter;
758 : }
759 :
760 : template <typename _realT, typename _detectorT>
761 0 : void directPhaseSensor<_realT, _detectorT>::beta_p( realT bp )
762 : {
763 0 : m_beta_p = bp;
764 0 : }
765 :
766 : template <typename realT, typename detectorT>
767 0 : realT directPhaseSensor<realT, detectorT>::beta_p()
768 : {
769 0 : return m_beta_p;
770 : }
771 :
772 : template <typename realT, typename detectorT>
773 0 : void directPhaseSensor<realT, detectorT>::pupil( pupilT *pupil )
774 : {
775 0 : m_pupil = pupil;
776 0 : }
777 :
778 : template <typename realT, typename detectorT>
779 0 : void directPhaseSensor<realT, detectorT>::pupil( pupilT &pupil )
780 : {
781 0 : m_pupil = &pupil;
782 0 : }
783 :
784 : template <typename realT, typename detectorT>
785 0 : typename directPhaseSensor<realT, detectorT>::pupilT *directPhaseSensor<realT, detectorT>::pupil()
786 : {
787 0 : return m_pupil;
788 : }
789 :
790 : } // namespace sim
791 :
792 : } // namespace AO
793 :
794 : } // namespace mx
795 :
796 : #endif // directPhaseSensor_hpp
|
