% -*- mode: Noweb; noweb-code-mode: python-mode -*-

\section{PXD file}
\label{sec:pxd-file}


<<imaging.pxd>>=
from source cimport source, Source
from utilities cimport Sensor, cuFloatArray, cuFloatComplexArray
cdef extern from "cuda.h":
    ctypedef struct float2:
        float x
        float y
cdef extern from "imaging.h":
    cdef cppclass imaging:
        int N_PX_PUPIL, N_DFT, N_SIDE_LENSLET, N_PX_CAMERA, N_SOURCE, N_FRAME
        float photoelectron_gain, N_PHOTON_PER_FRAME, N_PHOTON_TOTAL, pixel_scale
        float *d__frame
        float2 *d__wave_PUPIL
        char absolute_pointing
        void setup(int , int , int , int , int, int)
        void cleanup()
        void set_pointing_direction(float *, float *)
        void propagate(source *)
        void propagate_cpx(source *)
        void propagateThroughFieldStop(source *, float)
        void propagateThroughPyramid(source *, float)
        void propagateThroughModulatedPyramid(source *, float, int, float)
        void readout(float, float)
        void readout(float, float, float, float)
        void noiseless_readout(float)
        void reset()
        void reset_rng(int)
<<class definitions>>
@
\subsection{Class definitions}
\label{sec:class-definitions}

<<class definitions>>=
cdef class Imaging(Sensor):
    cdef:
        imaging *_c_imaging
        readonly int N_SIDE_LENSLET, N_PX_PUPIL, DFT_osf, \
                     N_PX_IMAGE, BIN_IMAGE, N_SOURCE, \
                     N_PX_FRAME, ghost_N_FRAME
        readonly cuFloatArray frame
        readonly cuFloatComplexArray cpx
        char isPtr
        public float exposureTime, readOutNoiseRms, nBackgroundPhoton, noiseFactor
        public float[:,:] ghost_frame
    cdef void set_ptr(Imaging, imaging *)

@
\section{PYX file}
\label{sec:pyx-file}

