Moffat gauss

spectractor/extractor/background.py

@@ -5,6 +5,7 @@
 import os
 import copy
 
+import photutils.utils.exceptions
 from spectractor import parameters
 from spectractor.tools import fit_poly1d_outlier_removal, fit_poly2d_outlier_removal, plot_image_simple
 
@@ -15,7 +16,8 @@
 from scipy.signal import medfilt2d
 from scipy import ndimage
 from scipy.interpolate import RegularGridInterpolator
-
+import warnings
+warnings.filterwarnings("ignore", category=photutils.utils.exceptions.NoDetectionsWarning)
 
 def _from_bkgd_interp_to_func(bgd_model_func_interp):
     def bgd_model_func(x, y):

spectractor/extractor/chromaticpsf.py

@@ -202,7 +202,9 @@ def crop_table(self, new_Nx):
     def set_polynomial_degrees(self, deg):
         self.deg = deg
         self.degrees = {key: deg for key in self.psf.params.labels}
-        self.degrees["x_c"] = max(self.degrees["x_c"], 1)
+        for key in self.psf.params.labels:
+            if "x_" in key:
+                self.degrees[key] = max(self.degrees[key], 1)
         self.degrees['saturation'] = 0
 
     def generate_test_poly_params(self):
@@ -222,9 +224,10 @@ def generate_test_poly_params(self):
         >>> s = ChromaticPSF(psf, Nx=5, Ny=4, deg=1, saturation=20000)
         >>> params = s.generate_test_poly_params()
 
-        ..  doctest::
+        .. doctest::
             :hide:
-            >>> assert(np.all(np.isclose(params,[10, 50, 100, 150, 200, 0, 0, 0, 0, 5, 0, 2, 0, -0.4, 0, 1, 0, 20000])))
+
+            >>> assert(np.all(np.isclose(params,[10, 50, 100, 150, 200, 0, 0, 0, 0, 5, 0, 2, 0, 0.4, 0, 0, 0, 4, 0, 1, 0, 20000])))
 
         """
         if not isinstance(self.psf, MoffatGauss) and not isinstance(self.psf, Moffat):
@@ -241,7 +244,9 @@ def generate_test_poly_params(self):
         params += [0.] * (self.degrees['gamma'] - 1) + [0, 5]  # gamma
         params += [0.] * (self.degrees['alpha'] - 1) + [0, 2]  # alpha
         if isinstance(self.psf, MoffatGauss):
-            params += [0.] * (self.degrees['eta_gauss'] - 1) + [0, -0.4]  # eta_gauss
+            params += [0.] * (self.degrees['eta_gauss'] - 1) + [0, 0.4]  # eta_gauss
+            params += [0.] * (self.degrees['x_g'] - 1) + [0, 0]  # x mean
+            params += [0.] * (self.degrees['y_g'] - 1) + [0, 4]  # y mean
             params += [0.] * (self.degrees['stddev'] - 1) + [0, 1]  # stddev
         params += [self.saturation]  # saturation
         poly_params = np.zeros_like(params)
@@ -1116,7 +1121,7 @@ def from_profile_params_to_poly_params(self, profile_params, indices=None):
 
         ..  doctest::
             :hide:
-            >>> assert(np.all(np.isclose(profile_params[0], [10, 0, 50, 5, 2, -0.4, 1, 8e3])))
+            >>> assert(np.all(np.isclose(profile_params[0], [10, 0, 50, 5, 2, 0.4, 0, 53, 1, 8e3])))
 
         From the profile parameters to the polynomial parameters:
 
@@ -1143,9 +1148,9 @@ def from_profile_params_to_poly_params(self, profile_params, indices=None):
             delta = 0
             if name != 'amplitude':
                 weights = np.copy(amplitude)[indices]
-                if name == 'x_c':
+                if 'x_' in name:
                     delta = self.x0
-                if name == 'y_c':
+                if 'y_' in name:
                     delta = self.y0
                 if parameters.PSF_POLY_TYPE == "legendre":
                     fit = np.polynomial.legendre.legfit(poly_x, profile_params[indices, k] - delta,
@@ -1255,24 +1260,24 @@ def from_poly_params_to_profile_params(self, poly_params, apply_bounds=False):
 
         ..  doctest::
             :hide:
-            >>> assert np.allclose(profile_params[0], [10, 0, 50, 5, 2, -0.4, 1, 8e3], rtol=1e-3, atol=1e-3)
+            >>> assert np.allclose(profile_params[0], [10, 0, 50, 5, 2, 0.4, 0, 54, 1, 8e3], rtol=1e-3, atol=1e-3)
 
         From the profile parameters to the polynomial parameters:
 
-        >>> profile_params = s.from_profile_params_to_poly_params(profile_params)
+        >>> poly_params = s.from_profile_params_to_poly_params(profile_params)
 
         ..  doctest::
             :hide:
 
-            >>> assert np.allclose(profile_params, poly_params_test)
+            >>> assert np.allclose(poly_params, poly_params_test)
 
         From the polynomial parameters to the profile parameters without Moffat amplitudes:
 
         >>> profile_params = s.from_poly_params_to_profile_params(poly_params_test[100:])
 
         ..  doctest::
             :hide:
-            >>> assert np.allclose(profile_params[0], [1, 0, 50, 5, 2, -0.4, 1, 8e3])
+            >>> assert np.allclose(profile_params[0], [1, 0, 50, 5, 2, 0.4, 0, 54, 1, 8e3])
 
         """
         length = len(self.table)
