pywfsSlopeReconstructor.hpp

#ifndef __pywfsSlopeReconstructor_hpp__
#define __pywfsSlopeReconstructor_hpp__
 
#include <mx/improc/eigenCube.hpp>
 
namespace mx
{
namespace AO
{
namespace sim
{
 
struct pywfsSlopeReconstructorSpec
{
    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 _floatT>
class pywfsSlopeReconstructor
{
  public:
    typedef _floatT floatT;
 
    /// The type of the measurement (i.e. the slope vector)
    // typedef Eigen::Array<floatT,-1,-1> measurementT;
 
    /// The type of the WFS image
    typedef Eigen::Array<floatT, -1, -1> imageT;
 
    /// The type of the response matrix
    typedef Eigen::Array<floatT, -1, -1> rmatT;
 
    typedef pywfsSlopeReconstructorSpec specT;
 
  protected:
    Eigen::Array<floatT, -1, -1> _recon; ///< The reconstructor matrix.
 
    int _maskType; /// 0 is centrally obscured circle, 1 is supplied by fits files
 
    floatT _maskRadius;      ///< The radius of the quadrant mask
    floatT _maskObscuration; ///< The central obscuration of the quadrant mask
 
    std::string _maskFile; ///< The name of the quadrant mask file
 
    bool _maskMade; ///< Whether or not the mask has been made
 
    Eigen::Array<floatT, -1, -1> _quadMask; ///< The quadrant mask
    void calcMask();                        ///< Calculates the quadrant mask
    int _measurementSize;                   ///< The number of slopes in the measurement
 
    floatT _calAmp; ///< The calibration amplitude used for response matrix acquisition
 
    int _nModes; ///< The number of modes to be reconstructed
 
    int _detRows; ///< The size of the WFS image, in rows
    int _detCols; ///< The size of the WFS image, in columns
 
    imageT _rMat; ///< The response matrix
    improc::eigenCube<floatT> _rImages;
 
  public:
    int _binFact; ///< The binning to apply before reconstructing.
 
    /// Default c'tor
    pywfsSlopeReconstructor();
 
    template <typename AOSysT>
    void initialize( AOSysT &AOSys, specT &spec )
    {
        std::string recMatrix = mx::AO::path::sys::cal::iMat(
            AOSys._sysName, spec.dmName, AOSys._wfsName, AOSys._pupilName, spec.basisName, spec.rMatId );
        loadRecon( recMatrix );
    }
 
    /// Get the quadrant mask radius (_maskRadius)
    floatT maskRadius();
 
    /// Set the quadrant mask radius (_maskRadius)
    /** Calling this will cause the quadrant mask to be recalculated next time it is needed.
     *
     * \param mr [in] the new mask radius
     */
    void maskRadius( floatT mr );
 
    /// Get the quadrant mask central obscuration ratio (_maskObscuration)
    floatT maskObscuration();
 
    /// Set the quadrant mask central obscuration ratio (_maskObscuration)
    /** Calling this will cause the quadrant mask to be recalculated next time it is needed.
     *
     * \param mo [in] the new central obscuration ratio
     */
    void maskObscuration( floatT mo );
 
    /// Get the quadrant mask file name (_maskFile)
    std::string maskFile();
 
    /// Set the quadrant mask file name (_maskFile)
    /** Calling this will cause the quadrant mask to be reloaded next time it is needed.
     *
     * \param mf [in] the new mask file
     */
    void maskFile( const std::string &mf );
 
    /// Get the calibration amplitude used in response matrix acquisition (_calAmp)
    floatT calAmp();
 
