directPhaseReconstructorD.hpp

#ifndef directPhaseReconstructor_hpp
#define directPhaseReconstructor_hpp
 
#pragma GCC system_header
#include <Eigen/Dense>
 
#include <mx/signalWindows.hpp>
 
#ifdef DEBUG
#define BREAD_CRUMB std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
#else
#define BREAD_CRUMB
#endif
 
namespace mx
{
 
namespace AO
{
 
namespace sim
{
 
struct slopeReconstructorSpec
{
    std::string dmName;
    std::string basisName;
 
    std::string rMatId;
};
 
/// A Pyramid Wavefront Sensor slope reconstructor.
/** Calculates slopes, normalized by total flux in the image.
 */
template <typename realT>
class directPhaseReconstructor
{
  public:
    /// The type of the measurement (i.e. the slope vector)
    // typedef Eigen::Array<realT,-1,-1> measurementT;
 
    /// The type of the WFS image
    typedef Eigen::Array<realT, -1, -1> imageT;
 
    /// The type of the response matrix
    typedef Eigen::Array<realT, -1, -1> rmatT;
 
    typedef slopeReconstructorSpec specT;
 
  protected:
    Eigen::Array<realT, -1, -1> _recon; ///< The reconstructor matrix.
 
    imageT _rMat; ///< The response matrix
    eigenCube<realT> _rImages;
 
    int _nModes{ 0 }; ///< The number of modes to be reconstructed
 
    int _detRows{ 0 }; ///< The size of the WFS image, in rows
    int _detCols{ 0 }; ///< The size of the WFS image, in columns
 
    int _measurementSize{ 0 }; ///< The number of values in the measurement
 
    // The mirror modes
    mx::eigenCube<realT> *_modes{ nullptr };
 
    imageT *_pupil{ nullptr };
 
    imageT _mask;
 
    realT norm{ 0 };
 
    ds9_interface ds9i;
 
  public:
    mx::eigenImaged _spectrum;
 
    imageT *_gains{ nullptr };
 
  public:
    /// Default c'tor
    directPhaseReconstructor();
 
    template <typename AOSysT>
    void initialize( AOSysT &AOSys, specT &spec )
    {
 
        _modes = &AOSys.dm._infF;
 
        _nModes = _modes->planes();
        _detRows = _modes->rows();
        _detCols = _modes->cols();
 
        _measurementSize = _detRows * _detCols;
 
        _pupil = &AOSys._pupil;
        _measurementSize = _pupil->sum();
 
        _mask.resize( _pupil->rows(), _pupil->cols() );
 
        //          tukey2d(realT *filt, int dim, realT N, realT alpha, realT xc, realT yc)
        // mx::tukey2d(_mask.data(), _mask.rows(), (realT) _mask.rows(), (realT) 0.0, (realT) 0.5*(_mask.rows()-1),
        // (realT) 0.5*(_mask.cols()-1));
        _mask = *_pupil;
 
        mx::fitsFile<realT> ff;
        ff.write( "dprMask.fits", _mask );
 
        std::string recMatrix = mx::AO::path::sys::cal::iMat(
            AOSys._sysName, spec.dmName, AOSys._wfsName, AOSys._pupilName, spec.basisName, spec.rMatId );
        loadRecon( recMatrix );
    }
 
    template <typename AOSysT>
    void linkSystem( AOSysT &AOSys );
 
    /// Get the calibration amplitude used in response matrix acquisition (_calAmp)
    realT calAmp();
 
    /// Set the calibration amplitude used in response matrix acquisition (_calAmp)
    /**
     * \param ca [in] the new calibration amplitude
     */
    void calAmp( realT ca );
 
    /// Get the number of modes (_nModes)
    int nModes();
 
    /// Get the number of detector rows (_detRows)
    int detRows();
 
    /// Set the number of detector rows (_detRows)
    void detRows( int dr );
 
    /// Get the number of detector columns (_detCols)
    int detCols();
 
    /// Set the number of detector columns (_detCols)
    void detCols( int dc );
 
    /// Load the reconstrutor from the specified FITS file
    /**
     * \param fname is the name of the FITS file, including path
     */
    void loadRecon( std::string fname );
 
    /// Return the size of the unbinned measurement
    int measurementSize();
 
    /// Calculate the slope measurement
    /**
     * \param slopes [out] a (_measurementSize X 2)  array of slopes
     * \param wfsImage [in] the WFS image from which to measure the slopes
     */
    template <typename measurementT, typename wfsImageT>
    void calcMeasurement( measurementT &slopes, wfsImageT &wfsImage );
 
    /// Reconstruct the wavefront from the input image, producing the modal amplitude vector
    template <typename measurementT, typename wfsImageT>
    void reconstruct( measurementT &commandVect, wfsImageT &wfsImage );
 
    /// Initialize the response matrix for acquisition
    /**
     * \param nmodes the number of modes
     * \param calamp the calibration amplitude
     * \param detrows the number of detector rows
     * \param detcols the number of detector columns
     */
    void initializeRMat( int nmodes, realT calamp, int detrows, int detcols );
 
