directPhaseSensor.hpp

Go to the documentation of this file.

/** \file directPhaseSensor.hpp
  * \author Jared R. Males (jaredmales@gmail.com)
  * \brief Declaration and definition of a direct phase sensor.
  * \ingroup mxAO_sim_files
  *
  */
 
#ifndef directPhaseSensor_hpp
#define directPhaseSensor_hpp
 
#include "../../math/randomT.hpp"
 
#include "../../ioutils/fits/fitsFile.hpp"
#include "../../improc/eigenCube.hpp"
 
//#include <mx/fraunhoferImager.hpp>
#include "../../sigproc/psdFilter.hpp"
 
#include "wavefront.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 _realT>
struct wfsImageT
{
   typedef _realT realT;
 
   realT iterNo;
 
   ///The wavefront sensor detector image type
   typedef Eigen::Array< realT, Eigen::Dynamic, Eigen::Dynamic> imageT;
 
   imageT image;
};
 
 
template<typename _realT, typename _detectorT>
class directPhaseSensor
{
public:
   typedef _realT realT;
 
   typedef std::complex<_realT> complexT;
 
   ///The wavefront data type
   typedef wavefront<_realT> wavefrontT;
 
   ///The pupil type
   typedef improc::eigenImage<realT> pupilT;
   
   ///The wavefront complex field type
   typedef mx::wfp::imagingArray<std::complex<_realT>,mx::wfp::fftwAllocator<std::complex<_realT> >, 0> complexFieldT;
 
   typedef _detectorT detectorT;
 
   
   typedef mx::math::randomT<realT, std::mt19937_64, std::normal_distribution<realT> > norm_distT;
 
   ///Alias for a poisson random variate
   typedef mx::math::randomT<long, std::mt19937_64, std::poisson_distribution<long> > poisson_distT;
 
 
protected:
 
   /* Standard WFS Interface: */
   int m_wfSz {0}; ///< Size of the wavefront in pixels
 
   int m_detRows {0}; ///<The number of rows of the WFS m_detector.  After forming the image the WFS detector plane is re-binned to this.
   int m_detCols {0}; ///<The number of columns of the WFS m_detector.  After forming the image the WFS detector plane is re-binned to this.
 
   realT m_lambda {0}; ///< Central wavelength, in meters
 
   int m_iTime {0}; ///<Integration time in loop steps
 
   int m_roTime {1}; ///<Readout time in loop steps
 
   realT m_simStep {0.001}; ///<The simulation stepsize in seconds.
 
 
   /* Direct Phase Sensor specific: */
 
   realT m_npix {0}; ///< The number of pixels assumed actually used for WFS in the S/N calc.
   realT m_Fbg {0.0}; ///< The background flux in photons/sec/pixel.
   
   int m_iTime_counter {0};
 
   int m_reading {0};
 
   int m_roTime_counter {0};
 
   bool m_firstRun {true};
 
   std::vector<wavefrontT> m_wavefronts;
 
   int m_lastWavefront;
 
   ///The image formed by the WFS
   wfsImageT<realT> m_wfsImage;
 
public:
   bool m_poissonNoise {true};
   
   ///Default c'tor
   directPhaseSensor();
 
   /* The standard WFS interface: */
 
   detectorT m_detector;
 
   ///The image on the detector, resized from wfsImage
   wfsImageT<realT> m_detectorImage;
  
   ///Get the wavefront size in pixels
   /**
     * \returns the wavefront size in pixels
     */
   int wfSz();
 
   ///Set the wavefront size in pixels.
   void wfSz( int sz /**< [in] The new size */);
 
   ///Get the detector rows  in pixels
   /**
     * \returns m_detRows
     */
   int detRows();
 
   ///Set the detector rows in pixels.
   /**
     * \param sz is the new size
     */
   void detRows(int sz);
 
   ///Get the detector Cols  in pixels
   /**
     * \returns m_detCols
     */
   int detCols();
 
   ///Set the detector columns in pixels.
   /**
     * \param sz is the new size
     */
   void detCols(int sz);
 
   ///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
     */
   realT lambda();
 