\index{imaging!python!Imaging}
<<imaging.pyx>>=
from scipy.optimize import bisect
import numpy as np
cdef class Imaging(Sensor):
    """
    Creates an imaging object.

    Parameters
    ----------
    N_PX_PUPIL : int
       The sampling in pixel of one lenslet.
    N_SIDE_LENSLET : int
       The linear size of the lenslet array (>=1).
    DFT_osf : int
       The oversampling factor for the Discrete Fourier Transform (>=1).
    N_PX_IMAGE : int
       The sampling in pixel of the imagelet (>=1).
    BIN_IMAGE : int
       The binning factor of the imagelet (>=1).
    N_SOURCE : int
       The number of guide stars (>=1).
    photoelectron_gain : float
        The photon to electron conversion gain of the camera; default: 1
    exposureTime : float
       The detector exposure time in second; default: 1s
    readOutNoiseRms : float
       The detector read-out noise rms; default: 0
    nBackgroundPhoton : float
       The number of background photon; default: 0
    noiseFactor : float
       The excess noise factor; default: 1

    Attributes
    ----------
    N_SIDE_LENSLET : int
       The linear size of the lenslet array (>=1).
    N_PX_PUPIL : int
       The sampling in pixel of one lenslet.
    DFT_osf : int
       The oversampling factor for the Discrete Fourier Transform
    N_PX_IMAGE : int
       The sampling in pixel of the imagelet
    BIN_IMAGE : int
       The binning factor of the imagelet
    N_GS : int
       The number of guide stars
    N_PX_FRAME : int
        The detector resolution
    frame : cuFloatArray
        The detector frame of size [`N_SOURCE`x`N_PX_CAMERA`,`N_PX_CAMERA`]
    photoelectron_gain : float
        The photon to electron conversion gain of the camera; default: 1
    exposureTime : float
       The detector exposure time in second; default: 1s
    readOutNoiseRms : float
       The detector read-out noise rms; default: 0
    nBackgroundPhoton : float
       The number of background photon; default: 0
    noiseFactor : float
       The excess noise factor; default: 1

    See also
    --------
    cuFloatArray : class acting as the interface between GPU host and device

    """

    def __cinit__(self, *args, **kwargs):
        self._c_imaging = new imaging()

    def __init__(self, int N_SIDE_LENSLET=1, int N_PX_PUPIL=128,
		 int DFT_osf=2, int N_PX_IMAGE=0, 
                 int BIN_IMAGE=1, int N_SOURCE=1,
                 float exposureTime=1.0, float readOutNoiseRms=0.0, 
                 float nBackgroundPhoton=0.0, float noiseFactor=1.0,
                 float photoElectronGain=1.0,
		 char isPtr=0, **kwargs):
        if N_PX_IMAGE==0:
            N_PX_IMAGE = N_PX_PUPIL;
        self.N_SIDE_LENSLET = N_SIDE_LENSLET
        self.N_PX_PUPIL     = N_PX_PUPIL
        self.DFT_osf        = DFT_osf
        self.N_PX_IMAGE     = N_PX_IMAGE
        self.BIN_IMAGE      = BIN_IMAGE
        self.N_PX_FRAME     = self.N_PX_IMAGE*self.N_SIDE_LENSLET/self.BIN_IMAGE
        self.N_SOURCE       = N_SOURCE
        """
        if wfs is not None:
            self.isPtr = 1
            self._c_imaging = &(wfs._c_shackHartmann.camera)
        elif sps is not None:
            self.isPtr = 1
            self._c_imaging = &(sps._c_segmentPistonSensor.camera)
	elif pym is not None:
            self.isPtr = 1
            self._c_imaging = &(pym._c_pyramid.camera)
        else:
            self.isPtr = 0
            if N_PX_IMAGE is None:
                N_PX_IMAGE = N_PX_PUPIL
            self._c_imaging.setup(N_PX_PUPIL, N_SIDE_LENSLET,
                                  DFT_osf, N_PX_IMAGE,
                                  BIN_IMAGE, N_SOURCE)
        """
        self.isPtr = isPtr
        self.frame = cuFloatArray(shape=(self.N_SOURCE*self.N_PX_FRAME,
				         self.N_PX_FRAME))
        self.cpx = cuFloatComplexArray(shape=(self.DFT_osf*(self.N_PX_PUPIL+1),
				                              self.DFT_osf*(self.N_PX_PUPIL+1)))
        if not self.isPtr:
            if N_PX_IMAGE is None:
                N_PX_IMAGE = N_PX_PUPIL
            self._c_imaging.setup(N_PX_PUPIL, N_SIDE_LENSLET,
                                  DFT_osf, N_PX_IMAGE,
                                  BIN_IMAGE, N_SOURCE)
            self.frame._c_gpu.dev_data = self._c_imaging.d__frame
            self.cpx._c_gpu.dev_data = self._c_imaging.d__wave_PUPIL
            self._c_imaging.photoelectron_gain = photoElectronGain

        self.exposureTime      = exposureTime
        self.readOutNoiseRms   = readOutNoiseRms
        self.nBackgroundPhoton = nBackgroundPhoton
        self.noiseFactor       = noiseFactor

    cdef void set_ptr(self, imaging *camera_ptr):
        self._c_imaging = camera_ptr
        self.frame._c_gpu.dev_data = self._c_imaging.d__frame
        
    def __dealloc__(self):
        if not self.isPtr:
            self._c_imaging.cleanup()

    def __invert__(self):
        """
        Reset the sensor
        """
        self.reset()

    def __pos__(self):
        """
        Read-out and reset the detector
        """
        self._c_imaging.readout(self.exposureTime,
                                self.readOutNoiseRms,
                                self.nBackgroundPhoton,
                                self.noiseFactor)
        self.reset()

    def calibrate(self, Source src, args=None):
        pass

    def process(self):
        pass

    def propagate(self, Source src):
        """
        Propgates a source through the lenset to detector in the focal plane

        Parameters
        ----------
        src : Source
            A source object

        See also
        --------
        Source : class modeling star objects
        """		       
        self._c_imaging.propagate(src._c_source)

    def propagate_cpx(self, Source src):
        self._c_imaging.propagate_cpx(src._c_source)

    def propagateThroughFieldStop(self, Source src, float field_stop_diam):
        """
        Propgates a source through the lenset then through a field stop in 
        the focal plane and to a pupil plane

        Parameters
        ----------
        src : Source
            A source object
        field_stop_diam : float
            The diameter of the field stop in units of :math:`\lambda/D`

        See also
        --------
        Source : class modeling star objects
        """		       
        self._c_imaging.propagateThroughFieldStop(src._c_source, field_stop_diam)

    def propagateThroughPyramid(self, Source src,
                                float modulation=0.0,
                                int modulation_sampling=0,
                                float alpha=0.5):
        """
        Propgates a source through a pyramid wavefront sensor to a detector in a pupil plane

        Parameters
        ----------
        src : Source
            A source object
        modulation : float
            The pyramid modulation amplitude in units of :math:`\lambda/D`; default: 0.0

        See also
        --------
        Source : class modeling star objects
        """		       
        if modulation>0.0:
            self._c_imaging.propagateThroughModulatedPyramid(src._c_source,
                                                             modulation,
                                                             modulation_sampling,
                                                             alpha)
        else:
            self._c_imaging.propagateThroughPyramid(src._c_source,alpha)

    def readOut(self,float exposureTime, float readOutNoiseRms, 
                float nBackgroundPhoton=0.0, float noiseFactor=1.0):
        """
        Reads-out the detector

        Parameters
        ----------
        exposureTime : float
            The exposure time of the camera
        readOutNoiseRms : float
            The read-out moise rms of the camera
        nBackgroundPhoton : float (optional)
            The number of background photon per second; default: 0
        noiseFactor : float (optional)
            The excess noise factor for EMCCD ($\sqrt{2}$); default: 1
        """
        self._c_imaging.readout(exposureTime,readOutNoiseRms,
                                nBackgroundPhoton, noiseFactor)

    def noiselessReadOut(self,float exposureTime):
        """
        Reads-out the detector

        Parameters
        ----------
        exposureTime : float
            The exposure time of the camera
        """
        self._c_imaging.noiseless_readout(exposureTime)

    def reset(self):
        """
        Resets the frame of the camera to 0
        """
        self.ghost_frame = self.frame.host()
        self.ghost_N_FRAME = self._c_imaging.N_FRAME
        self._c_imaging.reset()

    def reset_rng(self, int SEED):
        """
        Resets the random number generator
        """
        self._c_imaging.reset_rng(SEED)

    def pixelScale(self,Source src):
        """
        Returns the pixel scale in radian
        """
        self._c_imaging.pixel_scale = (src.wavelength/src.rays.L)*self.N_SIDE_LENSLET*self.BIN_IMAGE/self.DFT_osf
        return self._c_imaging.pixel_scale

    def pixelScaleArcsec(self,Source src):
        """
        Returns the pixel scale in arcsecond
        """
        return self.pixelScale(src)*180.0*3600.0/np.pi

    def pointing(self, float zen, float azi):
        """
        Points the imager away from the source direction

        Parameters
        ----------
        zen : float
            The zenith angle offset with respect to the source zenith angle
        azi : float
            The azimuth angle offset with respect to the source azimuth angle
        """
        assert self._c_imaging.pixel_scale>0, "Imager pixel scale is not set, set it with <Imaging object>.pixelScale(<Source object>)"
        self._c_imaging.set_pointing_direction(&zen,&azi)

    def ee80(self,double units=1.0, bint from_ghost=0):
        """
        Computes the 80% ensquared energy patch size
        
        Parameters
        ----------
        units : double
            The units to convert the 80% ensquared energy to

        Returns
        -------
        delta : double
            The 80% ensquared energy patch size in pixel (default) or units if given
        """
        cdef:
            int k, n, m, N_FRAME
            double[::1] q 
            double[:,::1] x, y, gate
            double[:] u, v
        if from_ghost:
            psf = np.asarray( self.ghost_frame )
            N_FRAME = self.ghost_N_FRAME
        else:
            psf = self.frame.host()  
            N_FRAME = self._c_imaging.N_FRAME
        n,m   = psf.shape
        n /= self.N_SOURCE
        #u     = np.linspace(-1.0,1.0,n)*(n/2.0)
        #v     = np.linspace(-1.0,1.0,m)*(m/2.0)
        u     = np.arange(n,dtype=np.float64)-(n-1)*0.5
        v     = np.arange(m,dtype=np.float64)-(m-1)*0.5
        x,y   = np.meshgrid(u,v)
        psf = psf.reshape(self.N_SOURCE,n,m)
        def ee80_fun(float ee_Delta, _psf_):
            _ee_Delta_ = ee_Delta/2.0
            #gate = np.logical_and(np.abs(x)<=_ee_Delta_,np.abs(y)<=_ee_Delta_)
            gate = np.hypot(x,y)<=_ee_Delta_
            return np.sum(_psf_*gate)/(self._c_imaging.N_PHOTON_PER_FRAME*\
                                       N_FRAME/self.N_SOURCE) - 0.8
        q = np.zeros(self.N_SOURCE)
        for k in range(self.N_SOURCE):
            try:
                q[k] = bisect(ee80_fun,3,n,args=(psf[k,:,:]))
            except ValueError:
                q[k] = np.float('inf')
        return np.asarray(q)*units


    def ee50(self,double units=1.0, bint from_ghost=0):
        """
        Computes the 50% ensquared energy patch size
        
        Parameters
        ----------
        units : double
            The units to convert the 80% ensquared energy to

        Returns
        -------
        delta : double
            The 80% ensquared energy patch size in pixel (default) or units if given
        """
        cdef:
            int k, n, m, N_FRAME
            double[::1] q 
            double[:,::1] x, y, gate
            double[:] u, v
        if from_ghost:
            psf = np.asarray( self.ghost_frame )
            N_FRAME = self.ghost_N_FRAME
        else:
            psf = self.frame.host()  
            N_FRAME = self._c_imaging.N_FRAME
        n,m   = psf.shape
        n /= self.N_SOURCE
        #u     = np.linspace(-1.0,1.0,n)*(n/2.0)
        #v     = np.linspace(-1.0,1.0,m)*(m/2.0)
        u     = np.arange(n,dtype=np.float64)-(n-1)*0.5
        v     = np.arange(m,dtype=np.float64)-(m-1)*0.5
        x,y   = np.meshgrid(u,v)
        psf = psf.reshape(self.N_SOURCE,n,m)
        def ee50_fun(float ee_Delta, _psf_):
            _ee_Delta_ = ee_Delta/2.0
            #gate = np.logical_and(np.abs(x)<=_ee_Delta_,np.abs(y)<=_ee_Delta_)
            gate = np.hypot(x,y)<=_ee_Delta_
            return np.sum(_psf_*gate)/(self._c_imaging.N_PHOTON_PER_FRAME*\
                                       N_FRAME/self.N_SOURCE) - 0.5
        q = np.zeros(self.N_SOURCE)
        for k in range(self.N_SOURCE):
            try:
                q[k] = bisect(ee50_fun,3,n,args=(psf[k,:,:]))
            except ValueError:
                q[k] = np.float('inf')
        return np.asarray(q)*units

    @property
    def pixel_scale(self):
        return self._c_imaging.pixel_scale

    @property
    def photoelectron_gain(self):
        return self._c_imaging.photoelectron_gain
    @photoelectron_gain.setter
    def photoelectron_gain(self,float value):
        self._c_imaging.photoelectron_gain = value

    @property
    def N_FRAME(self):
        return self._c_imaging.N_FRAME

    @property
    def N_PHOTON_PER_FRAME(self):
        return self._c_imaging.N_PHOTON_PER_FRAME

    @property
    def absolute_pointing(self):
        return self._c_imaging.absolute_pointing
    @absolute_pointing.setter
    def absolute_pointing(self, int value):
        self._c_imaging.absolute_pointing = value
@
\index{imaging!python!JImaging}
<<imaging.pyx>>=
from utilities import JSONAbstract
class JImaging(JSONAbstract,Imaging):
    """
    """
    def __init__(self, jprms = None, jsonfile = None):
        JSONAbstract.__init__(self,jprms=jprms, jsonfile=jsonfile)
	Imaging.__init__(self,1,self.jprms["pupil sampling"]-1,
			 DFT_osf    = 2*self.jprms["nyquist factor"],
			 N_PX_IMAGE = self.jprms["resolution"],
			 N_SOURCE   = self.jprms["guide star #"])