    /// Accumalte the next measurement in the response matrix
    /**
     * \param i the measurement index
     * \param measureVec is the i-th measurement vector
     */
    template <typename measurementT>
    void accumulateRMat( int i, measurementT &measureVec );
 
    template <typename measurementT, typename wfsImageT>
    void accumulateRMat( int i, measurementT &measureVec, wfsImageT &wfsImage );
 
    /// Write the accumulated response matrix to disk
    /**
     * \param fname the name, including path, of the response matrix
     */
    void saveRMat( std::string fname );
 
    void saveRImages( std::string fname );
};
 
template <typename realT>
directPhaseReconstructor<realT>::directPhaseReconstructor()
{
}
 
template <typename realT>
template <typename AOSysT>
void directPhaseReconstructor<realT>::linkSystem( AOSysT &AOSys )
{
    _modes = &AOSys.dm._modes;
 
    _nModes = _modes->planes();
    _detRows = _modes->rows();
    _detCols = _modes->cols();
 
    _measurementSize = _detRows * _detCols;
 
    _pupil = &AOSys._pupil;
    _measurementSize = _pupil->sum();
 
    _mask.resize( _pupil->rows(), _pupil->cols() );
    mx::tukey2d(
        _mask.data(), _mask.rows(), (realT)_mask.rows(), 0.0, 0.5 * ( _mask.rows() - 1 ), 0.5 * ( _mask.cols() - 1 ) );
 
    mx::fitsFile<realT> ff;
    ff.write( "dprMask.fits", _mask );
}
 
template <typename realT>
realT directPhaseReconstructor<realT>::calAmp()
{
    return 0.5 * 780.0e-9 / two_pi<realT>();
}
 
template <typename realT>
void directPhaseReconstructor<realT>::calAmp( realT ca )
{
    return;
}
 
template <typename realT>
int directPhaseReconstructor<realT>::nModes()
{
    return _nModes;
}
 
template <typename realT>
int directPhaseReconstructor<realT>::detRows()
{
    return _detRows;
}
 
template <typename realT>
int directPhaseReconstructor<realT>::detCols()
{
    return _detCols;
}
 
template <typename realT>
void directPhaseReconstructor<realT>::loadRecon( std::string fname )
{
    mx::fitsFile<realT> ff;
    mx::fitsHeader head;
 
    ff.read( fname, _recon, head );
}
 
template <typename realT>
int directPhaseReconstructor<realT>::measurementSize()
{
    return _measurementSize;
}
 
template <typename realT>
template <typename measurementT, typename wfsImageT>
void directPhaseReconstructor<realT>::calcMeasurement( measurementT &slopes, wfsImageT &wfsImage )
{
 
    slopes.measurement.resize( 1, _measurementSize );
 
    int k = 0;
 
    for( int i = 0; i < wfsImage.image.rows(); ++i )
    {
        for( int j = 0; j < wfsImage.image.cols(); ++j )
        {
            if( ( *_pupil )( i, j ) == 0 )
                continue;
 
            slopes.measurement( 0, k ) = wfsImage.image( i, j );
            ++k;
        }
    }
}
 
template <typename realT>
template <typename measurementT, typename wfsImageT>
void directPhaseReconstructor<realT>::reconstruct( measurementT &commandVect, wfsImageT &wfsImage )
{
 
    BREAD_CRUMB;
 
    measurementT slopes;
 
    calcMeasurement( slopes, wfsImage );
 
    BREAD_CRUMB;
    commandVect.measurement = slopes.measurement.matrix() * _recon.matrix();
 
    BREAD_CRUMB;
    commandVect.iterNo = wfsImage.iterNo;
}
 
template <typename realT>
void directPhaseReconstructor<realT>::initializeRMat( int nModes, realT calamp, int detRows, int detCols )
{
    _nModes = nModes;
 
    _detRows = detRows;
    _detCols = detCols;
 
    _rMat.resize( measurementSize(), nModes );
    _rMat.setZero();
 
    _rImages.resize( _detRows, _detCols, _nModes );
}
 
template <typename realT>
template <typename measurementT>
void directPhaseReconstructor<realT>::accumulateRMat( int i, measurementT &measureVec )
{
    _rMat.col( i ) = measureVec.measurement.row( 0 );
}
 
template <typename realT>
template <typename measurementT, typename wfsImageT>
void directPhaseReconstructor<realT>::accumulateRMat( int i, measurementT &measureVec, wfsImageT &wfsImage )
{
    BREAD_CRUMB;
    accumulateRMat( i, measureVec );
 
    BREAD_CRUMB;
 
    _rImages.image( i ) = wfsImage.image;
}
 
template <typename realT>
void directPhaseReconstructor<realT>::saveRMat( std::string fname )
{
    mx::fitsFile<realT> ff;
    mx::fitsHeader head;
 
    ff.write( fname, _rMat );
}
 
template <typename realT>
void directPhaseReconstructor<realT>::saveRImages( std::string fname )
{
    mx::fitsFile<realT> ff;
 
    ff.write( fname, _rImages );
}
 
} // namespace sim
} // namespace AO
} // namespace mx
 
#endif // directPhaseReconstructor_hpp