Commit 62a1c45f authored by rfetick's avatar rfetick
Browse files

remove some fftshift

parent 13c5d29e
......@@ -7,12 +7,12 @@ Created on Mon May 27 17:31:18 2019
"""
import numpy as _np
from numpy.fft import fft2, ifft2, fftshift
from numpy.fft import fft2, ifft2, fftshift, ifftshift
from scipy.optimize import least_squares
from astropy.io import fits
import time
import sys
from maoppy.utils import binning as _binning
from maoppy.utils import binning
__all__ = ["oversample", "lsq_flux_bck", "psffit", "reduced_coord",
......@@ -548,7 +548,7 @@ class ParametricPSFfromPSD(ParametricPSF):
@property
def otfDiffraction(self):
return self._otfDiffraction
return fftshift(self._otfDiffraction)
def strehlOTF(self,x0):
"""Compute the Strehl ratio based on the sum of the OTF"""
......@@ -572,7 +572,7 @@ class ParametricPSFfromPSD(ParametricPSF):
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):
def otf(self, x0, dx=0, dy=0, _caller='user', grad=False, shift=True):
"""
See __call__ for input arguments
Warning: result of otf will be unconsistent if undersampled!!!
......@@ -591,10 +591,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):
def _otfTurbulent(self, x0, grad=False):
L = self.system.D * self._samp_over
if grad:
PSD,integral,g,integral_g = self.psd(x0,grad=True)
......@@ -602,16 +605,16 @@ class ParametricPSFfromPSD(ParametricPSF):
PSD,integral = self.psd(x0,grad=False)
Bphi = _np.real(fft2(fftshift(PSD))) / L**2
#Bphi = fftshift(Bphi[0, 0] - Bphi) # normalized on the numerical FoV
Bphi = fftshift(integral - Bphi) # normalized up to infinity
#Bphi = Bphi[0, 0] - Bphi # normalized on the numerical FoV
Bphi = integral - Bphi # normalized up to infinity
otf = _np.exp(-Bphi)
if grad:
g2 = _np.zeros(g.shape,dtype=complex) # I cannot override 'g' here due to float to complex type
for i in range(len(x0)):
Bphi = _np.real(fft2(fftshift(g[i,...]))) / L**2
#Bphi = fftshift(Bphi[0, 0] - Bphi) # normalized on the numerical FoV
Bphi = fftshift(integral_g[i] - Bphi) # normalized up to infinity
#Bphi = Bphi[0, 0] - Bphi # normalized on the numerical FoV
Bphi = integral_g[i] - Bphi # normalized up to infinity
g2[i,...] = -Bphi*otf
return otf, g2
return otf
......@@ -622,9 +625,9 @@ class ParametricPSFfromPSD(ParametricPSF):
NpupY = _np.ceil(ny/self._samp_over)
tab = _np.zeros((nx, ny), dtype=complex)
tab[0:int(NpupX), 0:int(NpupY)] = self.system.pupil((NpupX,NpupY),samp=self._samp_over)
self._otfDiffraction = fftshift(_np.abs(ifft2(_np.abs(fft2(tab)) ** 2)) / _np.sum(tab))
self._otfDiffraction = _np.abs(ifft2(_np.abs(fft2(tab)) ** 2)) / _np.sum(tab)
def _otfShift(self,dx,dy):
def _otfShift(self, dx, dy):
Y, X = self._fxfy/self._pix2freq # so it is just an array of pixels
ny,nx = X.shape
# Compensate oversampling shift
......@@ -633,9 +636,9 @@ class ParametricPSFfromPSD(ParametricPSF):
# Compensate odd pixel shift
dx -= (self.npix[1]%2)/2
dy -= (self.npix[0]%2)/2
return _np.exp(-2j*_np.pi*self._k*(dx*X/nx + dy*Y/ny))
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):
def __call__(self, x0, dx=0, dy=0, grad=False):
"""
Parameters
----------
......@@ -648,13 +651,13 @@ class ParametricPSFfromPSD(ParametricPSF):
"""
if grad:
otf,g = self.otf(x0,dx=dx,dy=dy,_caller='self',grad=True)
otf,g = self.otf(x0, dx=dx, dy=dy, _caller='self', grad=True, shift=False)
for i in range(len(x0)):
g[i,...] = _np.real(fftshift(ifft2(fftshift(g[i,...]))))
g[i,...] = fftshift(_np.real(ifft2(g[i,...])))
else:
otf = self.otf(x0,dx=dx,dy=dy,_caller='self',grad=False)
otf = self.otf(x0, dx=dx, dy=dy, _caller='self', grad=False, shift=False)
psf = _np.real(fftshift(ifft2(fftshift(otf))))
psf = fftshift(_np.real(ifft2(otf)))
k = int(self._k)
if k==1:
......@@ -664,9 +667,9 @@ class ParametricPSFfromPSD(ParametricPSF):
if grad:
g2 = _np.zeros((len(x0),psf.shape[0]//k,psf.shape[1]//k))
for i in range(len(x0)):
g2[i,...] = _binning(g[i,...].astype(float),k)
return _binning(psf,k), g2
return _binning(psf,k) # undersample PSF if needed (if it was oversampled for computations)
g2[i,...] = binning(g[i,...].astype(float),k)
return binning(psf,k), g2
return binning(psf,k) # undersample PSF if needed (if it was oversampled for computations)
#%% TURBULENT PSF
......
......@@ -54,6 +54,9 @@ def binning(image, k):
>>> data_bin = binning(data,2)
"""
if k==1:
return _np.copy(image)
S = _np.shape(image)
S0 = int(S[0] / k)
S1 = int(S[1] / k)
......
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