Add possibility of fitting on a reduced FoV

maoppy/example/fit_muse_simulated_psf.py

@@ -3,6 +3,10 @@
 """
 Created on Fri May  8 18:45:28 2020
 
+-------------------
+SHOW FITTING POSSIBILITY WITH PSFAO and the MAOPPY LIBRARY
+-------------------
+
 Generate an observation of a star with MUSE-NFM, at a given wavelength
 Fit the simulated image with the Psfao model
 Plot results and display fitting accuracy

maoppy/example/fit_muse_simulated_psf_small_fov.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri May  8 21:07:53 2020
+
+--------------
+Show the possibility of fitting a PSF with Psfao on a small field of view
+and to retrieve the PSF shape on an increased field of view.
+--------------
+This script is a more complex version of 'fit_muse_simulated_psf.py'
+The other one should be run first.
+--------------
+
+Generate an observation of a star with MUSE-NFM, at a given wavelength
+Crop the image to simulate a limited field of view
+Fit the simulated image with the Psfao model
+Plot results and display fitting accuracy
+    Full images are shown
+    Red rectangles indicate the reduced field of view used for fitting
+
+@author: rfetick
+"""
+
+import matplotlib.pyplot as plt
+from matplotlib.patches import Rectangle
+import numpy as np
+
+from maoppy.psfmodel import Psfao, psffit, rmserror
+from maoppy.instrument import MUSE_NFM
+
+npix = 128 # pixel size of PSF
+wvl = 600*1e-9 # wavelength [m]
+
+flux = 1e6
+ron = MUSE_NFM.ron
+bckgd = 1.0
+
+#%% Initialize PSF model
+samp = MUSE_NFM.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
+Pmodel = Psfao((npix,npix),system=MUSE_NFM,symmetric=True,samp=samp)
+
+#%% Generate a MUSE-NFM PSF
+r0 = 0.15 # Fried parameter [m]
+b = 1e-7 # Phase PSD background [rad² m²]
+amp = 1.4 # Phase PSD Moffat amplitude [rad²]
+ax = 0.1 # Phase PSD Moffat alpha [1/m]
+beta = 1.6 # Phase PSD Moffat beta power law
+
+param = [r0,b,amp,ax,beta]
+
+psf = Pmodel(param,dx=0,dy=0)
+
+noise = ron*np.random.randn(npix,npix)
+image = flux*psf + bckgd + noise
+
+#%% Crop the image
+ncrop = 64
+im_crop = image[npix//2-ncrop//2:npix//2+ncrop//2,
+                npix//2-ncrop//2:npix//2+ncrop//2]
+
+#%% Fit the cropped image with Psfao
+guess = [0.145,2e-7,1.2,0.08,1.5]
+w = np.ones_like(im_crop)/ron**2.0
+
+out = psffit(im_crop,Psfao,guess,weights=w,npixfit=npix, # fit keywords
+             system=MUSE_NFM,samp=samp,symmetric=True) # Psfao keywords
+
+flux_fit, bck_fit = out.flux_bck
+fitao = flux_fit*out.psf + bck_fit
+
+#%% Plots
+print(31*'-')
+print("%9s %10s %10s"%('','MUSE-NFM','Psfao'))
+print(31*'-')
+print("%9s %10.4g %10.4g"%('Flux [e-]',flux,flux_fit))
+print("%9s %10.4g %10.4g"%('Bck [e-]',bckgd,bck_fit))
+print("%9s %10.2f %10.2f"%('r0 [cm]',r0*100,out.x[0]*100))
+print("%9s %10.2f %10.2f"%('A [rad²]',amp,out.x[2]))
+print(31*'-')
+print("shape + noise RMS error = %.2f%%"%(100*rmserror(fitao,image)))
+print("        noise RMS error = %.2f%%"%(np.sqrt(np.sum(noise**2))/np.sum(image)*100))
+print(31*'-')
+
+corner = (npix//2-ncrop//2,npix//2-ncrop//2)
+r1 = Rectangle(corner,ncrop,ncrop,linewidth=2,edgecolor='r',facecolor='none')
+r2 = Rectangle(corner,ncrop,ncrop,linewidth=2,edgecolor='r',facecolor='none')
+r3 = Rectangle(corner,ncrop,ncrop,linewidth=2,edgecolor='r',facecolor='none')
+
+vmin,vmax = np.arcsinh([image.min(),image.max()])
+
+plt.figure(1)
+plt.clf()
+
+ax1 = plt.subplot(131)
+plt.imshow(np.arcsinh(image),vmin=vmin,vmax=vmax)
+plt.axis('off')
+plt.title('MUSE-NFM simu')
+ax1.add_patch(r1) # actual field of view for fitting
+
+ax2 = plt.subplot(132)
+plt.imshow(np.arcsinh(fitao),vmin=vmin,vmax=vmax)
+plt.axis('off')
+plt.title('Psfao fit')
+ax2.add_patch(r2) # actual field of view for fitting
+
+ax3 = plt.subplot(133)
+plt.imshow(np.arcsinh(image-fitao))
+plt.axis('off')
+plt.title('residuals')
+ax3.add_patch(r3) # actual field of view for fitting

maoppy/psfmodel.py

@@ -63,7 +63,7 @@ def lsq_flux_bck(model, data, weights, background=True, positive_bck=False):
     return amp, bck
 
 def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
-           positive_bck=False,fixed=None,**kwargs):
+           positive_bck=False,fixed=None,npixfit=None,**kwargs):
     """Fit a PSF with a parametric model solving the least-square problem
        epsilon(x) = SUM_pixel { weights * (amp * Model(x) + bck - psf)² }
     
@@ -88,6 +88,8 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
         Force background to be positive or null
     fixed : numpy.ndarray
         Fix some parameters to their initial value (default: None)
+    npixfit : int (default=None)
+        Increased pixel size for improved fitting accuracy
     **kwargs :
         All keywords used to instantiate your `Model`
     
@@ -114,10 +116,23 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
        .cost : float
            Value of cost function at optimum
     """
-    Model_inst = Model(np.shape(psf),**kwargs)
+    
     if weights is None: weights = np.ones_like(psf)
     elif len(psf)!=len(weights): raise ValueError("Keyword `weights` must have same number of elements as `psf`")
+    
+    # Increase array size for improved fitting accuracy
+    if npixfit is not None:
+        sx,sy = np.shape(psf)
+        if (sx>npixfit) or (sy>npixfit): raise ValueError('npixfit must be greater or equal to both psf dimensions')
+        psfbig = np.zeros((npixfit,npixfit))
+        wbig = np.zeros((npixfit,npixfit))
+        psfbig[npixfit//2-sx//2:npixfit//2+sx//2,npixfit//2-sy//2:npixfit//2+sy//2] = psf
+        wbig[npixfit//2-sx//2:npixfit//2+sx//2,npixfit//2-sy//2:npixfit//2+sy//2] = weights
+        psf = psfbig
+        weights = wbig
+    
     sqW = np.sqrt(weights)
+    Model_inst = Model(np.shape(psf),**kwargs)
     
     class CostClass(object):
         def __init__(self):