ccdDetector.hpp

Go to the documentation of this file.

/** \file ccdDetector.hpp
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief Provides a class to simulate a CCD.
 * \ingroup mxAO_files
 *
 */
 
#ifndef ccdDetector_hpp
#define ccdDetector_hpp
 
// #include <Eigen/Dense>
 
#include "../../improc/imageTransforms.hpp"
#include "../../math/randomT.hpp"
 
#include "wavefront.hpp"
 
#include <fcntl.h>
#include <iostream>
 
namespace mx
{
namespace AO
{
namespace sim
{
 
/// A simulated CCD detector
/** A simulated CCD detector, including an optional EMCCD noise model.
 *
 * \ingroup mxAOSim
 */
template <typename _realT>
class ccdDetector
{
 
  public:
    typedef _realT realT;
 
    typedef wavefront<realT> wavefrontT;
 
    typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> imageT;
 
    typedef mx::math::randomT<realT, std::mt19937_64, std::normal_distribution<realT>> norm_distT;
 
    typedef mx::math::randomT<int, std::mt19937_64, std::poisson_distribution<int>> poisson_distT;
 
    typedef mx::math::randomT<realT, std::mt19937_64, std::gamma_distribution<realT>> gamma_distT;
 
    ccdDetector();
 
  protected:
    norm_distT m_normVar;       ///< Gets normal-distributed variates
    poisson_distT m_poissonVar; ///< Gets Poisson distributed variates
    gamma_distT m_gammaVar;     ///< Gets gamma distributed variates
 
    realT m_qe{ 1 }; ///< The quantum efficiency.
 
    realT m_darkCurrent{ 0 }; ///< The dark current, per pixel per second
    realT m_ron{ 0 };         ///< The readout noise, electrons per pixel per read #include "mx/randomT.hpp"
 
    realT m_cic{ 0 };  ///< EMCCD clock induced charge, electrons per pixel per read.
    realT m_gain{ 1 }; ///< Electron multiplication gain.  If >1, then EMCCD is modeled.
 
    realT m_expTime{ 1 }; ///< The exposure time, in seconds.
 
    int m_rows{ 0 }; ///< The detector size, in rows.
    int m_cols{ 0 }; ///< The detector size, in columns.
 
    bool m_noNoise{ false }; ///< If true no noise is added to the exposed image.
 
    // This is probably vestigial, commented out on 2018-07-23
    // Delete after testing with an aoSim compile.
    // long get_seed();
 
  public:
    /// Get the current value of qe.
    /**
     * \returns the current value of m_qe
     */
    realT qe();
 
    void qe( const realT &q );
 
    realT darkCurrent();
 
    void darkCurrent( const realT &dc );
 
    realT ron();
 
    void ron( const realT &r );
 
    realT cic();
 
    void cic( const realT &c );
 
    realT gain();
 
    void gain( const realT &g );
 
    realT expTime();
 
    void expTime( const realT &dt );
 
    int rows();
 
    void rows( const int &r );
 
    int cols();
 
    void cols( const int &c );
 
    void setSize( const int &r, const int &c );
 
    bool noNoise();
 
    void noNoise( bool nn );
 
