From 884af8ff2be5053ec7ba0e906bd225545a7f0df8 Mon Sep 17 00:00:00 2001 From: rcooke Date: Mon, 26 Aug 2024 20:19:45 +0100 Subject: [PATCH] revert to old masking --- deprecated/arc_old.py | 50 +++++++++++++++++++++++++++++++++ pypeit/core/arc.py | 51 ---------------------------------- pypeit/core/findobj_skymask.py | 10 +------ pypeit/tests/test_arc.py | 40 -------------------------- 4 files changed, 51 insertions(+), 100 deletions(-) diff --git a/deprecated/arc_old.py b/deprecated/arc_old.py index edecb94dd2..79b58ded49 100644 --- a/deprecated/arc_old.py +++ b/deprecated/arc_old.py @@ -855,3 +855,53 @@ def saturation_mask(a, satlevel): return mask.astype(int) +def mask_around_peaks(spec, inbpm): + """ + Find peaks in the input spectrum and mask pixels around them. + + All pixels to the left and right of a peak is masked until + a pixel has a lower value than the adjacent pixel. At this + point, we assume that spec has reached the noise level. + + Parameters + ---------- + spec: `numpy.ndarray`_ + Spectrum (1D array) in counts + inbpm: `numpy.ndarray`_ + Input bad pixel mask + + Returns + ------- + outbpm: `numpy.ndarray`_ + Bad pixel mask with pixels around peaks masked + """ + # Find the peak locations + pks = detect_peaks(spec) + + # Initialise some useful variables and the output bpm + xarray = np.arange(spec.size) + specdiff = np.append(np.diff(spec), 0.0) + outbpm = inbpm.copy() + + # Loop over the peaks and mask pixels around them + for i in range(len(pks)): + # Find all pixels to the left of the peak that are above the noise level + wl = np.where((xarray <= pks[i]) & (specdiff > 0.0))[0] + ww = (pks[i]-wl)[::-1] + # Find the first pixel to the left of the peak that is below the noise level + nmask = np.where(np.diff(ww) > 1)[0] + if nmask.size != 0 and nmask[0] > 5: + # Mask all pixels to the left of the peak + mini = max(0,wl.size-nmask[0]-1) + outbpm[wl[mini]:pks[i]] = True + # Find all pixels to the right of the peak that are above the noise level + ww = np.where((xarray >= pks[i]) & (specdiff < 0.0))[0] + # Find the first pixel to the right of the peak that is below the noise level + nmask = np.where(np.diff(ww) > 1)[0] + if nmask.size != 0 and nmask[0] > 5: + # Mask all pixels to the right of the peak + maxi = min(nmask[0], ww.size) + outbpm[pks[i]:ww[maxi]+2] = True + # Return the output bpm + return outbpm + diff --git a/pypeit/core/arc.py b/pypeit/core/arc.py index ca2df6ad71..e41fc7fc8b 100644 --- a/pypeit/core/arc.py +++ b/pypeit/core/arc.py @@ -1070,57 +1070,6 @@ def detect_lines(censpec, sigdetect=5.0, fwhm=4.0, fit_frac_fwhm=1.25, input_thr return tampl_true, tampl, tcent, twid, centerr, ww, arc, nsig -def mask_around_peaks(spec, inbpm): - """ - Find peaks in the input spectrum and mask pixels around them. - - All pixels to the left and right of a peak is masked until - a pixel has a lower value than the adjacent pixel. At this - point, we assume that spec has reached the noise level. - - Parameters - ---------- - spec: `numpy.ndarray`_ - Spectrum (1D array) in counts - inbpm: `numpy.ndarray`_ - Input bad pixel mask - - Returns - ------- - outbpm: `numpy.ndarray`_ - Bad pixel mask with pixels around peaks masked - """ - # Find the peak locations - pks = detect_peaks(spec) - - # Initialise some useful variables and the output bpm - xarray = np.arange(spec.size) - specdiff = np.append(np.diff(spec), 0.0) - outbpm = inbpm.copy() - - # Loop over the peaks and mask pixels around them - for i in range(len(pks)): - # Find all pixels to the left of the peak that are above the noise level - wl = np.where((xarray <= pks[i]) & (specdiff > 0.0))[0] - ww = (pks[i]-wl)[::-1] - # Find the first pixel to the left of the peak that is below the noise level - nmask = np.where(np.diff(ww) > 1)[0] - if nmask.size != 0 and nmask[0] > 5: - # Mask all pixels to the left of the peak - mini = max(0,wl.size-nmask[0]-1) - outbpm[wl[mini]:pks[i]] = True - # Find all pixels to the right of the peak that are above the noise level - ww = np.where((xarray >= pks[i]) & (specdiff < 0.0))[0] - # Find the first pixel to the right of the peak that is below the noise level - nmask = np.where(np.diff(ww) > 1)[0] - if nmask.size != 0 and nmask[0] > 5: - # Mask all pixels to the right of the peak - maxi = min(nmask[0], ww.size) - outbpm[pks[i]:ww[maxi]+2] = True - # Return the output bpm - return outbpm - - def fit_arcspec(xarray, yarray, pixt, fitp): """ Fit a series of pre-identified arc spectrum lines. diff --git a/pypeit/core/findobj_skymask.py b/pypeit/core/findobj_skymask.py index e100df877e..6e4698c844 100644 --- a/pypeit/core/findobj_skymask.py +++ b/pypeit/core/findobj_skymask.py @@ -1763,16 +1763,8 @@ def objs_in_slit(image, ivar, thismask, slit_left, slit_righ, gpm_smash = npix_smash > 0.3*nsmash flux_sum_smash = np.sum((image_rect*gpm_sigclip)[find_min_max_out[0]:find_min_max_out[1]], axis=0) flux_smash = flux_sum_smash*gpm_smash/(npix_smash + (npix_smash == 0.0)) - # Mask around the peaks - peak_bpm = arc.mask_around_peaks(flux_smash, np.logical_not(gpm_smash)) - # Check that the peak mask is not empty - if np.sum(np.logical_not(peak_bpm))/peak_bpm.size < 0.05: # 95% masked! - msgs.warn('Peak masking was too aggressive. Attempting to sigma clip.') - smash_mask = np.logical_not(gpm_smash) - else: - smash_mask = np.logical_not(gpm_smash) | peak_bpm flux_smash_mean, flux_smash_med, flux_smash_std = astropy.stats.sigma_clipped_stats( - flux_smash, mask=smash_mask, sigma_lower=3.0, sigma_upper=3.0) + flux_smash, mask=np.logical_not(gpm_smash), sigma_lower=3.0, sigma_upper=3.0) flux_smash_recen = flux_smash - flux_smash_med # Return if none found and no hand extraction diff --git a/pypeit/tests/test_arc.py b/pypeit/tests/test_arc.py index a112f30e4d..3cde322d37 100644 --- a/pypeit/tests/test_arc.py +++ b/pypeit/tests/test_arc.py @@ -15,45 +15,5 @@ def test_detect_lines(): = arc.detect_lines(arx_sky.flux.value) assert (len(arx_w) > 3275) - -def test_mask_around_peaks(): - # Generate a fake spectrum - np.random.seed(0) - npks = 5 # Number of peaks to generate - npix = 20 # Number of pixels in each segment - sigma = 2.0 # Fluctuations - linewid = 2.0 # Width of the peaks - noiselev = [np.random.normal(0.0, sigma, npix) for i in range(npks+1)] - ampl = np.random.randint(100*sigma, 200*sigma, npks) - spec = np.array([]) - outbpm = np.array([], dtype=bool) - for i in range(npks): - spec = np.append(spec, noiselev[i]) - outbpm = np.append(outbpm, np.zeros(npix, dtype=bool)) - gauss = ampl[i]*np.exp(-0.5*(np.arange(npix)-npix/2)**2/linewid**2) - spec = np.append(spec, np.random.normal(gauss, sigma)) - if i == 0 or i == 1: - tmpbpm = np.array(2*[False] + 16*[True] + 2*[False]) - elif i == 2: - tmpbpm = np.array(1*[False] + 18*[True] + 1*[False]) - elif i == 3: - tmpbpm = np.array(3*[False] + 15*[True] + 2*[False]) - elif i == 4: - tmpbpm = np.array(20*[True]) - else: - assert (False, "Input test data have changed") - outbpm = np.append(outbpm, tmpbpm.copy()) - outbpm = np.append(outbpm, np.zeros(npix, dtype=bool)) - spec = np.append(spec, noiselev[npks]) - inbpm = np.zeros(len(spec), dtype=bool) - # This is the expected bpm - bpm = arc.mask_around_peaks(spec, inbpm=inbpm) - # Check that all peaks were identified and masked - ind = arc.detect_peaks(bpm) - assert (len(ind) == npks) - # Check the BPM matches exactly - assert (np.array_equal(bpm, outbpm)) - - # Many more functions in pypeit.core.arc that need tests!