   ///Set the PyWFS central wavelength
   /**
     * \param d is the new central wavelength in meters
     */
   void lambda(realT 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.
   realT simStep();
 
   ///Set the simulation step-size, in seconds.
   void simStep(realT st);
 
   ///Get the number of pixels used for noise calculations
   /**
     * \returns the number of pixels used by the wfs, m_npix.
     */ 
   realT npix();
   
   ///Set the number of pixels used fo rnoise calculations
   /**
     */
   void npix( realT np /**< [in] the number of pixels to be used by the WFS*/);
 
   ///Get the background rate
   /**
     * \returns the background rate in photons/pixel/sec, the current value of m_Fbg.
     */ 
   realT Fbg();
   
   ///Set the background rate
   /**
     */
   void Fbg( realT bg /**< [in] the background rate in photons/pixel/sec.*/);
   
   template<typename AOSysT>
   void linkSystem(AOSysT & AOSys);
 
   ///Record a wavefront without sensing
   /** This merely inserts the wavefront in the circular buffer but does no other processing
     * Retruns 0 on success
     */
   int recordWavefront(wavefrontT & pupilPlane);
   
   ///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);
 
   void doSenseWavefront();
 
   /** \name Spatial Filtering
     * @{
     */
protected:
   bool m_applyFilter {false}; ///< Flag to control whether or not the spatial filter is applied.
   
   int m_filterWidth {0}; ///< The half-width of the filter, in Fourier domain pixels corresponding to 1/D spatial frequency units.
 
   sigproc::psdFilter<realT,2> m_filter; ///< The spatial filter class object.
   
public:
   
   /// Set the flag controlling whether or not the spatial filter is applied
   void applyFilter(bool af /**< [in] true to apply the filter, false to not apply*/); 
   
   /// Get the current value of the flag controlling whether or not the spatial filter is applied
   /**
     * \returns the current value of m_applyFilter
     */ 
   bool applyFilter();
   
   /// Set the spatial filter half-width
   /** This sets up the filter and loads it into m_filter.  The filter is a hard-edged square of width 2*width+1 pixels in the Fourier plane.
     * 
     * m_filterWidth must still be set to true after calling this function.
     */ 
   void filterWidth( int width /**< [in] the desired half-width of the filter*/);
   
   /// Get the current value of the filter half-width
   /**
     * \returns the current value of m_filterWidth
     */
   int filterWidth();
   
   /// Provides access to the filter for use by other parts of the simulation.
   /** 
     * \returns a const reference to m_filter.
     */ 
   const sigproc::psdFilter<realT, 2> & filter();
   
   ///@}
   
   /** \name Photon Noise
     * To turn on the addition of photon noise to the WFS image, you must set beta_p to be greater than 0, 
     * and must supply the pupil image.
     * 
     * @{
     */ 
protected:
   
   norm_distT m_normVar; ///< Gets normal-distributed variates
   poisson_distT m_poissonVar;
   
   
   /// The photon noise senstivity parameter.  
   /** If 0 (default) no noise is applied.  A value of 1 represents the ideal interferomter.
     * Note that this is constant vs. spatial frequency.  See Guyon (2005) \cite guyon_2005 for more information about \f$ \beta_p \f$.
     */ 
   realT m_beta_p {0};
   
   /// The pupil is needed to properly normalize poisson noise. 
   /** If null, then no noise will be applied.
     */
   pupilT * m_pupil {nullptr};
   
   /// Array used internally to calculated noise, global to avoid re-allocations.
   typename wfsImageT<realT>::imageT m_noiseIm;
   
public:
   
   /// Set the new value of the photon noise sensitivity parameter, beta_p.
   void beta_p(realT bp /**< [in] the new value of beta_p */);
   
   /// Get the current value of the photon noise sensitivity parameter.
   /**
     * \returns the current value of m_beta_p.
     */  
   realT beta_p();
   
   /// Set the pupil
   /** Sets the m_pupil pointer to the value provided.
     */
   void pupil(pupilT * pupil /**< [in] pointer to the pupil image */);
   
   /// Set the pupil
   /** Sets the m_pupil pointer to the address of the image provided.
     */
   void pupil(pupilT & pupil /**< [in] the pupil image */);
   
   /// Get the current pupil as a pointer
   /** 
     * \returns the current value of m_pupil
     */ 
   pupilT * pupil();
   
   ///@}
};
 
template<typename _realT, typename _detectorT>
directPhaseSensor<_realT, _detectorT>::directPhaseSensor()
{
   iTime(1);
 
   m_normVar.seed();
   m_poissonVar.seed();
}
 
 
template<typename _realT, typename _detectorT>
int directPhaseSensor<_realT, _detectorT>::wfSz()
{
   return m_wfSz;
}
 
template<typename _realT, typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::wfSz(int sz)
{
   if( m_wfSz == sz) return;
 
   m_wfSz = sz;
}
 
template<typename _realT, typename _detectorT>
int directPhaseSensor<_realT, _detectorT>::detRows()
{
   return m_detRows;
}
 
template<typename _realT, typename _detectorT>
int directPhaseSensor<_realT, _detectorT>::detCols()
{
   return m_detCols;
}
 
template<typename _realT,  typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::detSize(int nrows, int ncols)
{
   if( m_detRows == nrows && m_detCols == ncols) return;
 
   m_detRows = nrows;
   m_detCols = ncols;
 
   m_detector.setSize(m_detRows, m_detCols);
   m_detectorImage.image.resize(m_detRows, m_detCols);
 
 
}
 
template<typename _realT, typename _detectorT>
_realT directPhaseSensor<_realT, _detectorT>::lambda()
{
   return m_lambda;
}
 
template<typename _realT, typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::lambda(_realT l)
{
   m_lambda = l;
}
 
template<typename _realT,  typename _detectorT>
template<typename AOSysT>
void directPhaseSensor<_realT, _detectorT>::linkSystem(AOSysT & AOSys)
{
   static_cast<void>(AOSys);
}
 
template<typename _realT,  typename _detectorT>
int directPhaseSensor<_realT, _detectorT>::iTime()
{
   return m_iTime;
}
 
template<typename _realT,  typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::iTime(int it)
{
   if(it < 1)
   {
      std::cerr << "iTime must be >= 1.  Correcting\n";
      it = 1;
   }
 
   m_iTime = it;
 
   m_wavefronts.resize(m_iTime+2+100);
   m_lastWavefront = -1;
 
   m_detector.expTime(m_simStep*m_iTime);
 
}
 
template<typename _realT,  typename _detectorT>
int directPhaseSensor<_realT, _detectorT>::roTime()
{
   return roTime;
}
 
template<typename _realT,  typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::roTime(int rt)
{
   if(rt < 1)
   {
      std::cerr << "roTime must be >= 1.  Correcting\n";
      rt = 1;
   }
 
   m_roTime = rt;
 
 
}
 
template<typename _realT,  typename _detectorT>
_realT directPhaseSensor<_realT, _detectorT>::simStep()
{
   return simStep;
}
 
template<typename _realT,  typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::simStep(_realT st)
{
 
   m_simStep = st;
 
   m_detector.expTime(m_simStep*m_iTime);
 
 
}
 
template<typename realT,  typename _detectorT>
realT directPhaseSensor<realT, _detectorT>::npix()
{
   return m_npix;
}
   
template<typename realT,  typename _detectorT>
void directPhaseSensor<realT, _detectorT>::npix( realT np )
{
   m_npix = np;
}
 
template<typename realT,  typename _detectorT>
realT directPhaseSensor<realT, _detectorT>::Fbg()
{
   return m_Fbg;
}
   
template<typename realT,  typename _detectorT>
void directPhaseSensor<realT, _detectorT>::Fbg( realT bg)
{
   m_Fbg = bg;
}
 
template<typename _realT, typename _detectorT>
int directPhaseSensor<_realT, _detectorT>::recordWavefront(wavefrontT & pupilPlane)
{
   
   ++m_lastWavefront;
   if((size_t) m_lastWavefront >= m_wavefronts.size()) m_lastWavefront = 0;
 
 
   wavefrontT pPlane = pupilPlane;
 
   m_wavefronts[m_lastWavefront].amplitude = pPlane.amplitude;
   m_wavefronts[m_lastWavefront].phase = pPlane.phase;
   m_wavefronts[m_lastWavefront].iterNo = pPlane.iterNo;
 
   return 0;
}
 
 
template<typename _realT, typename _detectorT>
bool directPhaseSensor<_realT, _detectorT>::senseWavefront(wavefrontT & pupilPlane)
{
   using poisson_param_t = typename std::poisson_distribution<long>::param_type;
   
   ++m_lastWavefront;
   if((size_t) m_lastWavefront >= m_wavefronts.size()) m_lastWavefront = 0;
 
 
   wavefrontT pPlane = pupilPlane;
 
   m_wavefronts[m_lastWavefront].amplitude = pPlane.amplitude;
   m_wavefronts[m_lastWavefront].phase = pPlane.phase;
   m_wavefronts[m_lastWavefront].iterNo = pPlane.iterNo;
 
   //std::cerr << m_lastWavefront << " " << pPlane.iterNo << "\n";
   //Always skip the first one for averaging to center of iTime.
   if(m_firstRun)
   {
      m_firstRun = false;
      return false;
   }
 
   ++m_iTime_counter;
 
 
   bool rv = false;
 
   if(0)//m_reading)
   {
      ++m_roTime_counter;
 
      if(m_roTime_counter >= m_roTime)
      {
         m_detectorImage.image = m_wfsImage.image.block( 0.5*(m_wfsImage.image.rows()-1) - 0.5*(m_detectorImage.image.rows()-1), 0.5*(m_wfsImage.image.cols()-1) - 0.5*(m_detectorImage.image.cols()-1), m_detectorImage.image.rows(), m_detectorImage.image.cols());
 
         m_detectorImage.iterNo = m_wfsImage.iterNo;
 
         m_roTime_counter = 0;
         m_reading=0;
         rv = true;
      }
   }
 
   if( m_iTime_counter >= m_iTime)
   {
      doSenseWavefront();
 
      m_iTime_counter = 0;
 
      m_reading = 1;
      m_roTime_counter = 0;
 
      //Just do the read
      m_detectorImage.image = m_wfsImage.image.block( 0.5*(m_wfsImage.image.rows()-1) - 0.5*(m_detectorImage.image.rows()-1), 0.5*(m_wfsImage.image.cols()-1) - 0.5*(m_detectorImage.image.cols()-1), m_detectorImage.image.rows(), m_detectorImage.image.cols());
 
      //*** Spatial Filter:
      if(m_applyFilter)
      {
         //std::cerr << "Filtering . . . \n";
         m_filter.filter(m_detectorImage.image);
         if(m_pupil != nullptr) m_detectorImage.image *= *m_pupil; 
      }
 
      
      //*** Adding Noise:
      if(m_beta_p > 0 && m_pupil != nullptr)
      {
         m_noiseIm.resize(m_detectorImage.image.rows(), m_detectorImage.image.cols());
       
         //Calculate intensity at each pixel
         for(int c=0;c<m_noiseIm.cols();++c)
         {
            for(int r=0;r<m_noiseIm.rows();++r)
            {
               m_noiseIm(r,c) = pow(pupilPlane.amplitude(r,c),2)*(*m_pupil)(r,c) * m_detector.expTime(); 
            }
         }
         //std::cerr << "Total Phots: " << m_noiseIm.sum() << "\n";
         
         realT sqrtFbg = sqrt(m_Fbg*m_detector.expTime());
         
         //Add noise
         for(int c=0;c<m_noiseIm.cols();++c)
         {
            for(int r=0;r<m_noiseIm.rows();++r)
            {
               if(m_noiseIm(r,c) == 0) //assuming branch prediction makes this worth it
               {
                  continue;
               }
               else
               {
                  realT err = sqrt(m_noiseIm(r,c))*m_normVar;
                  if(sqrtFbg > 0) err += sqrtFbg*m_normVar;
                  if(m_detector.ron() > 0) err += m_detector.ron()*m_normVar;
                  
                  m_detectorImage.image(r,c) += m_beta_p * err/m_noiseIm(r,c); //Fractional noise in photons is noise in radians
               }
            }
         }
         
      } //if(m_beta_p > 0 && m_pupil != nullptr)
 
 
      m_detectorImage.iterNo = m_wfsImage.iterNo;
 
      m_roTime_counter = 0;
      m_reading=0;
      rv = true;
 
   }
 
   //std::cerr << "DPWFS: " << m_iTime_counter << " " << m_reading << " " << m_roTime_counter << " " << rv << "\n";
 
   return rv;
 
}
 
 
template<typename _realT,  typename _detectorT>
bool directPhaseSensor<_realT, _detectorT>::senseWavefrontCal(wavefrontT & pupilPlane)
{
 
   m_lastWavefront = 1;
 
   m_wavefronts[0].amplitude = pupilPlane.amplitude;
   m_wavefronts[0].phase = pupilPlane.phase;
 
   m_wavefronts[1].amplitude = pupilPlane.amplitude;
   m_wavefronts[1].phase = pupilPlane.phase;
 
   doSenseWavefront();
 
   m_detectorImage.image = m_wfsImage.image.block( 0.5*(m_wfsImage.image.rows()-1) - 0.5*(m_detectorImage.image.rows()-1), 0.5*(m_wfsImage.image.cols()-1) - 0.5*(m_detectorImage.image.cols()-1), m_detectorImage.image.rows(), m_detectorImage.image.cols());
 
   if(m_applyFilter)
   {
      m_filter.filter(m_detectorImage.image);
   }
 
 
 
   return true;
 
}
 
template<typename _realT,  typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::doSenseWavefront()
{
   wavefrontT pupilPlane;
 
   BREAD_CRUMB;
   
   /* Here make average wavefront for now */
   int _firstWavefront = m_lastWavefront - m_iTime;
   if(_firstWavefront < 0) _firstWavefront += m_wavefronts.size();
 
   //std::cerr << m_lastWavefront << " " << m_iTime_counter << " " << _firstWavefront << "\n";
   
   pupilPlane.amplitude = 0.5*m_wavefronts[_firstWavefront].amplitude;
   pupilPlane.phase = 0.5*m_wavefronts[_firstWavefront].phase;
 
   pupilPlane.iterNo = 0.5*m_wavefronts[_firstWavefront].iterNo;
 
 
   //std::cerr << "DPS Averaging: " << m_wavefronts[_firstWavefront].iterNo << " " ;
   BREAD_CRUMB;
 
   if(m_wavefronts[_firstWavefront].iterNo < m_iTime)
   {
      m_wfsImage.image = pupilPlane.phase;
      m_wfsImage.iterNo  = pupilPlane.iterNo;
      return;
   }
 
   for(int i=0; i < m_iTime - 1; ++i)
   {
      ++_firstWavefront;
      if( (size_t) _firstWavefront >= m_wavefronts.size()) _firstWavefront = 0;
 
      //std::cerr << m_lastWavefront << " " << m_iTime_counter << " " << _firstWavefront << "\n";
 
      pupilPlane.amplitude += m_wavefronts[_firstWavefront].amplitude;
      pupilPlane.phase += m_wavefronts[_firstWavefront].phase;
      pupilPlane.iterNo += m_wavefronts[_firstWavefront].iterNo;
 
   }
 
   ++_firstWavefront;
   if( (size_t) _firstWavefront >= m_wavefronts.size()) _firstWavefront = 0;
 
   //std::cerr << m_lastWavefront << " " << m_iTime_counter << " " << _firstWavefront << "\n";
 
   pupilPlane.amplitude += 0.5*m_wavefronts[_firstWavefront].amplitude;
   pupilPlane.phase += 0.5*m_wavefronts[_firstWavefront].phase;
   pupilPlane.iterNo += 0.5*m_wavefronts[_firstWavefront].iterNo;
      
   pupilPlane.amplitude /= (m_iTime);
   pupilPlane.phase /= (m_iTime);
   pupilPlane.iterNo /= (m_iTime);
 
   m_wfsImage.image = pupilPlane.phase;
   m_wfsImage.iterNo  = pupilPlane.iterNo;
 
}
 
template<typename _realT,  typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::applyFilter(bool af)
{
   m_applyFilter = af;
}
 
template<typename _realT,  typename _detectorT>
bool directPhaseSensor<_realT, _detectorT>::applyFilter()
{
   return m_applyFilter;
}
   
 
template<typename _realT,  typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::filterWidth( int width )
{
   int nr = m_detectorImage.image.rows();
   int nc = m_detectorImage.image.cols();
 
   typename wfsImageT<_realT>::imageT filterMask;
 
   filterMask.resize( nr, nc );
   filterMask.setZero();
 
   filterMask.block(0,0, width+1, width+1).setConstant(1.0);
   filterMask.block(0, nc - width, width+1, width).setConstant(1.0);
   filterMask.block(nr - width, 0, width, width+1).setConstant(1.0);
   filterMask.block(nr - width, nc - width, width, width).setConstant(1.0);
 
   m_filter.psdSqrt(filterMask,static_cast<realT>(1)/nr,static_cast<realT>(1)/nc);
 
   m_filterWidth = width;
 
}
 
template<typename _realT,  typename _detectorT>
int directPhaseSensor<_realT, _detectorT>::filterWidth()
{
   return m_filterWidth;
}
 
template<typename realT,  typename detectorT>
const sigproc::psdFilter<realT,2> & directPhaseSensor<realT, detectorT>::filter()
{
   return m_filter;
}
 
 
template<typename _realT,  typename _detectorT>
void directPhaseSensor<_realT, _detectorT>::beta_p(realT bp)
{
   m_beta_p = bp;
}
 
template<typename realT,  typename detectorT>
realT directPhaseSensor<realT, detectorT>::beta_p()
{
   return m_beta_p;
}
 
template<typename realT,  typename detectorT>
void directPhaseSensor<realT, detectorT>::pupil(pupilT * pupil)
{
   m_pupil = pupil;
}
 
template<typename realT,  typename detectorT>
void directPhaseSensor<realT, detectorT>::pupil(pupilT & pupil)
{
   m_pupil = &pupil;
}
   
template<typename realT,  typename detectorT>
typename directPhaseSensor<realT, detectorT>::pupilT * directPhaseSensor<realT, detectorT>::pupil()
{
   return m_pupil;
}
   
} //namespace sim
 
} //namespace AO
 
} //namespace mx
 
#endif //directPhaseSensor_hpp