From de5435f812bd8888e14cad2f22019e863f2ad23f Mon Sep 17 00:00:00 2001 From: Sean Lockwood Date: Wed, 29 Jan 2025 12:34:14 -0500 Subject: [PATCH 1/4] Version 3.0.0. Based on PR #37. Remove dependence on deprecated stsci.tools.gfit Gaussian fitting; replace with Astropy fitting. Allow for scipy API change for scipy.signal.boxcar to scipy.signal.windows.boxcar. Changed Gaussian fitting initial guess to peak. Aggregate Astropy fitting warnings (with ierr > 4 into a single warning for each call of Trace.gFitTrace(). Set data types to Float as they are created rather than after, avoiding conversions. --- stistools/mktrace.py | 76 +++++++++++++++++++++++++++----------------- 1 file changed, 46 insertions(+), 30 deletions(-) diff --git a/stistools/mktrace.py b/stistools/mktrace.py index f6dc1989..c64b0111 100644 --- a/stistools/mktrace.py +++ b/stistools/mktrace.py @@ -1,13 +1,18 @@ #!/usr/bin/env python -import numpy as np -from astropy.io import fits import os -import os.path import stat -from scipy import signal +import warnings +from collections import Counter +import numpy as np +from astropy.io import fits +from astropy.modeling import models, fitting from scipy import ndimage as ni +try: + from scipy.signal import boxcar +except ImportError: + from scipy.signal.windows import boxcar -from stsci.tools import gfit, linefit +from stsci.tools import linefit from stsci.tools import fileutil as fu __doc__ = """ @@ -39,9 +44,9 @@ - Python version: Nadia Dencheva """ +__version__ = '3.0.0' +__vdate__ = '2025-01-28' -__version__ = '2.0.0' -__vdate__ = '2017-03-20' def mktrace(fname, tracecen=0.0, weights=None): @@ -99,7 +104,7 @@ def mktrace(fname, tracecen=0.0, weights=None): w = np.zeros(1024) w[ind] = 1 - X = np.arange(1024).astype(np.float64) + X = np.arange(1024, dtype=float) sparams = linefit.linefit(X, trace1024, weights=w) rparams = linefit.linefit(X, interp_trace, weights=w) sciline = sparams[0] + sparams[1] * X @@ -119,10 +124,9 @@ def mktrace(fname, tracecen=0.0, weights=None): tr.writeTrace(fname, sciline, refline, interp_trace, trace1024, tr_ind, a2disp_ind) - #print 'time', time.time()-start #the minus sign is for consistency withthe way x2d reports the rotation - print("Traces were rotated by {:.10f} degrees \n".format((-(sparams[1]-rparams[1])*180 / np.pi))) - print('trace is centered on row {:.10f}'.format(tr._a2center)) + print("Traces were rotated by {:.10f} degrees \n".format(-(sparams[1] - rparams[1]) * 180 / np.pi)) + print('Trace is centered on row {:.10f}'.format(tr._a2center)) return tr @@ -140,10 +144,10 @@ def interp(y, n): """ m = float(len(y)) x = np.arange(m) - i = np.arange(n,dtype=np.float64) - xx = i * (m-1)/n - xind = np.searchsorted(x, xx)-1 - yy = y[xind]+(xx-x[xind])*(y[xind+1]-y[xind])/(x[xind+1]-x[xind]) + i = np.arange(n, dtype=float) + xx = i * (m - 1) / n + xind = np.searchsorted(x, xx) - 1 + yy = y[xind] + (xx - x[xind]) * (y[xind + 1] - y[xind]) / (x[xind + 1] - x[xind]) return yy @@ -266,8 +270,7 @@ def readTrace(self, tr_ind): return tr - def writeTrace(self, fname, sciline, refline, interp_trace, trace1024, - tr_ind, a2disp_ind): + def writeTrace(self, fname, sciline, refline, interp_trace, trace1024, tr_ind, a2disp_ind): """ The 'writeTrace' method performs the following steps: @@ -357,13 +360,13 @@ def generateTrace(self, data, kwinfo, tracecen=0.0, wind=None): sizex, sizey = data.shape subim_size = 40 - y1 = int(_tracecen - subim_size/2.) - y2 = int(_tracecen + subim_size/2.) + y1 = int(_tracecen - subim_size / 2.) + y2 = int(_tracecen + subim_size / 2.) if y1 < 0: y1 = 0 - if y2 > (sizex -1): + if y2 > (sizex - 1): y2 = sizex - 1 - specimage = data[y1:y2+1, :] + specimage = data[y1:y2 + 1, :] smoytrace = self.gFitTrace(specimage, y1, y2) yshift = int(np.median(smoytrace) - 20) y1 = y1 + yshift @@ -372,9 +375,9 @@ def generateTrace(self, data, kwinfo, tracecen=0.0, wind=None): y1 = 0 if y2 > sizex: y2 = sizex - specimage = data[y1:y2+1, :] + specimage = data[y1:y2 + 1, :] smoytrace = self.gFitTrace(specimage, y1, y2) - med11smoytrace = ni.median_filter(smoytrace, 11) + med11smoytrace = ni.median_filter(smoytrace, 11, output=float) med11smoytrace[0] = med11smoytrace[2] diffmed = abs(smoytrace - med11smoytrace) tolerence = 3 * np.median(abs(smoytrace[wind] - med11smoytrace[wind])) @@ -387,7 +390,7 @@ def generateTrace(self, data, kwinfo, tracecen=0.0, wind=None): # Convolve with a gaussian to smooth it. fwhm = 10. sigma = fwhm / 2.355 - gaussconvxsmoytrace = ni.gaussian_filter1d(smoytrace, sigma) + gaussconvxsmoytrace = ni.gaussian_filter1d(smoytrace, sigma, output=float) # Compute the trace center as the median of the pixels # with nonzero weights. @@ -406,14 +409,27 @@ def gFitTrace(self, specimage, y1, y2): """ sizex, sizey = specimage.shape - smoytrace = np.zeros(sizey).astype(np.float64) - boxcar_kernel = signal.boxcar(3) / 3.0 + smoytrace = np.zeros(sizey, dtype=float) + boxcar_kernel = boxcar(3) / 3.0 + fitter = fitting.LevMarLSQFitter() - for c in np.arange(sizey): + fit_info = [] + for c in range(sizey): col = specimage[:, c] col = col - np.median(col) - smcol = ni.convolve(col, boxcar_kernel).astype(np.float64) - fit = gfit.gfit1d(smcol, quiet=1, maxiter=15) - smoytrace[c] = fit.params[1] + smcol = ni.convolve(col, boxcar_kernel, output=float) + x = np.arange(len(smcol), dtype=float) + gauss = models.Gaussian1D(amplitude=smcol.max(), mean=x[smcol.argmax()]) + with warnings.catch_warnings(record=True) as w: + fit = fitter(gauss, x, smcol, maxiter=150) + if fitter.fit_info['ierr'] > 4: + fit_info.append(fitter.fit_info['message']) + smoytrace[c] = fit.mean.value + + if len(fit_info) > 0: + # Aggregate all fit_info['message'] into a single warning: + fit_info = sorted(Counter(fit_info).most_common(1024), key=lambda x: -x[1]) + fit_info = [f'{y[1]} x "{y[0]}"' for y in fit_info] + warnings.warn('Astropy fitting caused some warnings. Associated fit_info:\n' + '\n'.join(fit_info)) return np.array(smoytrace) From 3827d8097171c53846b21ae99c3d98b634eb097c Mon Sep 17 00:00:00 2001 From: Sean Lockwood Date: Wed, 29 Jan 2025 12:40:59 -0500 Subject: [PATCH 2/4] Updated stdout of test results. Changes was of lower significance. --- tests/test_mktrace.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_mktrace.py b/tests/test_mktrace.py index 5dd5a817..61cbe38b 100644 --- a/tests/test_mktrace.py +++ b/tests/test_mktrace.py @@ -28,9 +28,9 @@ def test_mktrace_t1(self, capsys): # Compare results captured = capsys.readouterr() - assert captured.out == "Traces were rotated by 0.0066273929 degrees \n" \ + assert captured.out == "Traces were rotated by 0.0066273851 degrees \n" \ "\n" \ - "trace is centered on row 509.2475494293\n" + "Trace is centered on row 509.2475494292\n" outputs = [("o45a03010_crj_1dt.fits", "o45a03010_crj_1dt_ref.fits"), ("o45a03010_crj_1dt_interp.fits", "o45a03010_crj_1dt_interp_ref.fits"), @@ -57,9 +57,9 @@ def test_mktrace_t2(self, capsys): # Compare results captured = capsys.readouterr() - assert captured.out == "Traces were rotated by 0.0317233207 degrees \n" \ + assert captured.out == "Traces were rotated by 0.0317182356 degrees \n" \ "\n" \ - "trace is centered on row 893.7059688464\n" + "Trace is centered on row 893.7059707738\n" outputs = [("o8pp31020_crj_1dt.fits", "o8pp31020_crj_1dt_ref.fits"), ("o8pp31020_crj_1dt_interp.fits", "o8pp31020_crj_1dt_interp_ref.fits"), From 2307165aebbc989194e89834008fd7f1d46aec4a Mon Sep 17 00:00:00 2001 From: Sean Lockwood Date: Mon, 17 Feb 2025 16:02:51 -0500 Subject: [PATCH 3/4] Weighting Gaussian1D fitting by Hann window to bias towards locations in the middle of the fit region. Increased maxiter to 250. Simplified scipy imports. Updated user warning stacklevel to be more useful. --- stistools/mktrace.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/stistools/mktrace.py b/stistools/mktrace.py index c64b0111..a684b93a 100644 --- a/stistools/mktrace.py +++ b/stistools/mktrace.py @@ -7,10 +7,7 @@ from astropy.io import fits from astropy.modeling import models, fitting from scipy import ndimage as ni -try: - from scipy.signal import boxcar -except ImportError: - from scipy.signal.windows import boxcar +from scipy.signal.windows import boxcar, hann from stsci.tools import linefit from stsci.tools import fileutil as fu @@ -413,6 +410,10 @@ def gFitTrace(self, specimage, y1, y2): boxcar_kernel = boxcar(3) / 3.0 fitter = fitting.LevMarLSQFitter() + # Bias 1D Gaussian fit toward the center of the range: + window = hann(specimage.shape[0], sym=True) + window /= window.sum() + fit_info = [] for c in range(sizey): col = specimage[:, c] @@ -421,7 +422,7 @@ def gFitTrace(self, specimage, y1, y2): x = np.arange(len(smcol), dtype=float) gauss = models.Gaussian1D(amplitude=smcol.max(), mean=x[smcol.argmax()]) with warnings.catch_warnings(record=True) as w: - fit = fitter(gauss, x, smcol, maxiter=150) + fit = fitter(gauss, x, smcol, weights=window, maxiter=250) if fitter.fit_info['ierr'] > 4: fit_info.append(fitter.fit_info['message']) smoytrace[c] = fit.mean.value @@ -430,6 +431,7 @@ def gFitTrace(self, specimage, y1, y2): # Aggregate all fit_info['message'] into a single warning: fit_info = sorted(Counter(fit_info).most_common(1024), key=lambda x: -x[1]) fit_info = [f'{y[1]} x "{y[0]}"' for y in fit_info] - warnings.warn('Astropy fitting caused some warnings. Associated fit_info:\n' + '\n'.join(fit_info)) + warnings.warn('Astropy fitting caused some warnings. Associated fit_info:\n' + '\n'.join(fit_info), + stacklevel=2) return np.array(smoytrace) From 49f42e9f0ecb40168457f4b451856a8ef49a4c13 Mon Sep 17 00:00:00 2001 From: Sean Lockwood Date: Mon, 17 Feb 2025 18:32:47 -0500 Subject: [PATCH 4/4] =?UTF-8?q?Make=20Hann=20window=20non-zero=20on=20edge?= =?UTF-8?q?s=20to=20avoid=20non-finite=20fitting=20errors.=20=20Perform=20?= =?UTF-8?q?a=20global=20initial=20guess=20for=20the=20trace=20amplitude=20?= =?UTF-8?q?(scalar=20median)=20and=20peak=20position=20(3=CF=83-clipped=20?= =?UTF-8?q?linear=20fit);=20default=20to=20scalar=20median=20peak=20positi?= =?UTF-8?q?on=20when=20this=20fails.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- stistools/mktrace.py | 26 +++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/stistools/mktrace.py b/stistools/mktrace.py index a684b93a..82be5ff8 100644 --- a/stistools/mktrace.py +++ b/stistools/mktrace.py @@ -6,6 +6,7 @@ import numpy as np from astropy.io import fits from astropy.modeling import models, fitting +from astropy.stats import sigma_clip from scipy import ndimage as ni from scipy.signal.windows import boxcar, hann @@ -410,8 +411,22 @@ def gFitTrace(self, specimage, y1, y2): boxcar_kernel = boxcar(3) / 3.0 fitter = fitting.LevMarLSQFitter() + linear_outlier_fitter = fitting.FittingWithOutlierRemoval( + fitting.LinearLSQFitter(), + sigma_clip, + sigma=3., + niter=5) + + # Some initial conditions common across all columns: + amplitude_guess = float(np.nanmedian(specimage.max(axis=0))) + mean_guess, _ = linear_outlier_fitter( + models.Linear1D(), + np.arange(specimage.shape[1], dtype=float), + specimage.argmax(axis=0)) + mean_guess_scalar = float(np.nanmedian(specimage.argmax(axis=0))) + # Bias 1D Gaussian fit toward the center of the range: - window = hann(specimage.shape[0], sym=True) + window = hann(specimage.shape[0] + 2, sym=True)[1:-1] # Avoid zeros on edges window /= window.sum() fit_info = [] @@ -420,9 +435,14 @@ def gFitTrace(self, specimage, y1, y2): col = col - np.median(col) smcol = ni.convolve(col, boxcar_kernel, output=float) x = np.arange(len(smcol), dtype=float) - gauss = models.Gaussian1D(amplitude=smcol.max(), mean=x[smcol.argmax()]) + gauss = models.Gaussian1D(amplitude=amplitude_guess, mean=mean_guess(c), stddev=1.) with warnings.catch_warnings(record=True) as w: - fit = fitter(gauss, x, smcol, weights=window, maxiter=250) + try: + fit = fitter(gauss, x, smcol, weights=window, maxiter=250) + except fitting.NonFiniteValueError: + # Slightly different mean guess: + gauss = models.Gaussian1D(amplitude=amplitude_guess, mean=mean_guess_scalar, stddev=1.) + fit = fitter(gauss, x, smcol, weights=window, maxiter=250) if fitter.fit_info['ierr'] > 4: fit_info.append(fitter.fit_info['message']) smoytrace[c] = fit.mean.value