directPhaseReconstructorOrtho.hpp

#ifndef directPhaseReconstructor_hpp
#define directPhaseReconstructor_hpp
 
 
#include "../../improc/eigenImage.hpp"
#include "../../improc/eigenCube.hpp"
#include "../../improc/fitsFile.hpp"
using namespace mx::improc;
 
#include "../../sigproc/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 directPhaseReconstructorSpec
{
   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 directPhaseReconstructorSpec specT;
      
protected:
    Eigen::Array<realT,-1,-1> _recon; ///< The reconstructor matrix.
    
   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
   
   int _measurementSize; ///<The number of values in the measurement
      
   realT _calAmp; ///<The calibration amplitude used for response matrix acquisition
   
   imageT _rMat; ///<The response matrix
      
   eigenCube<realT> _rImages;
   
   //The mirror modes
   improc::eigenCube<realT> *_modes;
   
   
   
   
   imageT * _pupil;
 
   
   
   realT norm;
   
public:
   
   //The orthogonalized basis
   improc::eigenCube<realT> ortho;
   
   //The orthogonal spectrum 
   improc::eigenImage<realT> spectrum;
   
   std::vector<realT> b;
   
   imageT _mask;
   
   improc::eigenImage<double> _spectrum;
   
   imageT * _gains; 
   
   int _npix;
   
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;
   
      _mask = *_pupil;
   }
   
   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(const 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()
{
   _nModes = 0;
   
   norm = 0;
   
   _gains =0;
   
   _npix = 0;
}
 
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;
   
   _mask = *_pupil;
}
 
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(const std::string & fname)
{
}
 
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)
{
 
}
     
template<typename realT> 
template<typename measurementT, typename wfsImageT>
void directPhaseReconstructor<realT>::reconstruct(measurementT & commandVect, wfsImageT & wfsImage)
{
   
   BREAD_CRUMB;
   
   if(_npix == 0)
   {
      _npix = _pupil->sum();
      
   }
   
   BREAD_CRUMB;   
 
   wfsImage.image *= _mask;
   
   BREAD_CRUMB;
   
 
   
   b.resize(_nModes);
   
   #pragma omp parallel for
   for(int j=0; j< _nModes; ++j)
   {
      b[j] = (wfsImage.image*ortho.image(j)).sum()/ _npix;
   }
 
   //commandVect.measurement.setZero();
      
   //#pragma omp parallel for
   for(int j=0; j< modes.planes(); ++j)
   {
      realT amp = b[0]*spectrum(0,j);
      
      for(int i=1;i<modes.planes();++i) 
      {
         amp += b[i]*spectrum(i,j);
      }
      commandVect.measurement(0,j) = amp;
   }
 
   //std::cerr << "--- " << b[0] << " " << b[0]*spectrum(0,0) << " " << commandVect.measurement(0,0) << "\n";
   for(int k=0; k< _nModes; ++k)
   {
      std::cerr << commandVect.measurement(0,k) << " ";
   }
   std::cerr << "\n";
   
   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)
{
   int l = 0;
   for(int j=0; j<measureVec.measurement.rows(); ++j)
   {
      for(int k=0; k<measureVec.measurement.cols(); ++k)
      {
         _rMat(l, i) = measureVec.measurement(j,k);
         ++l;
      }
   }
   
  //_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)
{
   accumulateRMat(i, measureVec);
   _rImages.image(i) = wfsImage.image;   
}
 
template<typename realT> 
void directPhaseReconstructor<realT>::saveRMat(std::string fname)
{
   fitsFile<realT> ff;
   fitsHeader head;
   
   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, head);
}
 
template<typename realT> 
void directPhaseReconstructor<realT>::saveRImages(std::string fname)
{
   fitsFile<realT> ff;
   fitsHeader head;
   
   
   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 //directPhaseReconstructor_hpp
 
 
//--- Code for mean and tip/tilt subtraction in reconstruct:
//    realT mean = (wfsImage.image * *_pupil).sum()/npix;
//    
//    wfsImage.image -= mean;
//    wfsImage.image *= *_pupil;
 
#if 0 
   
   BREAD_CRUMB;
   if(norm == 0)
   {
      for(int i = 0; i < wfsImage.rows(); ++i)
      {
         for(int j=0; j < wfsImage.cols(); ++j)
         {
            norm += pow((i - 0.5*(wfsImage.rows()-1.0))*(*_pupil)(i,j),2);
         }
      }
    BREAD_CRUMB;  
      norm = sqrt(norm);
      
      //std::cout << "NORM: === " << norm << "\n";
   }
   BREAD_CRUMB;
 
   //std::cout << "NORM: === " << norm << "\n";
   
   //Remove tip/tilt:
   realT ampx = 0, ampy = 0;
   BREAD_CRUMB;   
   for(int i = 0; i < wfsImage.rows(); ++i)
   {
      for(int j=0; j < wfsImage.cols(); ++j)
      {
         ampx += wfsImage(i,j) * (i - 0.5*(wfsImage.rows()-1.0))/norm;
         ampy += wfsImage(i,j) * (j - 0.5*(wfsImage.cols()-1.0))/norm;
      }
   }
   BREAD_CRUMB;
   for(int i = 0; i < wfsImage.rows(); ++i)
   {
      for(int j=0; j < wfsImage.cols(); ++j)
      {
         wfsImage(i,j) -= ampx * (i - 0.5*(wfsImage.rows()-1.0))/norm * (*_pupil)(i,j);
         wfsImage(i,j) -= ampy * (j - 0.5*(wfsImage.cols()-1.0))/norm * (*_pupil)(i,j);
      }
   }
#endif