@@ -1299,9 +1304,9 @@ def from_poly_params_to_profile_params(self, poly_params, apply_bounds=False):
                         elif parameters.PSF_POLY_TYPE == "polynomial":
                             profile_params[:, k] = np.polynomial.polynomial.polyval(poly_x, p)
                 shift += self.degrees[name] + 1
-                if name == 'x_c':
+                if 'x_' in name:
                     profile_params[:, k] += self.x0
-                if name == 'y_c':
+                if 'y_' in name:
                     profile_params[:, k] += self.y0
         if apply_bounds:
             for k, name in enumerate(self.psf.params.labels):
@@ -1450,9 +1455,9 @@ def plot_summary(self, truth=None):
         for i, name in enumerate(self.psf.params.labels):
             coeffs = [self.params.values[k] for k in range(self.params.ndim) if name in self.params.labels[k]]
             delta = 0
-            if name == 'x_c':
+            if "x_" in name:
                 delta = self.x0
-            if name == 'y_c':
+            if "y_" in name:
                 delta = self.y0
             if parameters.PSF_POLY_TYPE == "legendre":
                 PSF_models.append(np.polynomial.polynomial.legval(rescale_x_to_legendre(all_pixels), coeffs) + delta)
@@ -1546,6 +1551,8 @@ def fit_transverse_PSF1D_profile(self, data, err, w, ws, pixel_step=1, bgd_model
             gamma True
             alpha True
             eta_gauss True
+            x_g True
+            y_g True
             stddev True
             saturation True
             >>> assert(not np.any(np.isclose(s.table['flux_sum'][3:6], np.zeros(s.Nx)[3:6], rtol=1e-3)))
@@ -1573,8 +1580,6 @@ def fit_transverse_PSF1D_profile(self, data, err, w, ws, pixel_step=1, bgd_model
         # guess = [2 * np.nanmax(signal), middle, 0.5 * fwhm, 2, 0, 0.1 * fwhm, saturation]
         signal_sum = np.nanmax(signal)
         guess[0] = signal_sum
-        guess[1] = xmax_index
-        guess[2] = middle
         guess[-1] = saturation
         # bounds = [(0.1 * maxi, 10 * maxi), (middle - w, middle + w), (0.1, min(fwhm, Ny // 2)), (0.1, self.alpha_max),
         #           (-1, 0),
@@ -1584,8 +1589,13 @@ def fit_transverse_PSF1D_profile(self, data, err, w, ws, pixel_step=1, bgd_model
         initial_bounds = np.copy(psf.params.bounds)
         bounds = np.copy(psf.params.bounds)
         bounds[0] = (0.1 * signal_sum, 2 * signal_sum)
-        bounds[2] = (middle - w, middle + w)
         bounds[-1] = (0, 2 * saturation)
+        for k, label in enumerate(psf.params.labels):
+            if "x_" in label:
+                guess[k] = xmax_index
+            if "y_" in label:
+                guess[k] = middle
+                bounds[k] = (middle - w, middle + w)
         # moffat_guess = [2 * np.nanmax(signal), middle, 0.5 * fwhm, 2]
         # moffat_bounds = [(0.1 * maxi, 10 * maxi), (middle - w, middle + w), (0.1, min(fwhm, Ny // 2)), (0.1, 10)]
         # fit = fit_moffat1d_outlier_removal(index, signal, sigma=sigma, niter=2,
@@ -1630,7 +1640,9 @@ def fit_transverse_PSF1D_profile(self, data, err, w, ws, pixel_step=1, bgd_model
             # bounds[0] = (0.1 * np.nanstd(bgd), 2 * np.nanmax(y[middle - ws[0]:middle + ws[0]]))
             bounds[0] = (0.1 * signal_sum, 1.5 * signal_sum)
             guess[0] = signal_sum
-            guess[1] = x
+            for k, label in enumerate(psf.params.labels):
+                if "x_" in label:
+                    guess[k] = x
             # if guess[4] > -1:
             #    guess[0] = np.max(signal) / (1 + guess[4])
             # std = np.sqrt(np.nansum(pdf * (index - mean) ** 2))
@@ -1724,6 +1736,7 @@ def fit_chromatic_psf(self, data, mask=None, bgd_model_func=None, data_errors=No
         >>> parameters.PIXWIDTH_BACKGROUND = 10
         >>> parameters.PIXWIDTH_SIGNAL = 30
         >>> parameters.DEBUG = True
+        >>> parameters.VERBOSE = True
 
         Build a mock spectrogram with random Poisson noise using the full 2D PSF model:
 
@@ -1857,7 +1870,8 @@ def fit_chromatic_psf(self, data, mask=None, bgd_model_func=None, data_errors=No
         self.cov_matrix = np.copy(w.amplitude_cov_matrix)
 
         # add background crop to y_c
-        self.params.values[w.Nx + w.y_c_0_index] += w.bgd_width
+        for y0_index in w.y0_params:
+            self.params.values[self.Nx + y0_index] += w.bgd_width
 
         # fill results
         self.psf.apply_max_width_to_bounds(max_half_width=w.Ny + 2 * w.bgd_width)
@@ -1885,7 +1899,7 @@ def __init__(self, chromatic_psf, data, data_errors, mode, bgd_model_func=None,
                                axis_names=list(np.copy(chromatic_psf.params.axis_names[length:])), fixed=None,
                                truth=truth, filename=file_name)
         for k, par in enumerate(params.labels):
-            if "x_c" in par or "saturation" in par:
+            if "x_" in par or "saturation" in par:
                 params.fixed[k] = True
         FitWorkspace.__init__(self, params, file_name=file_name, verbose=verbose, plot=plot, live_fit=live_fit)
         self.my_logger = set_logger(self.__class__.__name__)
@@ -1895,11 +1909,10 @@ def __init__(self, chromatic_psf, data, data_errors, mode, bgd_model_func=None,
         self.bgd_model_func = bgd_model_func
         self.analytical = analytical
         self.poly_params = np.copy(self.chromatic_psf.params.values)
-        self.y_c_0_index = -1
+        self.y0_params = []
         for k, par in enumerate(params.labels):
-            if par == "y_c_0":
-                self.y_c_0_index = k
-                break
+            if "y_" in par and "_0" in par:
+                self.y0_params.append(k)
 
         # prepare the fit
         self.Ny, self.Nx = self.data.shape
@@ -1936,8 +1949,8 @@ def __init__(self, chromatic_psf, data, data_errors, mode, bgd_model_func=None,
             yy, xx = np.mgrid[:self.Ny, :self.Nx]
             self.pixels = np.asarray([xx, yy])
 
-
-        self.poly_params[self.Nx + self.y_c_0_index] -= self.bgd_width
+        for y0_index in self.y0_params:
+            self.poly_params[self.Nx + y0_index] -= self.bgd_width
         self.profile_params = self.chromatic_psf.from_poly_params_to_profile_params(self.poly_params)
         self.data_before_mask = np.copy(self.data)
         self.mask_before_mask = list(np.copy(self.mask))
@@ -2274,7 +2287,8 @@ def simulate(self, *shape_params):
         # linear regression for the amplitude parameters
         # prepare the vectors
         poly_params = np.concatenate([np.ones(self.Nx), shape_params])
-        poly_params[self.Nx + self.y_c_0_index] -= self.bgd_width
+        for y0_index in self.y0_params:
+            poly_params[self.Nx + y0_index] -= self.bgd_width
         profile_params = self.chromatic_psf.from_poly_params_to_profile_params(poly_params, apply_bounds=True)
         profile_params[:self.Nx, 0] = 1
         profile_params[:self.Nx, 1] = np.arange(self.Nx)
@@ -2348,7 +2362,8 @@ def simulate(self, *shape_params):
         if self.amplitude_priors_method == "fixed":
             self.model = self.chromatic_psf.evaluate(self.pixels, poly_params)
             # self.chromatic_psf.params.values is updated in evaluate(): reset to original values
-            self.chromatic_psf.params.values[self.Nx + self.y_c_0_index] += self.bgd_width
+            for y0_index in self.y0_params:
+                self.chromatic_psf.params.values[self.Nx + y0_index] += self.bgd_width
         self.poly_params = np.copy(poly_params)
         self.profile_params = np.copy(profile_params)
         self.model_err = np.zeros_like(self.model)

spectractor/extractor/extractor.py

@@ -477,8 +477,11 @@ def simulate(self, *params):
 
             # Fill spectrogram trace as a function of the pixel column x
             self.psf_profile_params[order][:, 0] = self.params[f"A{order}"] * self.tr[k](self.lambdas)
-            self.psf_profile_params[order][:, 1] = dispersion_law.real + self.spectrum.spectrogram_x0
-            self.psf_profile_params[order][:, 2] += dispersion_law.imag - self.bgd_width
+            for p, label in enumerate(self.spectrum.chromatic_psf.psf.params.labels):
+                if "x_" in label:
+                    self.psf_profile_params[order][:, p] = dispersion_law.real + self.spectrum.spectrogram_x0
+                if "y_" in label:
+                    self.psf_profile_params[order][:, p] += dispersion_law.imag - self.bgd_width
 
             # Matrix filling
             M_order = self.spectrum.chromatic_psf.build_sparse_M(self.pixels, self.psf_profile_params[order],

spectractor/extractor/images.py

@@ -421,10 +421,8 @@ def check_statistical_error(self):
             raise AttributeError(f"Noise map must be in ADU units to be plotted and analyzed. "
                                  f"Currently self.units={self.units}.")
         data = np.copy(self.data)
-        min_noz = np.min(data[data > 0])
-        data[data <= 0] = min_noz
-        y = self.err.flatten() ** 2
-        x = data.flatten()
+        y = self.err[data > 0].flatten() ** 2
+        x = data[data > 0].flatten()
         fit, cov, model = fit_poly1d(x, y, order=1)
         gain = 1 / fit[0]
         read_out = np.sqrt(fit[1]) * gain
@@ -911,7 +909,7 @@ def find_target(image, guess=None, rotated=False, widths=[parameters.XWINDOW, pa
         if parameters.PdfPages:
             parameters.PdfPages.savefig()
 
-    if parameters.SPECTRACTOR_FIT_TARGET_CENTROID == "fit" or rotated:
+    if parameters.SPECTRACTOR_FIT_TARGET_CENTROID == "fit":
         if target_pixcoords[0] == -1 and target_pixcoords[1] == -1:
             if guess is None:
                 raise ValueError(f"Guess parameter must be set if WCS solution is not found.")
@@ -961,7 +959,12 @@ def find_target(image, guess=None, rotated=False, widths=[parameters.XWINDOW, pa
             sub_image_y0 = target_pixcoords[1] - y0 + Dy
     elif parameters.SPECTRACTOR_FIT_TARGET_CENTROID == "guess":
         Dx, Dy = widths
-        sub_image_subtracted, x0, y0, Dx, Dy, sub_errors = find_target_init(image=image, guess=target_pixcoords,
+        if rotated:
+            guess2 = find_target_after_rotation(image)
+            x0 = int(guess2[0])
+            y0 = int(guess2[1])
+            guess = [x0, y0]
+        sub_image_subtracted, x0, y0, Dx, Dy, sub_errors = find_target_init(image=image, guess=guess,
                                                                             rotated=rotated, widths=(Dx, Dy))
         theX, theY = guess
         sub_image_x0 = theX - x0 + Dx
@@ -1186,7 +1189,7 @@ def find_target_Moffat2D(image, sub_image_subtracted, sub_errors=None):
     # fit a 2D star profile close to this position
     psf = Moffat(clip=True)
     total_flux = np.sum(sub_image_subtracted)
-    psf.params.values[:3] = [total_flux, avX, avY]
+    psf.params.values[:5] = [total_flux, avX, avY, np.mean([sigX, sigY]), 3]
     psf.params.values[-1] = image.saturation
     if image.target_star2D is not None:
         psf.params.values = image.target_star2D.params.values
@@ -1217,6 +1220,8 @@ def find_target_Moffat2D(image, sub_image_subtracted, sub_errors=None):
         star2D = psf.evaluate(pixels=np.array([X, Y]))
         plot_image_simple(ax1, data=sub_image_subtracted, scale="lin", title="", units=image.units,
                           target_pixcoords=[new_avX, new_avY], vmin=vmin, vmax=vmax)
+        plt.axvline(avX, color='r', linestyle='-', lw=2, label="profile**4 guess")
+        plt.axhline(avY, color='r', linestyle='-', lw=2)
         ax1.legend(loc=1)
 
         ax1.text(0.05, 0.05, f'Data', color="white",