From 97640621157e967ab3254f30d9e6e16d25b7c6e4 Mon Sep 17 00:00:00 2001
From: rfetick <r.fetick@gmail.com>
Date: Fri, 8 May 2020 18:37:29 +0200
Subject: [PATCH] Add small functions getgaussnoise and rmserror
---
maoppy/psfmodel.py | 38 +++++++++++++++++++-------------------
maoppy/utils.py | 11 +++++++++++
2 files changed, 30 insertions(+), 19 deletions(-)
diff --git a/maoppy/psfmodel.py b/maoppy/psfmodel.py
index a6f6141..4cd98df 100644
--- a/maoppy/psfmodel.py
+++ b/maoppy/psfmodel.py
@@ -12,9 +12,18 @@ from astropy.io import fits
import time
from numpy.fft import fft2, fftshift, ifft2
from maoppy.config import _EPSILON
-from maoppy.utils import binning
+from maoppy.utils import binning, getgaussnoise
#%% FITTING FUNCTION
+def rmserror(model,image,weights=1.0,background=True, positive_bck=False):
+ """Compute the RMS error between a PSF model and an image"""
+ amp,bck = lsq_flux_bck(model, image, weights=weights, background=background, positive_bck=positive_bck)
+ diff = amp*model+bck - image
+ noise_std = getgaussnoise(image)
+ err = np.sqrt(np.sum(diff**2.0 - noise_std**2.0))
+ err = err/np.sum(image)
+ return err
+
def lsq_flux_bck(model, data, weights, background=True, positive_bck=False):
"""Compute the analytical least-square solution for flux and background
LS = SUM_pixels { weights*(flux*model + bck - data)² }
@@ -104,20 +113,16 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
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`")
+ 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`")
sqW = np.sqrt(weights)
class CostClass(object):
def __init__(self):
self.iter = 0
def __call__(self,y):
- if (self.iter%3) == 0:
- print("-",end="")
+ if (self.iter%3)==0: print("-",end="")
self.iter += 1
-
x, dxdy = mini2input(y)
m = Model_inst(x,dx=dxdy[0],dy=dxdy[1])
amp, bck = lsq_flux_bck(m, psf, weights, background=flux_bck[1], positive_bck=positive_bck)
@@ -126,17 +131,14 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
cost = CostClass()
if fixed is not None:
- if len(fixed)!=len(x0):
- raise ValueError("When defined, `fixed` must be same size as `x0`")
+ if len(fixed)!=len(x0): raise ValueError("When defined, `fixed` must be same size as `x0`")
FREE = [not fixed[i] for i in range(len(fixed))]
INDFREE = np.where(FREE)[0]
def input2mini(x,dxdy):
# Transform user parameters to minimizer parameters
- if fixed is None:
- xfree = x
- else:
- xfree = np.take(x,INDFREE)
+ if fixed is None: xfree = x
+ else: xfree = np.take(x,INDFREE)
return np.concatenate((xfree,dxdy))
def mini2input(y):
@@ -151,12 +153,10 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
def get_bound(inst):
b_low = inst.bounds[0]
- if fixed is not None:
- b_low = np.take(b_low,INDFREE)
+ if fixed is not None: b_low = np.take(b_low,INDFREE)
b_low = np.concatenate((b_low,[-np.inf,-np.inf]))
b_up = inst.bounds[1]
- if fixed is not None:
- b_up = np.take(b_up,INDFREE)
+ if fixed is not None: b_up = np.take(b_up,INDFREE)
b_up = np.concatenate((b_up,[np.inf,np.inf]))
return (b_low,b_up)
@@ -173,7 +173,7 @@ def psffit(psf,Model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
result.psf = m
return result
-#%% BESIC FUNCTIONS
+#%% BASIC FUNCTIONS
def reduced_coord(XY,ax,ay,theta,cx,cy):
c = np.cos(theta)
s = np.sin(theta)
diff --git a/maoppy/utils.py b/maoppy/utils.py
index df80f0c..42d569a 100644
--- a/maoppy/utils.py
+++ b/maoppy/utils.py
@@ -8,12 +8,23 @@ Created on Mon May 27 17:27:44 2019
import numpy as np
from scipy.special import jv
+from numpy.fft import fft2, fftshift
#%% CONSTANTS
# Radian to arcsec
RAD2ARCSEC = 180./np.pi * 3600.
#%% IMAGE PROCESSING FUNCTIONS
+def getgaussnoise(image):
+ """Compute the std deviation of a gaussian noise in an image"""
+ sx, sy = np.shape(image)
+ ft_im2 = fftshift(np.abs(fft2(image))**2)
+ sn = np.sum(ft_im2[0,0:sy-1])+np.sum(ft_im2[0:sx-1,sy-1])+np.sum(ft_im2[sx-1,sy-1:0:-1])+np.sum(ft_im2[sx-1:0:-1,0])
+ sn = sn/(2*sx+2*sy-4)
+ sigma = np.sqrt(sn/(sx*sy))
+ return sigma
+
+
def circarr(Npix, center=(None, None)):
"""Compute array of radii to the center of array
--
GitLab