Merge pull request #31 from lsst/tickets/DM-38622

DM-38622: Replace deprecated make_source_mask function

requirements.txt

@@ -5,7 +5,7 @@ pandas
 llvmlite #==0.29.0
 numba #==0.44.1
 astropy>=3.2
-photutils>=0.7
+photutils>=1.7.0
 astroquery
 coloredlogs
 scikit-image

spectractor/extractor/background.py

@@ -10,8 +10,9 @@
 
 from astropy.stats import SigmaClip
 from photutils import Background2D, SExtractorBackground
-from photutils.segmentation import make_source_mask
+from photutils.segmentation import SegmentationImage, detect_threshold, detect_sources
 
+from scipy import ndimage
 from scipy.signal import medfilt2d
 from scipy.interpolate import RegularGridInterpolator
 
@@ -43,6 +44,58 @@ def remove_image_background_sextractor(data, sigma=3.0, box_size=(50, 50), filte
     return data_wo_bkg
 
 
+def make_source_mask(data, nsigma, npixels, mask=None, sigclip_sigma=3.0,
+                     sigclip_iters=5, dilate_size=11):
+    """
+    Make a source mask using source segmentation and binary dilation.
+
+    This is a slight stripped down version of the method which was removed from
+    photutils in 1.7.0.
+
+    Parameters
+    ----------
+    data : 2D `~numpy.ndarray`
+        The 2D array of the image.
+    nsigma : float
+        The number of standard deviations per pixel above the ``background``
+        for which to consider a pixel as possibly being part of a source.
+    npixels : int
+        The minimum number of connected pixels, each greater than
+        ``threshold``, that an object must have to be detected. ``npixels``
+        must be a positive integer.
+    mask : 2D bool `~numpy.ndarray`, optional
+        A boolean mask with the same shape as ``data``, where a `True` value
+        indicates the corresponding element of ``data`` is masked. Masked
+        pixels are ignored when computing the image background statistics.
+    sigclip_sigma : float, optional
+        The number of standard deviations to use as the clipping limit when
+        calculating the image background statistics.
+    sigclip_iters : int, optional
+        The maximum number of iterations to perform sigma clipping, or `None`
+        to clip until convergence is achieved (i.e., continue until the last
+        iteration clips nothing) when calculating the image background
+        statistics.
+    dilate_size : int, optional
+        The size of the square array used to dilate the segmentation image.
+
+    Returns
+    -------
+    mask : 2D bool `~numpy.ndarray`
+        A 2D boolean image containing the source mask.
+    """
+    sigma_clip = SigmaClip(sigma=sigclip_sigma, maxiters=sigclip_iters)
+    threshold = detect_threshold(data, nsigma, background=None, error=None,
+                                 mask=mask, sigma_clip=sigma_clip)
+
+    segm = detect_sources(data, threshold, npixels)
+    if segm is None:
+        return np.zeros(data.shape, dtype=bool)
+
+    footprint = np.ones((dilate_size, dilate_size))
+    # Replace with size= when photutils>=1.7 is enforced in rubin-env
+    return segm.make_source_mask(footprint=footprint)
+
+
 def extract_spectrogram_background_fit1D(data, err, deg=1, ws=(20, 30), pixel_step=1, sigma=5):
     """
     Fit a polynomial background slice per slice along the x axis,
@@ -200,7 +253,7 @@ def extract_spectrogram_background_sextractor(data, err, ws=(20, 30), mask_signa
     Ny, Nx = data.shape
     if Dy_disp_axis is None:
         Dy_disp_axis = np.ones(Nx) * (Ny // 2)
-    
+
     # mask sources
     mask = make_source_mask(data, nsigma=3, npixels=5, dilate_size=11)
     mask += data == 0

tests/test_utils.py

@@ -0,0 +1,24 @@
+import matplotlib as mpl
+mpl.use('Agg')  # must be run first! But therefore requires noqa E402 on all other imports
+
+import pytest  # noqa: E402
+import numpy as np  # noqa: E402
+from numpy.testing import run_module_suite  # noqa: E402
+from photutils.datasets import make_4gaussians_image  # noqa: E402
+from spectractor.extractor.background import make_source_mask  # noqa: E402
+from photutils.utils.exceptions import NoDetectionsWarning  # noqa: E402
+
+
+def test_make_source_mask():
+    data = make_4gaussians_image()
+    mask1 = make_source_mask(data, 5, 10)
+    mask2 = make_source_mask(data, 5, 10, dilate_size=20)
+    assert np.count_nonzero(mask2) > np.count_nonzero(mask1)
+
+    with pytest.warns(NoDetectionsWarning):
+        mask = make_source_mask(data, 100, 100)
+        assert np.count_nonzero(mask) == 0
+
+
+if __name__ == "__main__":
+    run_module_suite()