Line data Source code
1 : /** \file pywfsSlopeReconstructor.hpp
2 : * \brief Declares the pyramid-wavefront-sensor slope reconstructor.
3 : * \author Jared R. Males (jaredmales@gmail.com)
4 : */
5 :
6 : #ifndef __pywfsSlopeReconstructor_hpp__
7 : #define __pywfsSlopeReconstructor_hpp__
8 :
9 : #include <string>
10 :
11 : #include "../../ao/aoPaths.hpp"
12 : #include "../../improc/eigenCube.hpp"
13 : #include "../../ioutils/fits/fitsFile.hpp"
14 : #include "../../wfp/imagingUtils.hpp"
15 :
16 : namespace mx
17 : {
18 : namespace AO
19 : {
20 : namespace sim
21 : {
22 :
23 : struct pywfsSlopeReconstructorSpec
24 : {
25 : std::string dmName;
26 : std::string basisName;
27 :
28 : std::string rMatId;
29 : };
30 :
31 : /// A Pyramid Wavefront Sensor slope reconstructor.
32 : /** Calculates slopes, normalized by total flux in the image.
33 : */
34 : template <typename _floatT>
35 : class pywfsSlopeReconstructor
36 : {
37 : public:
38 : typedef _floatT floatT;
39 :
40 : /// The type of the measurement (i.e. the slope vector)
41 : // typedef Eigen::Array<floatT,-1,-1> measurementT;
42 :
43 : /// The type of the WFS image
44 : typedef Eigen::Array<floatT, -1, -1> imageT;
45 :
46 : /// The type of the response matrix
47 : typedef Eigen::Array<floatT, -1, -1> rmatT;
48 :
49 : typedef pywfsSlopeReconstructorSpec specT;
50 :
51 : protected:
52 : Eigen::Array<floatT, -1, -1> _recon; ///< The reconstructor matrix.
53 :
54 : int _maskType; /// 0 is centrally obscured circle, 1 is supplied by fits files
55 :
56 : floatT _maskRadius; ///< The radius of the quadrant mask
57 : floatT _maskObscuration; ///< The central obscuration of the quadrant mask
58 :
59 : std::string _maskFile; ///< The name of the quadrant mask file
60 :
61 : bool _maskMade; ///< Whether or not the mask has been made
62 :
63 : Eigen::Array<floatT, -1, -1> _quadMask; ///< The quadrant mask
64 : void calcMask(); ///< Calculates the quadrant mask
65 : int _measurementSize; ///< The number of slopes in the measurement
66 :
67 : floatT _calAmp; ///< The calibration amplitude used for response matrix acquisition
68 :
69 : int _nModes; ///< The number of modes to be reconstructed
70 :
71 : int _detRows; ///< The size of the WFS image, in rows
72 : int _detCols; ///< The size of the WFS image, in columns
73 :
74 : imageT _rMat; ///< The response matrix
75 : improc::eigenCube<floatT> _rImages;
76 :
77 : public:
78 : int _binFact; ///< The binning to apply before reconstructing.
79 :
80 : /// Default c'tor
81 : pywfsSlopeReconstructor();
82 :
83 : template <typename AOSysT>
84 : void initialize( AOSysT &AOSys, specT &spec )
85 : {
86 : std::string recMatrix = mx::AO::path::sys::cal::iMat(
87 : AOSys._sysName, spec.dmName, AOSys._wfsName, AOSys._pupilName, spec.basisName, spec.rMatId );
88 : loadRecon( recMatrix );
89 : }
90 :
91 : /// Get the quadrant mask radius (_maskRadius)
92 : floatT maskRadius();
93 :
94 : /// Set the quadrant mask radius (_maskRadius)
95 : /** Calling this will cause the quadrant mask to be recalculated next time it is needed.
96 : *
97 : * \param mr [in] the new mask radius
98 : */
99 : void maskRadius( floatT mr );
100 :
101 : /// Get the quadrant mask central obscuration ratio (_maskObscuration)
102 : floatT maskObscuration();
103 :
104 : /// Set the quadrant mask central obscuration ratio (_maskObscuration)
105 : /** Calling this will cause the quadrant mask to be recalculated next time it is needed.
106 : *
107 : * \param mo [in] the new central obscuration ratio
108 : */
109 : void maskObscuration( floatT mo );
110 :
111 : /// Get the quadrant mask file name (_maskFile)
112 : std::string maskFile();
113 :
114 : /// Set the quadrant mask file name (_maskFile)
115 : /** Calling this will cause the quadrant mask to be reloaded next time it is needed.
116 : *
117 : * \param mf [in] the new mask file
118 : */
119 : void maskFile( const std::string &mf );
120 :
121 : /// Get the calibration amplitude used in response matrix acquisition (_calAmp)
122 : floatT calAmp();
123 :
124 : /// Set the calibration amplitude used in response matrix acquisition (_calAmp)
125 : /**
126 : * \param ca [in] the new calibration amplitude
127 : */
128 : void calAmp( floatT ca );
129 :
130 : /// Get the number of modes (_nModes)
131 : int nModes();
132 :
133 : /// Get the number of detector rows (_detRows)
134 : int detRows();
135 :
136 : /// Set the number of detector rows (_detRows)
137 : void detRows( int dr );
138 :
139 : /// Get the number of detector columns (_detCols)
140 : int detCols();
141 :
142 : /// Set the number of detector columns (_detCols)
143 : void detCols( int dc );
144 :
145 : /// Load the reconstrutor from the specified FITS file
146 : /**
147 : * \param fname is the name of the FITS file, including path
148 : */
149 : void loadRecon( std::string fname );
150 :
151 : /// Return the size of the unbinned measurement
152 : int measurementSize();
153 :
154 : /// Calculate the slope measurement
155 : /**
156 : * \param slopes [out] a (_measurementSize X 2) array of slopes
157 : * \param wfsImage [in] the WFS image from which to measure the slopes
158 : */
159 : template <typename measurementT, typename wfsImageT>
160 : void calcMeasurement( measurementT &slopes, wfsImageT &wfsImage );
161 :
162 : /// Reconstruct the wavefront from the input image, producing the modal amplitude vector
163 : template <typename measurementT, typename wfsImageT>
164 : void reconstruct( measurementT &commandVect, wfsImageT &wfsImage );
165 :
166 : /// Initialize the response matrix for acquisition
167 : /**
168 : * \param nmodes the number of modes
169 : * \param calamp the calibration amplitude
170 : * \param detrows the number of detector rows
171 : * \param detcols the number of detector columns
172 : */
173 : void initializeRMat( int nmodes, floatT calamp, int detrows, int detcols );
174 :
175 : /// Accumalte the next measurement in the response matrix
176 : /**
177 : * \param i the measurement index
178 : * \param measureVec is the i-th measurement vector
179 : */
180 : template <typename measurementT>
181 : void accumulateRMat( int i, measurementT &measureVec );
182 :
183 : template <typename measurementT, typename wfsImageT>
184 : void accumulateRMat( int i, measurementT &measureVec, wfsImageT &wfsImage );
185 :
186 : /// Write the accumulated response matrix to disk
187 : /**
188 : * \param fname the name, including path, of the response matrix
189 : */
190 : void saveRMat( std::string fname );
191 :
192 : void saveRImages( std::string fname );
193 : };
194 :
195 : template <typename floatT>
196 1 : pywfsSlopeReconstructor<floatT>::pywfsSlopeReconstructor()
197 : {
198 1 : _maskRadius = 0;
199 1 : _maskObscuration = 0;
200 :
201 1 : _maskMade = false;
202 :
203 1 : _binFact = 1;
204 1 : }
205 :
206 : template <typename floatT>
207 0 : floatT pywfsSlopeReconstructor<floatT>::maskRadius()
208 : {
209 0 : return _maskRadius;
210 : }
211 :
212 : template <typename floatT>
213 0 : void pywfsSlopeReconstructor<floatT>::maskRadius( floatT mr )
214 : {
215 0 : _maskRadius = mr;
216 0 : _maskType = 0;
217 0 : _maskMade = false;
218 0 : }
219 :
220 : template <typename floatT>
221 0 : floatT pywfsSlopeReconstructor<floatT>::maskObscuration()
222 : {
223 0 : return _maskObscuration;
224 : }
225 :
226 : template <typename floatT>
227 0 : void pywfsSlopeReconstructor<floatT>::maskObscuration( floatT mo )
228 : {
229 0 : _maskObscuration = mo;
230 0 : _maskType = 0;
231 0 : _maskMade = false;
232 0 : }
233 :
234 : template <typename floatT>
235 0 : std::string pywfsSlopeReconstructor<floatT>::maskFile()
236 : {
237 0 : return _maskFile;
238 : }
239 :
240 : template <typename floatT>
241 1 : void pywfsSlopeReconstructor<floatT>::maskFile( const std::string &mf )
242 : {
243 1 : _maskFile = mf;
244 1 : _maskType = 1;
245 1 : _maskMade = false;
246 1 : }
247 :
248 : template <typename floatT>
249 0 : floatT pywfsSlopeReconstructor<floatT>::calAmp()
250 : {
251 0 : return _calAmp;
252 : }
253 :
254 : template <typename floatT>
255 0 : void pywfsSlopeReconstructor<floatT>::calAmp( floatT ca )
256 : {
257 0 : _calAmp = ca;
258 0 : }
259 :
260 : template <typename floatT>
261 0 : int pywfsSlopeReconstructor<floatT>::nModes()
262 : {
263 0 : return _nModes;
264 : }
265 :
266 : template <typename floatT>
267 0 : int pywfsSlopeReconstructor<floatT>::detRows()
268 : {
269 0 : return _detRows;
270 : }
271 :
272 : template <typename floatT>
273 1 : void pywfsSlopeReconstructor<floatT>::detRows( int dr )
274 : {
275 1 : _detRows = dr;
276 1 : _maskMade = false;
277 1 : }
278 :
279 : template <typename floatT>
280 0 : int pywfsSlopeReconstructor<floatT>::detCols()
281 : {
282 0 : return _detCols;
283 : }
284 :
285 : template <typename floatT>
286 1 : void pywfsSlopeReconstructor<floatT>::detCols( int dc )
287 : {
288 1 : _detCols = dc;
289 1 : _maskMade = false;
290 1 : }
291 :
292 : template <typename floatT>
293 0 : void pywfsSlopeReconstructor<floatT>::loadRecon( std::string fname )
294 : {
295 0 : fits::fitsFile<floatT> ff;
296 0 : fits::fitsHeader<> head;
297 :
298 0 : ff.read( _recon, head, fname );
299 :
300 0 : _maskRadius = head["MASKRAD"].value<floatT>();
301 0 : _maskObscuration = head["MASKOBS"].value<floatT>();
302 0 : _maskFile = head["MASKFILE"].value<std::string>();
303 :
304 0 : if( _maskFile != "" )
305 0 : _maskType = 1;
306 :
307 0 : _calAmp = head["CALAMP"].value<floatT>();
308 :
309 0 : _nModes = head["NMODES"].value<int>();
310 :
311 0 : _detRows = head["DETROWS"].value<int>();
312 0 : _detCols = head["DETCOLS"].value<int>();
313 :
314 0 : _maskMade = false;
315 0 : }
316 :
317 : template <typename floatT>
318 1 : void pywfsSlopeReconstructor<floatT>::calcMask()
319 : {
320 1 : if( _maskType == 1 )
321 : {
322 1 : fits::fitsFile<floatT> ff;
323 :
324 1 : std::cerr << "Loading Mask: " << _maskFile << "\n";
325 1 : ff.read( _quadMask, _maskFile );
326 1 : }
327 : else
328 : {
329 0 : _quadMask.resize( 0.5 * _detRows / _binFact, 0.5 * _detCols / _binFact );
330 0 : wfp::circularPupil( _quadMask, _maskObscuration, _maskRadius / _binFact );
331 : }
332 :
333 1 : _measurementSize = 2 * _quadMask.sum(); // + 64*64;
334 :
335 1 : _maskMade = true;
336 :
337 1 : fits::fitsFile<floatT> ff;
338 2 : ff.write( "quadMask.fits", _quadMask );
339 1 : }
340 :
341 : template <typename floatT>
342 1 : int pywfsSlopeReconstructor<floatT>::measurementSize()
343 : {
344 1 : if( !_maskMade )
345 1 : calcMask();
346 1 : return _measurementSize;
347 : }
348 :
349 : template <typename floatT>
350 : template <typename measurementT, typename wfsImageT>
351 : void pywfsSlopeReconstructor<floatT>::calcMeasurement( measurementT &slopes, wfsImageT &wfsImage )
352 : {
353 : if( !_maskMade )
354 : calcMask();
355 :
356 : imageT _wfsImage;
357 :
358 : if( _binFact > 1 )
359 : {
360 : std::cout << "rebinning" << "\n";
361 : _wfsImage.resize( wfsImage.image.rows() / _binFact, wfsImage.image.cols() / _binFact );
362 : imageDownSample( _wfsImage, wfsImage.image );
363 : }
364 : else
365 : {
366 : _wfsImage = wfsImage.image;
367 : }
368 :
369 : int nsz = _wfsImage.rows();
370 :
371 : int nPix = _quadMask.sum();
372 :
373 : std::vector<int> x( nPix ), y( nPix );
374 :
375 : int k = 0;
376 : for( int i = 0; i < _quadMask.rows(); ++i )
377 : {
378 : for( int j = 0; j < _quadMask.rows(); ++j )
379 : {
380 : if( _quadMask( i, j ) == 1 )
381 : {
382 : x[k] = i;
383 : y[k] = j;
384 : ++k;
385 : }
386 : }
387 : }
388 :
389 : slopes.measurement.resize( 1, 2. * nPix ); // + 64*64); //wfsImage.tipImage.rows()*wfsImage.tipImage.cols());
390 :
391 : floatT I0, I1, I2, I3;
392 :
393 : floatT norm = _wfsImage.sum();
394 :
395 : for( int i = 0; i < nPix; ++i )
396 : {
397 : I0 = _wfsImage( x[i], y[i] );
398 : I1 = _wfsImage( x[i] + 0.5 * nsz, y[i] );
399 : I2 = _wfsImage( x[i], y[i] + 0.5 * nsz );
400 : I3 = _wfsImage( x[i] + 0.5 * nsz, y[i] + 0.5 * nsz );
401 :
402 : if( norm == 0 )
403 : {
404 : slopes.measurement( 0, i ) = 0;
405 : ;
406 : slopes.measurement( 0, i + nPix ) = 0; //
407 : }
408 : else
409 : {
410 : slopes.measurement( 0, i ) = ( I0 + I1 - I2 - I3 ) / norm; //(I0+I1+I2+I3);
411 : slopes.measurement( 0, i + nPix ) = ( I0 + I2 - I1 - I3 ) / norm; //(I0+I1+I2+I3);
412 : }
413 : }
414 :
415 : /*
416 : for(int i=0; i< wfsImage.tipImage.rows(); ++i)
417 : {
418 : for(int j=0;j<wfsImage.tipImage.cols(); ++j)
419 : {
420 : slopes.measurement(0, 2*nPix + i*wfsImage.tipImage.cols() + j) = wfsImage.tipImage(i,j);
421 : }
422 : }*/
423 : }
424 :
425 : template <typename floatT>
426 : template <typename measurementT, typename wfsImageT>
427 : void pywfsSlopeReconstructor<floatT>::reconstruct( measurementT &commandVect, wfsImageT &wfsImage )
428 : {
429 : measurementT slopes;
430 :
431 : calcMeasurement( slopes, wfsImage );
432 :
433 : commandVect.measurement = slopes.measurement.matrix() * _recon.matrix();
434 :
435 : commandVect.iterNo = wfsImage.iterNo;
436 : }
437 :
438 : template <typename floatT>
439 0 : void pywfsSlopeReconstructor<floatT>::initializeRMat( int nModes, floatT calamp, int detRows, int detCols )
440 : {
441 :
442 0 : _nModes = nModes;
443 :
444 0 : calAmp( calamp );
445 :
446 0 : _detRows = detRows;
447 0 : _detCols = detCols;
448 :
449 0 : _rMat.resize( measurementSize(), nModes );
450 0 : _rMat.setZero();
451 :
452 0 : _rImages.resize( _detRows, _detCols, _nModes );
453 0 : }
454 :
455 : template <typename floatT>
456 : template <typename measurementT>
457 : void pywfsSlopeReconstructor<floatT>::accumulateRMat( int i, measurementT &measureVec )
458 : {
459 : _rMat.col( i ) = measureVec.measurement.row( 0 );
460 : }
461 :
462 : template <typename floatT>
463 : template <typename measurementT, typename wfsImageT>
464 : void pywfsSlopeReconstructor<floatT>::accumulateRMat( int i, measurementT &measureVec, wfsImageT &wfsImage )
465 : {
466 : accumulateRMat( i, measureVec );
467 :
468 : _rImages.image( i ) = wfsImage.image;
469 : }
470 :
471 : template <typename floatT>
472 0 : void pywfsSlopeReconstructor<floatT>::saveRMat( std::string fname )
473 : {
474 0 : fits::fitsFile<floatT> ff;
475 0 : fits::fitsHeader<> head;
476 :
477 0 : if( _maskType == 1 )
478 : {
479 0 : head.append( "MASKFILE", maskFile(), "Name of mask file" );
480 : }
481 : else
482 : {
483 0 : head.append( "MASKRAD", maskRadius(), "Mask radius, in pixels" );
484 0 : head.append( "MASKOBS", maskObscuration(), "Mask fractional central obscuration" );
485 : }
486 :
487 0 : head.append( "DETROWS", _detRows, "WFS detector rows" );
488 0 : head.append( "DETCOLS", _detCols, "WFS detector cols" );
489 0 : head.append( "CALAMP", _calAmp, "DM Calibration amplitude" );
490 0 : head.append( "NMODES", _nModes, "Number of modes included in the response matrix." );
491 :
492 : // ff.write(fname, _rMat.data(), _rMat.rows(), _rMat.cols(), 1, &head);
493 0 : ff.write( fname, _rMat, head );
494 0 : }
495 :
496 : template <typename floatT>
497 0 : void pywfsSlopeReconstructor<floatT>::saveRImages( std::string fname )
498 : {
499 0 : fits::fitsFile<floatT> ff;
500 0 : fits::fitsHeader<> head;
501 :
502 0 : if( _maskType == 1 )
503 : {
504 0 : head.append( "MASKFILE", maskFile(), "Name of mask file" );
505 : }
506 : else
507 : {
508 0 : head.append( "MASKRAD", maskRadius(), "Mask radius, in pixels" );
509 0 : head.append( "MASKOBS", maskObscuration(), "Mask fractional central obscuration" );
510 : }
511 :
512 0 : head.append( "DETROWS", _detRows, "WFS detector rows" );
513 0 : head.append( "DETCOLS", _detCols, "WFS detector cols" );
514 0 : head.append( "CALAMP", _calAmp, "DM Calibration amplitude" );
515 0 : head.append( "NMODES", _nModes, "Number of modes included in the response matrix." );
516 :
517 : // ff.write(fname, _rImages.data(), _rImages.rows(), _rImages.cols(), _rImages.planes(), &head);
518 0 : ff.write( fname, _rImages, head );
519 0 : }
520 :
521 : } // namespace sim
522 : } // namespace AO
523 : } // namespace mx
524 :
525 : #endif //__pywfsSlopeReconstructor_hpp__
|
