pyramidSensorSepQuad.hpp

Go to the documentation of this file.

/** \file
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief
 * \ingroup mxAO_sim_files
 *
 */
 
#ifndef __pyramidSensor_hpp__
#define __pyramidSensor_hpp__
 
#include "mx/imagingUtils.hpp"
#include "mx/fraunhoferImager.hpp"
#include "mx/timeUtils.hpp"
#include "mx/fitsFile.hpp"
// #include "mx/ds9_interface.h"
#include "../../math/constants.hpp"
 
#include "wavefront.hpp"
 
#include "mx/imageTransforms.hpp"
 
#ifdef DEBUG
#define BREAD_CRUMB std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
#else
#define BREAD_CRUMB
#endif
 
namespace mx
{
namespace AO
{
namespace sim
{
 
template <typename _floatT, typename _detectorT>
class pyramidSensor
{
  public:
    typedef _floatT floatT;
 
    typedef std::complex<floatT> complexT;
 
    /// The wavefront data type
    typedef wavefront<floatT> wavefrontT;
 
    /// The wavefront complex field type
    typedef mx::imagingArray<std::complex<floatT>, fftwAllocator<std::complex<floatT>>, 0> complexFieldT;
 
    /// The wavefront sensor detector image type
    typedef Eigen::Array<floatT, Eigen::Dynamic, Eigen::Dynamic> wfsImageT;
 
    typedef _detectorT detectorT;
 
  protected:
    /* Standard WFS Interface: */
    int _wfSz; ///< Size of the wavefront in pixels
 
    int _detRows; ///< The number of rows of the WFS detector.  After forming the image the WFS detector plane is
                  ///< re-binned to this.
    int _detCols; ///< The number of columns of the WFS detector.  After forming the image the WFS detector plane is
                  ///< re-binned to this.
 
    floatT _lambda; ///< Central wavelength, in meters
 
    int _iTime; ///< Integration time in loop steps
 
    int _roTime; ///< Readout time in loop steps
 
    floatT _simStep; ///< The simulation stepsize in seconds.
 
    /* PyWFS Specific: */
 
    floatT _wfPS; ///< Wavefront pixel scale, in meters/pixel
 
    floatT _D; ///< Telescope diameter, in meters
 
    int _modSteps; ///< Number of steps in the modulation simulation
 
    floatT _modRadius; ///< Radius of the modulation in pixels
 
    int _quadSz; ///< The size of the PyWFS quadrant
 
    mx::fraunhoferImager<complexFieldT> fi;
 
    bool tiltsMade;
    std::vector<complexFieldT> tilts;
 
    int _iTime_counter;
 
    int _reading;
 
    int _roTime_counter;
 
    std::vector<wavefrontT> _wavefronts;
 
    int _lastWavefront;
 
    ds9_interface ds9i;
 
  public:
    /// Default c'tor
    pyramidSensor();
 
    /* The standard WFS interface: */
 
    detectorT detector;
 
    /// The image on the detector, resized from wfsImage
    wfsImageT detectorImage;
 
    /// Get the wavefront size in pixels
    /**
     * \returns the wavefront size in pixels
     */
    int wfSz();
 
    /// Set the wavefront size in pixels.
    /**
     * \param sz is the new size
     */
    void wfSz( int sz );
 
    /// Get the detector rows  in pixels
    /**
     * \returns _detRows
     */
    int detRows();
 
    /// Get the detector columns  in pixels
    /**
     * \returns _detCols
     */
    int detCols();
 
    /// Set the detector columns in pixels.
    /**
     * \param sz is the new size
     */
    void detSize( int nrows, int ncols );
 
    /// Get the PyWFS central wavelength
    /**
     * \returns the central wavelength in meters
     */
    floatT lambda();
 
    /// Set the PyWFS central wavelength
    /**
     * \param d is the new central wavelength in meters
     */
    void lambda( floatT l );
 
    /// Get the PyWFS integration time, in time steps
    int iTime();
 
    /// Set the PyWFS integration time, in time steps
    void iTime( int it );
 
    /// Get the PyWFS detector readout time, in time steps
    int roTime();
 
    /// Set the PyWFS detector readout time, in time steps
    void roTime( int rt );
 
    /// Get the simulation step-size, in seconds.
    floatT simStep();
 
    /// Set the simulation step-size, in seconds.
    void simStep( floatT st );
 
    template <typename AOSysT>
    void linkSystem( AOSysT &AOSys );
 
