use broadest Ricker filter for pinpointing peaks in data

GeminiDRSoftware · Mar 19, 2021 · 11b43e3 · 11b43e3
1 parent 2e77a97
commit 11b43e3
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 6 deletions.
diff --git a/geminidr/core/primitives_spect.py b/geminidr/core/primitives_spect.py
@@ -603,7 +603,7 @@ def determineDistortion(self, adinputs=None, **params):
                                 format(direction, extract_slice.start + 1, extract_slice.stop))
 
                     # Find peaks; convert width FWHM to sigma
-                    widths = 0.42466 * fwidth * np.arange(0.8, 1.21, 0.05)  # TODO!
+                    widths = 0.42466 * fwidth * np.arange(0.75, 1.26, 0.05)  # TODO!
                     initial_peaks, _ = tracing.find_peaks(data, widths, mask=mask & DQ.not_signal,
                                                           variance=variance, min_snr=min_snr)
                     log.stdinfo("Found {} peaks".format(len(initial_peaks)))
@@ -1162,7 +1162,7 @@ def determineWavelengthSolution(self, adinputs=None, **params):
                     arc_lines *= 0.1
 
                 # Find peaks; convert width FWHM to sigma
-                widths = 0.42466 * fwidth * np.arange(0.7, 1.2, 0.05)  # TODO!
+                widths = 0.42466 * fwidth * np.arange(0.75, 1.26, 0.05)  # TODO!
                 peaks, peak_snrs = tracing.find_peaks(data, widths, mask=mask & DQ.not_signal,
                                                       variance=variance, min_snr=min_snr,
                                                       min_sep=min_sep, reject_bad=False)

diff --git a/gempy/library/tracing.py b/gempy/library/tracing.py
@@ -437,17 +437,22 @@ def find_peaks(data, widths, mask=None, variance=None, min_snr=1, min_sep=1,
     filtered = signal._peak_finding._filter_ridge_lines(
         wavelet_transformed_data, ridge_lines, window_size=window_size,
         min_length=int(min_frac * len(widths)), min_snr=min_snr)
-
     peaks = sorted([x[1][0] for x in filtered])
 
+    # We need to find accurate peak positions from convolved data so we're
+    # not affected by noise or problems with broad, flat-topped features.
+    # There appears to be a bias with narrow Ricker transforms, so we use the
+    # broadest one for this purpose.
+    pinpoint_data = wavelet_transformed_data[-1]
+
     # If no variance is supplied we estimate S/N from pixel-to-pixel variations
     if variance is not None:
-        snr = np.divide(wavelet_transformed_data[0], np.sqrt(variance),
+        snr = np.divide(pinpoint_data, np.sqrt(variance),
                         out=np.zeros_like(data, dtype=np.float32),
                         where=variance > 0)
     else:
         sigma = sigma_clip(data[~mask], masked=False).std() / np.sqrt(2)
-        snr = wavelet_transformed_data[0] / sigma
+        snr = pinpoint_data / sigma
     peaks = [x for x in peaks if snr[x] > min_snr]
 
     # remove adjacent points
@@ -472,7 +477,7 @@ def find_peaks(data, widths, mask=None, variance=None, min_snr=1, min_sep=1,
     peaks = [x for x in peaks if np.sum(mask[int(x-edge):int(x+edge+1)]) == 0]
 
     # Clip the really noisy parts of the data and get more accurate positions
-    pinpoint_data = np.where(snr < 0.5, 0, wavelet_transformed_data[0])
+    pinpoint_data[snr < 0.5] = 0
     peaks = pinpoint_peaks(pinpoint_data, mask, peaks,
                            halfwidth=int(np.median(widths)+0.5))