    /// Rebin and add noise to the input image, placing the result in the output image.
    /** The output image must be the same size or smaller than the input image.
     * The output image is resized only if necessary.
     * The input image is multiplied by expTime, so its flux should be in photons/sec.
     * Noise is modeled as follows:
     * -# The number of input photons for pixel (i,j) is \f$ n_{ie} = (in(i,j) + darkCurrent)*expTime + cic \f$
     * -# The Poisson distribution with mean and variance \f$ n_{ie} \f$ is sampled.
     * -# If gain > 1, the gamma distribution is sampled
     * -# Read noise is applied
     * -# Result is rounded to nearest whole integer.
     */
    void exposeImage( imageT &out, ///< [out] The output image, after all above steps applied.
                      imageT &in   ///< [in] the input image, in photons/sec flux units.
    );
};
 
template <typename realT>
ccdDetector<realT>::ccdDetector()
{
    m_normVar.seed();
    m_poissonVar.seed();
    m_gammaVar.seed();
}
 
/*template<typename realT>
long ccdDetector<realT>::get_seed()
{
   long int seedval;
 
   int fd;
 
   fd = open("/dev/urandom", O_RDONLY);
 
   int rv =read(fd, &seedval, sizeof(long int));
 
    close(fd);
 
 
   if(rv < 0)
   {
      std::cerr << "read from /dev/urandom returned error\n";
 
      return 0;
   }
 
   return seedval;
}
*/
 
template <typename realT>
realT ccdDetector<realT>::qe()
{
    return m_qe;
}
 
template <typename realT>
void ccdDetector<realT>::qe( const realT &q )
{
    m_qe = q;
}
 
template <typename realT>
realT ccdDetector<realT>::darkCurrent()
{
    return m_darkCurrent;
}
 
template <typename realT>
void ccdDetector<realT>::darkCurrent( const realT &dc )
{
    m_darkCurrent = dc;
}
 
template <typename realT>
realT ccdDetector<realT>::ron()
{
    return m_ron;
}
 
template <typename realT>
void ccdDetector<realT>::ron( const realT &r )
{
    m_ron = r;
}
 
template <typename realT>
realT ccdDetector<realT>::cic()
{
    return m_cic;
}
 
template <typename realT>
void ccdDetector<realT>::cic( const realT &c )
{
    m_cic = c;
}
 
template <typename realT>
realT ccdDetector<realT>::gain()
{
    return m_gain;
}
 
template <typename realT>
void ccdDetector<realT>::gain( const realT &g )
{
    m_gain = g;
}
 
template <typename realT>
realT ccdDetector<realT>::expTime()
{
    return m_expTime;
}
 
template <typename realT>
void ccdDetector<realT>::expTime( const realT &dt )
{
    m_expTime = dt;
}
 
template <typename realT>
int ccdDetector<realT>::rows()
{
    return m_rows;
}
 
template <typename realT>
void ccdDetector<realT>::rows( const int &r )
{
    m_rows = r;
}
 
template <typename realT>
int ccdDetector<realT>::cols()
{
    return m_cols;
}
 
template <typename realT>
void ccdDetector<realT>::cols( const int &c )
{
    m_cols = c;
}
 
template <typename realT>
void ccdDetector<realT>::setSize( const int &r, const int &c )
{
    m_rows = r;
    m_cols = c;
}
 
template <typename realT>
bool ccdDetector<realT>::noNoise()
{
    return m_noNoise;
}
 
template <typename realT>
void ccdDetector<realT>::noNoise( bool nn )
{
    m_noNoise = nn;
}
 
template <typename realT>
void ccdDetector<realT>::exposeImage( imageT &out, imageT &in )
{
 
    using poisson_param_t = typename std::poisson_distribution<int>::param_type;
    using gamma_param_t = typename std::gamma_distribution<realT>::param_type;
 
    BREAD_CRUMB
 
    out.resize( m_rows, m_cols );
 
    BREAD_CRUMB
 
    improc::imageDownSample( out, in );
 
    BREAD_CRUMB
 
    if( m_noNoise )
        return;
 
    for( int j = 0; j < m_cols; ++j )
    {
        for( int i = 0; i < m_rows; ++i )
        {
            realT charge = ( out( i, j ) * m_qe + m_darkCurrent ) * m_expTime + m_cic;
 
            if( charge > 1000 )
            {
                out( i, j ) = charge + m_normVar * sqrt( charge );
            }
            else
            {
                m_poissonVar.distribution.param( poisson_param_t{ charge } );
                out( i, j ) = m_poissonVar;
            }
 
            if( m_gain > 1 )
            {
                m_gammaVar.distribution.param( gamma_param_t{ out( i, j ), m_gain } );
 
                out( i, j ) = m_gammaVar;
            }
 
            if( m_ron > 0 )
            {
                out( i, j ) += m_normVar * m_ron;
            }
 
            if( m_gain > 1 )
            {
                out( i, j ) /= m_gain;
            }
 
            // out(i,j) = round(out(i,j));
        }
    }
}
 
} // namespace sim
} // namespace AO
} // namespace mx
 
#endif // ccdDetector_hpp