    /// Sense the wavefront aberrations
    /** Returns true if a new wavefront measurement is ready.
     * Retruns false if still integrating.
     */
    bool senseWavefront( wavefrontT &pupilPlane );
 
    /// Sense the wavefront aberrations in a calibration mode
    /** Allows for faster calibrations.
     */
    bool senseWavefrontCal( wavefrontT &pupilPlane );
 
    /* The PyWFS  Specific Interface */
 
    /// Get the wavefront pixel scale in meters per pixel
    /**
     * \returns the wavefront pixel scale in meters/pixel
     */
    int wfPS();
 
    /// Set the wavefront pixel scale in meters per pixel
    /**
     * \param ps is the pixel scale
     */
    void wfPS( floatT ps );
 
    /// Get the telescope diameter
    /**
     * \returns the telescope diameter in meters
     */
    floatT D();
 
    /// Set the telescope diameter
    /**
     * \param d is the new size in meters
     */
    void D( floatT d );
 
    /// Get the number of modulation steps
    /**
     * \returns _modSteps;
     */
    int modSteps();
 
    /// Set the number of modulation steps
    /**
     * \param mSt is the new number of modulation steps
     */
    void modSteps( int mSt );
 
    /// Get the radius of modulation
    /**
     * \returns _modRadius;
     */
    floatT modRadius();
 
    /// Set the modulation radius
    /**
     * \param mR is the new modulation radius in lambda/D
     */
    void modRadius( floatT mR );
 
    /// Get the quadrant size in pixels
    /** This is the size of the quadrant in un-binned wavefront space
     *
     *
     * \returns _quadSz
     */
    int quadSz();
 
    /// Set the quadrant size in pixels.
    /** This is the size of the quadrant in un-binned wavefront space
     * It should be at least the size of the Pupil.  Make larger than the pupil
     * to have smaller pupil images on the PyWFS detector.
     *
     * \param sz is the new size
     */
    void quadSz( int sz );
 
  protected:
    /// The image formed by the WFS
    wfsImageT wfsImage;
 
    void makeTilts();
 
    void doSenseWavefront();
    void doSenseWavefrontNoMod( wavefrontT & );
 
