Skip to content
GitLab
Commit 32e70419 authored by rfetick's avatar rfetick
Browse files

add a private '_otf' method to ParametricPSFfromPSD

parent 62a1c45f
Hide whitespace changes
Inline Side-by-side
maoppy/psfmodel.py
View file @ 32e70419
......@@ -469,9 +469,11 @@ class ParametricPSFfromPSD(ParametricPSF):
Attributes
----------
system : maoppy.Instrument
Do not touch after __init__.
Do not modify `system` after __init__.
samp : float
Required sampling (samp=2 for Shannon).
npix : (float,float)
Number of desired pixels for output arrays.
otfDiffraction: np.array
Beware that the array is numerically oversampled when samp<2. Do not use when samp<2 or be sure of what you do.
......@@ -484,7 +486,8 @@ class ParametricPSFfromPSD(ParametricPSF):
__call__ : return the PSF
tofits : write a FITS file of the PSF
"""
def __init__(self,nparam,npix,system=None,samp=None, fixed_k = None):
def __init__(self, nparam, npix, system=None, samp=None, fixed_k=None):
if not (type(npix) in [tuple,list,_np.ndarray]): raise ValueError("npix must be a tuple, list or numpy.ndarray")
if len(npix)!=2: raise ValueError("npix must be of length = 2")
if (npix[0]%2) or (npix[1]%2): raise ValueError("Each npix component must be even")
......@@ -502,16 +505,19 @@ class ParametricPSFfromPSD(ParametricPSF):
def npix(self):
return self._npix
@npix.setter
def npix(self,value):
def npix(self, value):
self._npix = value
self._computeXYarray()
self._computeOtfDiffraction()
@property
def samp(self):
return self._samp_over/self._k
@samp.setter
def samp(self, value):
# Manage cases of undersampling
......@@ -519,6 +525,7 @@ class ParametricPSFfromPSD(ParametricPSF):
self._computeXYarray()
self._computeOtfDiffraction()
def _computeXYarray(self):
nx,ny = self._nxnyOver
xyarray = _np.mgrid[0:nx, 0:ny].astype(float)
......@@ -527,43 +534,52 @@ class ParametricPSFfromPSD(ParametricPSF):
self._fxfy = xyarray * self._pix2freq
self._f2 = self._fxfy[0]**2. + self._fxfy[1]**2.
@property
def _pix2freq(self):
"""
pixel to frequency conversion in the PSD plane
"""
"""Pixel to frequency conversion in the PSD plane"""
return 1.0/(self.system.D*self._samp_over)
@property
def _nxnyOver(self):
"""
Return the number of pixels for the correctly sampled array (at least at Shannon-Nyquist)
"""
"""Return the number of pixels for the correctly sampled array (at least at Shannon-Nyquist)"""
return self.npix[0]*self._k, self.npix[1]*self._k
@property
def _nullFreqIndex(self):
nx,ny = self._nxnyOver
return nx//2, ny//2
@property
def otfDiffraction(self):
"""Return the diffraction OTF"""
if (self._k > 1):
raise ValueError("Cannot call `otfDiffraction(...)` when undersampled")
return fftshift(self._otfDiffraction)
def strehlOTF(self,x0):
def strehlOTF(self, x0):
"""Compute the Strehl ratio based on the sum of the OTF"""
otf = self.otf(x0,_caller='self')
return _np.real(_np.sum(otf)/_np.sum(self._otfDiffraction))
return _np.real(_np.sum(self._otf(x0))/_np.sum(self._otfDiffraction))
def strehlMarechal(self,x0):
"""Compute the Strehl ratio based on the Maréchal approximation"""
def strehlMarechal(self, x0):
"""
Compute the Strehl ratio based on the Maréchal approximation.
Note that `strehlOTF` might provide more accurate results.
"""
_,sig2 = self.psd(x0)
return _np.exp(-sig2)
def psd(self,x0,grad=False):
def psd(self, x0, grad=False):
raise ValueError("ParametricPSFfromPSD is not to be instantiated. the `psd` method must be override in the subclasses")
def check_parameters(self,x0):
def check_parameters(self, x0):
bd, bu = self.bounds
x0 = _np.array(x0)
bd = _np.array(bd)
......@@ -571,18 +587,31 @@ class ParametricPSFfromPSD(ParametricPSF):
if len(x0)!=self._nparam: raise ValueError('len(x0) is different from length of bounds')
if _np.any(x0<bd): raise ValueError('Lower bounds are not respected')
if _np.any(x0>bu): raise ValueError('Upper bounds are not respected')
def otf(self, x0, dx=0, dy=0, _caller='user', grad=False, shift=True):
def otf(self, x0, dx=0, dy=0, grad=False):
"""
See __call__ for input arguments
Warning: result of otf will be unconsistent if undersampled!!!
This issue is solved with oversampling + binning in __call__ but not here
For the moment, the `_caller` keyword prevents user to misuse otf
Public implemetation of the OTF.
Only available if samp>=2.
Null frequency is at the center of the array.
"""
if (self._k > 1):
raise ValueError("Cannot call `otf(...)` when undersampled")
if grad:
otf,g = self._otf(self, x0, dx=dx, dy=dy, grad=True)
return fftshift(otf), fftshift(g,axes=(1,2))
else:
otf = self._otf(self, x0, dx=dx, dy=dy, grad=False)
return fftshift(otf)
def _otf(self, x0, dx=0, dy=0, grad=False):
"""
Private implementation of the OTF.
It is correctly oversampled.
Null frequency is at position [0,0].
Computation should be slighly quicker than self.otf(...)
"""
if (self._k > 1) and (_caller != 'self'):
raise ValueError("Cannot call `otf(...)` when undersampled (functionality not implemented yet)")
otfShift = self._otfShift(dx,dy)
if grad:
otfTurb,g = self._otfTurbulent(x0,grad=True)
......@@ -591,13 +620,13 @@ class ParametricPSFfromPSD(ParametricPSF):
otfTurb = self._otfTurbulent(x0,grad=False)
otf = otfTurb * self._otfDiffraction * otfShift
if shift:
otf = fftshift(otf)
g = fftshift(g,axes=(1,2))
if grad: return otf,g
return otf
def _otfTurbulent(self, x0, grad=False):
"""Atmospheric part of the OTF"""
L = self.system.D * self._samp_over
if grad:
PSD,integral,g,integral_g = self.psd(x0,grad=True)
......@@ -619,7 +648,9 @@ class ParametricPSFfromPSD(ParametricPSF):
return otf, g2
return otf
def _computeOtfDiffraction(self):
"""Precompute the diffraction part of the OTF"""
nx, ny = self._nxnyOver
NpupX = _np.ceil(nx/self._samp_over)
NpupY = _np.ceil(ny/self._samp_over)
......@@ -627,7 +658,9 @@ class ParametricPSFfromPSD(ParametricPSF):
tab[0:int(NpupX), 0:int(NpupY)] = self.system.pupil((NpupX,NpupY),samp=self._samp_over)
self._otfDiffraction = _np.abs(ifft2(_np.abs(fft2(tab)) ** 2)) / _np.sum(tab)
def _otfShift(self, dx, dy):
"""Shifting part of the OTF"""
Y, X = self._fxfy/self._pix2freq # so it is just an array of pixels
ny,nx = X.shape
# Compensate oversampling shift
......@@ -638,24 +671,30 @@ class ParametricPSFfromPSD(ParametricPSF):
dy -= (self.npix[0]%2)/2
return ifftshift(_np.exp(-2j*_np.pi*self._k*(dx*X/nx + dy*Y/ny)))
def __call__(self, x0, dx=0, dy=0, grad=False):
"""
Parameters
----------
x0 : numpy.array (dim=1), tuple, list
Array of parameters for the PSD (see __doc__)
dx : float
PSF X shifting [pix] (default = 0)
dy : float
PSF Y shifting [pix] (default = 0)
Keywords
--------
dx : float (default = 0)
PSF X shifting [pix].
dy : float (default = 0)
PSF Y shifting [pix].
grad : boolean (default = False)
Also compute gradients.
"""
if grad:
otf,g = self.otf(x0, dx=dx, dy=dy, _caller='self', grad=True, shift=False)
otf,g = self._otf(x0, dx=dx, dy=dy, grad=True)
for i in range(len(x0)):
g[i,...] = fftshift(_np.real(ifft2(g[i,...])))
else:
otf = self.otf(x0, dx=dx, dy=dy, _caller='self', grad=False, shift=False)
otf = self._otf(x0, dx=dx, dy=dy, grad=False)
psf = fftshift(_np.real(ifft2(otf)))
......@@ -686,7 +725,7 @@ class Turbulent(ParametricPSFfromPSD):
x[0] - Fried parameter r0 [m]
x[1] - Von Karman external length [m]
"""
def __init__(self,npix,system=None,samp=None):
def __init__(self, npix, system=None, samp=None):
"""
Parameters
----------
......@@ -704,7 +743,7 @@ class Turbulent(ParametricPSFfromPSD):
bounds_up = [_np.inf, _np.inf]
self.bounds = (bounds_down,bounds_up)
def psd(self,x0,grad=False):
def psd(self, x0, grad=False):
"""Compute the PSD model from parameters
PSD is given in [rad²/f²] = [rad² m²]
......@@ -771,7 +810,8 @@ class Psfao(ParametricPSFfromPSD):
---------
Fétick et al., August 2019, A&A, Vol.628
"""
def __init__(self,npix,system=None,Lext=10.,samp=None, fixed_k = None):
def __init__(self, npix, system=None, Lext=10., samp=None, fixed_k=None):
"""
Parameters
----------
......@@ -802,16 +842,19 @@ class Psfao(ParametricPSFfromPSD):
bounds_up = [_np.inf]*4 + [1e2,_np.inf,5]
self.bounds = (bounds_down,bounds_up)
@property
def cutoffAOfreq(self):
return self.system.Nact/(2.0*self.system.D)
def _computeXYarray(self):
super()._computeXYarray()
self._maskin = (self._f2 < self.cutoffAOfreq**2.)
self._vk = (1-self._maskin) * 0.0229 * ((1. / self.Lext**2.) + self._f2)**(-11./6.)
def psd(self,x0,grad=False):
def psd(self, x0, grad=False):
"""Compute the PSD model from parameters
PSD is given in [rad²/f²] = [rad² m²]
......@@ -893,7 +936,8 @@ class Psfao(ParametricPSFfromPSD):
if grad: return PSD,integral,g,integral_g
return PSD, integral
def tofits(self,param,filename,*args,keys=None,overwrite=False,**kwargs):
def tofits(self, param, filename, *args, keys=None, overwrite=False, **kwargs):
if keys is None:
keys = ["R0","CST","SIGMA2","ALPHA","RATIO","THETA","BETA"]
keys_comment = ["Fried parameter [m]",
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment