PD on patches added to the main code

Author: kahil_dell

PD.py

@@ -1,18 +1,10 @@
 
-
-
 import numpy as np
-import matplotlib.pyplot as plt
 import scipy
 from scipy.fftpack import fftshift, ifftshift, fft2, ifft2
-from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
-from scipy.signal import correlate2d as correlate
-from scipy.signal import general_gaussian
-from astropy.io import fits
-from scipy import ndimage
 from functools import partial
-import time
-import imreg_dft
+
+
 import pyfits
 
 

README.md

@@ -46,13 +46,23 @@ The modules of this class:
   2. `plot_results`: takes in the returned Zernike coefficients and plots the wavefront error, the 2D MTF, the 1D MTF
   3. `restored_scene`: perform the deconvolution of a blurred image with an input of Zernike Coefficients. The user can use between two filter `Wiener` and `richardsonLucy`
 
+## The `patch_pd` class:
+For fitting the wavefront error is sub-regions of the full FOV of the PD dataset. The user can enter the size of the subregion (it has to be quadratic) plus the number of Zernike polynomials to be fit and the names of the output files (one for the wavefront error and one for the 2D MTF).
 # How to use the code?
-The code returns the best-fit zernike polynomials, a visualisation of the results (wavefront error+MTF), and the restored scene (in this case the focused image of the PD dataset). To get these specific results, type in the shell terminal:
+The `minimization` class returns the best-fit zernike polynomials, a visualisation of the results (wavefront error+MTF), and the restored scene (in this case the focused image of the PD dataset). To get these specific results, type in the shell terminal:
 ```
 python3 minimization.py -i 'path/input.fits' -s 150 -w 617.3e-6 -a 140 -f 4125.3 -p 0.5 -c 0.5 -r 1e-10 -ap 10 -x1 500 -x2 650 -y1 500 -y2 650 -z 10 -del 0.5 -o path/reduced.fits -fl 'Wiener'
 ```
 
 The specific description of the parsers can be found inside [the main code](https://github.com/fakahil/PyPD/blob/master/minimization.py). The values given above are for an example PD dataset taken by the PHI/HRT telescope. You can change the values according to your telescope.
 
+To use the `patch_pd` class:
+
+```
+python3 patch_pd.py -i 'path/input.fits' -z 10 -d 265 -ow 'path/output_wf.fits' -om 'path/output_wf.fits'
+
+```
+The parsers description can be found in the [main code](https://github.com/fakahil/PyPD/blob/master/patch_pd.py)
+
 
 

__pycache__/PD.cpython-38.pyc

@@ -1,17 +1,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-import scipy
-from scipy.fftpack import fftshift, ifftshift, fft2, ifft2
-from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
-from scipy.signal import correlate2d as correlate
-from scipy.signal import general_gaussian
-from astropy.io import fits
-from scipy import ndimage
-from functools import partial
-import time
-import imreg_dft
-import pyfits
+
 
 
 def mask_pupil(rpupil, size):

__pycache__/aperture.cpython-38.pyc

@@ -2,34 +2,30 @@
 import matplotlib.pyplot as plt
 import scipy
 from scipy.fftpack import fftshift, ifftshift, fft2, ifft2
-from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
-from scipy.signal import correlate2d as correlate
-from scipy.signal import general_gaussian
-from astropy.io import fits
-from scipy import ndimage
-from functools import partial
-#import time
-import imreg_dft
-import pyfits
+
+
 import tools
 import aperture
 import PD
 import noise
 import wavefront
+import telescope
+from telescope import *
+from noise import *
+from wavefront import *
+from tools import *
+from PD import *
 
 
-class fitting:
-   def __init__(self,size,cut_off,reg,ap,x1,x2,y1,y2,co_num,del_z,lam, diameter,focal_length,platescale):
+class cost_func:
+   def __init__(self,size,cut_off,reg,ap,co_num,del_z,lam, diameter,focal_length,platescale):
 
       self.cut_off = cut_off
       self.reg = reg
       self.size = size
 
       self.co_num = co_num
-      self.x1 = x1
-      self.x2 = x2
-      self.y1 = y1
-      self.y2 = y2
+    
       self.del_z = del_z
       self.lam = lam
       self.diameter = diameter
@@ -39,21 +35,20 @@ def __init__(self,size,cut_off,reg,ap,x1,x2,y1,y2,co_num,del_z,lam, diameter,foc
 
 
 
-  def Minimise(coefficients):
-	      A_f = wavefront.pupil_foc(coefficients,self.size,rpupil,self.co_num)
-	      A_def = wavefront.pupil_defocus(coefficients,size,del_z,rpupil,co_num)
 
-	      psf_foc = wavefront.PSF(Mask,A_f,False) 
-	      psf_defoc = wavefront.PSF(Mask,A_def,False)
-	      t0 = wavefront.OTF(psf_foc)
-	      tk = wavefront.OTF(psf_defoc)
 
-	      q,q2 = PD.Q_matrix(t0,tk,reg,gam)
-	  
-	 #noise_filter = sch_filter(N0, t0, tk, d0, dk,reg)
-	      F_m = PD.F_M(q2,d0, dk,t0,tk,noise_filter,gam)
+   def Minimize_res(self,coefficients):
+
+       noise_filter = fftshift(noise_mask_high(self.size,self.cut_off))
+       
+       ph = wavefront.phase_embed(coefficients,self.telescope.pupil_size(),self.size,self.co_num)
+       A_f = wavefront.pupil_foc(coefficients,self.size,self.telescope.pupil_size(),self.co_num)
+
+
+       psf_foc = wavefront.PSF( aperture.mask_pupil(self.telescope.pupil_size(),self.size),A_f,False) 
 
-	      E_metric = PD.Error_metric(t0,tk,d0,dk,q,noise_filter)
-	      L_m = PD.L_M(E_metric,size)
-	      return L_m
+       t0 = wavefront.OTF(psf_foc)
+       
+          
+       return t0,ph
 

__pycache__/cost_func.cpython-38.pyc

@@ -1,23 +1,16 @@
 
 import numpy as np
-import matplotlib.pyplot as plt
+
 import scipy
 from scipy.fftpack import fftshift, ifftshift, fft2, ifft2
-from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
-from scipy.signal import correlate2d as correlate
-from scipy.signal import general_gaussian
 from astropy.io import fits
-from scipy import ndimage
-from functools import partial
-import time
 import imreg_dft
 import pyfits
 import tools
 import aperture
 import PD
 import noise
 import wavefront
-from scipy.optimize import minimize, minimize_scalar
 import wavefront
 from wavefront import *
 import tools 

__pycache__/deconvolution.cpython-38.pyc

@@ -1,16 +1,10 @@
 import numpy as np
-import matplotlib.pyplot as plt
+
 import scipy
 from scipy.fftpack import fftshift, ifftshift, fft2, ifft2
-from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
-from scipy.signal import correlate2d as correlate
-from scipy.signal import general_gaussian
-from astropy.io import fits
-from scipy import ndimage
-from functools import partial
-import time
-import imreg_dft
-import pyfits
+
+
+
 
 
 # compute the gamma parameter

__pycache__/noise.cpython-38.pyc

@@ -0,0 +1,115 @@
+import time
+import imreg_dft
+import pyfits
+import tools
+import aperture
+import PD
+import noise
+import wavefront
+import deconvolution
+import cost_func
+from scipy.optimize import minimize, minimize_scalar
+import astropy 
+from astropy.io import fits
+import argparse
+import telescope
+from telescope import *
+
+from tools import imreg, apo2d
+from aperture import *
+from PD import *
+from noise import *
+from wavefront import *
+from deconvolution import *
+from cost_func import *
+
+
+class patch_pd(object):
+    def __init__(self,pd_data,Del,co_num,output_wf,output_mtf):
+      
+
+     self.data= fits.getdata(pd_data)
+     self.Del = Del
+     self.co_num = co_num
+     self.output_wf = output_wf
+     self.output_mtf = output_mtf
+
+
+     self.mean_im0 = self.data[0,500:1000,500:1000].mean()
+     self.mean_imk = self.data[1,500:1000,500:1000].mean()
+
+     self.Im0 = self.data[0,:,:]/self.mean_im0
+     self.Imk = self.data[1,:,:]/self.mean_imk
+
+     self.output_WF = np.zeros((2048,2048))
+     self.output_mtf = np.zeros((2048,2048))
+
+
+
+    def fit_patch(self):
+
+      upper = 1700
+      Nx = np.arange(300,upper,self.Del)
+      Ny = np.arange(300,upper,self.Del)
+      for n1 in Nx :
+        for n2 in Ny:
+
+            print('fitting the area in:'+str(n1)+','+str(n2))
+            im0 = self.Im0[n2:n2+self.Del,n1:n1+self.Del]
+            imk  = self.Imk[n2:n2+self.Del,n1:n1+self.Del]
+            im0 = im0/im0.mean()
+            imk = imk/imk.mean()
+            imk = imreg(im0,imk)
+
+
+           
+            im0 = apo2d(im0,10)
+            imk = apo2d(imk,10)
+        
+            d0,dk = FT(im0,imk)
+            gam =1# Gamma(d0,dk,M_gamma)
+
+            p0 =   np.zeros(self.co_num)
+            fit = cost_func(self.Del,0.5,1e-10,10,self.co_num,0.5,617.3e-6, 140,4125.3,0.5)
+            Mask = aperture.mask_pupil(fit.telescope.pupil_size(),fit.size)
+            noise_temp = noise_mask_high(fit.size,fit.cut_off)
+            noise_filter = fftshift(noise_temp)
+            def Minimise(coefficients):
+                  A_f = wavefront.pupil_foc(coefficients,fit.size,fit.telescope.pupil_size(),self.co_num)
+                  A_def = wavefront.pupil_defocus(coefficients,fit.size,fit.del_z,fit.telescope.pupil_size(),self.co_num)
+                  psf_foc = wavefront.PSF(Mask,A_f,False) 
+                  psf_defoc = wavefront.PSF(Mask,A_def,False)
+                  t0 = wavefront.OTF(psf_foc)
+                  tk = wavefront.OTF(psf_defoc) 
+
+                  q,q2 = PD.Q_matrix(t0,tk,fit.reg,gam)
+             
+
+                  F_m = PD.F_M(q2,d0, dk,t0,tk,noise_filter,gam)
+
+                  E_metric = PD.Error_metric(t0,tk,d0,dk,q,noise_filter)
+                  L_m = PD.L_M(E_metric,fit.size)
+                  return L_m
+
+       
+            Minimise_partial = partial(Minimise)
+            mini = scipy.optimize.minimize(Minimise_partial,p0,method= 'L-BFGS-B')
+            
+            result = fit.Minimize_res(mini.x)
+            self.output_WF[n2:n2+self.Del,n1:n1+self.Del] = result[1]
+            self.output_mtf[n2:n2+self.Del,n1:n1+self.Del] = MTF(fftshift(result[0]))
+            hdu = fits.PrimaryHDU(self.output_WF)
+            hdu.writeto(self.output_wf,overwrite=True)
+
+         
+if (__name__ == '__main__'):
+ parser = argparse.ArgumentParser(description='PD on sub-fields')
+ parser.add_argument('-i','--input', help='input')
+ parser.add_argument('-z','--Z', help='Zernike',default=10)
+ parser.add_argument('-d','--Del', help='Del',default=265)
+ parser.add_argument('-ow','--ow', help='output_WFE')
+ parser.add_argument('-om','--om', help='output MTF')
+ parsed = vars(parser.parse_args())
+ st = patch_pd(pd_data='{0}'.format(parsed['input']),Del=int(parsed['Del']),co_num=int(parsed['Z']),output_wf='{0}'.format(parsed['ow']),output_mtf='{0}'.format(parsed['om']))
+ st.fit_patch()
+   

__pycache__/telescope.cpython-38.pyc

@@ -1,17 +1,7 @@
 
 import numpy as np
-import matplotlib.pyplot as plt
-import scipy
-from scipy.fftpack import fftshift, ifftshift, fft2, ifft2
-from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
-from scipy.signal import correlate2d as correlate
-from scipy.signal import general_gaussian
-from astropy.io import fits
-from scipy import ndimage
-from functools import partial
-import time
-import imreg_dft
-import pyfits
+
+
 
 
 #from . import detector

__pycache__/tools.cpython-38.pyc

@@ -5,12 +5,10 @@
 from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
 from scipy.signal import correlate2d as correlate
 from scipy.signal import general_gaussian
-from astropy.io import fits
 from scipy import ndimage
-from functools import partial
-import time
-import imreg_dft
-import pyfits
+
+
+
 
 def GetPSD1D(psd2D):
     h  = psd2D.shape[0]

__pycache__/wavefront.cpython-38.pyc

@@ -1,16 +1,9 @@
 import numpy as np
-import matplotlib.pyplot as plt
+
 import scipy
 from scipy.fftpack import fftshift, ifftshift, fft2, ifft2
-from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
 from scipy.signal import correlate2d as correlate
-from scipy.signal import general_gaussian
-from astropy.io import fits
 from scipy import ndimage
-from functools import partial
-import time
-import imreg_dft
-import pyfits
 import zernike
 
 def phase(coefficients,rpupil, co_num):
@@ -42,7 +35,7 @@ def phase_aberr(del_z,rpupil):
 
 # Embed the phase in the center of a quadratic map of dimensions sizexsize
 def phase_embed(coefficients,rpupil,size,co_num):
-    ph =  phase(coefficients,rpupil)
+    ph =  phase(coefficients,rpupil,co_num)
     A = np.zeros([size,size])
     A[size//2-rpupil+1:size//2+rpupil+1,size//2-rpupil+1:size//2+rpupil+1]= phase(coefficients,rpupil,co_num)
     return A

__pycache__/zernike.cpython-38.pyc

@@ -1,20 +1,5 @@
 
 import numpy as np
-import matplotlib.pyplot as plt
-import scipy
-from scipy.fftpack import fftshift, ifftshift, fft2, ifft2
-from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
-from scipy.signal import correlate2d as correlate
-from scipy.signal import general_gaussian
-from astropy.io import fits
-from scipy import ndimage
-from functools import partial
-import time
-import imreg_dft
-import pyfits
-
-
-
 
 '''
 # Zernike polynomials in cartesian coordinates