    bool firstRun;
};
 
template <typename _floatT, typename _detectorT>
pyramidSensor<_floatT, _detectorT>::pyramidSensor()
{
    _wfSz = 0;
    _detRows = 0;
    _detCols = 0;
    _lambda = 0;
 
    _modSteps = 16;
    _modRadius = 16;
 
    iTime( 1 );
    _iTime_counter = 0;
 
    _reading = 0;
    _roTime = 1;
    _roTime_counter = 0;
 
    _simStep = 0.001;
 
    tiltsMade = false;
 
    firstRun = true;
 
    ds9_interface_init( &ds9i );
    ds9_interface_set_title( &ds9i, "PyWFS" );
}
 
template <typename _floatT, typename _detectorT>
int pyramidSensor<_floatT, _detectorT>::wfSz()
{
    return _wfSz;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::wfSz( int sz )
{
    if( _wfSz == sz )
        return;
 
    _wfSz = sz;
 
    fi.setWavefrontSizePixels( _wfSz );
 
    tiltsMade = false;
}
 
template <typename _floatT, typename _detectorT>
int pyramidSensor<_floatT, _detectorT>::detRows()
{
    return _detRows;
}
 
template <typename _floatT, typename _detectorT>
int pyramidSensor<_floatT, _detectorT>::detCols()
{
    return _detCols;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::detSize( int nrows, int ncols )
{
    if( _detRows == nrows && _detCols == ncols )
        return;
 
    _detRows = nrows;
    _detCols = ncols;
 
    detector.setSize( _detRows, _detCols );
    detectorImage.resize( _detRows, _detCols );
}
 
template <typename _floatT, typename _detectorT>
_floatT pyramidSensor<_floatT, _detectorT>::lambda()
{
    return _lambda;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::lambda( _floatT l )
{
    _lambda = l;
}
 
template <typename _floatT, typename _detectorT>
template <typename AOSysT>
void pyramidSensor<_floatT, _detectorT>::linkSystem( AOSysT &AOSys )
{
    AOSys.wfs.wfPS( AOSys._wfPS );
    AOSys.wfs.D( AOSys._D );
}
 
template <typename _floatT, typename _detectorT>
int pyramidSensor<_floatT, _detectorT>::wfPS()
{
    return _wfPS;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::wfPS( _floatT ps )
{
    _wfPS = ps;
}
 
template <typename _floatT, typename _detectorT>
_floatT pyramidSensor<_floatT, _detectorT>::D()
{
    return _D;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::D( _floatT d )
{
    _D = d;
}
 
template <typename _floatT, typename _detectorT>
int pyramidSensor<_floatT, _detectorT>::modSteps()
{
    return _modSteps;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::modSteps( int mSt )
{
    _modSteps = mSt;
 
    tiltsMade = false;
}
 
template <typename _floatT, typename _detectorT>
_floatT pyramidSensor<_floatT, _detectorT>::modRadius()
{
    return _modRadius; // * mx::fftPlateScale(_wfSz, _wfPS, _lambda)/(_lambda/D);
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::modRadius( _floatT mR )
{
    _modRadius = mR; // (_lambda/D)/mx::fftPlateScale(_wfSz, _wfPS, _lambda);
 
    tiltsMade = false;
}
 
template <typename _floatT, typename _detectorT>
int pyramidSensor<_floatT, _detectorT>::quadSz()
{
    return _quadSz;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::quadSz( int sz )
{
    if( _quadSz == sz )
        return;
 
    _quadSz = sz;
 
    wfsImage.resize( 2 * _quadSz, 2 * _quadSz );
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::makeTilts()
{
    constexpr _floatT pi = math::pi<_floatT>();
 
    _floatT dang = 2 * pi / _modSteps;
    _floatT dx, dy;
 
    tilts.resize( _modSteps );
 
    std::cout << "WF Size: " << _wfSz << "\n";
    std::cout << "WF PS:   " << _wfPS << "\n";
    std::cout << "Lambda:  " << _lambda << "\n";
 
    std::cout << "Pyr. PS: " << mx::fftPlateScale<floatT>( _wfSz, _wfPS, _lambda ) * 206265. << " (mas/pix)\n";
    std::cout << "Mod rad: " << _modRadius * ( _lambda / _D ) / mx::fftPlateScale<floatT>( _wfSz, _wfPS, _lambda )
              << " (pixels)\n";
 
    for( int i = 0; i < _modSteps; ++i )
    {
        dx = _modRadius * ( _lambda / _D ) / mx::fftPlateScale<floatT>( _wfSz, _wfPS, _lambda ) *
             cos( 0.5 * dang + dang * i );
        dy = _modRadius * ( _lambda / _D ) / mx::fftPlateScale<floatT>( _wfSz, _wfPS, _lambda ) *
             sin( 0.5 * dang + dang * i );
 
        tilts[i].resize( _wfSz, _wfSz );
        tilts[i].set( std::complex<_floatT>( 0, 1 ) );
 
        tiltWavefront( tilts[i], dx, dy );
    }
 
    tiltsMade = true;
}
 
template <typename _floatT, typename _detectorT>
int pyramidSensor<_floatT, _detectorT>::iTime()
{
    return _iTime;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::iTime( int it )
{
    if( it < 1 )
    {
        std::cerr << "iTime must be >= 1.  Correcting\n";
        it = 1;
    }
 
    _iTime = it;
 
    _wavefronts.resize( _iTime + 2 );
    _lastWavefront = -1;
 
    detector.expTime( _simStep * _iTime );
}
 
template <typename _floatT, typename _detectorT>
int pyramidSensor<_floatT, _detectorT>::roTime()
{
    return roTime;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::roTime( int rt )
{
    if( rt < 1 )
    {
        std::cerr << "roTime must be >= 1.  Correcting\n";
        rt = 1;
    }
 
    _roTime = rt;
}
 
template <typename _floatT, typename _detectorT>
_floatT pyramidSensor<_floatT, _detectorT>::simStep()
{
    return simStep;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::simStep( _floatT st )
{
 
    _simStep = st;
 
    detector.expTime( _simStep * _iTime );
}
 
template <typename _floatT, typename _detectorT>
bool pyramidSensor<_floatT, _detectorT>::senseWavefront( wavefrontT &pupilPlane )
{
 
    ++_lastWavefront;
    if( _lastWavefront >= _wavefronts.size() )
        _lastWavefront = 0;
    _wavefronts[_lastWavefront].amplitude = pupilPlane.amplitude;
    _wavefronts[_lastWavefront].phase = pupilPlane.phase;
 
    // Always skip the first one for averaging to center of iTime.
    if( firstRun )
    {
        firstRun = false;
        return false;
    }
 
    ++_iTime_counter;
 
    bool rv = false;
 
    if( _reading )
    {
        ++_roTime_counter;
 
        if( _roTime_counter >= _roTime )
        {
            std::cerr << "PyWFS: reading\n";
            detector.exposeImage( detectorImage, wfsImage );
            ds9_interface_display_raw( &ds9i,
                                       1,
                                       detectorImage.data(),
                                       detectorImage.rows(),
                                       detectorImage.cols(),
                                       1,
                                       mx::getFitsBITPIX<floatT>() );
 
            _roTime_counter = 0;
            _reading = 0;
            rv = true;
        }
    }
 
    if( _iTime_counter >= _iTime )
    {
        std::cerr << "PyWFS: sensing\n";
        doSenseWavefront();
        _iTime_counter = 0;
 
        _reading = 1;
        _roTime_counter = 0;
    }
 
    return rv;
}
 
template <typename _floatT, typename _detectorT>
bool pyramidSensor<_floatT, _detectorT>::senseWavefrontCal( wavefrontT &pupilPlane )
{
 
    _lastWavefront = 1;
 
    _wavefronts[0].amplitude = pupilPlane.amplitude;
    _wavefronts[0].phase = pupilPlane.phase;
 
    _wavefronts[1].amplitude = pupilPlane.amplitude;
    _wavefronts[1].phase = pupilPlane.phase;
 
    doSenseWavefront();
 
    detector.exposeImage( detectorImage, wfsImage );
 
    // ds9_interface_display_raw( &ds9i, 1, detectorImage.data(), detectorImage.rows(), detectorImage.cols(),1,
    // mx::getFitsBITPIX<floatT>());
 
    return true;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::doSenseWavefront()
{
    constexpr _floatT pi = math::pi<_floatT>();
 
    if( tiltsMade == false )
        makeTilts();
 
    // std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
    BREAD_CRUMB;
 
    wavefrontT pupilPlane;
 
    /* Here make average wavefront for now */
    int _firstWavefront = _lastWavefront - _iTime;
    if( _firstWavefront < 0 )
        _firstWavefront += _wavefronts.size();
 
    pupilPlane.amplitude = _wavefronts[_firstWavefront].amplitude;
    pupilPlane.phase = _wavefronts[_firstWavefront].phase;
 
    // std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
    BREAD_CRUMB;
 
    for( int i = 0; i < _iTime; ++i )
    {
        ++_firstWavefront;
        if( _firstWavefront >= _wavefronts.size() )
            _firstWavefront = 0;
 
        pupilPlane.amplitude += _wavefronts[_firstWavefront].amplitude;
        pupilPlane.phase += _wavefronts[_firstWavefront].phase;
    }
 
    pupilPlane.amplitude /= ( _iTime + 1 );
    pupilPlane.phase /= ( _iTime + 1 );
 
    /*=====================================*/
 
    if( _modRadius == 0 )
        return doSenseWavefrontNoMod( pupilPlane );
 
    // std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
    BREAD_CRUMB;
 
    wfsImage.resize( 2 * _quadSz, 2 * _quadSz );
    wfsImage.setZero();
 
    complexFieldT pupilPlaneCF;
    pupilPlane.getWavefront( pupilPlaneCF, _wfSz );
 
    complexFieldT focalPlaneCF;
    focalPlaneCF.resize( _wfSz, _wfSz );
 
#pragma omp parallel
    {
        complexFieldT tiltedPlane;
        complexFieldT focalPlane;
        complexFieldT sensorPlane;
        complexFieldT temp_ul;
        complexFieldT temp_ur;
        complexFieldT temp_ll;
        complexFieldT temp_lr;
 
        wfsImageT sensorImage;
 
        tiltedPlane.resize( _wfSz, _wfSz );
        focalPlane.resize( _wfSz, _wfSz );
        sensorPlane.resize( _wfSz, _wfSz );
 
        temp_ul.resize( _wfSz, _wfSz );
        temp_ul.set( (complexT)0 );
 
        temp_ur.resize( _wfSz, _wfSz );
        temp_ur.set( (complexT)0 );
 
        temp_ll.resize( _wfSz, _wfSz );
        temp_ll.set( (complexT)0 );
 
        temp_lr.resize( _wfSz, _wfSz );
        temp_lr.set( (complexT)0 );
 
        // sensorImage.resize(_wfSz, _wfSz);
        sensorImage.resize( _quadSz * 2, _quadSz * 2 );
        sensorImage.setZero();
 
        int dSz = ( 0.5 * _wfSz - 0.5 * _quadSz + 2 * _quadSz ) - _wfSz;
        if( dSz < 0 )
            dSz = 0;
 
#pragma omp for
        for( int i = 0; i < _modSteps; ++i )
        {
            int nelem = _wfSz * _wfSz;
 
            complexT *tp_data = tiltedPlane.data();
            complexT *pp_data = pupilPlaneCF.data();
            complexT *ti_data = tilts[i].data();
 
            for( int ii = 0; ii < nelem; ++ii )
            {
                tp_data[ii] = pp_data[ii] * ti_data[ii];
            }
 
            fi.propagatePupilToFocal( focalPlane, tiltedPlane );
            temp_ul.set( complexT( 0, 0 ) );
            extractBlock( temp_ul, 0, 0.5 * _wfSz, 0, 0.5 * _wfSz, focalPlane, 0, 0 );
 
            fi.propagateFocalToPupil( sensorPlane, temp_ul );
            extractIntensityImageAccum( sensorImage,
                                        0,
                                        2 * _quadSz - dSz,
                                        0,
                                        2 * _quadSz - dSz,
                                        sensorPlane,
                                        0.5 * _wfSz - 0.5 * _quadSz,
                                        0.5 * _wfSz - 0.5 * _quadSz );
 
            temp_ur.set( complexT( 0, 0 ) );
            extractBlock( temp_ur, 0.5 * _wfSz, 0.5 * _wfSz, 0, 0.5 * _wfSz, focalPlane, 0.5 * _wfSz, 0 );
            fi.propagateFocalToPupil( sensorPlane, temp_ur );
            extractIntensityImageAccum( sensorImage,
                                        dSz,
                                        2 * _quadSz - dSz,
                                        0,
                                        2 * _quadSz - dSz,
                                        sensorPlane,
                                        0.5 * _wfSz - 1.5 * _quadSz + dSz,
                                        0.5 * _wfSz - 0.5 * _quadSz );
 
            temp_ll.set( complexT( 0, 0 ) );
            extractBlock( temp_ll, 0, 0.5 * _wfSz, 0.5 * _wfSz, 0.5 * _wfSz, focalPlane, 0, 0.5 * _wfSz );
            fi.propagateFocalToPupil( sensorPlane, temp_ll );
            extractIntensityImageAccum( sensorImage,
                                        0,
                                        2 * _quadSz - dSz,
                                        dSz,
                                        2 * _quadSz - dSz,
                                        sensorPlane,
                                        0.5 * _wfSz - 0.5 * _quadSz,
                                        0.5 * _wfSz - 1.5 * _quadSz + dSz );
 
            temp_lr.set( complexT( 0, 0 ) );
            extractBlock(
                temp_lr, 0.5 * _wfSz, 0.5 * _wfSz, 0.5 * _wfSz, 0.5 * _wfSz, focalPlane, 0.5 * _wfSz, 0.5 * _wfSz );
            fi.propagateFocalToPupil( sensorPlane, temp_lr );
            extractIntensityImageAccum( sensorImage,
                                        dSz,
                                        2 * _quadSz - dSz,
                                        dSz,
                                        2 * _quadSz - dSz,
                                        sensorPlane,
                                        0.5 * _wfSz - 1.5 * _quadSz + dSz,
                                        0.5 * _wfSz - 1.5 * _quadSz + dSz );
 
        } // for
 
        BREAD_CRUMB;
 
        // std::cerr << wfsImage.rows() << " " << wfsImage.cols() << "\n";
        // std::cerr << sensorImage.rows() << " " << sensorImage.cols() << "\n";
 
#pragma omp critical
        {
            wfsImage += sensorImage;
        }
 
    } // #pragma omp parallel
 
    BREAD_CRUMB;
 
    wfsImage /= _modSteps;
}
 
template <typename _floatT, typename _detectorT>
void pyramidSensor<_floatT, _detectorT>::doSenseWavefrontNoMod( wavefrontT &pupilPlane )
{
    constexpr _floatT pi = math::pi<_floatT>();
 
    BREAD_CRUMB;
 
    wfsImage.resize( 2 * _quadSz, 2 * _quadSz );
    wfsImage.setZero();
 
    complexFieldT pupilPlaneCF;
    pupilPlane.getWavefront( pupilPlaneCF, _wfSz );
 
    complexFieldT focalPlaneCF;
    focalPlaneCF.resize( _wfSz, _wfSz );
 
    complexFieldT tiltedPlane;
    complexFieldT focalPlane;
    complexFieldT sensorPlane;
    complexFieldT temp_ul;
    complexFieldT temp_ur;
    complexFieldT temp_ll;
    complexFieldT temp_lr;
 
    wfsImageT sensorImage;
 
    tiltedPlane.resize( _wfSz, _wfSz );
    focalPlane.resize( _wfSz, _wfSz );
    sensorPlane.resize( _wfSz, _wfSz );
 
    temp_ul.resize( _wfSz, _wfSz );
    temp_ul.set( (complexT)0 );
 
    temp_ur.resize( _wfSz, _wfSz );
    temp_ur.set( (complexT)0 );
 
    temp_ll.resize( _wfSz, _wfSz );
    temp_ll.set( (complexT)0 );
 
    temp_lr.resize( _wfSz, _wfSz );
    temp_lr.set( (complexT)0 );
 
    sensorImage.resize( _quadSz * 2, _quadSz * 2 );
    sensorImage.setZero();
 
    int dSz = ( 0.5 * _wfSz - 0.5 * _quadSz + 2 * _quadSz ) - _wfSz;
    if( dSz < 0 )
        dSz = 0;
    int nelem = _wfSz * _wfSz;
 
    complexT *tp_data = tiltedPlane.data();
    complexT *pp_data = pupilPlaneCF.data();
 
    BREAD_CRUMB;
    for( int ii = 0; ii < nelem; ++ii )
    {
        tp_data[ii] = pp_data[ii]; //*ti_data[ii];
    }
 
    BREAD_CRUMB;
 
    fi.propagatePupilToFocal( focalPlane, tiltedPlane );
    temp_ul.set( complexT( 0, 0 ) );
    extractBlock( temp_ul, 0, 0.5 * _wfSz, 0, 0.5 * _wfSz, focalPlane, 0, 0 );
 
    fi.propagateFocalToPupil( sensorPlane, temp_ul );
    extractIntensityImageAccum( sensorImage,
                                0,
                                2 * _quadSz - dSz,
                                0,
                                2 * _quadSz - dSz,
                                sensorPlane,
                                0.5 * _wfSz - 0.5 * _quadSz,
                                0.5 * _wfSz - 0.5 * _quadSz );
 
    temp_ur.set( complexT( 0, 0 ) );
    extractBlock( temp_ur, 0.5 * _wfSz, 0.5 * _wfSz, 0, 0.5 * _wfSz, focalPlane, 0.5 * _wfSz, 0 );
    fi.propagateFocalToPupil( sensorPlane, temp_ur );
    extractIntensityImageAccum( sensorImage,
                                dSz,
                                2 * _quadSz - dSz,
                                0,
                                2 * _quadSz - dSz,
                                sensorPlane,
                                0.5 * _wfSz - 1.5 * _quadSz + dSz,
                                0.5 * _wfSz - 0.5 * _quadSz );
 
    temp_ll.set( complexT( 0, 0 ) );
    extractBlock( temp_ll, 0, 0.5 * _wfSz, 0.5 * _wfSz, 0.5 * _wfSz, focalPlane, 0, 0.5 * _wfSz );
    fi.propagateFocalToPupil( sensorPlane, temp_ll );
    extractIntensityImageAccum( sensorImage,
                                0,
                                2 * _quadSz - dSz,
                                dSz,
                                2 * _quadSz - dSz,
                                sensorPlane,
                                0.5 * _wfSz - 0.5 * _quadSz,
                                0.5 * _wfSz - 1.5 * _quadSz + dSz );
 
    temp_lr.set( complexT( 0, 0 ) );
    extractBlock( temp_lr, 0.5 * _wfSz, 0.5 * _wfSz, 0.5 * _wfSz, 0.5 * _wfSz, focalPlane, 0.5 * _wfSz, 0.5 * _wfSz );
    fi.propagateFocalToPupil( sensorPlane, temp_lr );
    extractIntensityImageAccum( sensorImage,
                                dSz,
                                2 * _quadSz - dSz,
                                dSz,
                                2 * _quadSz - dSz,
                                sensorPlane,
                                0.5 * _wfSz - 1.5 * _quadSz + dSz,
                                0.5 * _wfSz - 1.5 * _quadSz + dSz );
 
    BREAD_CRUMB;
 
    wfsImage += sensorImage;
 
    BREAD_CRUMB;
}
 
} // namespace sim
} // namespace AO
} // namespace mx
 
#endif //__pyramidSensor_hpp__