Commit c5726862 authored by rfetick's avatar rfetick
Browse files

psffit takes an instance instead of a class

parent a5b9650d
......@@ -30,7 +30,7 @@ bckgd = 1.0
#%% Initialize PSF model
samp = muse_nfm.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
Pmodel = Psfao((npix,npix),system=muse_nfm,samp=samp)
Pmodel = Psfao((npix,npix), system=muse_nfm, samp=samp)
#%% Generate a MUSE-NFM PSF
r0 = 0.15 # Fried parameter [m]
......@@ -53,7 +53,7 @@ guess = [0.145,2e-2,1.2,0.08,ratio,theta,1.5]
w = np.ones_like(image)/ron**2.0
fixed = [False,False,False,False,True,True,False]
out = psffit(image,Psfao,guess,weights=w,system=muse_nfm,samp=samp,fixed=fixed)
out = psffit(image, Pmodel, guess, weights=w, fixed=fixed)
flux_fit, bck_fit = out.flux_bck
fitao = flux_fit*out.psf + bck_fit
......
......@@ -38,7 +38,7 @@ bckgd = 1.0
#%% Initialize PSF model
samp = muse_nfm.samp(wvl) # sampling (2.0 for Shannon-Nyquist)
Pmodel = Psfao((npix,npix),system=muse_nfm,samp=samp)
Pmodel = Psfao((npix,npix), system=muse_nfm, samp=samp)
#%% Generate a MUSE-NFM PSF
r0 = 0.15 # Fried parameter [m]
......@@ -66,8 +66,7 @@ guess = [0.145,2e-2,1.2,0.08,ratio,theta,1.5]
w = np.ones_like(im_crop)/ron**2.0
fixed = [False,False,False,False,True,True,False]
out = psffit(im_crop,Psfao,guess,weights=w,fixed=fixed,npixfit=npix, # fit keywords
system=muse_nfm,samp=samp) # Psfao keywords
out = psffit(im_crop, Pmodel, guess, weights=w, fixed=fixed, npixfit=npix)
flux_fit, bck_fit = out.flux_bck
fitao = flux_fit*out.psf + bck_fit
......
......@@ -21,7 +21,7 @@ from maoppy.instrument import zimpol
# %% PARAMETERS TO MODIFY
psf_path = "path/to/your/zimpol_psf.fits" # set the path to your PSF *.fits file
filt = zimpol.filters["N_R"] # ZIMPOL filter
Npix = 512 # size of PSF for fitting [pixels]. The larger the better
npix = 512 # size of PSF for fitting [pixels]. The larger the better
r0 = 0.18 # Fried parameter [m]
moff_C = 1e-2 # PSD constant [rad²m²]
......@@ -36,7 +36,7 @@ fitsPSF = fits.open(psf_path)
hdr = fitsPSF[0].header
psf = fitsPSF[0].data
psf = imcenter(psf, (Npix, Npix), maxi=True) * zimpol.gain
psf = imcenter(psf, (npix,npix), maxi=True) * zimpol.gain
# %% FIT PSFAO model
wvl = filt[0]
......@@ -48,7 +48,8 @@ fixed = [False,False,False,False,True,True,False]
weights = np.ones_like(psf)
print("PSFAO fitting (please wait)")
out = psffit(psf,Psfao,x0,weights=weights,system=zimpol,samp=samp,fixed=fixed)
model = Psfao((npix,npix), system=zimpol, samp=samp)
out = psffit(psf, model, x0, weights=weights, fixed=fixed)
out.x[5] = out.x[5]%np.pi # make angle between 0 and PI
......@@ -111,4 +112,4 @@ plt.legend()
plt.grid()
plt.xlabel("Radius [pix]")
plt.ylabel("PSF [photon]")
plt.xlim(-Npix//2,Npix//2)
\ No newline at end of file
plt.xlim(-npix//2,npix//2)
\ No newline at end of file
......@@ -50,7 +50,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),
def psffit(psf,model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
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)² }
......@@ -59,8 +59,8 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
----------
psf : numpy.ndarray
The experimental image to be fitted
Model : class
The class representing the fitting model
model : instance of ParametricPSF (or its subclasses)
The model to fit
x0 : tuple, list, numpy.ndarray
Initial guess for parameters
weights : numpy.ndarray
......@@ -79,7 +79,7 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
npixfit : int (default=None)
Increased pixel size for improved fitting accuracy
**kwargs :
All keywords used to instantiate your `Model`
All keywords used to call your `model`
Returns
-------
......@@ -123,7 +123,6 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
weights = wbig
sqw = np.sqrt(weights)
model_inst = Model(np.shape(psf),**kwargs)
class CostClass:
"""A cost function that can print iterations"""
......@@ -134,9 +133,9 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
print("-",end="")
self.iter += 1
x, dxdy = mini2input(y)
mm = model_inst(x,dx=dxdy[0],dy=dxdy[1])
mm = model(x, dx=dxdy[0], dy=dxdy[1], **kwargs)
amp, bck = lsq_flux_bck(mm, psf, weights, background=flux_bck[1], positive_bck=positive_bck)
return np.reshape(sqw * (amp * mm + bck - psf), np.size(psf))
return np.reshape(sqw * (amp*mm + bck - psf), np.size(psf))
cost = CostClass()
......@@ -174,13 +173,13 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
b_up = np.concatenate((b_up,[np.inf,np.inf]))
return (b_low,b_up)
result = least_squares(cost, input2mini(x0,dxdy), bounds=get_bound(model_inst))
result = least_squares(cost, input2mini(x0,dxdy), bounds=get_bound(model))
print("") #finish line of "-"
result.x, result.dxdy = mini2input(result.x) #split output between x and dxdy
m = model_inst(result.x,dx=result.dxdy[0],dy=result.dxdy[1])
m = model(result.x,dx=result.dxdy[0],dy=result.dxdy[1])
amp, bck = lsq_flux_bck(m, psf, weights, background=flux_bck[1], positive_bck=positive_bck)
result.flux_bck = (amp,bck)
......
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