    /// Set the calibration amplitude used in response matrix acquisition (_calAmp)
    /**
     * \param ca [in] the new calibration amplitude
     */
    void calAmp( floatT 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, floatT 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 floatT>
pywfsSlopeReconstructor<floatT>::pywfsSlopeReconstructor()
{
    _maskRadius = 0;
    _maskObscuration = 0;
 
    _maskMade = false;
 
    _binFact = 1;
}
 
template <typename floatT>
floatT pywfsSlopeReconstructor<floatT>::maskRadius()
{
    return _maskRadius;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::maskRadius( floatT mr )
{
    _maskRadius = mr;
    _maskType = 0;
    _maskMade = false;
}
 
template <typename floatT>
floatT pywfsSlopeReconstructor<floatT>::maskObscuration()
{
    return _maskObscuration;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::maskObscuration( floatT mo )
{
    _maskObscuration = mo;
    _maskType = 0;
    _maskMade = false;
}
 
template <typename floatT>
std::string pywfsSlopeReconstructor<floatT>::maskFile()
{
    return _maskFile;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::maskFile( const std::string &mf )
{
    _maskFile = mf;
    _maskType = 1;
    _maskMade = false;
}
 
template <typename floatT>
floatT pywfsSlopeReconstructor<floatT>::calAmp()
{
    return _calAmp;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::calAmp( floatT ca )
{
    _calAmp = ca;
}
 
template <typename floatT>
int pywfsSlopeReconstructor<floatT>::nModes()
{
    return _nModes;
}
 
template <typename floatT>
int pywfsSlopeReconstructor<floatT>::detRows()
{
    return _detRows;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::detRows( int dr )
{
    _detRows = dr;
    _maskMade = false;
}
 
template <typename floatT>
int pywfsSlopeReconstructor<floatT>::detCols()
{
    return _detCols;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::detCols( int dc )
{
    _detCols = dc;
    _maskMade = false;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::loadRecon( std::string fname )
{
    improc::fitsFile<floatT> ff;
    improc::fitsHeader head;
 
    ff.read( fname, _recon, head );
 
    _maskRadius = head["MASKRAD"].Value<floatT>();
    _maskObscuration = head["MASKOBS"].Value<floatT>();
    _maskFile = head["MASKFILE"].String();
 
    if( _maskFile != "" )
        _maskType = 1;
 
    _calAmp = head["CALAMP"].Value<floatT>();
 
    _nModes = head["NMODES"].Value<int>();
 
    _detRows = head["DETROWS"].Int();
    _detCols = head["DETCOLS"].Int();
 
    _maskMade = false;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::calcMask()
{
    if( _maskType == 1 )
    {
        improc::fitsFile<floatT> ff;
 
        std::cerr << "Loading Mask: " << _maskFile << "\n";
        ff.read( _maskFile, _quadMask );
    }
    {
        _quadMask.resize( 0.5 * _detRows / _binFact, 0.5 * _detCols / _binFact );
        wfp::circularPupil( _quadMask, _maskObscuration, _maskRadius / _binFact );
    }
 
    _measurementSize = 2 * _quadMask.sum(); // + 64*64;
 
    _maskMade = true;
 
    improc::fitsFile<floatT> ff;
    ff.write( "quadMask.fits", _quadMask );
}
 
template <typename floatT>
int pywfsSlopeReconstructor<floatT>::measurementSize()
{
    if( !_maskMade )
        calcMask();
    return _measurementSize;
}
 
template <typename floatT>
template <typename measurementT, typename wfsImageT>
void pywfsSlopeReconstructor<floatT>::calcMeasurement( measurementT &slopes, wfsImageT &wfsImage )
{
    if( !_maskMade )
        calcMask();
 
    imageT _wfsImage;
 
    if( _binFact > 1 )
    {
        std::cout << "rebinning" << "\n";
        _wfsImage.resize( wfsImage.image.rows() / _binFact, wfsImage.image.cols() / _binFact );
        imageDownSample( _wfsImage, wfsImage.image );
    }
    else
    {
        _wfsImage = wfsImage.image;
    }
 
    int nsz = _wfsImage.rows();
 
    int nPix = _quadMask.sum();
 
    std::vector<int> x( nPix ), y( nPix );
 
    int k = 0;
    for( int i = 0; i < _quadMask.rows(); ++i )
    {
        for( int j = 0; j < _quadMask.rows(); ++j )
        {
            if( _quadMask( i, j ) == 1 )
            {
                x[k] = i;
                y[k] = j;
                ++k;
            }
        }
    }
 
    slopes.measurement.resize( 1, 2. * nPix ); // + 64*64); //wfsImage.tipImage.rows()*wfsImage.tipImage.cols());
 
    floatT I0, I1, I2, I3;
 
    floatT norm = _wfsImage.sum();
 
    for( int i = 0; i < nPix; ++i )
    {
        I0 = _wfsImage( x[i], y[i] );
        I1 = _wfsImage( x[i] + 0.5 * nsz, y[i] );
        I2 = _wfsImage( x[i], y[i] + 0.5 * nsz );
        I3 = _wfsImage( x[i] + 0.5 * nsz, y[i] + 0.5 * nsz );
 
        if( norm == 0 )
        {
            slopes.measurement( 0, i ) = 0;
            ;
            slopes.measurement( 0, i + nPix ) = 0; //
        }
        else
        {
            slopes.measurement( 0, i ) = ( I0 + I1 - I2 - I3 ) / norm;        //(I0+I1+I2+I3);
            slopes.measurement( 0, i + nPix ) = ( I0 + I2 - I1 - I3 ) / norm; //(I0+I1+I2+I3);
        }
    }
 
    /*
    for(int i=0; i< wfsImage.tipImage.rows(); ++i)
    {
       for(int j=0;j<wfsImage.tipImage.cols(); ++j)
       {
          slopes.measurement(0, 2*nPix + i*wfsImage.tipImage.cols() + j) = wfsImage.tipImage(i,j);
       }
    }*/
}
 
template <typename floatT>
template <typename measurementT, typename wfsImageT>
void pywfsSlopeReconstructor<floatT>::reconstruct( measurementT &commandVect, wfsImageT &wfsImage )
{
    measurementT slopes;
 
    calcMeasurement( slopes, wfsImage );
 
    commandVect.measurement = slopes.measurement.matrix() * _recon.matrix();
 
    commandVect.iterNo = wfsImage.iterNo;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::initializeRMat( int nModes, floatT calamp, int detRows, int detCols )
{
 
    _nModes = nModes;
 
    calAmp( calamp );
 
    _detRows = detRows;
    _detCols = detCols;
 
    _rMat.resize( measurementSize(), nModes );
    _rMat.setZero();
 
    _rImages.resize( _detRows, _detCols, _nModes );
}
 
template <typename floatT>
template <typename measurementT>
void pywfsSlopeReconstructor<floatT>::accumulateRMat( int i, measurementT &measureVec )
{
    _rMat.col( i ) = measureVec.measurement.row( 0 );
}
 
template <typename floatT>
template <typename measurementT, typename wfsImageT>
void pywfsSlopeReconstructor<floatT>::accumulateRMat( int i, measurementT &measureVec, wfsImageT &wfsImage )
{
    accumulateRMat( i, measureVec );
 
    _rImages.image( i ) = wfsImage.image;
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::saveRMat( std::string fname )
{
    improc::fitsFile<floatT> ff;
    improc::fitsHeader head;
 
    if( _maskType == 1 )
    {
        head.append( "MASKFILE", maskFile(), "Name of mask file" );
    }
    else
    {
        head.append( "MASKRAD", maskRadius(), "Mask radius, in pixels" );
        head.append( "MASKOBS", maskObscuration(), "Mask fractional central obscuration" );
    }
 
    head.append( "DETROWS", _detRows, "WFS detector rows" );
    head.append( "DETCOLS", _detCols, "WFS detector cols" );
    head.append( "CALAMP", _calAmp, "DM Calibration amplitude" );
    head.append( "NMODES", _nModes, "Number of modes included in the response matrix." );
 
    // ff.write(fname, _rMat.data(), _rMat.rows(), _rMat.cols(), 1, &head);
    ff.write( fname, _rMat, head );
}
 
template <typename floatT>
void pywfsSlopeReconstructor<floatT>::saveRImages( std::string fname )
{
    improc::fitsFile<floatT> ff;
    improc::fitsHeader head;
 
    if( _maskType == 1 )
    {
        head.append( "MASKFILE", maskFile(), "Name of mask file" );
    }
    else
    {
        head.append( "MASKRAD", maskRadius(), "Mask radius, in pixels" );
        head.append( "MASKOBS", maskObscuration(), "Mask fractional central obscuration" );
    }
 
    head.append( "DETROWS", _detRows, "WFS detector rows" );
    head.append( "DETCOLS", _detCols, "WFS detector cols" );
    head.append( "CALAMP", _calAmp, "DM Calibration amplitude" );
    head.append( "NMODES", _nModes, "Number of modes included in the response matrix." );
 
    // ff.write(fname, _rImages.data(), _rImages.rows(), _rImages.cols(), _rImages.planes(), &head);
    ff.write( fname, _rImages, head );
}
 
} // namespace sim
} // namespace AO
} // namespace mx
 
#endif //__pywfsSlopeReconstructor_hpp__