diff --git a/CHANGES.rst b/CHANGES.rst index 13bce6bd0e..8969e0514c 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -294,6 +294,13 @@ New Features are divided into chunks that are processed concurrently, producing results identical to the single-threaded computation. [#2407] + - Significantly improved the performance of source deblending in + ``deblend_sources`` and ``SourceFinder``, producing identical results. + The multithreshold watershed markers are now built by compiled code. + Deblending is typically ~7-30 times faster, both for fields of many + small blended sources and for large segments with many markers. + [#2408] + - ``photutils.utils`` - Added a new ``DeblendWarning`` class, a subclass of astropy's @@ -1019,6 +1026,15 @@ API Changes are now computed for all sources at once in compiled code, so no progress bar is displayed and the keyword has no effect. [#2406] + - The ``deblend_sources`` and ``SourceFinder`` ``progress_bar`` + and ``n_processes`` keywords are now deprecated and will be removed + in version 4.0, and both keywords no longer have any effect. + Deblending is now dominated by compiled code, so no progress bar is + displayed, and the multiprocessing implementation has been removed + because its process startup and data-pickling overheads made it + slower than the serial implementation at any number of sources. + [#2408] + - ``photutils.utils`` - The ``ShepardIDWInterpolator`` ``ncoords`` attribute has been diff --git a/benchmarks/bench_deblend.py b/benchmarks/bench_deblend.py new file mode 100755 index 0000000000..c4a0b592c0 --- /dev/null +++ b/benchmarks/bench_deblend.py @@ -0,0 +1,550 @@ +#!/usr/bin/env python3 +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +Benchmarks for source deblending (deblend_sources). + +The benchmarks cover the two axes along which deblending is known to +be slow: + +* many sources: an image with a grid of blended Gaussian-source + pairs, sweeping the number of sources (per-source Python and numpy + call overhead dominates) + +* large sources: a single connected segment made of a broad Gaussian + envelope with superposed peaks, sweeping the segment size and the + number of peaks (per-level full-cutout work and the iterative + watershed contrast loop dominate) + +A per-stage breakdown of the single-source deblending pipeline +(multithresholding, marker building, watershed, contrast loop) and a +cProfile mode are also provided for bottleneck analysis. + +Run ``python benchmarks/bench_deblend.py --help`` to see the +available options. +""" + +import argparse +import cProfile +import pstats +import warnings +from functools import partial + +import numpy as np +from astropy.modeling.models import Gaussian2D +from astropy.stats import gaussian_fwhm_to_sigma +from bench_helpers import print_environment, time_best +from bench_segmentation import N_PIXELS, THRESHOLD, make_inputs + +from photutils.segmentation import detect_sources +from photutils.segmentation._deblend_reference import _SingleSourceDeblender +from photutils.segmentation.deblend import _DeblendParams, deblend_sources +from photutils.segmentation.utils import _make_binary_structure +from photutils.utils.exceptions import DeblendWarning + +BLEND_THRESHOLD = 0.5 +ENVELOPE_AMPLITUDE = 5.0 +PEAK_FWHM = 8.0 + + +def make_blended_image(size, n_peaks, *, amp_range=(3.0, 100.0), seed=0): + """ + Return an image containing a single large blended source. + + The source is a broad Gaussian envelope with ``n_peaks`` compact + Gaussian peaks superposed on it, so that detection at + ``BLEND_THRESHOLD`` yields one large connected segment that + deblending must split. + + Parameters + ---------- + size : int + The image size; the image is ``(size, size)``. The envelope + sigma is ``size / 6``, so the segment area scales with the + image area. + + n_peaks : int + The number of compact peaks placed within one envelope sigma + of the center. + + amp_range : tuple of float, optional + The (min, max) peak amplitudes. The amplitudes are + logarithmically spaced over this range. + + seed : int, optional + The random number generator seed. + + Returns + ------- + data : 2D `~numpy.ndarray` + The image. + """ + rng = np.random.default_rng(seed) + data = rng.normal(0.0, 0.01, (size, size)) + + cen = (size - 1) / 2.0 + sigma_env = size / 6.0 + yy, xx = np.mgrid[0:size, 0:size] + envelope = Gaussian2D(ENVELOPE_AMPLITUDE, cen, cen, sigma_env, + sigma_env) + data += envelope(xx, yy) + + sigma = PEAK_FWHM * gaussian_fwhm_to_sigma + half = int(np.ceil(4.0 * sigma)) + yy_cut, xx_cut = np.mgrid[0:2 * half + 1, 0:2 * half + 1] + amplitudes = np.geomspace(amp_range[1], amp_range[0], n_peaks) + radii = sigma_env * np.sqrt(rng.uniform(0.0, 1.0, n_peaks)) + angles = rng.uniform(0.0, 2.0 * np.pi, n_peaks) + for amplitude, radius, angle in zip(amplitudes, radii, angles, + strict=True): + xc = cen + radius * np.cos(angle) + yc = cen + radius * np.sin(angle) + x0 = int(xc) - half + y0 = int(yc) - half + model = Gaussian2D(amplitude, xc - x0, yc - y0, sigma, sigma) + data[y0:y0 + 2 * half + 1, + x0:x0 + 2 * half + 1] += model(xx_cut, yy_cut) + + return data + + +def make_blended_inputs(size, n_peaks, *, amp_range=(3.0, 100.0), seed=0): + """ + Return the image and segmentation image for a single large + blended source. + + Parameters + ---------- + size : int + The image size; the image is ``(size, size)``. + + n_peaks : int + The number of compact peaks. + + amp_range : tuple of float, optional + The (min, max) peak amplitudes. + + seed : int, optional + The random number generator seed. + + Returns + ------- + data : 2D `~numpy.ndarray` + The image. + + segm : `~photutils.segmentation.SegmentationImage` + The segmentation image containing a single label. + """ + data = make_blended_image(size, n_peaks, amp_range=amp_range, + seed=seed) + segm = detect_sources(data, BLEND_THRESHOLD, N_PIXELS) + if segm.n_labels != 1: + msg = (f'expected a single blended segment, got ' + f'{segm.n_labels}') + raise ValueError(msg) + return data, segm + + +def n_fallbacks(segm): + """ + Return the number of deblending mode fallbacks recorded in a + deblended segmentation image. + + Parameters + ---------- + segm : `~photutils.segmentation.SegmentationImage` + The deblended segmentation image. + + Returns + ------- + result : int + The total number of input labels whose deblending mode fell + back to linear. + """ + return (len(segm.info.get('nonposmin_labels', ())) + + len(segm.info.get('n_markers_labels', ()))) + + +def bench_many_sources(*, n_sources_sweep=(500, 1000, 2000, 4000), + repeats=3, seed=0): + """ + Benchmark deblend_sources versus the number of sources. + + Each image contains a grid of blended Gaussian-source pairs, so + the number of segments is half the number of sources and every + segment deblends into two sources. The detect_sources time is + included as a reference point. + + Parameters + ---------- + n_sources_sweep : tuple of int, optional + The numbers of Gaussian sources. + + repeats : int, optional + The number of repeats for each timing (best time is kept). + + seed : int, optional + The random number generator seed. + """ + print('\n== deblend_sources: many small sources ==') + header = (f'{"benchmark":>24}{"segments":>10}{"time":>12}' + f'{"ms/segment":>12}{"n_labels":>10}') + + for n_sources in n_sources_sweep: + _, data, segm = make_inputs(n_sources, seed=seed) + + print(f'\n-- {n_sources} sources, {segm.n_labels} segments, ' + f'{data.shape[0]}x{data.shape[1]} image --') + print(header) + + bench = partial(detect_sources, data, THRESHOLD, N_PIXELS) + t_best = time_best(bench, repeats=repeats) + per_segment = 1000.0 * t_best / segm.n_labels + print(f'{"detect_sources (ref)":>24}{segm.n_labels:>10}' + f'{f"{t_best:.4f}s":>12}{per_segment:>12.3f}' + f'{segm.n_labels:>10}') + + for mode in ('linear', 'exponential', 'sinh'): + bench = partial(deblend_sources, data, segm, N_PIXELS, + mode=mode) + segm_deblended = bench() + t_best = time_best(bench, repeats=repeats) + per_segment = 1000.0 * t_best / segm.n_labels + name = f'mode={mode}' + print(f'{name:>24}{segm.n_labels:>10}{f"{t_best:.4f}s":>12}' + f'{per_segment:>12.3f}{segm_deblended.n_labels:>10}') + + +def bench_large_source(*, size_sweep=(250, 500, 1000, 2000), n_peaks=8, + repeats=3, seed=0): + """ + Benchmark deblend_sources versus the size of a single segment. + + Each image contains one connected segment (a Gaussian envelope + with bright peaks) whose area scales with the image area. The + contrast is set to 0 so that no markers are removed and the sweep + measures the pure multithreshold plus watershed scaling. + + Parameters + ---------- + size_sweep : tuple of int, optional + The image sizes. + + n_peaks : int, optional + The number of compact peaks in the segment. + + repeats : int, optional + The number of repeats for each timing (best time is kept). + + seed : int, optional + The random number generator seed. + """ + print(f'\n== deblend_sources: single large segment ({n_peaks} peaks, ' + 'contrast=0) ==') + print(f'{"benchmark":>28}{"seg area":>10}{"time":>12}{"n_labels":>10}') + + for size in size_sweep: + data, segm = make_blended_inputs(size, n_peaks, + amp_range=(50.0, 100.0), + seed=seed) + area = int(segm.areas[0]) + for mode in ('linear', 'exponential'): + bench = partial(deblend_sources, data, segm, N_PIXELS, + mode=mode, contrast=0.0) + segm_deblended = bench() + t_best = time_best(bench, repeats=repeats) + name = f'size={size}, mode={mode}' + print(f'{name:>28}{area:>10}{f"{t_best:.4f}s":>12}' + f'{segm_deblended.n_labels:>10}') + + +def bench_many_peaks(*, size=1000, n_peaks_sweep=(10, 25, 50, 100), + contrast_sweep=(0.0, 0.001, 0.01, 0.03), + repeats=3, seed=0): + """ + Benchmark deblend_sources versus the number of peaks in one + segment. + + The contrast criterion applies to the flux in each watershed + basin (which includes a share of the envelope flux), so larger + contrast values remove more markers. Each removal iteration + re-runs the watershed over the full segment, so comparing the + contrast=0 row (a single watershed call) to the larger-contrast + rows isolates the cost of the iterative marker-removal loop. The + number of watershed calls is at most the difference between the + contrast=0 n_labels and the row's n_labels, plus one (batched + removal may use fewer). + + Parameters + ---------- + size : int, optional + The image size. + + n_peaks_sweep : tuple of int, optional + The numbers of peaks. + + contrast_sweep : tuple of float, optional + The contrast values. + + repeats : int, optional + The number of repeats for each timing (best time is kept). + + seed : int, optional + The random number generator seed. + """ + print(f'\n== deblend_sources: many peaks in one segment ' + f'({size}x{size} image, mode=exponential) ==') + print(f'{"benchmark":>32}{"time":>12}{"n_labels":>10}' + f'{"fallbacks":>10}') + + for n_peaks in n_peaks_sweep: + data, segm = make_blended_inputs(size, n_peaks, seed=seed) + for contrast in contrast_sweep: + bench = partial(deblend_sources, data, segm, N_PIXELS, + contrast=contrast) + segm_deblended = bench() + t_best = time_best(bench, repeats=repeats) + name = f'n_peaks={n_peaks}, contrast={contrast}' + print(f'{name:>32}{f"{t_best:.4f}s":>12}' + f'{segm_deblended.n_labels:>10}' + f'{n_fallbacks(segm_deblended):>10}') + + +def bench_stages(*, size=1000, n_peaks=25, mode='exponential', + contrast=0.001, n_levels=32, repeats=3, seed=0): + """ + Benchmark the stages of the single-source deblending pipeline. + + The stages are timed on the segment cutout of a single large + blended source, using the private _SingleSourceDeblender class: + + * constructor: the segment mask and min/max/sum reductions + * multithreshold: the ``n_levels`` per-level detection passes of + the reference marker construction + * make_markers: the multithreshold levels, the level quantization, + and the compiled component-tree kernel + * watershed: a single watershed call over the cutout + * apply_watershed: the watershed contrast loop (one watershed + call per removed marker) + * deblend_source: the full pipeline + + Parameters + ---------- + size : int, optional + The image size. + + n_peaks : int, optional + The number of compact peaks in the segment. + + mode : str, optional + The mode used for spacing the multithreshold levels. + + contrast : float, optional + The deblending contrast criterion. + + n_levels : int, optional + The number of multithreshold levels. + + repeats : int, optional + The number of repeats for each timing (best time is kept). + + seed : int, optional + The random number generator seed. + """ + data, segm = make_blended_inputs(size, n_peaks, seed=seed) + label = segm.labels[0] + slc = segm.slices[0] + cutout = data[slc] + segment_cutout = segm.data[slc] + footprint = _make_binary_structure(2, 8) + params = _DeblendParams(N_PIXELS, footprint, n_levels, contrast, mode) + + def _make_deblender(): + return _SingleSourceDeblender(cutout, segment_cutout, label, + params) + + deblender = _make_deblender() + markers = deblender.make_markers() + n_markers = len(np.unique(markers[markers > 0])) + final = deblender.apply_watershed(markers) + n_final = len(np.unique(final[final > 0])) + n_watershed = n_markers - n_final + 1 + + from photutils.segmentation._deblend_watershed import deblend_watershed + + data_neg = np.ascontiguousarray(-cutout, dtype=np.float64) + connectivity = 8 if footprint[0, 0] else 4 + + def _run_constructor(): + _make_deblender() + + def _run_multithreshold(): + _make_deblender().multithreshold() + + def _run_make_markers(): + _make_deblender().make_markers() + + def _run_watershed(): + deblend_watershed(data_neg, markers, deblender.segment_mask, + connectivity) + + def _run_apply_watershed(): + deblender.apply_watershed(markers) + + def _run_deblend_source(): + _make_deblender().deblend_source() + + print(f'\n== single-source pipeline stages ({size}x{size} image, ' + f'cutout {cutout.shape[0]}x{cutout.shape[1]}, ' + f'{n_peaks} peaks, mode={mode}, contrast={contrast}, ' + f'n_levels={n_levels}) ==') + print(f'{n_markers} markers, {n_final} final labels, ' + f'<={n_watershed} watershed calls') + print(f'{"stage":>36}{"time":>12}') + + benchmarks = [ + ('constructor (mask + min/max/sum)', _run_constructor), + (f'multithreshold ({n_levels} per-level detects)', + _run_multithreshold), + ('make_markers (component tree)', _run_make_markers), + ('watershed (single call)', _run_watershed), + ('apply_watershed (contrast loop)', _run_apply_watershed), + ('deblend_source (full)', _run_deblend_source), + ] + for name, func in benchmarks: + t_best = time_best(func, repeats=repeats) + print(f'{name:>36}{f"{t_best:.4f}s":>12}') + + +def profile_case(name, func, *, limit=20): + """ + Profile a callable with cProfile and print the top functions. + + Parameters + ---------- + name : str + The name of the profiled case. + + func : callable + The zero-argument callable to profile. + + limit : int, optional + The number of functions to print. + """ + profiler = cProfile.Profile() + profiler.enable() + func() + profiler.disable() + + print(f'\n== cProfile: {name} ==') + stats = pstats.Stats(profiler) + stats.sort_stats('tottime') + stats.print_stats(limit) + + +def bench_profile(*, n_sources=2000, size=1000, n_peaks=25, seed=0): + """ + Profile deblend_sources for the many-source and large-source + scenarios. + + Parameters + ---------- + n_sources : int, optional + The number of sources for the many-source case. + + size : int, optional + The image size for the large-source case. + + n_peaks : int, optional + The number of peaks for the large-source case. + + seed : int, optional + The random number generator seed. + """ + _, data, segm = make_inputs(n_sources, seed=seed) + profile_case( + f'many small sources ({segm.n_labels} segments)', + partial(deblend_sources, data, segm, N_PIXELS)) + + data, segm = make_blended_inputs(size, n_peaks, seed=seed) + profile_case( + f'single large segment ({size}x{size}, {n_peaks} peaks)', + partial(deblend_sources, data, segm, N_PIXELS)) + + +def parse_int_list(text): + """ + Parse a comma-separated list of positive integers. + + Parameters + ---------- + text : str + The comma-separated integers (e.g., ``'250,500,1000'``). + + Returns + ------- + result : list of int + The parsed integers. + """ + values = [int(item) for item in text.split(',')] + if any(value < 1 for value in values): + msg = 'values must be positive integers' + raise ValueError(msg) + return values + + +def main(): + """ + Run the source deblending benchmarks. + """ + parser = argparse.ArgumentParser( + description='Benchmarks for source deblending.') + parser.add_argument('--n-sources', type=parse_int_list, + default=[500, 1000, 2000, 4000], + help='comma-separated source counts for the ' + 'many-source benchmark ' + '(default: 500,1000,2000,4000)') + parser.add_argument('--sizes', type=parse_int_list, + default=[250, 500, 1000, 2000], + help='comma-separated image sizes for the ' + 'large-source benchmark ' + '(default: 250,500,1000,2000)') + parser.add_argument('--n-peaks', type=parse_int_list, + default=[10, 25, 50, 100], + help='comma-separated peak counts for the ' + 'many-peak benchmark ' + '(default: 10,25,50,100)') + parser.add_argument('--repeats', type=int, default=3, + help='number of repeats per timing; the best ' + 'time is reported (default: %(default)s)') + parser.add_argument('--seed', type=int, default=0, + help='random number generator seed ' + '(default: %(default)s)') + parser.add_argument('--which', default='all', + choices=['all', 'many', 'large', 'peaks', + 'stages', 'profile'], + help='which benchmark to run ' + '(default: %(default)s)') + args = parser.parse_args() + + warnings.filterwarnings('ignore', category=DeblendWarning) + print_environment() + + if args.which in ('all', 'many'): + bench_many_sources(n_sources_sweep=args.n_sources, + repeats=args.repeats, seed=args.seed) + if args.which in ('all', 'large'): + bench_large_source(size_sweep=args.sizes, repeats=args.repeats, + seed=args.seed) + if args.which in ('all', 'peaks'): + bench_many_peaks(n_peaks_sweep=args.n_peaks, + repeats=args.repeats, seed=args.seed) + if args.which in ('all', 'stages'): + bench_stages(repeats=args.repeats, seed=args.seed) + bench_stages(contrast=0.03, repeats=args.repeats, + seed=args.seed) + if args.which == 'profile': + bench_profile(seed=args.seed) + + +if __name__ == '__main__': + main() diff --git a/benchmarks/bench_segmentation.py b/benchmarks/bench_segmentation.py index 6a1a0eef4d..76a3eb36cb 100755 --- a/benchmarks/bench_segmentation.py +++ b/benchmarks/bench_segmentation.py @@ -34,8 +34,7 @@ detect_sources, detect_threshold, make_2dgaussian_kernel) from photutils.utils import circular_footprint -from photutils.utils._optional_deps import (HAS_RASTERIO, HAS_SHAPELY, - HAS_SKIMAGE) +from photutils.utils._optional_deps import HAS_RASTERIO, HAS_SHAPELY FWHM = 4.0 THRESHOLD = 5.0 @@ -189,10 +188,6 @@ def bench_deblend(*, n_sources=1000, process_counts=(1, 4), repeats=3, seed : int, optional The random number generator seed. """ - if not HAS_SKIMAGE: - print('\n== deblend_sources: skipped (scikit-image required) ==') - return - data, convolved_data, segm = make_inputs(n_sources, seed=seed) print(f'\n== deblend_sources ({n_sources} sources, ' @@ -202,9 +197,9 @@ def bench_deblend(*, n_sources=1000, process_counts=(1, 4), repeats=3, for mode in ('linear', 'exponential', 'sinh'): segm_deblended = deblend_sources(convolved_data, segm, N_PIXELS, - mode=mode, progress_bar=False) + mode=mode) bench = partial(deblend_sources, convolved_data, segm, N_PIXELS, - mode=mode, progress_bar=False) + mode=mode) t_best = time_best(bench, repeats=repeats) name = f'mode={mode}' print(f'{name:>36}{f"{t_best:.4f}s":>12}' @@ -214,8 +209,7 @@ def bench_deblend(*, n_sources=1000, process_counts=(1, 4), repeats=3, if n_processes == 1: continue bench = partial(deblend_sources, convolved_data, segm, N_PIXELS, - mode='exponential', n_processes=n_processes, - progress_bar=False) + mode='exponential', n_processes=n_processes) t_best = time_best(bench, repeats=repeats) name = f'mode=exponential, n_processes={n_processes}' print(f'{name:>36}{f"{t_best:.4f}s":>12}{"":>10}') @@ -244,12 +238,8 @@ def bench_finder(*, n_sources=1000, repeats=3, seed=0): f'{data.shape[0]}x{data.shape[1]} image) ==') print(f'{"benchmark":>24}{"time":>12}{"n_labels":>10}') - deblend_options = [False] - if HAS_SKIMAGE: - deblend_options.append(True) - for deblend in deblend_options: - finder = SourceFinder(n_pixels=N_PIXELS, deblend=deblend, - progress_bar=False) + for deblend in (False, True): + finder = SourceFinder(n_pixels=N_PIXELS, deblend=deblend) segm = finder(convolved_data, THRESHOLD) bench = partial(finder, convolved_data, THRESHOLD) t_best = time_best(bench, repeats=repeats) diff --git a/docs/user_guide/segmentation.rst b/docs/user_guide/segmentation.rst index 7062d5d58f..a8280b2333 100644 --- a/docs/user_guide/segmentation.rst +++ b/docs/user_guide/segmentation.rst @@ -148,14 +148,11 @@ and ``contrast``. ``n_levels`` is the number of multi-thresholding levels to use. ``contrast`` is the fraction of the total source flux that a local peak must have to be considered as a separate object. -Here's a simple example of source deblending: - -.. doctest-requires:: skimage +Here's a simple example of source deblending:: >>> from photutils.segmentation import deblend_sources >>> segment_map2 = deblend_sources(convolved_data, segment_map, - ... n_pixels=10, n_levels=32, contrast=0.001, - ... progress_bar=False) + ... n_pixels=10, n_levels=32, contrast=0.001) where ``segment_map`` is the :class:`~photutils.segmentation.SegmentationImage` that was @@ -190,8 +187,7 @@ deblended segmentation image: n_pixels = 10 segment_map = detect_sources(convolved_data, threshold, n_pixels=n_pixels) deblended_segment_map = deblend_sources(convolved_data, segment_map, - n_pixels=n_pixels, - progress_bar=False) + n_pixels=n_pixels) fig, ax = plt.subplots(figsize=(10, 6.5)) deblended_segment_map.imshow(ax=ax) @@ -224,8 +220,7 @@ Let's plot one of the deblended sources: n_pixels = 10 segment_map = detect_sources(convolved_data, threshold, n_pixels=n_pixels) deblended_segment_map = deblend_sources(convolved_data, segment_map, - n_pixels=n_pixels, - progress_bar=False) + n_pixels=n_pixels) fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(10, 4)) slc = (slice(273, 297), slice(425, 444)) @@ -250,12 +245,10 @@ is a convenience class that combines the functionality of `~photutils.segmentation.detect_sources` and `~photutils.segmentation.deblend_sources`. After defining the object with the desired detection and deblending parameters, you call it with -the background-subtracted (convolved) image and threshold: - -.. doctest-requires:: skimage +the background-subtracted (convolved) image and threshold:: >>> from photutils.segmentation import SourceFinder - >>> finder = SourceFinder(n_pixels=10, progress_bar=False) + >>> finder = SourceFinder(n_pixels=10) >>> segment_map = finder(convolved_data, threshold) >>> print(segment_map) @@ -292,9 +285,7 @@ measuring source photometry and other source properties, including: Remove labeled segments located within a masked region. Here's a simple example of removing border labels and relabeling the -result: - -.. doctest-requires:: skimage +result:: >>> segment_map3 = segment_map.copy() >>> segment_map3.remove_border_labels(border_width=10, relabel=True) @@ -314,16 +305,12 @@ image. The source mask can be used, for example, to mask sources when estimating the background level. The source mask can optionally be dilated using the ``size`` or ``footprint`` keyword to mask a larger area around each source. Dilating the source mask is useful for -excluding the faint wings of sources when estimating the background: - -.. doctest-requires:: skimage +excluding the faint wings of sources when estimating the background:: >>> mask = segment_map.make_source_mask() >>> dilated_mask = segment_map.make_source_mask(size=11) -A circular footprint can also be used to dilate the source mask: - -.. doctest-requires:: skimage +A circular footprint can also be used to dilate the source mask:: >>> from photutils.utils import circular_footprint >>> footprint = circular_footprint(radius=5) @@ -408,7 +395,7 @@ segmentation image and the science image: kernel = make_2dgaussian_kernel(3.0, size=5) convolved_data = convolve(data, kernel) - finder = SourceFinder(n_pixels=10, progress_bar=False) + finder = SourceFinder(n_pixels=10) segment_map = finder(convolved_data, threshold) fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(10, 12.5)) @@ -491,9 +478,7 @@ from which the source centroids and shape/morphological properties are measured (if not input, the unconvolved image is used instead). Let's continue our example from above and measure the properties of the -detected sources: - -.. doctest-requires:: skimage +detected sources:: >>> from photutils.segmentation import SourceCatalog >>> cat = SourceCatalog(data, segment_map, convolved_data=convolved_data) @@ -516,9 +501,7 @@ generate a `~astropy.table.QTable` of source properties. Each row in the table represents a source. The columns represent the calculated source properties. The ``label`` column corresponds to the label value in the input segmentation image. Note that only a small subset of the source -properties are shown below: - -.. doctest-requires:: skimage +properties are shown below:: >>> tbl = cat.to_table() >>> tbl['x_centroid'].info.format = '.2f' # optional format @@ -593,7 +576,7 @@ of each source) on the data: convolved_data = convolve(data, kernel) n_pixels = 10 - finder = SourceFinder(n_pixels=n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=n_pixels) segment_map = finder(convolved_data, threshold) cat = SourceCatalog(data, segment_map, convolved_data=convolved_data) @@ -610,9 +593,7 @@ of each source) on the data: We can also create a `~photutils.segmentation.SourceCatalog` object containing only a specific subset of sources, defined by their -label numbers in the segmentation image: - -.. doctest-requires:: skimage +label numbers in the segmentation image:: >>> cat = SourceCatalog(data, segment_map, convolved_data=convolved_data) >>> labels = [1, 5, 20, 50, 75, 80] @@ -635,9 +616,7 @@ label numbers in the segmentation image: By default, the :meth:`~photutils.segmentation.SourceCatalog.to_table` includes only a small subset of source properties. The output table properties can be customized in the `~astropy.table.QTable` using the -``columns`` keyword: - -.. doctest-requires:: skimage +``columns`` keyword:: >>> cat = SourceCatalog(data, segment_map, convolved_data=convolved_data) >>> labels = [1, 5, 20, 50, 75, 80] @@ -672,9 +651,7 @@ array that is input to :class:`~photutils.segmentation.SourceCatalog` should be background subtracted. If you input the background image that was subtracted from the data into the ``background`` keyword of :class:`~photutils.segmentation.SourceCatalog`, the background -properties for each source will also be calculated: - -.. doctest-requires:: skimage +properties for each source will also be calculated:: >>> cat = SourceCatalog(data, segment_map, background=bkg.background) >>> labels = [1, 5, 20, 50, 75, 80] @@ -725,9 +702,7 @@ class. When a total ``error`` is input, the `~photutils.segmentation.SourceCatalog.kron_flux_err` properties are calculated. `~photutils.segmentation.SourceCatalog.segment_flux` and `~photutils.segmentation.SourceCatalog.segment_flux_err` are the -instrumental flux and propagated flux error within the source segments: - -.. doctest-requires:: skimage +instrumental flux and propagated flux error within the source segments:: >>> from photutils.utils import calc_total_error >>> effective_gain = 500.0 @@ -766,9 +741,7 @@ properties and their per-axis `~photutils.segmentation.SourceCatalog.y_centroid_win_err`, `~photutils.segmentation.SourceCatalog.x_centroid_quad_err`, and `~photutils.segmentation.SourceCatalog.y_centroid_quad_err` -equivalents: - -.. doctest-requires:: skimage +equivalents:: >>> columns = ['label', 'x_centroid', 'x_centroid_err', 'y_centroid', ... 'y_centroid_err'] @@ -840,9 +813,7 @@ The ``wcs``, ``aperture_mask_method``, and ``kron_params`` keywords are inherited from the ``detection_catalog`` and are therefore ignored when ``detection_catalog`` is input. Note that the segmentation image used to create the detection catalog must be the same one input to the -measurement catalog: - -.. doctest-requires:: skimage +measurement catalog:: >>> det_cat = SourceCatalog(data, segment_map, ... convolved_data=convolved_data) @@ -868,7 +839,7 @@ undefined or degenerate shape properties, windowed or quadratic centroid failures, and Kron-aperture issues (``kron_``-prefixed flags). Accessing ``flags`` computes the flagged quantities (moments, windowed and quadratic centroids, and Kron photometry) if they have -not already been computed; the results are cached and shared with +not already been computed. The results are cached and shared with the corresponding source properties. Flags that describe the same condition as an `~photutils.aperture.ApertureStats` flag, with the source segment as the region, use the same flag name. However, the diff --git a/docs/whats_new/3.1.rst b/docs/whats_new/3.1.rst index 11db702a80..be81992a36 100644 --- a/docs/whats_new/3.1.rst +++ b/docs/whats_new/3.1.rst @@ -578,7 +578,7 @@ quadratic-fit centroid >>> error = np.full(data.shape, 2.0) >>> kernel = make_2dgaussian_kernel(3.0, size=5) >>> convolved_data = convolve(data, kernel) - >>> finder = SourceFinder(n_pixels=10, progress_bar=False) + >>> finder = SourceFinder(n_pixels=10) >>> segment_map = finder(convolved_data, 4.5) >>> cat = SourceCatalog(data, segment_map, ... convolved_data=convolved_data, error=error) @@ -649,6 +649,20 @@ sources in Python any more, no progress bar is displayed and the `). +Source Deblending Performance Improvements +========================================== + +Source deblending with :func:`~photutils.segmentation.deblend_sources` +and `~photutils.segmentation.SourceFinder` is now typically 7--30 times +faster, producing identical results. The whole watershed step now runs +in compiled code. + +The compiled watershed kernel produces results identical to +``skimage.segmentation.watershed``, so scikit-image is no longer +required for source deblending (it remains an optional dependency +for `~photutils.utils.ImageDepth`). + + ``GriddedPSFModel`` Performance Improvements ============================================ @@ -675,6 +689,18 @@ The ``PixelAperture.do_photometry`` method is now deprecated and will be removed in a future version. Use the new `~photutils.aperture.AperturePhotometry` class instead. +Deblending ``progress_bar`` and ``n_processes`` Deprecated +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The ``progress_bar`` and ``n_processes`` keywords of +:func:`~photutils.segmentation.deblend_sources` and +`~photutils.segmentation.SourceFinder` are deprecated and will be +removed in version 4.0, and both keywords no longer have any effect. +Deblending is now dominated by compiled code, so no progress bar is +displayed, and the multiprocessing implementation has been removed +because its process startup and data-pickling overheads made it slower +than the serial implementation. + ``EPSFBuildResult`` Renamed to ``EPSFBuildResults`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/photutils/aperture/flags.py b/photutils/aperture/flags.py index 1a1540e746..97dfba7397 100644 --- a/photutils/aperture/flags.py +++ b/photutils/aperture/flags.py @@ -152,7 +152,7 @@ class _ApertureFlags(FlagRegistry): FlagDefinition( bit_value=4096, name='undefined_shape', - description='non-positive net flux; shape properties undefined', + description='non-positive net flux (shape properties undefined)', detailed_description=('The net flux within the aperture ' '(the zeroth image moment of the ' 'unmasked "center"-method pixels) is ' diff --git a/photutils/segmentation/_batch_catalog.pyx b/photutils/segmentation/_batch_catalog.pyx index f8e0587479..c4c10956d0 100644 --- a/photutils/segmentation/_batch_catalog.pyx +++ b/photutils/segmentation/_batch_catalog.pyx @@ -139,7 +139,7 @@ cdef void _centroid_win_source(const double *data, const double *error, """ Compute the raw windowed-centroid quantities for a single source. - Replicates the previous per-source Python implementation: an + Replicates the previous per-source Python implementation, an iterative Gaussian-weighted centroid within a binary circular window of radius ``4 * sigma``, with masked pixels contributing zero and neighbor-source pixels excluded or mirror-corrected per @@ -236,7 +236,7 @@ cdef void _centroid_win_source(const double *data, const double *error, weighted_flux = sumw # 0/0 yields NaN (cdivision), matching the suppressed NumPy - # RuntimeWarning path; a NaN dcen ends the loop + # RuntimeWarning path. A NaN dcen ends the loop dx_mom = sumwx / sumw dy_mom = sumwy / sumw dcen = sqrt(dx_mom * dx_mom + dy_mom * dy_mom) @@ -265,7 +265,8 @@ cdef void _centroid_win_source(const double *data, const double *error, return # Final pass over the last completed window using the pre-update - # center: windowed central 2nd-order moments and raw error sums. + # center, computing the windowed central 2nd-order moments and + # raw error sums. sxx = 0.0 syy = 0.0 sxy = 0.0 @@ -477,7 +478,7 @@ cdef void _kron_radius_source(const double *data, Accumulate the Kron radius numerator and denominator for a single source. - Replicates the previous per-source Python implementation: the + Replicates the previous per-source Python implementation, the sums of ``data * r`` and ``data`` over the pixels whose centers fall inside the ellipse of elliptical radius ``scale`` (or the circle of radius ``min_circ_radius`` when both axes are zero), @@ -736,7 +737,7 @@ cdef void _flux_radius_cutout(const double *data, Fill the cleaned, background-subtracted cutout of a single source for the flux-radius root-find. - Replicates the previous per-source Python preparation: masked and + Replicates the previous per-source Python preparation. Masked and non-finite pixels are zero, neighbor-source pixels are zeroed or mirror-corrected per ``seg_method`` (an uncorrectable neighbor pixel is zero), and every other pixel is ``data - local_bkg``. @@ -1019,7 +1020,7 @@ cdef void _bin_flux_radius_pixels(_FluxRadiusArgs *p, Py_ssize_t *order, ---------- p : _FluxRadiusArgs * The per-source arguments. On input the grid and data members - are set; on output the ``order``, ``bin_starts``, and + are set. On output the ``order``, ``bin_starts``, and ``bin_cumsum`` members point at the filled arrays. order : Py_ssize_t * @@ -1089,11 +1090,11 @@ cdef double _flux_radius_objective(double r, The flux is ``sum(data * overlap)`` with the same per-pixel overlap arithmetic as ``circular_overlap_grid`` (grid edges relative to the source centroid). The pixels are visited through - the radial bins of ``_bin_flux_radius_pixels``: bins whose pixels + the radial bins of ``_bin_flux_radius_pixels``. Bins whose pixels all lie within ``r - pixel_radius`` of the centroid are fully enclosed (overlap fraction exactly 1, as in the interior fast path - of ``circle_frac_from_d2``) and are added from the prefix sums; - bins beyond ``r + pixel_radius`` have zero overlap; the pixels of + of ``circle_frac_from_d2``) and are added from the prefix sums. + Bins beyond ``r + pixel_radius`` have zero overlap. The pixels of the remaining bins are evaluated individually. """ cdef _FluxRadiusArgs *p = <_FluxRadiusArgs *>args @@ -1154,7 +1155,7 @@ def batch_flux_radius_solve(const double[::1] values, *, evaluates the overlap only for the pixels near the circle boundary (see ``_flux_radius_objective``). Every pixel receives the same overlap fraction as before, but the flux is summed in a different - order, which perturbs the Brent iteration path: the roots agree + order, which perturbs the Brent iteration path. The roots agree with the previous implementation to within the root-finder's absolute tolerance (``xtol`` = 2e-12 pixels), not to rounding. @@ -1313,7 +1314,7 @@ def batch_flux_radius_solve(const double[::1] values, *, rmax = max_radius[i] # Radial bins of a quarter pixel diagonal out to the - # farthest cutout corner; the boundary band of a circle + # farthest cutout corner. The boundary band of a circle # then spans the pixel diagonal plus two bins p.bin_width = 0.5 * p.pixel_radius max_extent = sqrt(fmax(p.xmin_e * p.xmin_e, xmax_e * xmax_e) @@ -1407,7 +1408,7 @@ def batch_perimeter(const unsigned char[:, ::1] mask, *, cross footprint, and the border image is convolved with the ``[[10, 2, 10], [2, 1, 2], [10, 2, 10]]`` kernel at every bounding-box pixel. The histogram of the convolved values below 34 - is returned; the caller applies the perimeter weights of the + is returned. The caller applies the perimeter weights of the estimator of Benkrid et al. (2000) to it. This replicates the previous per-source implementation exactly. @@ -1553,7 +1554,7 @@ def batch_quad_boxes(const double[:, ::1] data, *, peak : 2D ndarray of intp The cutout-frame ``(x, y)`` index of the peak pixel, with shape - ``(n_sources, 2)``. Valid for status 0 and 3; ``(-1, -1)`` + ``(n_sources, 2)``. Valid for status 0 and 3 and ``(-1, -1)`` otherwise. boxes : 2D ndarray of float64 @@ -1712,7 +1713,7 @@ def batch_segment_gather(const double[:, ::1] values, *, offsets : 1D ndarray of intp The start offset of each source in ``packed``, with shape - ``(n_sources + 1,)``; the values of source ``i`` are + ``(n_sources + 1,)``. The values of source ``i`` are ``packed[offsets[i]:offsets[i + 1]]``. counts : 1D ndarray of intp @@ -1924,8 +1925,8 @@ def batch_raw_moments(const double[:, ::1] convdata, *, mask : 2D ndarray of uint8 (C-contiguous) A mask array where bit 1 (value 1) marks input-masked pixels and bit 2 (value 2) marks non-finite data pixels folded into - the mask by the caller. Only bit 1 excludes a pixel here; - non-finite convolved values are excluded by their own test. + the mask by the caller. Only bit 1 excludes a pixel here. + Non-finite convolved values are excluded by their own test. Must have the same shape as ``convdata``. segm : 2D ndarray of intp (C-contiguous) @@ -2039,8 +2040,8 @@ def batch_central_moments(const double[:, ::1] convdata, *, mask : 2D ndarray of uint8 (C-contiguous) A mask array where bit 1 (value 1) marks input-masked pixels and bit 2 (value 2) marks non-finite data pixels folded into - the mask by the caller. Only bit 1 excludes a pixel here; - non-finite convolved values are excluded by their own test. + the mask by the caller. Only bit 1 excludes a pixel here. + Non-finite convolved values are excluded by their own test. Must have the same shape as ``convdata``. segm : 2D ndarray of intp (C-contiguous) diff --git a/photutils/segmentation/_deblend_markers.pyx b/photutils/segmentation/_deblend_markers.pyx new file mode 100644 index 0000000000..be7bfdc0cf --- /dev/null +++ b/photutils/segmentation/_deblend_markers.pyx @@ -0,0 +1,806 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True +# cython: freethreading_compatible=True +""" +Cython kernel that builds the deblending watershed markers from a +level-quantized source cutout in a single component-tree pass. + +The multithreshold marker construction is defined level by level. +At each threshold level, the connected components (with fewer than +``n_pixels`` pixels removed) of the pixels above the threshold are the +candidate sources, and a marker is replaced by its components at a +higher level whenever it contains at least two of them. Computing this +directly requires one full labeling pass per level. + +This kernel instead builds the quantized component tree once, by +adding pixels in decreasing level order to a union-find structure and +snapshotting the components at each populated level (levels between +populated values have identical components and are provably no-ops in +the per-level construction, so they need no snapshots). The marker set +is then derived by descending the tree. Starting from the components +of the lowest level that has at least two components with ``n_pixels`` +or more pixels, each marker is replaced by its sufficiently large +components at the first higher level that contains at least two of them. +The resulting markers are identical to the per-level construction, +including the raster-scan ordering of the marker labels. + +The multithreshold levels themselves are computed by the caller +(vectorized in NumPy over all the sources of a chunk, see +``photutils.segmentation.deblend``) and passed in, so that they are +bitwise identical to the pure-Python reference implementation on every +platform. The kernels here compute the per-source data extrema that the +levels depend on, quantize each cutout against its levels, and build the +markers. + +The kernels run without the GIL and use no global mutable state, so +this module is safe to use from multiple threads, including on +free-threaded Python builds. +""" + +import numpy as np + +from libc.math cimport NAN, isnan +from libc.stdlib cimport free, malloc, realloc + +__all__ = ['deblend_markers_chunk', 'deblend_source_extrema', + 'make_deblend_markers'] + +ctypedef fused data_t: + float + double + +ctypedef fused segm_t: + int + long long + + +cdef struct _Nodes: + # Growable per-node storage for the component tree + int* level + int* repr_pix + int* child + int* sib + unsigned char* kept + Py_ssize_t n + Py_ssize_t cap + + +cdef inline bint _nodes_grow(_Nodes* nodes) noexcept nogil: + """ + Double the node storage capacity and return False on failure. + """ + cdef Py_ssize_t cap = nodes.cap * 2 + cdef int* level = realloc(nodes.level, cap * sizeof(int)) + if level == NULL: + return False + nodes.level = level + cdef int* repr_pix = realloc(nodes.repr_pix, + cap * sizeof(int)) + if repr_pix == NULL: + return False + nodes.repr_pix = repr_pix + cdef int* child = realloc(nodes.child, cap * sizeof(int)) + if child == NULL: + return False + nodes.child = child + cdef int* sib = realloc(nodes.sib, cap * sizeof(int)) + if sib == NULL: + return False + nodes.sib = sib + cdef unsigned char* kept = realloc( + nodes.kept, cap * sizeof(unsigned char)) + if kept == NULL: + return False + nodes.kept = kept + nodes.cap = cap + return True + + +cdef inline Py_ssize_t _find(int* parent, Py_ssize_t x) noexcept nogil: + """ + Return the union-find root of ``x``, compressing the path. + """ + cdef Py_ssize_t root = x + cdef Py_ssize_t nxt + while parent[root] != root: + root = parent[root] + while parent[x] != root: + nxt = parent[x] + parent[x] = root + x = nxt + return root + + +cdef inline void _union(int* parent, int* size, Py_ssize_t a, + Py_ssize_t b) noexcept nogil: + """ + Union the components containing ``a`` and ``b`` by size. + """ + cdef Py_ssize_t ra = _find(parent, a) + cdef Py_ssize_t rb = _find(parent, b) + cdef Py_ssize_t tmp + if ra == rb: + return + if size[ra] < size[rb]: + tmp = ra + ra = rb + rb = tmp + parent[rb] = ra + size[ra] += size[rb] + + +cdef Py_ssize_t _markers_core(const int* qflat, Py_ssize_t ny, + Py_ssize_t nx, int n_pixels, bint conn8, + int* parent, int* size, + unsigned char* added, int* stamp, + int* node_of_root, int* order, + int* flood, int* markers) noexcept nogil: + """ + Build the deblending markers for one level-quantized cutout. + + The caller must provide the pixel-sized workspace arrays with + ``added`` all zero and ``stamp`` all -1. Both invariants are + restored before returning. The marker labels are written into + ``markers`` at the pixels of the cutout that have a nonzero + quantized level (the caller is responsible for treating the + other ``markers`` entries as zero). + + Returns the number of markers (0 if no threshold level has at + least two sufficiently large components), or -1 if a memory + allocation failed. + """ + cdef Py_ssize_t n_tot = ny * nx + cdef Py_ssize_t result = 0 + cdef Py_ssize_t p, i, j, s, m, batch_start, nid, pid, root, nb + cdef Py_ssize_t prev_lo, prev_hi, n_snaps, kept_cnt, base + cdef Py_ssize_t node, chain, child, first_kept, n_kept_children + cdef Py_ssize_t n_stack, n_final, fs, vthr, mn, rank + cdef Py_ssize_t py, px, dy, dx + cdef int v, lbl + + # Count the active (nonzero) pixels and the maximum level + cdef int qmax = 0 + cdef Py_ssize_t n_active = 0 + for p in range(n_tot): + if qflat[p] > 0: + n_active += 1 + if qflat[p] > qmax: + qmax = qflat[p] + if n_active == 0: + return 0 + + cdef _Nodes nodes + nodes.n = 0 + nodes.cap = 256 + nodes.level = malloc(nodes.cap * sizeof(int)) + nodes.repr_pix = malloc(nodes.cap * sizeof(int)) + nodes.child = malloc(nodes.cap * sizeof(int)) + nodes.sib = malloc(nodes.cap * sizeof(int)) + nodes.kept = malloc(nodes.cap + * sizeof(unsigned char)) + + # Counting-sort bookkeeping and per-snapshot bookkeeping + cdef Py_ssize_t* counts = malloc( + (qmax + 1) * sizeof(Py_ssize_t)) + cdef Py_ssize_t* fill = malloc( + (qmax + 1) * sizeof(Py_ssize_t)) + cdef Py_ssize_t* snap_lo = malloc( + qmax * sizeof(Py_ssize_t)) + cdef Py_ssize_t* snap_hi = malloc( + qmax * sizeof(Py_ssize_t)) + cdef Py_ssize_t* snap_kept = malloc( + qmax * sizeof(Py_ssize_t)) + cdef Py_ssize_t* stack = NULL + cdef Py_ssize_t* final = NULL + cdef Py_ssize_t* min_pix = NULL + cdef Py_ssize_t* order_rank = NULL + cdef int* lut = NULL + + if (nodes.level == NULL or nodes.repr_pix == NULL + or nodes.child == NULL or nodes.sib == NULL + or nodes.kept == NULL or counts == NULL or fill == NULL + or snap_lo == NULL or snap_hi == NULL + or snap_kept == NULL): + result = -1 + + if result == 0: + # Sort the active pixels by decreasing level with a counting + # sort. Pixels within a level stay in raster order + for v in range(qmax + 1): + counts[v] = 0 + for p in range(n_tot): + if qflat[p] > 0: + counts[qflat[p]] += 1 + i = 0 + for v in range(qmax, 0, -1): + fill[v] = i + i += counts[v] + for p in range(n_tot): + v = qflat[p] + if v > 0: + order[fill[v]] = p + fill[v] += 1 + + # The caller only guarantees zeros at inactive pixels + for i in range(n_active): + markers[order[i]] = 0 + + # Build the component tree by adding pixels in decreasing + # level order, snapshotting components at populated levels + i = 0 + prev_lo = 0 + prev_hi = 0 + n_snaps = 0 + while i < n_active: + v = qflat[order[i]] + batch_start = i + + # Add pixels at this level and union with added neighbors + while i < n_active and qflat[order[i]] == v: + p = order[i] + parent[p] = p + size[p] = 1 + added[p] = 1 + py = p // nx + px = p % nx + for dy in range(-1, 2): + if py + dy < 0 or py + dy >= ny: + continue + for dx in range(-1, 2): + if px + dx < 0 or px + dx >= nx: + continue + if dy == 0 and dx == 0: + continue + if not conn8 and dy != 0 and dx != 0: + continue + nb = p + dy * nx + dx + if added[nb]: + _union(parent, size, p, nb) + i += 1 + + # Snapshot the components of the level set with threshold + # index v - 1 and attach the previous (higher) level + # components as children + kept_cnt = 0 + for nid in range(prev_lo, prev_hi): + root = _find(parent, nodes.repr_pix[nid]) + if stamp[root] != n_snaps: + stamp[root] = n_snaps + if nodes.n == nodes.cap and not _nodes_grow(&nodes): + result = -1 + break + nodes.level[nodes.n] = v - 1 + nodes.repr_pix[nodes.n] = root + nodes.child[nodes.n] = -1 + nodes.kept[nodes.n] = size[root] >= n_pixels + kept_cnt += nodes.kept[nodes.n] + node_of_root[root] = nodes.n + nodes.n += 1 + pid = node_of_root[root] + nodes.sib[nid] = nodes.child[pid] + nodes.child[pid] = nid + if result != 0: + break + + # Create nodes for components new at this level + for j in range(batch_start, i): + root = _find(parent, order[j]) + if stamp[root] != n_snaps: + stamp[root] = n_snaps + if nodes.n == nodes.cap and not _nodes_grow(&nodes): + result = -1 + break + nodes.level[nodes.n] = v - 1 + nodes.repr_pix[nodes.n] = root + nodes.child[nodes.n] = -1 + nodes.kept[nodes.n] = size[root] >= n_pixels + kept_cnt += nodes.kept[nodes.n] + node_of_root[root] = nodes.n + nodes.n += 1 + if result != 0: + break + + snap_lo[n_snaps] = prev_hi + snap_hi[n_snaps] = nodes.n + snap_kept[n_snaps] = kept_cnt + prev_lo = prev_hi + prev_hi = nodes.n + n_snaps += 1 + + if result == 0: + # Find the base snapshot: the lowest level with at least two + # sufficiently large components + base = -1 + for s in range(n_snaps - 1, -1, -1): + if snap_kept[s] >= 2: + base = s + break + + if base >= 0: + # Descend the tree: replace each marker by its + # sufficiently large components at the first higher + # level containing at least two of them + stack = malloc(nodes.n * sizeof(Py_ssize_t)) + final = malloc(nodes.n * sizeof(Py_ssize_t)) + if stack == NULL or final == NULL: + result = -1 + else: + n_stack = 0 + n_final = 0 + for nid in range(snap_lo[base], snap_hi[base]): + if nodes.kept[nid]: + stack[n_stack] = nid + n_stack += 1 + while n_stack > 0: + n_stack -= 1 + node = stack[n_stack] + chain = node + while True: + n_kept_children = 0 + first_kept = -1 + child = nodes.child[chain] + while child != -1: + if nodes.kept[child]: + n_kept_children += 1 + if n_kept_children == 1: + first_kept = child + elif n_kept_children == 2: + stack[n_stack] = first_kept + n_stack += 1 + stack[n_stack] = child + n_stack += 1 + else: + stack[n_stack] = child + n_stack += 1 + child = nodes.sib[child] + if n_kept_children >= 2: + break + if n_kept_children == 1: + chain = first_kept + continue + final[n_final] = node + n_final += 1 + break + + # Paint each final marker by flood filling its + # component from the recorded representative pixel; + # the marker regions are disjoint, so each pixel is + # visited at most once + min_pix = malloc(n_final + * sizeof(Py_ssize_t)) + order_rank = malloc(n_final + * sizeof(Py_ssize_t)) + lut = malloc((n_final + 1) * sizeof(int)) + if min_pix == NULL or order_rank == NULL or lut == NULL: + result = -1 + + if result == 0: + for m in range(n_final): + nid = final[m] + vthr = nodes.level[nid] + 1 + lbl = (m + 1) + p = nodes.repr_pix[nid] + markers[p] = lbl + flood[0] = p + fs = 1 + mn = p + while fs > 0: + fs -= 1 + p = flood[fs] + if p < mn: + mn = p + py = p // nx + px = p % nx + for dy in range(-1, 2): + if py + dy < 0 or py + dy >= ny: + continue + for dx in range(-1, 2): + if px + dx < 0 or px + dx >= nx: + continue + if dy == 0 and dx == 0: + continue + if not conn8 and dy != 0 and dx != 0: + continue + nb = p + dy * nx + dx + if markers[nb] == 0 and qflat[nb] >= vthr: + markers[nb] = lbl + flood[fs] = nb + fs += 1 + min_pix[m] = mn + + # Relabel the markers in raster-scan order of their + # first pixels (an insertion sort, as the first pixels + # are distinct), matching the ordering that per-level + # labeling would produce + for m in range(n_final): + order_rank[m] = m + for m in range(1, n_final): + j = m + while (j > 0 and min_pix[order_rank[j - 1]] + > min_pix[order_rank[j]]): + rank = order_rank[j - 1] + order_rank[j - 1] = order_rank[j] + order_rank[j] = rank + j -= 1 + lut[0] = 0 + for m in range(n_final): + lut[order_rank[m] + 1] = (m + 1) + for i in range(n_active): + p = order[i] + if markers[p] != 0: + markers[p] = lut[markers[p]] + + result = n_final + + # Restore the workspace invariants for the added and stamp + # arrays. The active pixels are found from qflat rather than from + # order, which is not filled if an early allocation failed + for p in range(n_tot): + if qflat[p] > 0: + added[p] = 0 + stamp[p] = -1 + + free(nodes.level) + free(nodes.repr_pix) + free(nodes.child) + free(nodes.sib) + free(nodes.kept) + free(counts) + free(fill) + free(snap_lo) + free(snap_hi) + free(snap_kept) + free(stack) + free(final) + free(min_pix) + free(order_rank) + free(lut) + + return result + + +def make_deblend_markers(const int[:, ::1] quantized, int n_pixels, + int connectivity): + """ + Build the deblending watershed markers for a single source. + + This single-source entry point is exported only for the pure-Python + reference implementation in ``_deblend_reference`` and the + cross-implementation tests. Production deblending goes through + ``deblend_markers_chunk``. + + Parameters + ---------- + quantized : 2D int `~numpy.ndarray` + The level-quantized source cutout. Each pixel value is the + number of multithreshold levels below the pixel data value + (i.e., the pixel is above threshold level ``i`` if ``i < + quantized``). Pixels outside the source segment and NaN + pixels must be 0. + + n_pixels : int + The minimum number of connected pixels an above-threshold + component must have to be considered a source. + + connectivity : {8, 4} + The pixel connectivity. + + Returns + ------- + markers : 2D int `~numpy.ndarray` + The marker image, with markers labeled consecutively from 1 + in raster-scan order. All values are 0 if no markers were + found. + + n_markers : int + The number of markers. Zero means no threshold level had at + least two components with ``n_pixels`` or more pixels. + """ + cdef Py_ssize_t ny = quantized.shape[0] + cdef Py_ssize_t nx = quantized.shape[1] + cdef Py_ssize_t n_tot = ny * nx + + markers_arr = np.zeros((ny, nx), dtype=np.int32) + parent_arr = np.empty(n_tot, dtype=np.int32) + size_arr = np.zeros(n_tot, dtype=np.int32) + added_arr = np.zeros(n_tot, dtype=np.uint8) + stamp_arr = np.full(n_tot, -1, dtype=np.int32) + node_of_root_arr = np.zeros(n_tot, dtype=np.int32) + order_arr = np.empty(n_tot, dtype=np.int32) + flood_arr = np.empty(n_tot, dtype=np.int32) + + cdef const int* qflat = &quantized[0, 0] + cdef int[:, ::1] markers_mv = markers_arr + cdef int[::1] parent_mv = parent_arr + cdef int[::1] size_mv = size_arr + cdef unsigned char[::1] added_mv = added_arr + cdef int[::1] stamp_mv = stamp_arr + cdef int[::1] node_of_root_mv = node_of_root_arr + cdef int[::1] order_mv = order_arr + cdef int[::1] flood_mv = flood_arr + + cdef Py_ssize_t n_markers + with nogil: + n_markers = _markers_core(qflat, ny, nx, n_pixels, + connectivity == 8, + &parent_mv[0], &size_mv[0], + &added_mv[0], &stamp_mv[0], + &node_of_root_mv[0], &order_mv[0], + &flood_mv[0], &markers_mv[0, 0]) + if n_markers < 0: + raise MemoryError + + return markers_arr, int(n_markers) + + +cdef inline int _count_below(const double* thresholds, int n_levels, + double value) noexcept nogil: + """ + Return the number of thresholds strictly below ``value``. + + This matches ``np.searchsorted(thresholds, value, side='left')`` + for non-NaN values. + """ + cdef int lo = 0 + cdef int hi = n_levels + cdef int mid + while lo < hi: + mid = (lo + hi) // 2 + if thresholds[mid] < value: + lo = mid + 1 + else: + hi = mid + return lo + + +cdef void _source_extrema(const data_t* data, const segm_t* segm, + Py_ssize_t img_nx, long long label, + Py_ssize_t y0, Py_ssize_t y1, Py_ssize_t x0, + Py_ssize_t x1, double* smin, + double* smax) noexcept nogil: + """ + Compute the minimum and maximum data value of one source segment. + + NaN pixels are excluded. Both outputs are NaN if the segment has no + finite pixel. + """ + cdef Py_ssize_t iy, ix, idx + cdef double value + cdef bint has_value = False + + smin[0] = NAN + smax[0] = NAN + for iy in range(y1 - y0): + for ix in range(x1 - x0): + idx = (y0 + iy) * img_nx + x0 + ix + if segm[idx] != label: + continue + value = data[idx] + if isnan(value): + continue + if not has_value: + smin[0] = value + smax[0] = value + has_value = True + elif value < smin[0]: + smin[0] = value + elif value > smax[0]: + smax[0] = value + + +def deblend_source_extrema(const data_t[:, ::1] data, + const segm_t[:, ::1] segm_data, + const long long[::1] labels, + const long long[::1] y0, + const long long[::1] y1, + const long long[::1] x0, + const long long[::1] x1): + """ + Compute the minimum and maximum data value of each source segment. + + NaN pixels are excluded, so the results are identical to the + ``nanmin`` and ``nanmax`` reductions over the segment pixels. + + Parameters + ---------- + data : 2D float `~numpy.ndarray` + The full data array. + + segm_data : 2D int `~numpy.ndarray` + The full segmentation array. + + labels : 1D int64 `~numpy.ndarray` + The label of each source. + + y0, y1, x0, x1 : 1D int64 `~numpy.ndarray` + The bounding-box slice bounds of each source. + + Returns + ------- + source_min, source_max : 1D float64 `~numpy.ndarray` + The minimum and maximum data value of each source segment. Both + are NaN for a segment without any finite pixel. + """ + cdef Py_ssize_t n_src = labels.shape[0] + cdef Py_ssize_t img_nx = data.shape[1] + cdef Py_ssize_t isrc + + smin_arr = np.empty(n_src, dtype=np.float64) + smax_arr = np.empty(n_src, dtype=np.float64) + cdef double[::1] smin_mv = smin_arr + cdef double[::1] smax_mv = smax_arr + + with nogil: + for isrc in range(n_src): + _source_extrema(&data[0, 0], &segm_data[0, 0], img_nx, + labels[isrc], y0[isrc], y1[isrc], x0[isrc], + x1[isrc], &smin_mv[isrc], &smax_mv[isrc]) + + return smin_arr, smax_arr + + +cdef Py_ssize_t _source_markers(const data_t* data, const segm_t* segm, + Py_ssize_t img_nx, long long label, + Py_ssize_t y0, Py_ssize_t y1, + Py_ssize_t x0, Py_ssize_t x1, + int n_pixels, bint conn8, int n_levels, + const double* thresholds, int* q, + int* parent, int* size, + unsigned char* added, int* stamp, + int* node_of_root, int* order, + int* flood, int* markers) noexcept nogil: + """ + Build the deblending markers for one source of a chunk. + + Quantizes the cutout against the multithreshold levels of the + source (the number of levels strictly below each pixel value, with + the pixels outside the segment and the NaN pixels at 0) and runs + the component-tree core. Returns the number of markers (0 means the + source does not split) or -1 if a memory allocation failed. + """ + cdef Py_ssize_t ny_c = y1 - y0 + cdef Py_ssize_t nx_c = x1 - x0 + cdef Py_ssize_t iy, ix, idx + cdef double value + + for iy in range(ny_c): + for ix in range(nx_c): + idx = (y0 + iy) * img_nx + x0 + ix + if segm[idx] != label: + q[iy * nx_c + ix] = 0 + continue + value = data[idx] + if isnan(value): + q[iy * nx_c + ix] = 0 + continue + q[iy * nx_c + ix] = _count_below(thresholds, n_levels, + value) + + return _markers_core(q, ny_c, nx_c, n_pixels, conn8, parent, size, + added, stamp, node_of_root, order, flood, + markers) + + +def deblend_markers_chunk(const data_t[:, ::1] data, + const segm_t[:, ::1] segm_data, + const long long[::1] labels, + const long long[::1] y0, + const long long[::1] y1, + const long long[::1] x0, + const long long[::1] x1, + const double[:, ::1] thresholds, *, + int n_pixels, int connectivity, + int max_markers): + """ + Build the deblending watershed markers for a chunk of sources. + + Each cutout is quantized against its own multithreshold levels and + its markers are built by the component-tree kernel, in compiled + code that releases the GIL and reuses one workspace sized to the + largest cutout in the chunk. + + Parameters + ---------- + data : 2D float `~numpy.ndarray` + The full data array. + + segm_data : 2D int `~numpy.ndarray` + The full segmentation array. + + labels : 1D int64 `~numpy.ndarray` + The label of each source in the chunk. + + y0, y1, x0, x1 : 1D int64 `~numpy.ndarray` + The bounding-box slice bounds of each source. + + thresholds : 2D float64 `~numpy.ndarray` + The multithreshold levels of each source, with shape + ``(n_sources, n_levels)`` and ascending along the second axis. + + n_pixels : int + The minimum number of connected pixels an above-threshold + component must have to be considered a source. + + connectivity : {8, 4} + The pixel connectivity. + + max_markers : int + The number of markers above which the marker image of a source + is not built. Its marker count is still returned, so that the + caller can retry the source with other levels. A negative value + disables the limit. + + Returns + ------- + markers_list : list of 2D int `~numpy.ndarray` or `None` + The cutout marker image of each source that splits into two or + more markers (and no more than ``max_markers``), otherwise + `None`. + + n_markers : 1D intp `~numpy.ndarray` + The number of markers found for each source. + """ + cdef Py_ssize_t n_src = labels.shape[0] + cdef Py_ssize_t img_nx = data.shape[1] + cdef int n_levels = thresholds.shape[1] + cdef Py_ssize_t isrc, n_tot, max_ntot, ny_c, nx_c + cdef Py_ssize_t n_markers + + if thresholds.shape[0] != n_src: + msg = 'thresholds must have one row per source' + raise ValueError(msg) + + max_ntot = 1 + for isrc in range(n_src): + n_tot = (y1[isrc] - y0[isrc]) * (x1[isrc] - x0[isrc]) + if n_tot > max_ntot: + max_ntot = n_tot + + q_arr = np.empty(max_ntot, dtype=np.int32) + parent_arr = np.empty(max_ntot, dtype=np.int32) + size_arr = np.zeros(max_ntot, dtype=np.int32) + added_arr = np.zeros(max_ntot, dtype=np.uint8) + stamp_arr = np.full(max_ntot, -1, dtype=np.int32) + node_of_root_arr = np.zeros(max_ntot, dtype=np.int32) + order_arr = np.empty(max_ntot, dtype=np.int32) + flood_arr = np.empty(max_ntot, dtype=np.int32) + markers_arr = np.zeros(max_ntot, dtype=np.int32) + n_markers_arr = np.zeros(n_src, dtype=np.intp) + + cdef int[::1] q_mv = q_arr + cdef int[::1] parent_mv = parent_arr + cdef int[::1] size_mv = size_arr + cdef unsigned char[::1] added_mv = added_arr + cdef int[::1] stamp_mv = stamp_arr + cdef int[::1] node_of_root_mv = node_of_root_arr + cdef int[::1] order_mv = order_arr + cdef int[::1] flood_mv = flood_arr + cdef int[::1] markers_mv = markers_arr + cdef Py_ssize_t[::1] n_markers_mv = n_markers_arr + + markers_list = [] + for isrc in range(n_src): + ny_c = y1[isrc] - y0[isrc] + nx_c = x1[isrc] - x0[isrc] + with nogil: + n_markers = _source_markers( + &data[0, 0], &segm_data[0, 0], img_nx, labels[isrc], + y0[isrc], y1[isrc], x0[isrc], x1[isrc], n_pixels, + connectivity == 8, n_levels, &thresholds[isrc, 0], + &q_mv[0], &parent_mv[0], &size_mv[0], &added_mv[0], + &stamp_mv[0], &node_of_root_mv[0], &order_mv[0], + &flood_mv[0], &markers_mv[0]) + if n_markers < 0: + raise MemoryError + n_markers_mv[isrc] = n_markers + if n_markers >= 2 and (max_markers < 0 + or n_markers <= max_markers): + n_tot = ny_c * nx_c + quantized = q_arr[:n_tot].reshape(ny_c, nx_c) + markers = markers_arr[:n_tot].reshape(ny_c, nx_c) + markers_list.append(np.where(quantized > 0, markers, + np.int32(0))) + else: + markers_list.append(None) + + return markers_list, n_markers_arr diff --git a/photutils/segmentation/_deblend_reference.py b/photutils/segmentation/_deblend_reference.py new file mode 100644 index 0000000000..6d4d370e67 --- /dev/null +++ b/photutils/segmentation/_deblend_reference.py @@ -0,0 +1,504 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +Pure-Python reference implementation of the deblending pipeline. + +This module is a pure-Python mirror of the compiled deblending +pipeline used by :func:`~photutils.segmentation.deblend_sources`, +kept for verification and debugging. It is not used in +production. It must track the compiled (Cython) semantics +exactly, which is enforced by the cross-implementation tests in +``photutils/segmentation/tests/test_deblend.py``. +""" + +from functools import cached_property + +import numpy as np +from scipy.ndimage import label as ndi_label +from scipy.ndimage import sum_labels + +from photutils.segmentation._deblend_markers import make_deblend_markers +from photutils.segmentation._deblend_watershed import deblend_watershed +from photutils.segmentation.core import _get_labels +from photutils.segmentation.deblend import _MAX_MARKERS, _create_relabel_map +from photutils.utils._stats import nanmax, nanmin, nansum + + +def _detect_sources_deblend(data, threshold, n_pixels, *, footprint, + segment_mask): + """ + Detect sources for a single multithreshold level during deblending. + + This is the deblending analogue of + `photutils.segmentation.detect._detect_sources`. It differs in + that the detected segments keep their (possibly non-consecutive) + label numbers from `~scipy.ndimage.label`, the small segments are + removed with a bincount-based area filter (the per-label cutout + loop used by ``_detect_sources`` has a fixed per-label overhead + that dominates for the small cutouts and the many calls made + during deblending), and `None` is returned when fewer than two + segments are found. + + Parameters + ---------- + data : 2D `~numpy.ndarray` + The cutout data array for a single source. + + threshold : float + The data value to be used for the detection threshold. + + n_pixels : int + The minimum number of connected pixels, each greater than + ``threshold``, that an object must have to be detected. + + footprint : array_like + A footprint that defines feature connections. + + segment_mask : 2D bool `~numpy.ndarray` + A boolean mask of the source segment, with the same shape as + ``data``. Pixels outside the segment will not be included in + any source. + + Returns + ------- + segment_img : 2D int `~numpy.ndarray` or `None` + A 2D segmentation image, with the same shape as ``data``, + where sources are marked by different positive integer + values. A value of zero is reserved for the background. If + fewer than two sources are found then `None` is returned. + """ + # NaN values compare as False, so NaN pixels are never included + # in any source. The comparison is never empty because the + # deblending thresholds are strictly below the source maximum. + segment_img = data > threshold + segment_img &= segment_mask + + segment_img, n_labels = ndi_label(segment_img, structure=footprint) + + # Remove objects with less than n_pixels + areas = np.bincount(segment_img.ravel()) + keep = areas >= n_pixels + keep[0] = False + n_keep = np.count_nonzero(keep) + if n_keep <= 1: + return None + + if n_keep < n_labels: + label_map = np.where( + keep, np.arange(areas.size, dtype=segment_img.dtype), 0) + segment_img = label_map[segment_img] + + return segment_img + + +class _SingleSourceDeblender: + """ + Class to deblend a single labeled source. + + Parameters + ---------- + data : 2D `~numpy.ndarray` + The cutout data array for a single source. ``data`` should + also already be smoothed by the same filter used in + :func:`~photutils.segmentation.detect_sources`, if applicable. + + segment_data : 2D int `~numpy.ndarray` + The cutout segmentation image for a single source. Must have the + same shape as ``data``. + + label : int + The label of the source to deblend. This is needed because there + may be more than one source label within the cutout. + + deblend_params : `~photutils.segmentation.deblend._DeblendParams` + The parameters for deblending the source. + """ + + def __init__(self, data, segment_data, label, deblend_params): + self.data = data + self.segment_data = segment_data + self.label = label + self.n_pixels = deblend_params.n_pixels + self.footprint = deblend_params.footprint + self.n_levels = deblend_params.n_levels + self.contrast = deblend_params.contrast + self.mode = deblend_params.mode + + self.segment_mask = segment_data == label + data_values = data[self.segment_mask] + self.source_min = nanmin(data_values) + self.source_max = nanmax(data_values) + self.source_sum = nansum(data_values) + self.warnings = {} + + @cached_property + def linear_thresholds(self): + """ + Linearly spaced thresholds between the source minimum and + maximum (inclusive). + + The source min/max are excluded later, giving n_levels + thresholds between min and max (noninclusive). + """ + return np.linspace(self.source_min, self.source_max, self.n_levels + 2) + + @cached_property + def normalized_thresholds(self): + """ + Normalized thresholds (from 0 to 1) between the source minimum + and maximum (inclusive). + """ + return ((self.linear_thresholds - self.source_min) + / (self.source_max - self.source_min)) + + def compute_thresholds(self): + """ + Compute the multi-level detection thresholds for the source. + + Note that this method has side effects. When the mode is + "exponential" and the source minimum is non-positive, it + changes ``self.mode`` to "linear" and records the fallback in + ``self.warnings``. Later calls (e.g., from ``make_markers``) + therefore use the fallback mode, mirroring the sticky per-source + mode fallback in the compiled pipeline. + + Returns + ------- + thresholds : 1D `~numpy.ndarray` + The multi-level detection thresholds for the source. + """ + if self.mode == 'exponential' and self.source_min <= 0: + self.warnings['nonposmin'] = 'non-positive minimum' + self.mode = 'linear' + + if self.mode == 'linear': + thresholds = self.linear_thresholds + elif self.mode == 'sinh': + a = 0.25 + minval = self.source_min + maxval = self.source_max + thresholds = self.normalized_thresholds + thresholds = np.sinh(thresholds / a) / np.sinh(1.0 / a) + thresholds *= (maxval - minval) + thresholds += minval + elif self.mode == 'exponential': + minval = self.source_min + maxval = self.source_max + thresholds = self.normalized_thresholds + thresholds = minval * (maxval / minval) ** thresholds + + return thresholds[1:-1] # do not include source min and max + + def multithreshold(self): + """ + Perform multithreshold detection for each source. + + This method is useful for debugging and testing. + + Returns + ------- + segments : list of 2D `~numpy.ndarray` or `None` + A list of segmentation images, one for each threshold. + `None` is returned for thresholds that do not have more than + one label. + """ + thresholds = self.compute_thresholds() + segms = [] + for threshold in thresholds: + segm = _detect_sources_deblend(self.data, threshold, + self.n_pixels, + footprint=self.footprint, + segment_mask=self.segment_mask) + segms.append(segm) + return segms + + def make_markers(self): + """ + Make markers (possible sources) for the watershed algorithm. + + The markers are built from a single component-tree pass over + the level-quantized cutout (see + `~photutils.segmentation._deblend_markers.make_deblend_markers`), + which produces markers identical to the per-level + multithreshold construction of ``make_markers_per_level``. + + Returns + ------- + markers : 2D `~numpy.ndarray` or `None` + A segmentation image that contains markers for possible + sources. `None` is returned if there is only one source + at every threshold. + """ + thresholds = self.compute_thresholds() + + # A pixel is above threshold level i if i < quantized. NaN + # pixels compare as False against every threshold. + quantized = np.searchsorted(thresholds, self.data.ravel(), + side='left') + quantized = quantized.reshape(self.data.shape).astype(np.int32) + quantized[~self.segment_mask | np.isnan(self.data)] = 0 + connectivity = 8 if self.footprint[0, 0] else 4 + markers, n_markers = make_deblend_markers(quantized, + self.n_pixels, + connectivity) + if n_markers == 0: + return None + return markers + + def make_markers_per_level(self): + """ + Make markers with the per-level multithreshold construction. + + The markers are refined level by level, replacing every + marker that splits at the next-higher threshold by its + children (see ``make_marker_segment``). The last list + element contains the final markers, identical to the + ``make_markers`` result. This method is useful for debugging + and testing. + + Returns + ------- + markers : list of (2D `~numpy.ndarray` or `None`) + The segmentation marker image after each threshold level + that had two or more segments. + """ + thresholds = self.compute_thresholds() + segm_lower = _detect_sources_deblend( + self.data, thresholds[0], self.n_pixels, + footprint=self.footprint, segment_mask=self.segment_mask) + all_segms = [segm_lower] + for threshold in thresholds[1:]: + segm_upper = _detect_sources_deblend( + self.data, threshold, self.n_pixels, + footprint=self.footprint, + segment_mask=self.segment_mask) + if segm_upper is None: # 0 or 1 labels + continue + segm_lower = self.make_marker_segment(segm_lower, + segm_upper) + all_segms.append(segm_lower) + return all_segms + + def make_marker_segment(self, segment_lower, segment_upper): + """ + Make markers (possible sources) for the watershed algorithm. + + Parameters + ---------- + segment_lower : 2D `~numpy.ndarray` + The "lower" threshold level segmentation image. + + segment_upper : 2D `~numpy.ndarray` + The next-highest threshold level segmentation image. + + Returns + ------- + markers : 2D `~numpy.ndarray` + A segmentation image that contain markers for possible + sources. + + Notes + ----- + For a given label in the lower level, find the labels in the + upper level (higher threshold value) that are its children + (i.e., the labels within the same mask as the lower level). If + there are multiple children, then the lower-level parent label + is replaced by its children. Parent labels that do not have + multiple children in the upper level are kept as is (maximizing + the marker size). + """ + if segment_lower is None: + return segment_upper + + # Count the upper-level children of each lower-level label from + # the unique (lower, upper) label pairs, encoded as a combined + # integer key. + both = (segment_lower > 0) & (segment_upper > 0) + stride = np.int64(np.max(segment_upper)) + 1 + keys = segment_lower[both].astype(np.int64) * stride + keys += segment_upper[both] + parents, n_children = np.unique(np.unique(keys) // stride, + return_counts=True) + multi_parents = parents[n_children >= 2] + if multi_parents.size == 0: + return segment_lower + + # Replace each multi-child parent by its children. Pixels of + # the parent mask that are below the upper threshold are unset. + # Single-child parents are kept as is (maximizing the marker + # size). + replace_lut = np.zeros(np.max(segment_lower) + 1, dtype=bool) + replace_lut[multi_parents] = True + replace = replace_lut[segment_lower] + markers = segment_lower > 0 + markers[replace] = segment_upper[replace] > 0 + + # Convert bool markers to integer labels + return ndi_label(markers, structure=self.footprint)[0] + + def apply_watershed(self, markers): + """ + Apply the watershed algorithm to the source markers. + + Parameters + ---------- + markers : 2D int `~numpy.ndarray` + The marker image for the source. + + Returns + ------- + segment_data : 2D int `~numpy.ndarray` + A 2D int array containing the deblended source labels. Note + that the source labels may not be consecutive if a label was + removed. + """ + # Deblend using watershed. If any source does not meet the + # contrast criterion, then remove the faintest such source(s) + # and repeat until all sources meet the contrast criterion. + data_neg = np.ascontiguousarray(-self.data, dtype=np.float64) + # NaN pixels are flooded after all finite pixels, as in the + # compiled contrast loop + data_neg[np.isnan(data_neg)] = np.inf + connectivity = 8 if self.footprint[0, 0] else 4 + remove_marker = True + while remove_marker: + markers = deblend_watershed(data_neg, markers, + self.segment_mask, connectivity) + + labels = _get_labels(markers) + if labels.size == 1: # only 1 source left + remove_marker = False + else: + flux_frac = (sum_labels(self.data, markers, index=labels) + / self.source_sum) + remove_marker = any(flux_frac < self.contrast) + + if remove_marker: + self._remove_faint_markers(markers, labels, flux_frac) + + return markers + + def _remove_faint_markers(self, markers, labels, flux_frac): + """ + Remove the faintest below-contrast marker(s) in place. + + The faintest marker is always removed. When the source data + values are all nonnegative, the largest batch of the faintest + markers is removed whose total flux fraction is below both + the contrast and the next-faintest marker flux fraction. + Removing such a batch in one step is equivalent to removing + its markers one at a time. Every batch member stays below the + contrast no matter how the flux of the other removed members is + redistributed (their total is below the contrast), and markers + outside the batch can only become brighter, so the faintest + below-contrast marker always lies inside the batch until the + batch is exhausted. + + Faint markers that do not fit in such a batch are removed one at + a time because several faint sources could combine to meet the + contrast criterion. + + Parameters + ---------- + markers : 2D int `~numpy.ndarray` + The watershed-labeled marker image, modified in place. + + labels : 1D `~numpy.ndarray` + The sorted marker labels. + + flux_frac : 1D `~numpy.ndarray` + The flux fraction in each marker basin, in the same + order as ``labels``. + """ + if self.source_min >= 0 and labels.size > 2: + order = np.argsort(flux_frac) + sorted_frac = flux_frac[order] + csum = np.cumsum(sorted_frac) + # A batch of the n faintest markers (2 <= n < N) is valid + # if its total flux fraction is below both the contrast + # and the next-faintest marker flux fraction. + batch_ok = ((csum[1:-1] < self.contrast) + & (csum[1:-1] < sorted_frac[2:])) + valid = np.nonzero(batch_ok)[0] + if valid.size > 0: + n_remove = int(valid[-1]) + 2 + remove_lut = np.zeros(int(labels[-1]) + 1, dtype=bool) + remove_lut[labels[order[:n_remove]]] = True + markers[remove_lut[markers]] = 0 + return + + markers[markers == labels[np.argmin(flux_frac)]] = 0 + + def deblend_from_markers(self, markers): + """ + Deblend the source given its precomputed watershed markers. + + Parameters + ---------- + markers : 2D int `~numpy.ndarray` + The marker image for the source. + + Returns + ------- + segment_data : 2D int `~numpy.ndarray` or `None` + A 2D int array containing the deblended source labels. + The source labels are consecutive starting at 1. `None` + is returned if only one source remains after applying + the contrast criterion. + """ + # Deblend using the watershed algorithm using the markers as seeds + markers = self.apply_watershed(markers) + + if not np.array_equal(self.segment_mask, markers.astype(bool)): + msg = (f'Deblending failed for source {self.label!r}. ' + 'Please ensure you used the same pixel connectivity ' + 'in detect_sources and deblend_sources.') + raise ValueError(msg) + + if len(_get_labels(markers)) == 1: # no deblending + return None + + # Markers may not be consecutive if a label was removed due to + # the contrast criterion + relabel_map = _create_relabel_map(markers, start_label=1) + if relabel_map is not None: + markers = relabel_map[markers] + return markers + + def deblend_source(self): + """ + Deblend a single labeled source. + + This method computes the markers and the watershed steps + entirely in Python, mirroring what ``deblend_sources`` + computes through the compiled chunk driver. It is useful for + debugging and testing. + + Returns + ------- + segment_data : 2D int `~numpy.ndarray` or `None` + A 2D int array containing the deblended source labels. The + source labels are consecutive starting at 1. + """ + if self.source_min == self.source_max: # no deblending + return None + + # Define the markers (possible sources) for the watershed algorithm + markers = self.make_markers() + if markers is None: + return None + + # If there are too many markers (e.g., due to low threshold + # and/or small n_pixels), the watershed step can be very slow. + # This mostly affects the "exponential" mode, where there are + # many levels at low thresholds, so here we try again with + # "linear" mode. + n_labels = len(_get_labels(markers)) + if self.mode != 'linear' and n_labels > _MAX_MARKERS: + del markers # free memory + self.warnings['n_markers'] = 'too many markers' + self.mode = 'linear' + markers = self.make_markers() + if markers is None: + return None + + return self.deblend_from_markers(markers) diff --git a/photutils/segmentation/_deblend_watershed.pyx b/photutils/segmentation/_deblend_watershed.pyx new file mode 100644 index 0000000000..49ea95ea9e --- /dev/null +++ b/photutils/segmentation/_deblend_watershed.pyx @@ -0,0 +1,570 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True +# cython: freethreading_compatible=True +""" +Cython marker-based watershed kernel for source deblending. + +This implements the classic priority-flood watershed (Soille 1990) +used for deblending. Pixels are flooded from the markers in order of +increasing image value, with the queue-entry age breaking ties so +that plateaus are split between the markers that reach them first. +The algorithm, the neighbor ordering (orthogonal neighbors before +diagonal ones, each group in raster order), and the tie-breaking +match ``skimage.segmentation.watershed`` (with ``compactness=0`` +and ``watershed_line=False``), so the results are identical, but +without the per-call validation, padding, and cropping overhead of the +general-purpose function, which dominates for small cutouts. + +The flood order is only defined for ordered image values, so the +deblending entry point maps NaN data pixels to a +inf flooding cost. +Such pixels are assigned to the basins that reach them after all the +finite pixels have been assigned. + +The flood core runs without the GIL and uses no global mutable state, +so this module is safe to use from multiple threads, including on +free-threaded Python builds. +""" + +import numpy as np + +from libc.math cimport INFINITY, isnan +from libc.stdlib cimport free, malloc + +__all__ = ['deblend_source_contrast', 'deblend_watershed'] + +ctypedef fused data_t: + float + double + +ctypedef fused segm_t: + int + long long + + +cdef struct _Heap: + # Binary min-heap on (value, age) + double* value + long long* age + int* index + Py_ssize_t n + + +cdef inline bint _heap_less(_Heap* heap, Py_ssize_t a, + Py_ssize_t b) noexcept nogil: + """ + Compare two heap slots by (value, age). + """ + if heap.value[a] != heap.value[b]: + return heap.value[a] < heap.value[b] + return heap.age[a] < heap.age[b] + + +cdef inline void _heap_swap(_Heap* heap, Py_ssize_t a, + Py_ssize_t b) noexcept nogil: + """ + Swap two heap slots. + """ + cdef double value = heap.value[a] + cdef long long age = heap.age[a] + cdef int index = heap.index[a] + heap.value[a] = heap.value[b] + heap.age[a] = heap.age[b] + heap.index[a] = heap.index[b] + heap.value[b] = value + heap.age[b] = age + heap.index[b] = index + + +cdef inline void _heap_push(_Heap* heap, double value, long long age, + Py_ssize_t index) noexcept nogil: + """ + Push an item onto the heap. + """ + cdef Py_ssize_t pos = heap.n + cdef Py_ssize_t parent + heap.value[pos] = value + heap.age[pos] = age + heap.index[pos] = index + heap.n += 1 + while pos > 0: + parent = (pos - 1) // 2 + if _heap_less(heap, pos, parent): + _heap_swap(heap, pos, parent) + pos = parent + else: + break + + +cdef inline Py_ssize_t _heap_pop(_Heap* heap, + double* value) noexcept nogil: + """ + Pop the smallest item, returning its pixel index and value. + """ + cdef Py_ssize_t result = heap.index[0] + cdef Py_ssize_t pos = 0 + cdef Py_ssize_t child + value[0] = heap.value[0] + heap.n -= 1 + if heap.n > 0: + _heap_swap(heap, 0, heap.n) + while True: + child = 2 * pos + 1 + if child >= heap.n: + break + if (child + 1 < heap.n + and _heap_less(heap, child + 1, child)): + child += 1 + if _heap_less(heap, child, pos): + _heap_swap(heap, pos, child) + pos = child + else: + break + return result + + +cdef int _watershed_core(const double* image, unsigned char* mask, + int* output, Py_ssize_t ny, Py_ssize_t nx, + bint conn8) noexcept nogil: + """ + Flood the masked pixels of ``output`` from its nonzero markers. + + Returns 0 on success or -1 if a memory allocation failed. + """ + # Neighbor offsets: orthogonal neighbors before diagonal ones, + # each group in raster order (the stable distance ordering) + cdef Py_ssize_t[8] off_y + cdef Py_ssize_t[8] off_x + cdef Py_ssize_t n_off + off_y[0] = -1 + off_x[0] = 0 + off_y[1] = 0 + off_x[1] = -1 + off_y[2] = 0 + off_x[2] = 1 + off_y[3] = 1 + off_x[3] = 0 + if conn8: + off_y[4] = -1 + off_x[4] = -1 + off_y[5] = -1 + off_x[5] = 1 + off_y[6] = 1 + off_x[6] = -1 + off_y[7] = 1 + off_x[7] = 1 + n_off = 8 + else: + n_off = 4 + + # Each masked pixel enters the queue at most once + cdef Py_ssize_t n_tot = ny * nx + cdef Py_ssize_t cap = 0 + cdef Py_ssize_t p + for p in range(n_tot): + if mask[p]: + cap += 1 + if cap == 0: + return 0 + + cdef _Heap heap + heap.n = 0 + heap.value = malloc(cap * sizeof(double)) + heap.age = malloc(cap * sizeof(long long)) + heap.index = malloc(cap * sizeof(int)) + if heap.value == NULL or heap.age == NULL or heap.index == NULL: + free(heap.value) + free(heap.age) + free(heap.index) + return -1 + + # Push the marker pixels in raster order with age 0 + cdef long long age = 0 + for p in range(n_tot): + if output[p] != 0 and mask[p]: + _heap_push(&heap, image[p], 0, p) + + cdef Py_ssize_t py, px, ny_i, nx_i, nb, i + cdef double pop_value, push_value + while heap.n > 0: + p = _heap_pop(&heap, &pop_value) + py = p // nx + px = p % nx + for i in range(n_off): + ny_i = py + off_y[i] + nx_i = px + off_x[i] + if ny_i < 0 or ny_i >= ny or nx_i < 0 or nx_i >= nx: + continue + nb = ny_i * nx + nx_i + if not mask[nb]: + continue + if output[nb] != 0: + continue + age += 1 + output[nb] = output[p] + # The flooding cost of a pixel is at least the cost of + # the pixel it was reached from, so that plateaus and + # basins below the current flood level are distributed + # by queue-entry age between contesting markers + push_value = image[nb] + if push_value < pop_value: + push_value = pop_value + _heap_push(&heap, push_value, age, nb) + + free(heap.value) + free(heap.age) + free(heap.index) + return 0 + + +def deblend_watershed(image, markers, mask, connectivity): + """ + Compute the marker-based watershed of an image. + + This is equivalent to ``skimage.segmentation.watershed(image, + markers, mask=mask, connectivity=footprint)`` for the deblending + use case (markers inside the mask, ``compactness=0``, and + ``watershed_line=False``), but avoids the per-call validation, + padding, and cropping overhead. + + This entry point is exported only for the pure-Python + reference implementation in ``_deblend_reference`` and the + cross-implementation tests. The production contrast loop calls the + watershed core directly through ``deblend_source_contrast``. + + Parameters + ---------- + image : 2D `~numpy.ndarray` + The image to flood (the lowest values are flooded first). It + must not contain NaN values. The deblending callers map NaN + data pixels to +inf so that they are flooded last. + + markers : 2D int `~numpy.ndarray` + The marker image. Zero means not a marker. All markers must + lie inside the mask. + + mask : 2D bool `~numpy.ndarray` + Only pixels where the mask is `True` are labeled. + + connectivity : {8, 4} + The pixel connectivity. + + Returns + ------- + output : 2D int `~numpy.ndarray` + The labeled basins, with the same shape as ``image``. + """ + image_arr = np.ascontiguousarray(image, dtype=np.float64) + output_arr = np.array(markers, dtype=np.int32, copy=True, order='C') + mask_arr = np.ascontiguousarray(mask, dtype=np.uint8) + + cdef const double[:, ::1] image_mv = image_arr + cdef int[:, ::1] output_mv = output_arr + cdef unsigned char[:, ::1] mask_mv = mask_arr + cdef Py_ssize_t ny = image_mv.shape[0] + cdef Py_ssize_t nx = image_mv.shape[1] + + cdef bint conn8 = connectivity == 8 + cdef int status + with nogil: + status = _watershed_core(&image_mv[0, 0], &mask_mv[0, 0], + &output_mv[0, 0], ny, nx, conn8) + if status < 0: + raise MemoryError + + return output_arr + + +cdef int _contrast_core(const double* posimg, const double* negimg, + unsigned char* mask, int* output, + Py_ssize_t ny, Py_ssize_t nx, bint conn8, + double contrast, double source_sum, + double source_min, + Py_ssize_t n_max_labels) noexcept nogil: + """ + Run the watershed contrast loop for one source in place. + + ``output`` holds the markers on input and the final relabeled + (consecutive from 1) basins on output. Returns the number of + final labels, or -1 if a memory allocation failed, or -2 if the + flooded basins do not cover the segment mask (a connectivity + mismatch). + + The loop replicates the NumPy implementation operation for + operation. The basin fluxes are accumulated in raster order in + float64 (as np.bincount does), the below-contrast markers are + removed one at a time or in the largest provably equivalent + batch of the faintest markers, and NaN basin fluxes compare + false against every threshold, with np.argmin's first-NaN + behavior for the single-marker removal. The only divergence is + the order of bitwise-equal basin fluxes in the batch sort (NumPy + uses an unstable sort there, while this uses a stable one). + """ + cdef Py_ssize_t n_tot = ny * nx + cdef Py_ssize_t n_cap = n_max_labels + 1 + cdef long long* counts = malloc( + n_cap * sizeof(long long)) + cdef double* flux = malloc(n_cap * sizeof(double)) + cdef int* lab = malloc(n_cap * sizeof(int)) + cdef double* frac = malloc(n_cap * sizeof(double)) + cdef Py_ssize_t* order = malloc( + n_cap * sizeof(Py_ssize_t)) + cdef double* csum = malloc(n_cap * sizeof(double)) + cdef unsigned char* removed = malloc( + n_cap * sizeof(unsigned char)) + cdef int* lut = malloc(n_cap * sizeof(int)) + + cdef Py_ssize_t p, i, j, k, n_labels, n_remove, min_idx, last_ok + cdef Py_ssize_t pos + cdef int status, current + cdef bint remove_marker, a_nan, b_nan + cdef double value_a, value_b + + if (counts == NULL or flux == NULL or lab == NULL or frac == NULL + or order == NULL or csum == NULL or removed == NULL + or lut == NULL): + status = -1 + else: + status = 0 + + n_labels = 0 + while status == 0: + status = _watershed_core(negimg, mask, output, ny, nx, conn8) + if status != 0: + break + + # Present labels (ascending) and their fluxes, accumulated + # in raster order in float64 as np.bincount does. + for i in range(n_cap): + counts[i] = 0 + flux[i] = 0.0 + for p in range(n_tot): + current = output[p] + if current != 0: + counts[current] += 1 + flux[current] += posimg[p] + n_labels = 0 + for i in range(1, n_cap): + if counts[i] > 0: + lab[n_labels] = i + frac[n_labels] = flux[i] / source_sum + n_labels += 1 + + if n_labels == 1: # only 1 source left + break + + remove_marker = False + for i in range(n_labels): + if frac[i] < contrast: + remove_marker = True + break + if not remove_marker: + break + + # Remove the faintest below-contrast marker(s). See + # _SingleSourceDeblender._remove_faint_markers. + n_remove = 1 + if source_min >= 0 and n_labels > 2: + # Stable sort of the label positions by flux fraction, + # with NaN values sorting to the end. + for i in range(n_labels): + order[i] = i + for i in range(1, n_labels): + j = i + while j > 0: + value_a = frac[order[j]] + value_b = frac[order[j - 1]] + a_nan = isnan(value_a) + b_nan = isnan(value_b) + if ((not a_nan and b_nan) + or (not a_nan and value_a < value_b)): + pos = order[j - 1] + order[j - 1] = order[j] + order[j] = pos + j -= 1 + else: + break + csum[0] = frac[order[0]] + for i in range(1, n_labels): + csum[i] = csum[i - 1] + frac[order[i]] + # A batch of the n faintest markers (2 <= n < N) is + # valid if its total flux fraction is below both the + # contrast and the next-faintest marker flux fraction. + last_ok = -1 + for k in range(1, n_labels - 1): + if (csum[k] < contrast + and csum[k] < frac[order[k + 1]]): + last_ok = k + if last_ok >= 0: + n_remove = last_ok + 1 + + for i in range(n_cap): + removed[i] = 0 + if n_remove == 1: + # np.argmin: the first NaN if any, else the first minimum + min_idx = -1 + for i in range(n_labels): + if isnan(frac[i]): + min_idx = i + break + if min_idx == -1: + min_idx = 0 + for i in range(1, n_labels): + if frac[i] < frac[min_idx]: + min_idx = i + removed[lab[min_idx]] = 1 + else: + for j in range(n_remove): + removed[lab[order[j]]] = 1 + for p in range(n_tot): + if output[p] != 0 and removed[output[p]]: + output[p] = 0 + + if status == 0: + # The flooded basins must cover the segment mask exactly + # (they cannot with mismatched detection and deblending + # connectivities). + for p in range(n_tot): + if mask[p] and output[p] == 0: + status = -2 + break + + if status == 0: + # Relabel the surviving labels consecutively from 1 in + # ascending label order. + for i in range(n_labels): + lut[lab[i]] = (i + 1) + for p in range(n_tot): + if output[p] != 0: + output[p] = lut[output[p]] + status = n_labels + + free(counts) + free(flux) + free(lab) + free(frac) + free(order) + free(csum) + free(removed) + free(lut) + + return status + + +def deblend_source_contrast(const data_t[:, ::1] data, + const segm_t[:, ::1] segm_data, + long long label, Py_ssize_t y0, + Py_ssize_t y1, Py_ssize_t x0, + Py_ssize_t x1, markers, *, + int connectivity, double contrast, + double source_sum, double source_min): + """ + Apply the watershed contrast loop to one source's markers. + + The watershed flooding, the basin flux measurements, the + below-contrast marker removal (single or batched), and the final + consecutive relabeling all run in compiled code that releases + the GIL, producing results identical to the per-step NumPy + implementation in ``_SingleSourceDeblender``. + + Parameters + ---------- + data : 2D float `~numpy.ndarray` + The full data array. NaN pixels within the segment are flooded + after all the finite pixels, so they are assigned to a + neighboring basin. They contribute NaN to the flux of that + basin, as in the NumPy implementation. + + segm_data : 2D int `~numpy.ndarray` + The full segmentation array. + + label : int + The label of the source segment. + + y0, y1, x0, x1 : int + The bounding-box slice bounds of the source. + + markers : 2D int `~numpy.ndarray` + The marker image cutout, with markers labeled from 1. It is + not modified. + + connectivity : {8, 4} + The pixel connectivity. + + contrast : float + The contrast criterion (the minimum fraction of the total + source flux that a watershed basin must contain). + + source_sum : float + The total flux of the source segment (NaN pixels excluded). + + source_min : float + The minimum data value of the source segment (NaN pixels + excluded). + + Returns + ------- + output : 2D int `~numpy.ndarray` or `None` + The deblended cutout with consecutive labels starting from + 1, or `None` if only one basin remains after applying the + contrast criterion. + + Raises + ------ + ValueError + If the flooded basins do not cover the segment mask, which + happens when the detection and deblending connectivities + differ. + """ + cdef Py_ssize_t ny_c = y1 - y0 + cdef Py_ssize_t nx_c = x1 - x0 + cdef Py_ssize_t img_nx = data.shape[1] + + posimg_arr = np.empty((ny_c, nx_c), dtype=np.float64) + negimg_arr = np.empty((ny_c, nx_c), dtype=np.float64) + mask_arr = np.empty((ny_c, nx_c), dtype=np.uint8) + output_arr = np.array(markers, dtype=np.int32, copy=True, order='C') + + cdef double[:, ::1] posimg_mv = posimg_arr + cdef double[:, ::1] negimg_mv = negimg_arr + cdef unsigned char[:, ::1] mask_mv = mask_arr + cdef int[:, ::1] output_mv = output_arr + cdef double* posimg = &posimg_mv[0, 0] + cdef double* negimg = &negimg_mv[0, 0] + cdef unsigned char* mask = &mask_mv[0, 0] + cdef int* output = &output_mv[0, 0] + cdef const data_t* data_ptr = &data[0, 0] + cdef const segm_t* segm_ptr = &segm_data[0, 0] + + cdef bint conn8 = connectivity == 8 + cdef Py_ssize_t iy, ix, idx, p + cdef Py_ssize_t n_max_labels = 0 + cdef int status + + with nogil: + for iy in range(ny_c): + for ix in range(nx_c): + idx = (y0 + iy) * img_nx + x0 + ix + p = iy * nx_c + ix + posimg[p] = data_ptr[idx] + if isnan(posimg[p]): + # NaN pixels are flooded after all finite pixels + negimg[p] = INFINITY + else: + negimg[p] = -posimg[p] + mask[p] = segm_ptr[idx] == label + if output[p] > n_max_labels: + n_max_labels = output[p] + status = _contrast_core(posimg, negimg, mask, output, ny_c, + nx_c, conn8, contrast, source_sum, + source_min, n_max_labels) + + if status == -1: + raise MemoryError + if status == -2: + msg = (f'Deblending failed for source {int(label)!r}. ' + 'Please ensure you used the same pixel connectivity ' + 'in detect_sources and deblend_sources.') + raise ValueError(msg) + if status == 1: # no deblending + return None + return output_arr diff --git a/photutils/segmentation/catalog.py b/photutils/segmentation/catalog.py index 5d242ab6b7..5ed92a6dae 100644 --- a/photutils/segmentation/catalog.py +++ b/photutils/segmentation/catalog.py @@ -272,9 +272,10 @@ def _batch_gini(values): sources. This is the vectorized form of `photutils.morphology.gini` applied - to each source's unmasked (finite) absolute pixel values: NaN for a - source with no pixels, 0.0 for a single pixel or a zero mean, and - otherwise the Lotz et al. (2004) sum over the sorted values. + to each source's unmasked (finite) absolute pixel values. The + result is NaN for a source with no pixels, 0.0 for a single pixel + or a zero mean, and otherwise the Lotz et al. (2004) sum over the + sorted values. The values of each source are sorted with NumPy in a Python loop rather than with the compiled ``batch_gini`` kernel shared with @@ -664,9 +665,9 @@ class SourceCatalog: **Scalar vs. Multi-source Catalogs** A `SourceCatalog` can represent a single source or multiple - sources. Most properties adapt their return type accordingly: for + sources. Most properties adapt their return type accordingly. For a multi-source catalog, properties return arrays or lists (one - element per source); for a single-source (scalar) catalog, the + element per source). For a single-source (scalar) catalog, the same properties return a scalar value or a single object. For example, `kron_aperture` returns a list of aperture objects for a multi-source catalog, but a single aperture object for a scalar @@ -790,7 +791,7 @@ def _validate_array(self, array, name, *, shape=True): array = None if array is not None: # UFuncTypeError is raised when subtracting float - # local_background from int data; convert to float + # local_background from int data, so convert to float array = np.asanyarray(array) if array.ndim != 2: msg = f'{name} must be a 2D array' @@ -1059,7 +1060,7 @@ def __getattr__(self, name): since='3.0', until='4.0') def __getattribute__(self, name): - # Centralized scalar collapse: for a scalar (single-source) + # Centralized scalar collapse. For a scalar (single-source) # catalog, public properties that return a length-1 # array/list/tuple (or a scalar SkyCoord) are returned as a # scalar value. Private ('_'-prefixed) attributes are always @@ -2322,7 +2323,7 @@ def _centroid_err_cov(self): pixel_var = 1.0 / 12.0 # Ignore divide-by-zero and invalid-value RuntimeWarnings for - # sources with non-positive or non-finite total flux; those + # sources with non-positive or non-finite total flux. Those # values are replaced by NaN below. with warnings.catch_warnings(): warnings.simplefilter('ignore', RuntimeWarning) @@ -2785,8 +2786,8 @@ def _centroid_win_err_cov(self): windowed centroid, ``[[var_x, cov_xy], [cov_xy, var_y]]``. Fallback sources hold the isophotal covariance (see - ``_centroid_err_cov``); matrices are all-NaN where errors are - unavailable. + ``_centroid_err_cov``). The matrices are all-NaN where errors + are unavailable. """ results = self._centroid_win_results cov = np.empty((len(results), 2, 2)) @@ -3886,7 +3887,7 @@ def _singular_covariance_flag_mask(self): This is the mask used for the ``'singular_covariance'`` flag. It matches the equivalent aperture flag (see - `~photutils.aperture.decode_aperture_flags`): in addition to + `~photutils.aperture.decode_aperture_flags`). In addition to the determinant test used by ``_singular_covariance_mask``, a source is flagged when its minor-axis variance (the smaller eigenvalue of the raw covariance matrix) is less than ``1 / @@ -4998,7 +4999,7 @@ def _calc_circular_photometry(self, radius): fcounts = result.flag_counts # Membership matches the previous cutout-based - # ``in_aperture & ~mask`` rule: under 'correct', uncorrectable + # ``in_aperture & ~mask`` rule. Under 'correct', uncorrectable # neighbor pixels stay members (with value zero) members = fcounts[:, FLAG_COL_VALID].copy() if seg_code == 3: @@ -5043,8 +5044,9 @@ def _kron_aperture_params(self, kron_params): defined : 1D `~numpy.ndarray` (bool) `False` where the aperture is undefined (`None` in - `kron_aperture`): where the source is completely masked or - its centroid or elliptical shape parameters are not finite. + `kron_aperture`). The aperture is undefined where the + source is completely masked or its centroid or elliptical + shape parameters are not finite. """ # NOTE: if kron_radius = NaN, scale = NaN and the aperture is # undefined @@ -5171,7 +5173,7 @@ def _calc_kron_photometry(self, *, kron_params=None): clipped = fcounts[:, FLAG_COL_BBOX_CLIPPED].astype(bool) weights_out = weights_out.astype(bool) # Membership matches the previous cutout-based - # ``in_aperture & ~mask`` rule: under 'correct', + # ``in_aperture & ~mask`` rule. Under 'correct', # uncorrectable neighbor pixels stay members (with value # zero), so they keep the flux at 0.0 rather than NaN members = fcounts[:, FLAG_COL_VALID].copy() diff --git a/photutils/segmentation/core.py b/photutils/segmentation/core.py index 524e9e136b..360968fce5 100644 --- a/photutils/segmentation/core.py +++ b/photutils/segmentation/core.py @@ -240,7 +240,7 @@ def __getitem__(self, key): and all(isinstance(key[i], slice) for i in (0, 1))): result = self.data[key] if result.size == 0: - msg = ('The sliced result is empty; cannot create ' + msg = ('The sliced result is empty. Cannot create ' 'a SegmentationImage with zero size') raise ValueError(msg) return SegmentationImage(result) @@ -1800,7 +1800,7 @@ def _geojson_polygons(self): polygon_dict[int(label)].append(polygon) # Check that the polygon labels match the segmentation image - # labels; this is a sanity check to ensure that the rasterio + # labels. This is a sanity check to ensure that the rasterio # library is working correctly. # Note that polygons have been sorted by label. if not np.all(np.array(list(polygon_dict.keys())) == self.labels): diff --git a/photutils/segmentation/deblend.py b/photutils/segmentation/deblend.py index 14edd5d4bd..8bd28066de 100644 --- a/photutils/segmentation/deblend.py +++ b/photutils/segmentation/deblend.py @@ -5,28 +5,29 @@ """ import warnings -from concurrent.futures import ProcessPoolExecutor, as_completed from dataclasses import dataclass -from functools import cached_property, partial -from multiprocessing import cpu_count, get_context import numpy as np from astropy.units import Quantity -from scipy.ndimage import label as ndi_label -from scipy.ndimage import sum_labels +from photutils.segmentation._deblend_markers import (deblend_markers_chunk, + deblend_source_extrema) +from photutils.segmentation._deblend_watershed import deblend_source_contrast from photutils.segmentation.core import (SegmentationImage, _get_labels, _remap_deblend_label_map) -from photutils.segmentation.detect import _detect_sources from photutils.segmentation.flags import SEGMENTATION_FLAGS from photutils.segmentation.utils import _make_binary_structure from photutils.utils._deprecation import deprecated_renamed_argument -from photutils.utils._progress_bars import add_progress_bar, tqdm -from photutils.utils._stats import nanmax, nanmin, nansum +from photutils.utils._stats import nansum from photutils.utils.exceptions import DeblendWarning __all__ = ['deblend_sources'] +# The number of markers above which a source is deblended again with +# linearly spaced threshold levels, which have fewer levels at low +# thresholds. The value is arbitrary but works well in practice +_MAX_MARKERS = 200 + @dataclass class _DeblendParams: @@ -41,11 +42,15 @@ class _DeblendParams: until='4.0') @deprecated_renamed_argument('npixels', 'n_pixels', '3.0', until='4.0') @deprecated_renamed_argument('nlevels', 'n_levels', '3.0', until='4.0') -@deprecated_renamed_argument('nproc', 'n_processes', '3.0', until='4.0') +@deprecated_renamed_argument('nproc', None, '3.0', until='4.0') +@deprecated_renamed_argument('n_processes', None, '3.1', until='4.0') +@deprecated_renamed_argument('progress_bar', None, '3.1', until='4.0') def deblend_sources(data, segmentation_image, n_pixels, *, labels=None, n_levels=32, contrast=0.001, mode='exponential', - connectivity=8, relabel=True, n_processes=1, - progress_bar=True): + connectivity=8, relabel=True, + nproc=1, # noqa: ARG001 + n_processes=1, # noqa: ARG001 + progress_bar=True): # noqa: ARG001 """ Deblend overlapping sources labeled in a segmentation image. @@ -60,7 +65,10 @@ def deblend_sources(data, segmentation_image, n_pixels, *, labels=None, data : 2D `~numpy.ndarray` The 2D array of the image. If filtering is desired, please input a convolved image here. This array should be the same array used - in `~photutils.segmentation.detect_sources`. + in `~photutils.segmentation.detect_sources`. NaN pixels within + a source segment are excluded from the multithreshold levels + and the source flux, and are assigned to a neighboring + deblended source after all finite pixels have been assigned. segmentation_image : `~photutils.segmentation.SegmentationImage` The segmentation image to deblend. @@ -98,8 +106,8 @@ def deblend_sources(data, segmentation_image, n_pixels, *, labels=None, the threshold levels between the source minimum and maximum. The ``'exponential'`` and ``'sinh'`` modes differ in that the ``'exponential'`` levels are dependent on the source - maximum/minimum ratio (smaller ratios are more linear; larger - ratios are more exponential), while the ``'sinh'`` levels + maximum/minimum ratio (smaller ratios are more linear and + larger ratios are more exponential), while the ``'sinh'`` levels are not. Also, the ``'exponential'`` mode will be changed to ``'linear'`` for sources with non-positive minimum data values. @@ -116,24 +124,32 @@ def deblend_sources(data, segmentation_image, n_pixels, *, labels=None, relabeled such that the labels are in consecutive order starting from 1. + nproc : int, optional + This keyword is deprecated and has no effect. It was the name of + the ``n_processes`` keyword before version 3.0. + + .. deprecated:: 3.0 + The ``nproc`` keyword is deprecated and will be removed in + version 4.0. + n_processes : int, optional - The number of processes to use for multiprocessing (if larger - than 1). If set to 1, then a serial implementation is used - instead of a parallel one. If `None`, then the number of - processes will be set to the number of CPUs detected on the - machine. Please note that due to overheads, multiprocessing may - be slower than serial processing if only a small number of - sources are to be deblended. The benefits of multiprocessing - require ~1000 or more sources to deblend, with larger gains as - the number of sources increase. + This keyword is deprecated and has no effect. Multiprocessing + no longer provides any benefit. The deblending computation is + now dominated by compiled code, and the process startup and + data-pickling overheads of multiprocessing made it slower than + the serial implementation. + + .. deprecated:: 3.1 + The ``n_processes`` keyword is deprecated and will be + removed in version 4.0. progress_bar : bool, optional - Whether to display a progress bar. If ``n_processes = 1``, then the - ID shown after the progress bar is the source label being - deblended. If multiprocessing is used (``n_processes > 1``), the ID - shown is the last source label that was deblended. The progress - bar requires that the `tqdm `_ optional - dependency be installed. + This keyword is deprecated and has no effect. Deblending no + longer displays a progress bar. + + .. deprecated:: 3.1 + The ``progress_bar`` keyword is deprecated and will be + removed in version 4.0. Returns ------- @@ -181,8 +197,8 @@ def deblend_sources(data, segmentation_image, n_pixels, *, labels=None, msg = f'n_pixels must be a positive integer, got {n_pixels!r}' raise ValueError(msg) - if n_levels < 1: - msg = 'n_levels must be >= 1' + if (n_levels < 1) or (int(n_levels) != n_levels): + msg = f'n_levels must be a positive integer, got {n_levels!r}' raise ValueError(msg) if contrast < 0 or contrast > 1: msg = 'contrast must be >= 0 and <= 1' @@ -204,8 +220,8 @@ def deblend_sources(data, segmentation_image, n_pixels, *, labels=None, labels = np.atleast_1d(labels) segmentation_image.check_labels(labels) - # Include only sources that have at least (2 * n_pixels); - # this is required for a source to be deblended into multiple + # Include only sources that have at least (2 * n_pixels). + # This is required for a source to be deblended into multiple # sources, each with a minimum of n_pixels mask = (segmentation_image.areas[ segmentation_image.get_indices(labels)] @@ -218,111 +234,48 @@ def deblend_sources(data, segmentation_image, n_pixels, *, labels=None, segm_deblended = segmentation_image.data.copy() label_indices = segmentation_image.get_indices(labels) + all_slices = [segmentation_image.slices[idx] for idx in label_indices] + + # Contiguous, native-byte-order views of the inputs for the + # compiled kernels. Casting the non-float data dtypes to float64 is + # exact for the threshold comparisons, matching the NumPy promotion + # rules + if data.dtype.type in (np.float32, np.float64): + driver_dtype = data.dtype.newbyteorder('=') + else: + driver_dtype = np.float64 + driver_data = np.ascontiguousarray(data, dtype=driver_dtype) + segm_data = segmentation_image.data + if segm_data.dtype.type in (np.int32, np.int64): + driver_dtype = segm_data.dtype.newbyteorder('=') + else: + driver_dtype = np.int64 + driver_segm = np.ascontiguousarray(segm_data, dtype=driver_dtype) - if n_processes is None: - n_processes = cpu_count() + results = _deblend_sources_chunk(data, segm_data, driver_data, + driver_segm, labels, all_slices, + deblend_params) deblend_label_map = {} max_label = segmentation_image.max_label - if n_processes == 1: - if progress_bar: - desc = 'Deblending' - label_indices = add_progress_bar(label_indices, desc=desc) - - nonposmin_labels = [] - n_markers_labels = [] - for label, label_idx in zip(labels, label_indices, strict=True): - if not isinstance(label_indices, np.ndarray): - label_indices.set_postfix_str(f'ID: {label}') - source_slice = segmentation_image.slices[label_idx] - source_data = data[source_slice] - source_segment = segmentation_image.data[source_slice] - source_deblended, warns = _deblend_source(source_data, - source_segment, - label, - deblend_params) - - if warns: - if 'nonposmin' in warns: - nonposmin_labels.append(label) - if 'n_markers' in warns: - n_markers_labels.append(label) - - if source_deblended is not None: - source_mask = source_deblended > 0 - new_segm = source_deblended[source_mask] # min label = 1 - segm_deblended[source_slice][source_mask] = ( - new_segm + max_label) - new_labels = _get_labels(new_segm) + max_label - deblend_label_map[label] = new_labels - max_label += len(new_labels) - - else: - # Use multiprocessing to deblend sources - - # Prepare the arguments for the worker function - all_source_data = [] - all_source_segments = [] - all_source_slices = [] - for label_idx in label_indices: - source_slice = segmentation_image.slices[label_idx] - source_data = data[source_slice] - source_segment = segmentation_image.data[source_slice] - all_source_data.append(source_data) - all_source_segments.append(source_segment) - all_source_slices.append(source_slice) - - args_all = zip(all_source_data, all_source_segments, labels, - strict=True) - - # Create a partial function to pass the deblend_params to the - # worker function - worker = partial(_deblend_source, deblend_params=deblend_params) - - # Prepare to store futures and results to preserve the input - # order of the labels when using as_completed() - futures_dict = {} - results = [None] * len(labels) - - disable_pbar = not progress_bar - mp_context = get_context('spawn') - with ProcessPoolExecutor(mp_context=mp_context, - max_workers=n_processes) as executor: - # Submit all jobs at once - for index, args in enumerate(args_all): - futures_dict[executor.submit(worker, *args)] = index - - with tqdm(total=len(labels), desc='Deblending', - disable=disable_pbar) as pbar: - # Process the results as they are completed - for future in as_completed(futures_dict): - pbar.update(1) - idx = futures_dict[future] - pbar.set_postfix_str(f'ID: {labels[idx]}') - results[idx] = future.result() - - # Process the results - nonposmin_labels = [] - n_markers_labels = [] - for label, source_slice, source_deblended in zip(labels, - all_source_slices, - results, strict=True): - source_deblended, warns = source_deblended - - if warns: - if 'nonposmin' in warns: - nonposmin_labels.append(label) - if 'n_markers' in warns: - n_markers_labels.append(label) - - if source_deblended is not None: - source_mask = source_deblended > 0 - new_segm = source_deblended[source_mask] # min label = 1 - segm_deblended[source_slice][source_mask] = ( - new_segm + max_label) - new_labels = _get_labels(new_segm) + max_label - deblend_label_map[label] = new_labels - max_label += len(new_labels) + nonposmin_labels = [] + n_markers_labels = [] + for label, source_slice, (source_deblended, warns) in zip( + labels, all_slices, results, strict=True): + if warns: + if 'nonposmin' in warns: + nonposmin_labels.append(label) + if 'n_markers' in warns: + n_markers_labels.append(label) + + if source_deblended is not None: + source_mask = source_deblended > 0 + new_segm = source_deblended[source_mask] # min label = 1 + segm_deblended[source_slice][source_mask] = ( + new_segm + max_label) + new_labels = _get_labels(new_segm) + max_label + deblend_label_map[label] = new_labels + max_label += len(new_labels) if nonposmin_labels or n_markers_labels: msg = ('The deblending mode of one or more source labels from the ' @@ -355,315 +308,228 @@ def deblend_sources(data, segmentation_image, n_pixels, *, labels=None, return segm_img -def _deblend_source(data, segment_data, label, deblend_params): +def _linspace_rows(start, stop, num): """ - Convenience function to deblend a single labeled source. + Evaluate ``np.linspace(start[i], stop[i], num)`` for every row. + + This replicates the NumPy implementation operation for operation + (the evaluation dtype, the step computation, and the zero-step + special case), so that every row is bitwise identical to the + scalar ``np.linspace`` call of the pure-Python reference + implementation. + + Parameters + ---------- + start, stop : 1D `~numpy.ndarray` + The start and stop values of each row. + + num : int + The number of samples per row. + + Returns + ------- + result : 2D `~numpy.ndarray` + The samples, with shape ``(len(start), num)``. """ - deblender = _SingleSourceDeblender(data, segment_data, label, - deblend_params) - return deblender.deblend_source(), deblender.warnings + dtype = np.result_type(start, stop) + if not np.issubdtype(dtype, np.inexact): + dtype = np.float64 + div = num - 1 + delta = np.subtract(stop, start, dtype=dtype) + samples = np.arange(0, num, dtype=dtype) + step = delta / div + result = samples[None, :] * step[:, None] + zero_step = step == 0 + if np.any(zero_step): + # The np.linspace special case for subnormal steps + result[zero_step] = (samples / div)[None, :] * delta[zero_step, None] + result += start[:, None] + result[:, -1] = stop + return result + + +def _compute_thresholds(source_min, source_max, n_levels, mode): + """ + Compute the multithreshold levels of a set of sources. + This is the vectorized form of the per-source threshold + computation of the pure-Python reference implementation + (``_SingleSourceDeblender.compute_thresholds``). It performs the + same NumPy operations in the same dtypes, so the levels are bitwise + identical to the reference implementation on every platform. -class _SingleSourceDeblender: + Parameters + ---------- + source_min, source_max : 1D `~numpy.ndarray` + The minimum and maximum data value of each source segment, in + the data dtype. Each maximum must be larger than its minimum. + + n_levels : int + The number of levels per source. + + mode : {'exponential', 'linear', 'sinh'} + The level spacing. For the ``'exponential'`` mode, the sources + with a non-positive minimum use the ``'linear'`` spacing. + + Returns + ------- + thresholds : 2D float64 `~numpy.ndarray` + The levels of each source, with shape ``(n_sources, + n_levels)`` and ascending along the second axis. The source + minimum and maximum are excluded. + + nonposmin : 1D bool `~numpy.ndarray` + Whether each source fell back to the ``'linear'`` spacing + because of a non-positive minimum. + """ + source_min = np.asarray(source_min) + source_max = np.asarray(source_max) + nonposmin = np.zeros(source_min.shape, dtype=bool) + if mode == 'exponential': + nonposmin = source_min <= 0 + + thresholds = _linspace_rows(source_min, source_max, n_levels + 2) + if mode != 'linear': + delta = source_max - source_min + normalized = (thresholds - source_min[:, None]) / delta[:, None] + if mode == 'sinh': + a = 0.25 + thresholds = np.sinh(normalized / a) / np.sinh(1.0 / a) + thresholds *= delta[:, None] + thresholds += source_min[:, None] + else: + keep = ~nonposmin + ratio = source_max[keep] / source_min[keep] + thresholds[keep] = (source_min[keep, None] + * ratio[:, None] ** normalized[keep]) + + # Do not include the source minimum and maximum + return (np.ascontiguousarray(thresholds[:, 1:-1], dtype=np.float64), + nonposmin) + + +def _deblend_sources_chunk(data, segm_data, driver_data, driver_segm, + labels, slices, deblend_params): """ - Class to deblend a single labeled source. + Deblend a chunk of labeled sources. + + The data extrema of every source in the chunk are computed by a + compiled kernel, the multithreshold levels and the mode fallbacks + are computed for all sources at once in NumPy, the level + quantization and the marker construction run in the compiled chunk + driver (which releases the GIL and reuses one workspace across the + chunk), and the watershed and contrast steps then run for the + sources that split into two or more markers. Parameters ---------- data : 2D `~numpy.ndarray` - The cutout data array for a single source. ``data`` should - also already be smoothed by the same filter used in - :func:`~photutils.segmentation.detect_sources`, if applicable. + The data array. - segment_data : 2D int `~numpy.ndarray` - The cutout segmentation image for a single source. Must have the - same shape as ``data``. + segm_data : 2D int `~numpy.ndarray` + The segmentation array. - label : int - The label of the source to deblend. This is needed because there - may be more than one source label within the cutout. + driver_data, driver_segm : 2D `~numpy.ndarray` + Contiguous views (or exact casts) of ``data`` and + ``segm_data`` with dtypes supported by the compiled kernels. + + labels : 1D `~numpy.ndarray` + The labels of the sources in the chunk. + + slices : list of tuple of slice + The bounding-box slices of the sources in the chunk. deblend_params : `_DeblendParams` - The parameters for deblending the source. - """ + The parameters for deblending the sources. - def __init__(self, data, segment_data, label, deblend_params): - self.data = data - self.segment_data = segment_data - self.label = label - self.n_pixels = deblend_params.n_pixels - self.footprint = deblend_params.footprint - self.n_levels = deblend_params.n_levels - self.contrast = deblend_params.contrast - self.mode = deblend_params.mode - - self.segment_mask = segment_data == label - data_values = data[self.segment_mask] - self.source_min = nanmin(data_values) - self.source_max = nanmax(data_values) - self.source_sum = nansum(data_values) - self.warnings = {} - - @cached_property - def linear_thresholds(self): - """ - Linearly spaced thresholds between the source minimum and - maximum (inclusive). - - The source min/max are excluded later, giving n_levels - thresholds between min and max (noninclusive). - """ - return np.linspace(self.source_min, self.source_max, self.n_levels + 2) - - @cached_property - def normalized_thresholds(self): - """ - Normalized thresholds (from 0 to 1) between the source minimum - and maximum (inclusive). - """ - return ((self.linear_thresholds - self.source_min) - / (self.source_max - self.source_min)) - - def compute_thresholds(self): - """ - Compute the multi-level detection thresholds for the source. - - Returns - ------- - thresholds : 1D `~numpy.ndarray` - The multi-level detection thresholds for the source. - """ - if self.mode == 'exponential' and self.source_min <= 0: - self.warnings['nonposmin'] = 'non-positive minimum' - self.mode = 'linear' - - if self.mode == 'linear': - thresholds = self.linear_thresholds - elif self.mode == 'sinh': - a = 0.25 - minval = self.source_min - maxval = self.source_max - thresholds = self.normalized_thresholds - thresholds = np.sinh(thresholds / a) / np.sinh(1.0 / a) - thresholds *= (maxval - minval) - thresholds += minval - elif self.mode == 'exponential': - minval = self.source_min - maxval = self.source_max - thresholds = self.normalized_thresholds - thresholds = minval * (maxval / minval) ** thresholds - - return thresholds[1:-1] # do not include source min and max - - def multithreshold(self): - """ - Perform multithreshold detection for each source. - - This method is useful for debugging and testing. - - Returns - ------- - segments : list of 2D `~numpy.ndarray` or `None` - A list of segmentation images, one for each threshold. - `None` is returned for thresholds that do not have more than - one label. - """ - thresholds = self.compute_thresholds() - segms = [] - for threshold in thresholds: - segm = _detect_sources(self.data, threshold, self.n_pixels, - self.footprint, self.segment_mask, - relabel=False, return_segmimg=False) - segms.append(segm) - return segms - - def make_markers(self, *, return_all=False): - """ - Make markers (possible sources) for the watershed algorithm. - - Parameters - ---------- - return_all : bool, optional - If `False` then return only the final segmentation marker - image. If `True` then return all segmentation marker images. - This keyword is useful for debugging and testing. - - Returns - ------- - markers : 2D `~numpy.ndarray` or list of 2D `~numpy.ndarray` - A segmentation image that contain markers for possible - sources. If ``return_all=True`` then a list of all - segmentation marker images is returned. `None` is returned - if there is only one source at every threshold. - """ - thresholds = self.compute_thresholds() - segm_lower = _detect_sources(self.data, thresholds[0], self.n_pixels, - self.footprint, self.segment_mask, - relabel=False, return_segmimg=False) - - if return_all: - all_segms = [segm_lower] - - for threshold in thresholds[1:]: - segm_upper = _detect_sources(self.data, threshold, self.n_pixels, - self.footprint, self.segment_mask, - relabel=False, return_segmimg=False) - if segm_upper is None: # 0 or 1 labels - continue - - segm_lower = self.make_marker_segment(segm_lower, segm_upper) - - if return_all: - all_segms.append(segm_lower) - - if return_all: - return all_segms - - return segm_lower - - def make_marker_segment(self, segment_lower, segment_upper): - """ - Make markers (possible sources) for the watershed algorithm. - - Parameters - ---------- - segment_lower : 2D `~numpy.ndarray` - The "lower" threshold level segmentation image. - - segment_upper : 2D `~numpy.ndarray` - The next-highest threshold level segmentation image. - - Returns - ------- - markers : 2D `~numpy.ndarray` - A segmentation image that contain markers for possible - sources. - - Notes - ----- - For a given label in the lower level, find the labels in the - upper level (higher threshold value) that are its children - (i.e., the labels within the same mask as the lower level). If - there are multiple children, then the lower-level parent label - is replaced by its children. Parent labels that do not have - multiple children in the upper level are kept as is (maximizing - the marker size). - """ - if segment_lower is None: - return segment_upper - - labels = _get_labels(segment_lower) - new_markers = False - markers = segment_lower.astype(bool) - for label in labels: - mask = (segment_lower == label) - # Find label mapping from the lower to upper level - upper_labels = _get_labels(segment_upper[mask]) - if upper_labels.size >= 2: # new child markers found - new_markers = True - markers[mask] = segment_upper[mask].astype(bool) - - if new_markers: - # Convert bool markers to integer labels - return ndi_label(markers, structure=self.footprint)[0] - - return segment_lower - - def apply_watershed(self, markers): - """ - Apply the watershed algorithm to the source markers. - - Parameters - ---------- - markers : list of `~photutils.segmentation.SegmentationImage` - A list of segmentation images that contain possible sources - as markers. The last list element contains all the potential - source markers. - - Returns - ------- - segment_data : 2D int `~numpy.ndarray` - A 2D int array containing the deblended source labels. Note - that the source labels may not be consecutive if a label was - removed. - """ - from skimage.segmentation import watershed - - # Deblend using watershed. If any source does not meet the contrast - # criterion, then remove the faintest such source and repeat until - # all sources meet the contrast criterion. - remove_marker = True - while remove_marker: - markers = watershed(-self.data, markers, mask=self.segment_mask, - connectivity=self.footprint) - - labels = _get_labels(markers) - if labels.size == 1: # only 1 source left - remove_marker = False - else: - flux_frac = (sum_labels(self.data, markers, index=labels) - / self.source_sum) - remove_marker = any(flux_frac < self.contrast) - - if remove_marker: - # Remove only the faintest source (one at a time) - # because several faint sources could combine to meet - # the contrast criterion - markers[markers == labels[np.argmin(flux_frac)]] = 0.0 - - return markers - - def deblend_source(self): - """ - Deblend a single labeled source. - - Returns - ------- - segment_data : 2D int `~numpy.ndarray` - A 2D int array containing the deblended source labels. The - source labels are consecutive starting at 1. - """ - if self.source_min == self.source_max: # no deblending - return None - - # Define the markers (possible sources) for the watershed algorithm - markers = self.make_markers() + Returns + ------- + results : list of (2D `~numpy.ndarray` or `None`, dict) + For each source, the deblended cutout (`None` if the source + did not deblend) and the mode-fallback warnings dictionary. + """ + labels = np.asarray(labels, dtype=np.int64) + y0 = np.array([slc[0].start for slc in slices], dtype=np.int64) + y1 = np.array([slc[0].stop for slc in slices], dtype=np.int64) + x0 = np.array([slc[1].start for slc in slices], dtype=np.int64) + x1 = np.array([slc[1].stop for slc in slices], dtype=np.int64) + connectivity = 8 if deblend_params.footprint[0, 0] else 4 + n_levels = int(deblend_params.n_levels) + mode = deblend_params.mode + chunk_kwargs = {'n_pixels': int(deblend_params.n_pixels), + 'connectivity': connectivity} + + source_min, source_max = deblend_source_extrema( + driver_data, driver_segm, labels, y0, y1, x0, x1) + + # Constant (or all-NaN) sources do not deblend. The multithreshold + # levels of the other sources are computed in the data dtype, as + # the pure-Python reference implementation does + active = np.flatnonzero(source_min < source_max) + markers_list = [None] * len(labels) + nonposmin = np.zeros(len(labels), dtype=bool) + n_markers_fallback = np.zeros(len(labels), dtype=bool) + if active.size > 0: + values_dtype = data.dtype.newbyteorder('=') + smin = source_min[active].astype(values_dtype) + smax = source_max[active].astype(values_dtype) + thresholds, fallback = _compute_thresholds(smin, smax, n_levels, + mode) + nonposmin[active] = fallback + markers, n_markers = deblend_markers_chunk( + driver_data, driver_segm, labels[active], y0[active], + y1[active], x0[active], x1[active], thresholds, + max_markers=_MAX_MARKERS, **chunk_kwargs) + for index, source_markers in zip(active, markers, strict=True): + markers_list[index] = source_markers + + # Too many markers make the watershed step very slow, so such + # sources are deblended again with linearly spaced levels. A + # source that already fell back to the linear spacing keeps its + # markers + if mode != 'linear': + retry = np.flatnonzero((n_markers > _MAX_MARKERS) & ~fallback) + if retry.size > 0: + thresholds, _ = _compute_thresholds( + smin[retry], smax[retry], n_levels, 'linear') + retry = active[retry] + markers, _ = deblend_markers_chunk( + driver_data, driver_segm, labels[retry], y0[retry], + y1[retry], x0[retry], x1[retry], thresholds, + max_markers=-1, **chunk_kwargs) + for index, source_markers in zip(retry, markers, + strict=True): + markers_list[index] = source_markers + n_markers_fallback[retry] = True + + results = [] + for index, (label, slc, markers) in enumerate( + zip(labels, slices, markers_list, strict=True)): + warns = {} + if nonposmin[index]: + warns['nonposmin'] = 'non-positive minimum' + if n_markers_fallback[index]: + warns['n_markers'] = 'too many markers' if markers is None: - return None - - # If there are too many markers (e.g., due to low threshold - # and/or small n_pixels), the watershed step can be very slow - # (the threshold of 200 is arbitrary, but seems to work well). - # This mostly affects the "exponential" mode, where there are - # many levels at low thresholds, so here we try again with - # "linear" mode. - n_labels = len(_get_labels(markers)) - if self.mode != 'linear' and n_labels > 200: - del markers # free memory - self.warnings['n_markers'] = 'too many markers' - self.mode = 'linear' - markers = self.make_markers() - if markers is None: - return None - - # Deblend using the watershed algorithm using the markers as seeds - markers = self.apply_watershed(markers) - - if not np.array_equal(self.segment_mask, markers.astype(bool)): - msg = (f'Deblending failed for source {self.label!r}. ' - 'Please ensure you used the same pixel connectivity ' - 'in detect_sources and deblend_sources.') - raise ValueError(msg) - - if len(_get_labels(markers)) == 1: # no deblending - return None - - # Markers may not be consecutive if a label was removed due to - # the contrast criterion - relabel_map = _create_relabel_map(markers, start_label=1) - if relabel_map is not None: - markers = relabel_map[markers] - return markers + results.append((None, warns)) + continue + + # The total source flux is computed with the same reduction as + # the per-source Python path (its rounding depends on the + # summation order) and the source minimum comes from the + # compiled extrema. Both are passed to the compiled contrast + # loop. + values = data[slc][segm_data[slc] == label] + source_deblended = deblend_source_contrast( + driver_data, driver_segm, int(label), int(y0[index]), + int(y1[index]), int(x0[index]), int(x1[index]), markers, + connectivity=connectivity, + contrast=float(deblend_params.contrast), + source_sum=float(nansum(values)), + source_min=float(source_min[index])) + results.append((source_deblended, warns)) + + return results def _make_flags_map(deblend_label_map, nonposmin_labels, n_markers_labels, diff --git a/photutils/segmentation/detect.py b/photutils/segmentation/detect.py index a01186223b..523d66ec35 100644 --- a/photutils/segmentation/detect.py +++ b/photutils/segmentation/detect.py @@ -137,8 +137,7 @@ def detect_threshold(data, n_sigma, *, background=None, error=None, mask=None, return threshold -def _detect_sources(data, threshold, n_pixels, footprint, inverse_mask, *, - relabel=True, return_segmimg=True): +def _detect_sources(data, threshold, n_pixels, footprint, inverse_mask): """ Detect sources above a specified threshold value in an image. @@ -150,8 +149,7 @@ def _detect_sources(data, threshold, n_pixels, footprint, inverse_mask, *, `detect_sources` in that it does not perform any boilerplate checks, it accepts a ``footprint`` argument instead of a ``connectivity`` argument, and it accepts an ``inverse_mask`` argument instead of a - ``mask`` argument. It is also used by the source deblending function - for multithresholding. + ``mask`` argument. Parameters ---------- @@ -180,31 +178,17 @@ def _detect_sources(data, threshold, n_pixels, footprint, inverse_mask, *, `False` values indicate masked pixels (the inverse of usual pixel masks). Masked pixels will not be included in any source. - relabel : bool, optional - If `True`, relabel the segmentation image with consecutive - numbers. - - return_segmimg : bool, optional - If `True`, return a `~photutils.segmentation.SegmentationImage` - object. If `False`, return a 2D `~numpy.ndarray` segmentation - image. The latter is used by the source deblending function. - In that case, if only one source is found, then `None` is - returned. - Returns ------- - segment_image : `~photutils.segmentation.SegmentationImage`, \ - 2D `~numpy.ndarray`, or `None` + segment_image : `~photutils.segmentation.SegmentationImage` or `None` A 2D segmentation image, with the same shape as ``data``, where sources are marked by different positive integer values. A value - of zero is reserved for the background. If ``return_segmimg`` - is `False`, then a 2D `~numpy.ndarray` segmentation image is - returned. If no sources are found then `None` is returned. + of zero is reserved for the background. If no sources are found + then `None` is returned. """ - # Ignore RuntimeWarning caused by > comparison when data contains NaNs - with warnings.catch_warnings(): - warnings.simplefilter('ignore', category=RuntimeWarning) - segment_img = data > threshold + # NaN values compare as False, so NaN pixels are never included + # in any source + segment_img = data > threshold if inverse_mask is not None: segment_img &= inverse_mask @@ -239,35 +223,23 @@ def _detect_sources(data, threshold, n_pixels, footprint, inverse_mask, *, segm_slices.append(slc) segm_areas.append(area) - if np.count_nonzero(segment_img) == 0: + if not segm_labels: return None - if relabel: - # Relabel the segmentation image with consecutive numbers; - # ndimage.label returns segment_img with dtype = np.int32 - # unless the input array has more than 2**31 - 1 pixels - n_labels = len(segm_labels) - if len(labels) != n_labels: - label_map = np.zeros(np.max(labels) + 1, - dtype=segment_img.dtype) - labels = np.arange(n_labels, dtype=segment_img.dtype) + 1 - label_map[segm_labels] = labels - segment_img = label_map[segment_img] - else: - # Use an ndarray so that seeded labels are always an array, - # matching the relabel path - labels = np.asarray(segm_labels) - - if return_segmimg: - return SegmentationImage._from_data(segment_img, labels=labels, - areas=np.array(segm_areas), - slices=segm_slices) - - # This is used by deblend_sources - if len(labels) == 1: - return None + # Relabel the segmentation image with consecutive numbers. + # ndimage.label returns segment_img with dtype = np.int32 + # unless the input array has more than 2**31 - 1 pixels + n_labels = len(segm_labels) + if len(labels) != n_labels: + label_map = np.zeros(np.max(labels) + 1, + dtype=segment_img.dtype) + labels = np.arange(n_labels, dtype=segment_img.dtype) + 1 + label_map[segm_labels] = labels + segment_img = label_map[segment_img] - return segment_img + return SegmentationImage._from_data(segment_img, labels=labels, + areas=np.array(segm_areas), + slices=segm_slices) @deprecated_renamed_argument('npixels', 'n_pixels', '3.0', until='4.0') @@ -390,7 +362,7 @@ class to detect and deblend sources in a single step. footprint = _make_binary_structure(data.ndim, connectivity) segm = _detect_sources(data, threshold, n_pixels, footprint, - inverse_mask, relabel=True, return_segmimg=True) + inverse_mask) if segm is None: msg = ('No sources were found. Try lowering the threshold or ' diff --git a/photutils/segmentation/finder.py b/photutils/segmentation/finder.py index 00fffc9156..9f7d21b218 100644 --- a/photutils/segmentation/finder.py +++ b/photutils/segmentation/finder.py @@ -83,8 +83,8 @@ class SourceFinder: the threshold levels between the source minimum and maximum. The ``'exponential'`` and ``'sinh'`` modes differ in that the ``'exponential'`` levels are dependent on the source - maximum/minimum ratio (smaller ratios are more linear; larger - ratios are more exponential), while the ``'sinh'`` levels + maximum/minimum ratio (smaller ratios are more linear and + larger ratios are more exponential), while the ``'sinh'`` levels are not. Also, the ``'exponential'`` mode will be changed to ``'linear'`` for sources with non-positive minimum data values. This keyword is ignored unless ``deblend=True``. @@ -95,25 +95,29 @@ class SourceFinder: consecutive order starting from 1. This keyword is ignored unless ``deblend=True``. + nproc : int, optional + This keyword is deprecated and has no effect. It was the name of + the ``n_processes`` keyword before version 3.0. + + .. deprecated:: 3.0 + The ``nproc`` keyword is deprecated and will be removed in + version 4.0. + n_processes : int, optional - The number of processes to use for source deblending. If set to - 1, then a serial implementation is used instead of a parallel - one. If `None`, then the number of processes will be set to the - number of CPUs detected on the machine. Please note that due to - overheads, multiprocessing may be slower than serial processing - if only a small number of sources are to be deblended. The - benefits of multiprocessing require ~1000 or more sources to - deblend, with larger gains as the number of sources increase. - This keyword is ignored unless ``deblend=True``. + This keyword is deprecated and has no effect. Multiprocessing + no longer provides any benefit for source deblending. + + .. deprecated:: 3.1 + The ``n_processes`` keyword is deprecated and will be + removed in version 4.0. progress_bar : bool, optional - Whether to display a progress bar. If ``n_processes = 1``, then the - ID shown after the progress bar is the source label being - deblended. If multiprocessing is used (``n_processes > 1``), the - ID shown is the last source label that was deblended. The - progress bar requires that the `tqdm `_ - optional dependency be installed. This keyword is ignored unless - ``deblend=True``. + This keyword is deprecated and has no effect. Deblending no + longer displays a progress bar. + + .. deprecated:: 3.1 + The ``progress_bar`` keyword is deprecated and will be + removed in version 4.0. See Also -------- @@ -148,7 +152,7 @@ class SourceFinder: # Detect the sources threshold = 1.5 * bkg.background_rms # per-pixel detection threshold - finder = SourceFinder(n_pixels=10, progress_bar=False) + finder = SourceFinder(n_pixels=10) segment_map = finder(convolved_data, threshold) # Plot the image and the segmentation image @@ -161,9 +165,12 @@ class SourceFinder: @deprecated_renamed_argument('npixels', 'n_pixels', '3.0', until='4.0') @deprecated_renamed_argument('nlevels', 'n_levels', '3.0', until='4.0') - @deprecated_renamed_argument('nproc', 'n_processes', '3.0', until='4.0') + @deprecated_renamed_argument('nproc', None, '3.0', until='4.0') + @deprecated_renamed_argument('n_processes', None, '3.1', until='4.0') + @deprecated_renamed_argument('progress_bar', None, '3.1', until='4.0') def __init__(self, n_pixels, *, connectivity=8, deblend=True, n_levels=32, contrast=0.001, mode='exponential', relabel=True, + nproc=1, # noqa: ARG002 n_processes=1, progress_bar=True): self.n_pixels = as_pair('n_pixels', n_pixels, check_odd=False) self.deblend = deblend @@ -225,15 +232,12 @@ def __call__(self, data, threshold, mask=None): if segment_img is None: return None - # Source deblending requires scikit-image if self.deblend: segment_img = deblend_sources(data, segment_img, self.n_pixels[1], n_levels=self.n_levels, contrast=self.contrast, mode=self.mode, connectivity=self.connectivity, - relabel=self.relabel, - n_processes=self.n_processes, - progress_bar=self.progress_bar) + relabel=self.relabel) return segment_img diff --git a/photutils/segmentation/flags.py b/photutils/segmentation/flags.py index 29e28b2097..37d3a9ac67 100644 --- a/photutils/segmentation/flags.py +++ b/photutils/segmentation/flags.py @@ -123,8 +123,8 @@ class _SegmentationFlags(FlagRegistry): FlagDefinition( bit_value=256, name='undefined_shape', - description=('non-positive net flux; shape properties ' - 'undefined'), + description=('non-positive net flux (shape properties ' + 'undefined)'), detailed_description=('The net source flux (the zeroth ' 'image moment over the source ' 'segment) is not positive, so the ' @@ -160,10 +160,10 @@ class _SegmentationFlags(FlagRegistry): 'the isophotal centroid'), detailed_description=('The windowed centroid fell ' 'outside the 1-sigma moment ' - 'ellipse; the windowed flux was ' - 'non-positive; the windowed ' + 'ellipse, the windowed flux was ' + 'non-positive, the windowed ' '2nd-order moments or covariance ' - 'determinant were negative; or ' + 'determinant were negative, or ' 'the iterated centroid was NaN. ' 'In each of these cases, the ' 'isophotal ``centroid`` value ' diff --git a/photutils/segmentation/tests/test_batch_centroid_win.py b/photutils/segmentation/tests/test_batch_centroid_win.py index 943d5d670c..84a6f26856 100644 --- a/photutils/segmentation/tests/test_batch_centroid_win.py +++ b/photutils/segmentation/tests/test_batch_centroid_win.py @@ -256,8 +256,8 @@ def test_matches_reference(scene, method, with_error, with_mask): ref = _reference_iterate_centroid_win( label, inp['xcen0'][i], inp['ycen0'][i], inp['radius_hl'][i], inp['nan_hl'][i], **ref_kwargs) - # atol covers cancellation-limited near-zero central moments; - # positions and flux are far from zero and are effectively + # atol covers cancellation-limited near-zero central moments. + # Positions and flux are far from zero and are effectively # checked at rtol assert_allclose(result[i], np.array(ref), rtol=1e-12, atol=1e-10, equal_nan=True) diff --git a/photutils/segmentation/tests/test_batch_circular.py b/photutils/segmentation/tests/test_batch_circular.py index 926867712c..a9b1d35f86 100644 --- a/photutils/segmentation/tests/test_batch_circular.py +++ b/photutils/segmentation/tests/test_batch_circular.py @@ -118,7 +118,7 @@ def test_local_background(scene): def test_all_masked_and_nonfinite_centroid(scene): - # A fully masked source has no centroid and a NaN flux; a source + # A fully masked source has no centroid and a NaN flux. A source # whose centroid is non-finite is also NaN data = scene['data'].copy() mask = scene['mask'].copy() diff --git a/photutils/segmentation/tests/test_batch_flux_radius.py b/photutils/segmentation/tests/test_batch_flux_radius.py index 9dfa33c2b8..5483541bee 100644 --- a/photutils/segmentation/tests/test_batch_flux_radius.py +++ b/photutils/segmentation/tests/test_batch_flux_radius.py @@ -315,7 +315,7 @@ def test_all_skipped(scene): def test_bracket_shrink_and_no_solution(scene): - # Force the shrink path: negative data beyond a ring makes the + # Force the shrink path. Negative data beyond a ring makes the # enclosed flux non-monotonic, so the initial bracket has equal # signs at both ends cat = make_catalog(scene) diff --git a/photutils/segmentation/tests/test_batch_kron.py b/photutils/segmentation/tests/test_batch_kron.py index ff7b2cec76..e702c1dcfb 100644 --- a/photutils/segmentation/tests/test_batch_kron.py +++ b/photutils/segmentation/tests/test_batch_kron.py @@ -204,7 +204,7 @@ def test_custom_kron_params(scene): def test_mixed_aperture_types(scene): # A minimum circular radius makes the sources whose scaled Kron # ellipse is smaller than that radius circular, while the larger - # sources stay elliptical; then both driver groups run in one + # sources stay elliptical. Both driver groups then run in one # _calc_kron_photometry call cat = SourceCatalog(scene['data'], scene['segm'], error=scene['error'], mask=scene['mask'], diff --git a/photutils/segmentation/tests/test_batch_moments.py b/photutils/segmentation/tests/test_batch_moments.py index 3dabc83257..f34a43af70 100644 --- a/photutils/segmentation/tests/test_batch_moments.py +++ b/photutils/segmentation/tests/test_batch_moments.py @@ -149,7 +149,7 @@ def test_central_moments(scene): def test_central_moments_nan_centroid(scene): - # A fully-masked source has zero total flux -> NaN centroid; the + # A fully-masked source has zero total flux -> NaN centroid. The # previous implementation then produced NaN everywhere except # [0, 0], which holds the (zero) flux sum cat = make_catalog(scene) @@ -204,8 +204,8 @@ def _edge_catalog(): convdata[4, 4] = np.nan # non-finite convolved value convdata[8, 8] = np.inf convdata[9, 9] = 0.0 # zero flux weight, excluded from the errors - # Non-finite data with a finite convolved value: included in the - # moments, excluded from the errors + # Non-finite data with a finite convolved value is included in + # the moments and excluded from the errors data[5, 5] = np.nan data[9, 7] = np.nan mask = np.zeros((ny, nx), dtype=bool) @@ -299,8 +299,8 @@ def test_catalog_moments(scene): def test_catalog_centroid_err(scene): - # _centroid_err_cov = normalized batch_moment_err accumulators; - # rebuild it from the reference accumulation of test_moment_err + # _centroid_err_cov = normalized batch_moment_err accumulators. + # Rebuild it from the reference accumulation of test_moment_err cat = make_catalog(scene) cov = cat._centroid_err_cov assert cov.shape == (cat.n_labels, 2, 2) diff --git a/photutils/segmentation/tests/test_batch_segment_gather.py b/photutils/segmentation/tests/test_batch_segment_gather.py index 01039d35b4..4203a2767e 100644 --- a/photutils/segmentation/tests/test_batch_segment_gather.py +++ b/photutils/segmentation/tests/test_batch_segment_gather.py @@ -81,7 +81,7 @@ def test_matches_reference(scene, with_error, with_mask, with_background): def test_all_masked_source(scene): # A completely masked source gathers a single NaN and is flagged as - # all masked; a source with a single unmasked pixel is not + # all masked. A source with a single unmasked pixel is not mask = scene['mask'].copy() segm = scene['segm'] slc = segm.slices[0] diff --git a/photutils/segmentation/tests/test_catalog.py b/photutils/segmentation/tests/test_catalog.py index 62bf498858..640d3e0f91 100644 --- a/photutils/segmentation/tests/test_catalog.py +++ b/photutils/segmentation/tests/test_catalog.py @@ -31,8 +31,7 @@ from photutils.segmentation.finder import SourceFinder from photutils.segmentation.flags import SEGMENTATION_FLAGS from photutils.segmentation.utils import make_2dgaussian_kernel -from photutils.utils._optional_deps import (HAS_GWCS, HAS_MATPLOTLIB, - HAS_SKIMAGE) +from photutils.utils._optional_deps import HAS_GWCS, HAS_MATPLOTLIB from photutils.utils._wcs_helpers import compute_pixel_to_sky_jacobians from photutils.utils.cutouts import CutoutImage @@ -96,7 +95,7 @@ def centroid_win_data(): kernel = make_2dgaussian_kernel(3.0, size=5) convolved_data = convolve(data, kernel) n_pixels = 10 - finder = SourceFinder(n_pixels=n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=n_pixels) threshold = 107.9 segment_map = finder(convolved_data, threshold) return data, segment_map, convolved_data @@ -1035,7 +1034,7 @@ def test_custom_properties(self, scalar): # Built-in cached property cat.add_property('area', segment_snr) with pytest.raises(ValueError, match=match): - # Built-in method; must raise even with overwrite=True + # Built-in method, which must raise even with overwrite=True cat.add_property('to_table', segment_snr, overwrite=True) cat.add_property('segment_snr', segment_snr) @@ -1215,7 +1214,6 @@ def test_make_kron_apertures(self): aper = obj.make_kron_apertures() assert isinstance(aper, EllipticalAperture) - @pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_make_cutouts(self): """ Test make cutouts. @@ -1229,7 +1227,7 @@ def test_make_cutouts(self): kernel = make_2dgaussian_kernel(3.0, size=5) convolved_data = convolve(data, kernel) n_pixels = 10 - finder = SourceFinder(n_pixels=n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=n_pixels) segment_map = finder(convolved_data, threshold) cat = SourceCatalog(data, segment_map, convolved_data=convolved_data) @@ -1607,9 +1605,8 @@ def _make_deblended_pair(): data = (Gaussian2D(100, 25, 30, 4, 4)(xx, yy) + Gaussian2D(100, 38, 30, 4, 4)(xx, yy)) segm = detect_sources(data, 10, 5) - return data, deblend_sources(data, segm, 5, progress_bar=False) + return data, deblend_sources(data, segm, 5) - @pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_kron_neighbor_pixels(self): """ Test the kron_neighbor_pixels flag when a neighbor segment falls @@ -1620,7 +1617,6 @@ def test_kron_neighbor_pixels(self): assert np.all(cat.flags & SEGMENTATION_FLAGS.KRON_NEIGHBOR_PIXELS) - @pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_kron_uncorrected_pixels(self): """ Test the kron_uncorrected_pixels flag for neighbor pixels within @@ -1823,7 +1819,6 @@ def test_kron_minimum_radius_two_element_params(self): assert np.all(cat.flags & SEGMENTATION_FLAGS.KRON_MINIMUM_RADIUS) - @pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_provenance_from_deblending(self): """ Test that deblending provenance flags propagate into the catalog @@ -1833,7 +1828,7 @@ def test_provenance_from_deblending(self): data = (Gaussian2D(100, 50, 50, 5, 5)(xx, yy) + Gaussian2D(100, 35, 50, 5, 5)(xx, yy)) segm = detect_sources(data, 10, 5) - segm2 = deblend_sources(data, segm, 5, progress_bar=False) + segm2 = deblend_sources(data, segm, 5) cat = SourceCatalog(data, segm2) assert np.all(cat.flags & SEGMENTATION_FLAGS.DEBLENDED) @@ -2029,7 +2024,6 @@ def worker(index): assert getattr(obj, f'prop{index}') == float(index) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_kron_params(): """ Test kron params. @@ -2046,7 +2040,7 @@ def test_kron_params(): convolved_data = convolve(data, kernel) n_pixels = 10 - finder = SourceFinder(n_pixels=n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=n_pixels) segm = finder(convolved_data, threshold) minrad = 1.4 @@ -2086,7 +2080,6 @@ def test_kron_params(): assert isinstance(cat.kron_aperture[0], CircularAperture) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_centroid_win(centroid_win_data): """ Test centroid win. @@ -2576,7 +2569,6 @@ def test_centroid_win_oom_guard(gauss_101_catalog): assert_allclose(cwin[:, 1], cat.y_centroid) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_centroid_win_aperture_mask_mask(centroid_win_data): """ Test centroid_win with aperture_mask_method='mask' to cover the @@ -2595,7 +2587,7 @@ def test_centroid_win_aperture_mask_mask(centroid_win_data): def test_make_scalar(single_source_catalog): """ - Test the scalar collapse of method results: a length-1 sequence + Test the scalar collapse of method results. A length-1 sequence is collapsed for a scalar catalog, a longer sequence is returned unchanged, and a multi-source catalog never collapses. """ @@ -2757,7 +2749,6 @@ def test_measured_kron_radius_circular_no_min_radius(gauss_101_data): assert np.all(np.isnan(kr)) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_centroid_win_err(): """ Test that centroid_win_err returns finite 1-sigma position @@ -2775,7 +2766,7 @@ def test_centroid_win_err(): kernel = make_2dgaussian_kernel(3.0, size=5) convolved_data = convolve(data, kernel) n_pixels = 10 - finder = SourceFinder(n_pixels=n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=n_pixels) threshold = 107.9 segment_map = finder(convolved_data, threshold) cat = SourceCatalog(data, segment_map, convolved_data=convolved_data, @@ -2783,7 +2774,7 @@ def test_centroid_win_err(): errors = cat.centroid_win_err assert errors.shape == (cat.n_labels, 2) - # Source 0 converged; errors should be finite and positive + # Source 0 converged, so errors should be finite and positive assert np.all(np.isfinite(errors[0])) assert np.all(errors[0] > 0) # Source 1 fell back to the isophotal centroid, so its errors @@ -2840,7 +2831,6 @@ def test_centroid_win_err_scalar(): assert single.y_centroid_win_err == errors[1] -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_centroid_win_err_sliced(): """ Test that centroid_win_err works on a sliced catalog after @@ -2858,7 +2848,7 @@ def test_centroid_win_err_sliced(): error = np.full(data.shape, 65.0) kernel = make_2dgaussian_kernel(3.0, size=5) convolved_data = convolve(data, kernel) - finder = SourceFinder(n_pixels=10, progress_bar=False) + finder = SourceFinder(n_pixels=10) segment_map = finder(convolved_data, 107.9) cat = SourceCatalog(data, segment_map, convolved_data=convolved_data, error=error, aperture_mask_method='none') @@ -2903,12 +2893,11 @@ def test_centroid_win_err_singularity(): assert np.all(errors > 0) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_centroid_win_err_cov(): """ - Test the windowed pixel error covariance: symmetric for every - source, near-zero off-diagonal for a circular source with uniform - errors, and consistent with centroid_win_err. + Test the windowed pixel error covariance. It is symmetric for + every source, near-zero off-diagonal for a circular source with + uniform errors, and consistent with centroid_win_err. """ yy, xx = np.mgrid[0:31, 0:31] data = Gaussian2D(500.0, 15.2, 15.6, 2.5, 2.5)(xx, yy) @@ -2928,7 +2917,6 @@ def test_centroid_win_err_cov(): np.sqrt((cov[0, 0, 0], cov[0, 1, 1]))) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_centroid_win_err_cov_fallback(): """ Test that fallback sources use the full isophotal covariance. @@ -2942,7 +2930,7 @@ def test_centroid_win_err_cov_fallback(): error = np.full(data.shape, 65.0) kernel = make_2dgaussian_kernel(3.0, size=5) convolved_data = convolve(data, kernel) - finder = SourceFinder(n_pixels=10, progress_bar=False) + finder = SourceFinder(n_pixels=10) segment_map = finder(convolved_data, 107.9) cat = SourceCatalog(data, segment_map, convolved_data=convolved_data, error=error, @@ -2955,7 +2943,6 @@ def test_centroid_win_err_cov_fallback(): cat._centroid_err_cov[1]) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_centroid_err(): """ Test that centroid_err returns finite 1-sigma position errors @@ -2971,7 +2958,7 @@ def test_centroid_err(): kernel = make_2dgaussian_kernel(3.0, size=5) convolved_data = convolve(data, kernel) n_pixels = 10 - finder = SourceFinder(n_pixels=n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=n_pixels) threshold = 107.9 segment_map = finder(convolved_data, threshold) cat = SourceCatalog(data, segment_map, convolved_data=convolved_data, @@ -3116,7 +3103,6 @@ def test_centroid_err_zero_flux(): assert np.all(np.isnan(errors)) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_centroid_err_columns(): """ Test the x/y centroid error properties and their use as to_table @@ -3131,7 +3117,7 @@ def test_centroid_err_columns(): error = np.full(data.shape, 65.0) kernel = make_2dgaussian_kernel(3.0, size=5) convolved_data = convolve(data, kernel) - finder = SourceFinder(n_pixels=10, progress_bar=False) + finder = SourceFinder(n_pixels=10) segment_map = finder(convolved_data, 107.9) cat = SourceCatalog(data, segment_map, convolved_data=convolved_data, error=error, aperture_mask_method='none') @@ -3308,7 +3294,6 @@ def test_centroid_quad_err_peak_at_edge(): assert np.all(np.isnan(cat.centroid_quad_err)) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_centroid_quad_err_columns(): """ Test the quadratic centroid error properties as to_table columns @@ -3321,7 +3306,7 @@ def test_centroid_quad_err_columns(): error = np.full(data.shape, 65.0) kernel = make_2dgaussian_kernel(3.0, size=5) convolved_data = convolve(data, kernel) - finder = SourceFinder(n_pixels=10, progress_bar=False) + finder = SourceFinder(n_pixels=10) segment_map = finder(convolved_data, 107.9) cat = SourceCatalog(data, segment_map, convolved_data=convolved_data, error=error, aperture_mask_method='none') diff --git a/photutils/segmentation/tests/test_core.py b/photutils/segmentation/tests/test_core.py index 9c45bfaac8..9a38bbb81c 100644 --- a/photutils/segmentation/tests/test_core.py +++ b/photutils/segmentation/tests/test_core.py @@ -603,7 +603,7 @@ def test_segment_no_full_array_reference(self): assert segment._segment_data_cutout.shape == (10, 10) assert segment._segment_data_shape == (1000, 1000) - # Delete the SegmentationImage; the segment should not keep + # Delete the SegmentationImage. The segment should not keep # the full array alive full_refcount = sys.getrefcount(large) del segm @@ -1240,8 +1240,8 @@ def test_polygons_complex(self): assert isinstance(region, (Regions, PolygonPixelRegion)) assert isinstance(regions[2], Regions) - # Combine all segments into a single segment; - # now have multipolygon objects, some with holes + # Combine all segments into a single segment. + # Now have multipolygon objects, some with holes segm.reassign_labels(segm.labels, new_label=4) polygons = segm.polygons assert len(polygons) == 1 @@ -1570,7 +1570,7 @@ class TestGetLabelMapping: def setup(self, segm_data): self.parent = SegmentationImage(segm_data) # A deblend-like version of segm_data with the identical - # non-zero footprint: label 5 is split into labels 8 and 9, + # non-zero footprint. Label 5 is split into labels 8 and 9, # and label 7 is split into labels 10 and 11 child_data = np.array([[1, 1, 0, 0, 4, 4], [0, 0, 0, 0, 0, 4], @@ -2231,7 +2231,7 @@ def test_concurrent_cached_properties_introspection(self): Test concurrent first access of the class-level _cached_properties introspection cache on a fresh subclass. - The lazy class-level cache is written without a lock; the race + The lazy class-level cache is written without a lock. The race is benign because the computed list is identical for every thread. """ diff --git a/photutils/segmentation/tests/test_deblend.py b/photutils/segmentation/tests/test_deblend.py index 45be941a9d..51049dad4b 100644 --- a/photutils/segmentation/tests/test_deblend.py +++ b/photutils/segmentation/tests/test_deblend.py @@ -3,24 +3,31 @@ Tests for the deblend module. """ +import warnings from unittest.mock import patch import numpy as np import pytest from astropy.modeling.models import Gaussian2D -from astropy.utils.exceptions import AstropyUserWarning +from astropy.utils.exceptions import (AstropyDeprecationWarning, + AstropyUserWarning) from numpy.testing import assert_allclose, assert_equal - -from photutils.segmentation import (SegmentationImage, deblend_sources, - detect_sources) -from photutils.segmentation.deblend import (_DeblendParams, - _SingleSourceDeblender) +from scipy import ndimage as ndi + +from photutils.segmentation import SegmentationImage +from photutils.segmentation import deblend as deblend_module +from photutils.segmentation import deblend_sources, detect_sources +from photutils.segmentation._deblend_reference import _SingleSourceDeblender +from photutils.segmentation._deblend_watershed import deblend_watershed +from photutils.segmentation.deblend import (_compute_thresholds, + _create_relabel_map, + _DeblendParams) from photutils.segmentation.flags import SEGMENTATION_FLAGS +from photutils.segmentation.utils import _make_binary_structure from photutils.utils._optional_deps import HAS_SKIMAGE from photutils.utils.exceptions import DeblendWarning -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') class TestDeblendSources: @pytest.fixture(autouse=True) def setup(self): @@ -43,19 +50,8 @@ def test_deblend_sources(self, mode): Test deblend sources. """ result = deblend_sources(self.data, self.segm, self.n_pixels, - mode=mode, progress_bar=False) + mode=mode) assert result.data.dtype == self.segm.data.dtype - - if mode == 'linear': - # Test multiprocessing - result2 = deblend_sources(self.data, self.segm, - self.n_pixels, - mode=mode, - progress_bar=False, - n_processes=2) - assert_equal(result.data, result2.data) - assert result2.data.dtype == self.segm.data.dtype - assert result.n_labels == 2 assert result.n_labels == len(result.slices) mask1 = (result.data == 1) @@ -77,7 +73,7 @@ def test_deblend_multiple_sources(self): y = self.y data = self.data + g4(x, y) + g5(x, y) + g6(x, y) + g7(x, y) segm = detect_sources(data, self.threshold, self.n_pixels) - result = deblend_sources(data, segm, self.n_pixels, progress_bar=False) + result = deblend_sources(data, segm, self.n_pixels) assert result.n_labels == 6 assert result.n_labels == len(result.slices) assert result.areas[0] == result.areas[1] @@ -98,7 +94,7 @@ def test_deblend_multiple_sources_with_neighbor(self): y = self.y data = (g1 + g2 + g3)(x, y) segm = detect_sources(data, self.threshold, self.n_pixels) - result = deblend_sources(data, segm, self.n_pixels, progress_bar=False) + result = deblend_sources(data, segm, self.n_pixels) assert result.n_labels == 3 def test_deblend_labels(self): @@ -112,8 +108,7 @@ def test_deblend_labels(self): y = self.y data = (g1 + g2 + g3)(x, y) segm = detect_sources(data, self.threshold, self.n_pixels) - result = deblend_sources(data, segm, self.n_pixels, labels=1, - progress_bar=False) + result = deblend_sources(data, segm, self.n_pixels, labels=1) assert result.n_labels == 2 @pytest.mark.parametrize(('contrast', 'n_labels'), @@ -135,8 +130,7 @@ def test_deblend_contrast(self, contrast, n_labels): n_pixels = 5 segm = detect_sources(data, 1.0, n_pixels) segm2 = deblend_sources(data, segm, n_pixels, mode='linear', - n_levels=32, contrast=contrast, - progress_bar=False) + n_levels=32, contrast=contrast) assert segm2.n_labels == n_labels def test_deblend_contrast_levels(self): @@ -158,8 +152,7 @@ def test_deblend_contrast_levels(self): segm = detect_sources(data, 1.0, n_pixels) for contrast in np.arange(1, 11) / 10.0: segm3 = deblend_sources(data, segm, n_pixels, mode='linear', - n_levels=32, contrast=contrast, - progress_bar=False) + n_levels=32, contrast=contrast) assert segm3.n_labels >= 1 def test_deblend_connectivity(self): @@ -179,20 +172,17 @@ def test_deblend_connectivity(self): segm = detect_sources(data, 0.1, 1, connectivity=4) assert segm.n_labels == 9 - segm2 = deblend_sources(data, segm, 1, mode='linear', connectivity=4, - progress_bar=False) + segm2 = deblend_sources(data, segm, 1, mode='linear', connectivity=4) assert segm2.n_labels == 9 segm = detect_sources(data, 0.1, 1, connectivity=8) assert segm.n_labels == 1 - segm2 = deblend_sources(data, segm, 1, mode='linear', connectivity=8, - progress_bar=False) + segm2 = deblend_sources(data, segm, 1, mode='linear', connectivity=8) assert segm2.n_labels == 3 match = 'Deblending failed for source' with pytest.raises(ValueError, match=match): - deblend_sources(data, segm, 1, mode='linear', connectivity=4, - progress_bar=False) + deblend_sources(data, segm, 1, mode='linear', connectivity=4) def test_deblend_label_assignment(self): """ @@ -216,7 +206,7 @@ def test_deblend_label_assignment(self): n_pixels = 5 segm1 = detect_sources(data, 5.0, n_pixels) segm2 = deblend_sources(data, segm1, n_pixels, mode='linear', - n_levels=32, contrast=0.3, progress_bar=False) + n_levels=32, contrast=0.3) assert segm2.n_labels == 4 @pytest.mark.parametrize('mode', ['exponential', 'linear']) @@ -225,7 +215,7 @@ def test_deblend_sources_norelabel(self, mode): Test deblend sources norelabel. """ result = deblend_sources(self.data, self.segm, self.n_pixels, - mode=mode, relabel=False, progress_bar=False) + mode=mode, relabel=False) assert result.n_labels == 2 assert_equal(result.labels, [2, 3]) assert_equal(result.parent_to_deblended_labels, {1: [2, 3]}) @@ -239,7 +229,7 @@ def test_deblend_three_sources(self, mode): Test deblend three sources. """ result = deblend_sources(self.data3, self.segm3, self.n_pixels, - mode=mode, progress_bar=False) + mode=mode) assert result.n_labels == 3 assert_allclose(np.nonzero(self.segm3), np.nonzero(result)) @@ -250,14 +240,12 @@ def test_segmentation_image(self): segm_wrong = np.ones((2, 2), dtype=int) # ndarray match = 'segmentation_image must be a SegmentationImage' with pytest.raises(TypeError, match=match): - deblend_sources(self.data, segm_wrong, self.n_pixels, - progress_bar=False) + deblend_sources(self.data, segm_wrong, self.n_pixels) segm_wrong = SegmentationImage(segm_wrong) # wrong shape match = 'segmentation_image must have the same shape as data' with pytest.raises(ValueError, match=match): - deblend_sources(self.data, segm_wrong, self.n_pixels, - progress_bar=False) + deblend_sources(self.data, segm_wrong, self.n_pixels) @pytest.mark.parametrize('relabel', [False, True]) def test_contrast_one_relabel(self, relabel): @@ -268,8 +256,7 @@ def test_contrast_one_relabel(self, relabel): segm = self.segm.copy() segm.reassign_label(1, 1000) result = deblend_sources(self.data, segm, self.n_pixels, - contrast=1, relabel=relabel, - progress_bar=False) + contrast=1, relabel=relabel) expected = [1] if relabel else [1000] assert_equal(result.labels, expected) @@ -281,8 +268,7 @@ def test_empty_segmentation_image(self): segm = SegmentationImage(np.zeros(self.data.shape, dtype=int)) match = 'segmentation_image must have at least one non-zero label' with pytest.raises(ValueError, match=match): - deblend_sources(self.data, segm, self.n_pixels, - progress_bar=False) + deblend_sources(self.data, segm, self.n_pixels) @pytest.mark.parametrize('n_pixels', [0, -5, 2.5]) def test_invalid_n_pixels(self, n_pixels): @@ -291,17 +277,17 @@ def test_invalid_n_pixels(self, n_pixels): """ match = 'n_pixels must be a positive integer' with pytest.raises(ValueError, match=match): - deblend_sources(self.data, self.segm, n_pixels, - progress_bar=False) + deblend_sources(self.data, self.segm, n_pixels) - def test_invalid_n_levels(self): + @pytest.mark.parametrize('n_levels', [0, -3, 2.7]) + def test_invalid_n_levels(self, n_levels): """ - Test invalid n_levels. + Test that invalid n_levels values raise a ValueError. """ - match = 'n_levels must be >= 1' + match = 'n_levels must be a positive integer' with pytest.raises(ValueError, match=match): - deblend_sources(self.data, self.segm, self.n_pixels, n_levels=0, - progress_bar=False) + deblend_sources(self.data, self.segm, self.n_pixels, + n_levels=n_levels) def test_invalid_contrast(self): """ @@ -309,8 +295,7 @@ def test_invalid_contrast(self): """ match = 'contrast must be >= 0 and <= 1' with pytest.raises(ValueError, match=match): - deblend_sources(self.data, self.segm, self.n_pixels, contrast=-1, - progress_bar=False) + deblend_sources(self.data, self.segm, self.n_pixels, contrast=-1) def test_invalid_mode(self): """ @@ -319,7 +304,7 @@ def test_invalid_mode(self): match = "mode must be 'exponential', 'linear', or 'sinh'" with pytest.raises(ValueError, match=match): deblend_sources(self.data, self.segm, self.n_pixels, - mode='invalid', progress_bar=False) + mode='invalid') def test_invalid_connectivity(self): """ @@ -328,7 +313,7 @@ def test_invalid_connectivity(self): match = 'Invalid connectivity' with pytest.raises(ValueError, match=match): deblend_sources(self.data, self.segm, self.n_pixels, - connectivity='invalid', progress_bar=False) + connectivity='invalid') def test_constant_source(self): """ @@ -336,8 +321,7 @@ def test_constant_source(self): """ data = self.data.copy() data[data.nonzero()] = 1.0 - result = deblend_sources(data, self.segm, self.n_pixels, - progress_bar=False) + result = deblend_sources(data, self.segm, self.n_pixels) assert_allclose(result, self.segm) def test_source_with_negval(self): @@ -348,8 +332,7 @@ def test_source_with_negval(self): data -= 20 match = 'The deblending mode of one or more source labels from the' with pytest.warns(DeblendWarning, match=match): - segm = deblend_sources(data, self.segm, self.n_pixels, - progress_bar=False) + segm = deblend_sources(data, self.segm, self.n_pixels) assert list(segm.info) == ['nonposmin_labels'] assert_equal(segm.info['nonposmin_labels'], [1]) @@ -358,8 +341,7 @@ def test_flags_deblended(self): Test that deblended children carry the deblended flag and nothing else when no mode fallback occurred. """ - result = deblend_sources(self.data, self.segm, self.n_pixels, - progress_bar=False) + result = deblend_sources(self.data, self.segm, self.n_pixels) assert_equal(result.flags, np.full(result.n_labels, SEGMENTATION_FLAGS.DEBLENDED)) @@ -374,8 +356,7 @@ def test_flags_nonposmin_children(self): data -= 20 match = 'The deblending mode of one or more source labels' with pytest.warns(DeblendWarning, match=match): - segm = deblend_sources(data, self.segm, self.n_pixels, - progress_bar=False) + segm = deblend_sources(data, self.segm, self.n_pixels) expected = (SEGMENTATION_FLAGS.DEBLENDED | SEGMENTATION_FLAGS.DEBLEND_NONPOSMIN) for label in segm.parent_to_deblended_labels[1]: @@ -397,8 +378,7 @@ def test_source_zero_min(self): data -= data[self.segm.data > 0].min() match = 'The deblending mode of one or more source labels from the' with pytest.warns(DeblendWarning, match=match): - segm = deblend_sources(data, self.segm, self.n_pixels, - progress_bar=False) + segm = deblend_sources(data, self.segm, self.n_pixels) assert_equal(segm.info['nonposmin_labels'], [1]) def test_connectivity(self): @@ -415,13 +395,11 @@ def test_connectivity(self): segm[data.nonzero()] = 1 segm = SegmentationImage(segm) data = data * 100.0 - segm_deblend = deblend_sources(data, segm, n_pixels=1, connectivity=8, - progress_bar=False) + segm_deblend = deblend_sources(data, segm, n_pixels=1, connectivity=8) assert segm_deblend.n_labels == 1 match = 'Deblending failed for source' with pytest.raises(ValueError, match=match): - deblend_sources(data, segm, n_pixels=1, connectivity=4, - progress_bar=False) + deblend_sources(data, segm, n_pixels=1, connectivity=4) def test_data_nan(self): """ @@ -432,8 +410,10 @@ def test_data_nan(self): """ data = self.data.copy() data[50, 50] = np.nan - segm2 = deblend_sources(data, self.segm, 5, progress_bar=False) + segm2 = deblend_sources(data, self.segm, 5) assert segm2.n_labels == 2 + # The NaN pixel is assigned to the source surrounding it + assert segm2.data[50, 50] == segm2.data[50, 51] != 0 def test_watershed(self): """ @@ -445,8 +425,7 @@ def test_watershed(self): """ segm = self.segm.copy() segm.reassign_label(1, 512) - result = deblend_sources(self.data, segm, self.n_pixels, - progress_bar=False) + result = deblend_sources(self.data, segm, self.n_pixels) assert result.n_labels == 2 def test_nondetection(self): @@ -461,7 +440,7 @@ def test_nondetection(self): data[50, 50] = 1000.0 data[50, 70] = 500.0 self.segm = detect_sources(data, self.threshold, self.n_pixels) - deblend_sources(data, self.segm, self.n_pixels, progress_bar=False) + deblend_sources(data, self.segm, self.n_pixels) def test_nonconsecutive_labels(self): """ @@ -469,8 +448,7 @@ def test_nonconsecutive_labels(self): """ segm = self.segm.copy() segm.reassign_label(1, 1000) - result = deblend_sources(self.data, segm, self.n_pixels, - progress_bar=False) + result = deblend_sources(self.data, segm, self.n_pixels) assert result.n_labels == 2 def test_single_source_methods(self): @@ -492,7 +470,7 @@ def test_single_source_methods(self): segms = single_debl.multithreshold() assert len(segms) == 32 - markers = single_debl.make_markers(return_all=True) + markers = single_debl.make_markers_per_level() assert len(markers) == 19 def test_info_empty_without_warnings(self): @@ -501,32 +479,465 @@ def test_info_empty_without_warnings(self): attribute, which is an empty dict when no deblending warnings occurred. """ - result = deblend_sources(self.data, self.segm, self.n_pixels, - progress_bar=False) + result = deblend_sources(self.data, self.segm, self.n_pixels) assert result.info == {} # detect_sources output must also have an info attribute assert self.segm.info == {} - def test_deblend_progress_bar(self): + @pytest.mark.parametrize('kwargs', [{'progress_bar': False}, + {'n_processes': 2}, {'nproc': 2}]) + def test_deprecated_keywords(self, kwargs): """ - Test deblend_sources with progress_bar=True (serial). + Test that the progress_bar, n_processes, and nproc keywords + each emit a single deprecation warning and have no effect on + the results. """ - result = deblend_sources(self.data, self.segm, self.n_pixels, - mode='linear', progress_bar=True) + name = next(iter(kwargs)) + with pytest.warns(AstropyDeprecationWarning, match=name) as record: + result = deblend_sources(self.data, self.segm, + self.n_pixels, mode='linear', + **kwargs) + assert len(record) == 1 assert result.n_labels == 2 - def test_deblend_nproc_none(self): - """ - Test deblend_sources with n_processes=None (auto-detect CPU count). - """ - result = deblend_sources(self.data, self.segm, self.n_pixels, - mode='linear', progress_bar=False, - n_processes=None) - assert result.n_labels == 2 + +def make_marker_test_image(kind): + """ + Return an image containing a single connected source for the + marker-path equivalence tests. + + Parameters + ---------- + kind : {'blend', 'quantized', 'plateau'} + The image type. 'blend' is a Gaussian envelope with compact + peaks spanning a wide amplitude range, 'quantized' is the + same image coarsely quantized (duplicate data values produce + empty multithreshold levels), and 'plateau' contains flat + stepped square annuli with two embedded peaks. + + Returns + ------- + data : 2D `~numpy.ndarray` + The image. + """ + if kind == 'plateau': + data = np.zeros((101, 101)) + data[10:90, 10:90] = 1.0 + data[30:70, 30:70] = 2.0 + data[35:45, 35:45] = 5.0 + data[55:65, 55:65] = 4.0 + return data + + rng = np.random.default_rng(0) + y, x = np.mgrid[0:151, 0:151] + data = Gaussian2D(20, 75, 75, 40, 40)(x, y) + amplitudes = np.geomspace(3.0, 100.0, 12) + radii = 40 * np.sqrt(rng.uniform(0.0, 1.0, 12)) + angles = rng.uniform(0.0, 2.0 * np.pi, 12) + for amplitude, radius, angle in zip(amplitudes, radii, angles, + strict=True): + xc = 75 + radius * np.cos(angle) + yc = 75 + radius * np.sin(angle) + data += Gaussian2D(amplitude, xc, yc, 2.5, 2.5)(x, y) + if kind == 'quantized': + data = np.round(data * 2.0) / 2.0 + return data + + +def normalize_markers(markers): + """ + Relabel a marker image with consecutive raster-ordered labels. + + Parameters + ---------- + markers : 2D int `~numpy.ndarray` + The marker image. + + Returns + ------- + result : 2D int `~numpy.ndarray` + The relabeled marker image. + """ + relabel_map = _create_relabel_map(markers) + if relabel_map is not None: + markers = relabel_map[markers] + return markers + + +@pytest.mark.parametrize('kind', ['blend', 'quantized', 'plateau']) +@pytest.mark.parametrize('mode', ['exponential', 'linear', 'sinh']) +@pytest.mark.parametrize('connectivity', [8, 4]) +def test_make_markers_matches_legacy(kind, mode, connectivity): + """ + Test that make_markers produces the same markers as the legacy + per-level path (the last image of the make_markers_per_level + chain). + + The markers must contain the same regions with the same + raster-scan label ordering, since the ordering determines the + final deblended label assignment. + """ + data = make_marker_test_image(kind) + segm = detect_sources(data, 0.5, 5, connectivity=connectivity) + footprint = _make_binary_structure(2, connectivity) + n_seen = 0 + for label, slc in zip(segm.labels, segm.slices, strict=True): + params = _DeblendParams(5, footprint, 32, 0.001, mode) + deblender = _SingleSourceDeblender(data[slc], segm.data[slc], + label, params) + markers = deblender.make_markers() + + params = _DeblendParams(5, footprint, 32, 0.001, mode) + deblender = _SingleSourceDeblender(data[slc], segm.data[slc], + label, params) + legacy = deblender.make_markers_per_level() + legacy = None if legacy is None else legacy[-1] + + if markers is None or legacy is None: + assert markers is None + assert legacy is None + else: + n_seen += 1 + assert_equal(normalize_markers(markers), + normalize_markers(legacy)) + assert n_seen >= 1 @pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') +@pytest.mark.parametrize('connectivity', [8, 4]) +def test_watershed_matches_skimage(connectivity): + """ + Test that the deblending watershed kernel produces results identical + to skimage.segmentation.watershed over randomized images, including + integer-valued and constant images whose plateaus exercise the + queue-age tie-breaking. + """ + from skimage.segmentation import watershed + + footprint = _make_binary_structure(2, connectivity) + rng = np.random.default_rng(987) + n_run = 0 + for trial in range(150): + ny, nx = rng.integers(5, 35, 2) + kind = trial % 4 + if kind == 0: + image = rng.normal(0.0, 1.0, (ny, nx)) + elif kind == 1: + image = rng.integers(0, 4, (ny, nx)).astype(float) + elif kind == 2: + image = np.zeros((ny, nx)) + else: + image = np.round(rng.normal(0.0, 1.0, (ny, nx)), 1) + mask = rng.random((ny, nx)) < 0.8 + indices = np.flatnonzero(mask) + if indices.size < 4: + continue + seeds = np.zeros((ny, nx), dtype=bool) + pick = rng.choice(indices, size=min(6, indices.size), + replace=False) + seeds.ravel()[pick] = True + markers = ndi.label(seeds, + structure=footprint)[0].astype(np.int32) + expected = watershed(image, markers, mask=mask, + connectivity=footprint) + result = deblend_watershed(image, markers, mask, connectivity) + assert_equal(result, expected) + n_run += 1 + assert n_run > 100 + + +def python_deblend_chunk(data, segm_data, driver_data, # noqa: ARG001 + driver_segm, # noqa: ARG001 + labels, slices, deblend_params): + """ + Deblend a chunk of sources with the pure-Python reference path. + + This mirrors the compiled chunk driver used by deblend_sources, + computing the markers and fallbacks per source in Python. + + Parameters + ---------- + data, segm_data, driver_data, driver_segm : 2D `~numpy.ndarray` + The (driver-compatible) data and segmentation arrays. + + labels : 1D `~numpy.ndarray` + The labels of the sources in the chunk. + + slices : list of tuple of slice + The bounding-box slices of the sources in the chunk. + + deblend_params : `_DeblendParams` + The parameters for deblending the sources. + + Returns + ------- + results : list of (2D `~numpy.ndarray` or `None`, dict) + The deblended cutout and warnings for each source. + """ + results = [] + for label, slc in zip(labels, slices, strict=True): + deblender = _SingleSourceDeblender(data[slc], segm_data[slc], + label, deblend_params) + results.append((deblender.deblend_source(), + deblender.warnings)) + return results + + +@pytest.mark.parametrize('dtype', ['float64', 'float32', '>f4', 'int32']) +@pytest.mark.parametrize('scene', ['blend', 'negmin', 'checkerboard', + 'gaussian', 'flat', + 'contrast-batch', 'contrast-single', + 'contrast-all', 'contrast-negmin', + 'neighbors', 'nan']) +def test_chunk_driver_matches_python_path(dtype, scene): + """ + Test that the compiled chunk driver and contrast loop produce + results identical to the pure-Python per-source path, including + the threshold computation, the mode fallbacks, the below-contrast + marker removal, and the recorded warnings, for float64, float32, and + integer data. + + The scenes cover deblending sources, the non-positive-minimum and + too-many-markers mode fallbacks, a source that does not split, + a constant source, and contrast values that trigger the batched + removal, the one-at-a-time removal, the removal of all but one + basin, and the removal path for sources with a negative minimum + (which always removes one marker at a time), a scene of several + segments with overlapping bounding boxes, and NaN pixels within a + segment. + """ + contrast = 0.001 + if scene == 'blend': + data = make_marker_test_image('blend') + threshold, n_pixels = 0.5, 5 + elif scene == 'negmin': + data = make_marker_test_image('blend') - 15.0 + threshold, n_pixels = -14.5, 5 + elif scene.startswith('contrast'): + data, _ = make_multipeak_source() + threshold, n_pixels = 0.5, 5 + contrast = {'contrast-batch': 0.15, 'contrast-single': 0.07, + 'contrast-all': 0.35, + 'contrast-negmin': 0.15}[scene] + if scene == 'contrast-negmin': + data = data - 5.0 + threshold = -4.5 + elif scene == 'gaussian': + y, x = np.mgrid[0:51, 0:51] + data = Gaussian2D(10, 25, 25, 5, 5)(x, y) + threshold, n_pixels = 0.5, 5 + elif scene == 'flat': + data = np.zeros((51, 51)) + data[20:40, 20:40] = 5.0 + threshold, n_pixels = 0.5, 5 + elif scene == 'neighbors': + # Blended pairs whose bounding boxes contain other segments + y, x = np.mgrid[0:121, 0:161] + data = (Gaussian2D(100, 50, 60, 6, 6)(x, y) + + Gaussian2D(90, 68, 60, 6, 6)(x, y) + + Gaussian2D(60, 82, 40, 3, 3)(x, y) + + Gaussian2D(80, 120, 70, 5, 5)(x, y) + + Gaussian2D(70, 135, 70, 5, 5)(x, y) + + Gaussian2D(50, 118, 95, 3, 3)(x, y)) + threshold, n_pixels = 5.0, 5 + contrast = 0.01 + elif scene == 'nan': + if not np.issubdtype(np.dtype(dtype), np.floating): + pytest.skip('NaN requires a floating-point dtype') + data = make_marker_test_image('blend') + threshold, n_pixels = 0.5, 5 + else: + # The n_markers fallback checkerboard scene + size = 51 + data1 = np.resize([0, 0, 1, 1], size) + data1 = np.abs(data1 - np.atleast_2d(data1).T) + 2.0 + for i in range(size): + if i % 2 == 0: + data1[i, :] = 1 + data1[:, i] = 1 + data = np.zeros((101, 101)) + data[25:25 + size, 25:25 + size] = data1 + data[50:60, 50:60] = 10.0 + threshold, n_pixels = 0.01, 1 + + data = data.astype(dtype) + segm = detect_sources(data, threshold, n_pixels) + if scene == 'nan': + # NaN pixels within the segment, set after the detection + rng = np.random.default_rng(11) + ys, xs = np.nonzero(segm.data > 0) + pick = rng.choice(ys.size, size=40, replace=False) + data[ys[pick], xs[pick]] = np.nan + + with warnings.catch_warnings(): + warnings.simplefilter('ignore', DeblendWarning) + result = deblend_sources(data, segm, n_pixels, + contrast=contrast) + with patch.object(deblend_module, '_deblend_sources_chunk', + python_deblend_chunk): + expected = deblend_sources(data, segm, n_pixels, + contrast=contrast) + + assert_equal(result.data, expected.data) + assert result.info.keys() == expected.info.keys() + for key in expected.info: + assert_equal(result.info[key], expected.info[key]) + assert result._flags_map == expected._flags_map + + +def test_deblend_segm_dtype(): + """ + Test that deblending a segmentation image with a non-native + integer dtype gives the same result as the int32 one. + """ + data, segm = make_multipeak_source() + expected = deblend_sources(data, segm, 5) + segm16 = SegmentationImage(segm.data.astype(np.int16)) + result = deblend_sources(data, segm16, 5) + assert_equal(result.data, expected.data) + + +def test_deblend_byte_order(): + """ + Test that non-native byte order data and segmentation images give + results identical to the native ones. + """ + data, segm = make_multipeak_source() + for dtype in ('f4', 'f8'): + expected = deblend_sources(data.astype(f'<{dtype}'), segm, 5, + contrast=0.01) + result = deblend_sources(data.astype(f'>{dtype}'), segm, 5, + contrast=0.01) + assert_equal(result.data, expected.data) + + expected = deblend_sources(data, segm, 5, contrast=0.01) + segm_be = SegmentationImage(segm.data.astype('>i4')) + result = deblend_sources(data, segm_be, 5, contrast=0.01) + assert_equal(result.data, expected.data) + + +@pytest.mark.parametrize('dtype', ['float64', 'float32', 'int32']) +@pytest.mark.parametrize('mode', ['exponential', 'linear', 'sinh']) +@pytest.mark.parametrize('n_levels', [1, 7, 32]) +def test_compute_thresholds_matches_reference(dtype, mode, n_levels): + """ + Test that the vectorized multithreshold levels are bitwise identical + to the per-source levels of the reference implementation, including + the exponential-mode fallback for non-positive minima and the + zero-step special case of np.linspace. + """ + rng = np.random.default_rng(42) + n_src = 60 + if dtype == 'int32': + smin = rng.integers(-50, 50, n_src) + smax = smin + rng.integers(1, 1000, n_src) + else: + smin = rng.uniform(-5.0, 5.0, n_src) + smax = smin + 10.0 ** rng.uniform(-3.0, 3.0, n_src) + smin[1] = 0.0 + smax[1] = np.finfo(dtype).smallest_subnormal + smin[0] = 0 + smin = smin.astype(dtype) + smax = smax.astype(dtype) + assert np.all(smax > smin) + + thresholds, nonposmin = _compute_thresholds(smin, smax, n_levels, + mode) + assert thresholds.shape == (n_src, n_levels) + assert thresholds.dtype == np.float64 + assert thresholds.flags.c_contiguous + + params = _DeblendParams(1, np.ones((3, 3)), n_levels, 0.001, mode) + segment = np.ones((1, 2), dtype=int) + for i in range(n_src): + data = np.array([[smin[i], smax[i]]], dtype=dtype) + deblender = _SingleSourceDeblender(data, segment, 1, params) + expected = np.asarray(deblender.compute_thresholds(), + dtype=np.float64) + # Compare the bit patterns, not just the values + assert_equal(thresholds[i].view(np.int64), + expected.view(np.int64)) + assert nonposmin[i] == ('nonposmin' in deblender.warnings) + + +def test_python_path_connectivity_mismatch(): + """ + Test that the pure-Python reference path raises the same error + as the compiled contrast loop when the detection and deblending + connectivities differ. + """ + data = np.zeros((51, 51)) + data[15:36, 15:36] = 10.0 + data[14, 36] = 1.0 + data[13, 37] = 10 + data[14, 14] = 5.0 + data[13, 13] = 10.0 + data[36, 14] = 10.0 + data[37, 13] = 10.0 + data[36, 36] = 10.0 + data[37, 37] = 10.0 + segm = detect_sources(data, 0.1, 1, connectivity=8) + match = 'Deblending failed for source' + with (patch.object(deblend_module, '_deblend_sources_chunk', + python_deblend_chunk), + pytest.raises(ValueError, match=match)): + deblend_sources(data, segm, 1, mode='linear', connectivity=4) + + +def make_multipeak_source(): + """ + Return the image and segmentation image for a single connected + source with peaks spanning a wide range of basin fluxes. + + The watershed basin flux fractions are approximately 0.061, 0.062, + 0.225, and 0.652, so the contrast keyword controls how many of the + faintest basins fail the contrast criterion. + + Returns + ------- + data : 2D `~numpy.ndarray` + The image. + + segm : `~photutils.segmentation.SegmentationImage` + The segmentation image containing a single label. + """ + y, x = np.mgrid[0:101, 0:101] + envelope = Gaussian2D(1.0, 50, 50, 30, 30) + g_bright = Gaussian2D(100, 40, 50, 3, 3) + g_medium = Gaussian2D(30, 68, 50, 3, 3) + g_faint1 = Gaussian2D(3.0, 50, 30, 3, 3) + g_faint2 = Gaussian2D(4.0, 50, 72, 3, 3) + data = (envelope(x, y) + g_bright(x, y) + g_medium(x, y) + + g_faint1(x, y) + g_faint2(x, y)) + segm = detect_sources(data, 0.5, 5) + return data, segm + + +@pytest.mark.parametrize(('contrast', 'n_labels'), + [(0.0, 4), (0.07, 2), (0.15, 2), (0.35, 1)]) +def test_contrast_removal(contrast, n_labels): + """ + Test the below-contrast marker removal over a range of contrasts. + + The contrast=0.15 case removes the two faintest basins together + (their total flux fraction is below both the contrast and the + next-faintest basin flux), the contrast=0.07 case removes the same + two basins one at a time (their total is above the contrast), + and the contrast=0.35 case removes all but one basin so that no + deblending occurs. + """ + data, segm = make_multipeak_source() + result = deblend_sources(data, segm, 5, contrast=contrast) + assert result.n_labels == n_labels + assert_equal(np.nonzero(segm.data), np.nonzero(result.data)) + if n_labels > 1: + assert_equal(result.parent_to_deblended_labels, + {1: list(range(1, n_labels + 1))}) + else: + assert_equal(result.parent_to_deblended_labels, {}) + + def test_n_markers_fallback(): """ Test that if there are too many markers, a warning is raised. @@ -551,7 +962,6 @@ def test_n_markers_fallback(): assert segm2.info['n_markers_labels'][0] == 1 -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_flags_n_markers_fallback(): """ Test that the n_markers fallback flag is set on the output @@ -586,7 +996,6 @@ def test_flags_n_markers_fallback(): assert_equal(flagged, [1]) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') @pytest.mark.parametrize('relabel', [True, False]) def test_flags_fallback_without_deblending(relabel): """ @@ -611,55 +1020,23 @@ def test_flags_fallback_without_deblending(relabel): match = 'The deblending mode of one or more source labels' with pytest.warns(DeblendWarning, match=match): - segm2 = deblend_sources(data, segm, 5, progress_bar=False, + segm2 = deblend_sources(data, segm, 5, relabel=relabel) bit = SEGMENTATION_FLAGS.DEBLEND_NONPOSMIN - # Both remaining sources fell back (both have negative minima); - # neither splits, so each output label carries the fallback bit + # Both remaining sources fell back (both have negative minima). + # Neither splits, so each output label carries the fallback bit # but not the deblended bit assert_equal(segm2.flags & bit, [bit] * segm2.n_labels) assert_equal(segm2.flags & SEGMENTATION_FLAGS.DEBLENDED, [0] * segm2.n_labels) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') -def test_n_markers_fallback_multiproc(): - """ - Test the n_markers fallback warning via multiprocessing - (n_processes=2). This covers the multiprocessing result-processing - block for n_markers. - """ - size = 51 - data1 = np.resize([0, 0, 1, 1], size) - data1 = np.abs(data1 - np.atleast_2d(data1).T) + 2 - - for i in range(size): - if i % 2 == 0: - data1[i, :] = 1 - data1[:, i] = 1 - - data = np.zeros((101, 101)) - data[25:25 + size, 25:25 + size] = data1 - data[50:60, 50:60] = 10.0 - - segm = detect_sources(data, 0.01, 10) - match = 'The deblending mode of one or more source labels from the' - with pytest.warns(DeblendWarning, match=match): - segm2 = deblend_sources(data, segm, 1, mode='exponential', - n_processes=2) - assert segm2.info['n_markers_labels'][0] == 1 - - -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') -def test_nonposmin_multiproc(): +def test_nonposmin_astropy_user_warning(): """ - Test nonposmin warning via multiprocessing (n_processes=2). - - This covers the multiprocessing result-processing block for - nonposmin. The warning is caught as an AstropyUserWarning to check - that DeblendWarning is a subclass of it, so existing warning filters - continue to work. + Test that the nonposmin warning is caught as an + AstropyUserWarning, checking that DeblendWarning is a subclass of + it so that existing warning filters continue to work. """ g1 = Gaussian2D(100, 50, 50, 8, 8) g2 = Gaussian2D(100, 35, 50, 8, 8) @@ -669,13 +1046,11 @@ def test_nonposmin_multiproc(): segm = detect_sources(data + 20, 10, 5) # detect sources on positive data match = 'The deblending mode of one or more source labels from the' with pytest.warns(AstropyUserWarning, match=match): - segm2 = deblend_sources(data, segm, 5, progress_bar=False, - n_processes=2) + segm2 = deblend_sources(data, segm, 5) assert 'nonposmin_labels' in segm2.info assert np.all(segm2.flags & SEGMENTATION_FLAGS.DEBLEND_NONPOSMIN) -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_n_markers_fallback_returns_none(): """ Test that deblend_source returns None when make_markers returns @@ -696,7 +1071,7 @@ def test_n_markers_fallback_returns_none(): call_count = [0] - def mock_make_markers(*, _return_all=False): + def mock_make_markers(): call_count[0] += 1 if call_count[0] == 1: # First call: return markers with > 200 labels diff --git a/photutils/segmentation/tests/test_detect.py b/photutils/segmentation/tests/test_detect.py index 0bf08ee942..7fc90445f0 100644 --- a/photutils/segmentation/tests/test_detect.py +++ b/photutils/segmentation/tests/test_detect.py @@ -233,6 +233,19 @@ def test_small_sources(self): with pytest.warns(NoDetectionsWarning, match=match): detect_sources(self.data, threshold=0.9, n_pixels=5) + def test_nan_data(self): + """ + Test that NaN pixels are excluded from sources and that the + threshold comparison with NaN values does not emit warnings. + """ + data = np.ones((5, 5)) + data[2, 2] = np.nan + segm = detect_sources(data, threshold=0.5, n_pixels=5) + assert segm.n_labels == 1 + expected = np.ones((5, 5), dtype=np.int32) + expected[2, 2] = 0 + assert_equal(segm.data, expected) + def test_n_pixels(self): """ Test removal of sources whose size is less than n_pixels. diff --git a/photutils/segmentation/tests/test_finder.py b/photutils/segmentation/tests/test_finder.py index 80bcd6ecf6..04801ce0d6 100644 --- a/photutils/segmentation/tests/test_finder.py +++ b/photutils/segmentation/tests/test_finder.py @@ -14,7 +14,6 @@ from photutils.segmentation.finder import SourceFinder from photutils.segmentation.flags import SEGMENTATION_FLAGS from photutils.segmentation.utils import make_2dgaussian_kernel -from photutils.utils._optional_deps import HAS_SKIMAGE from photutils.utils.exceptions import NoDetectionsWarning @@ -25,12 +24,11 @@ class TestSourceFinder: threshold = 1.5 * 2.0 n_pixels = 10 - @pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_deblend(self): """ Test deblend. """ - finder = SourceFinder(n_pixels=self.n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=self.n_pixels) segm1 = finder(self.convolved_data, self.threshold) assert segm1.n_labels == 94 @@ -42,7 +40,7 @@ def test_invalid_units(self): """ Test invalid units. """ - finder = SourceFinder(n_pixels=self.n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=self.n_pixels) match = 'must all have the same units' with pytest.raises(ValueError, match=match): finder(self.convolved_data << u.uJy, self.threshold) @@ -55,8 +53,7 @@ def test_no_deblend(self): """ Test no deblend. """ - finder = SourceFinder(n_pixels=self.n_pixels, deblend=False, - progress_bar=False) + finder = SourceFinder(n_pixels=self.n_pixels, deblend=False) segm = finder(self.convolved_data, self.threshold) assert segm.n_labels == 87 @@ -64,15 +61,13 @@ def test_no_sources(self): """ Test no sources. """ - finder = SourceFinder(n_pixels=self.n_pixels, deblend=True, - progress_bar=False) + finder = SourceFinder(n_pixels=self.n_pixels, deblend=True) match = 'No sources were found' with pytest.warns(NoDetectionsWarning, match=match): segm = finder(self.convolved_data, 1000) assert segm is None - @pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_n_pixels_tuple(self): """ Test n_pixels tuple. @@ -95,22 +90,29 @@ def test_repr(self): """ Test repr. """ - finder = SourceFinder(n_pixels=self.n_pixels, deblend=False, - progress_bar=False) + finder = SourceFinder(n_pixels=self.n_pixels, deblend=False) cls_repr = repr(finder) assert cls_repr.startswith(finder.__class__.__name__) def test_finder_deprecations(): - finder = SourceFinder(n_pixels=10, progress_bar=False) + finder = SourceFinder(n_pixels=10) match = 'attribute was deprecated' with pytest.warns(AstropyDeprecationWarning, match=match): _ = finder.npixels with pytest.warns(AstropyDeprecationWarning, match=match): _ = finder.nlevels + # Each deprecated keyword emits exactly one warning + for kwargs in ({'progress_bar': False}, {'n_processes': 2}, + {'nproc': 2}): + name = next(iter(kwargs)) + with pytest.warns(AstropyDeprecationWarning, + match=name) as record: + SourceFinder(n_pixels=10, **kwargs) + assert len(record) == 1 + -@pytest.mark.skipif(not HAS_SKIMAGE, reason='skimage is required') def test_finder_flags_passthrough(): """ Test that deblending provenance flags survive the SourceFinder @@ -119,6 +121,6 @@ def test_finder_flags_passthrough(): yy, xx = np.mgrid[0:101, 0:101] data = (Gaussian2D(100, 50, 50, 5, 5)(xx, yy) + Gaussian2D(100, 35, 50, 5, 5)(xx, yy)) - finder = SourceFinder(n_pixels=5, progress_bar=False) + finder = SourceFinder(n_pixels=5) segm = finder(data, 10) assert np.any(segm.flags & SEGMENTATION_FLAGS.DEBLENDED) diff --git a/photutils/segmentation/utils.py b/photutils/segmentation/utils.py index 71d2379abc..fb6daa7ab1 100644 --- a/photutils/segmentation/utils.py +++ b/photutils/segmentation/utils.py @@ -98,7 +98,7 @@ def _make_binary_structure(ndim, connectivity): elif connectivity == 8: footprint = np.ones((3, 3), dtype=int) else: - msg = f'Invalid connectivity={connectivity} -- options are 4 or 8' + msg = f'Invalid connectivity={connectivity}. Options are 4 or 8' raise ValueError(msg) else: footprint = generate_binary_structure(ndim, 1) @@ -115,6 +115,12 @@ def _mask_to_mirrored_value(data, replace_mask, xycenter, *, mask=None, If the mirror pixel is unavailable (i.e., it is outside the image or masked), then the masked pixel value is set to zero. + Note that this function is not used in production. The + `~photutils.segmentation.SourceCatalog` batch kernels apply the + mirroring correction in compiled code. It is kept as the pure-Python + mirror of that correction and must track the compiled semantics + exactly, which is enforced by the batch cross-implementation tests. + Parameters ---------- data : 2D `~numpy.ndarray` diff --git a/photutils/utils/depths.py b/photutils/utils/depths.py index 1773156373..46fbca90c2 100644 --- a/photutils/utils/depths.py +++ b/photutils/utils/depths.py @@ -178,7 +178,7 @@ class ImageDepth: >>> convolved_data = convolve(data, kernel) >>> n_pixels = 10 >>> threshold = 3.2 - >>> finder = SourceFinder(n_pixels=n_pixels, progress_bar=False) + >>> finder = SourceFinder(n_pixels=n_pixels) >>> segment_map = finder(convolved_data, threshold) >>> mask = segment_map.make_source_mask() >>> radius = 4 @@ -208,7 +208,7 @@ class ImageDepth: convolved_data = convolve(data, kernel) n_pixels = 10 threshold = 3.2 - finder = SourceFinder(n_pixels=n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=n_pixels) segment_map = finder(convolved_data, threshold) mask = segment_map.make_source_mask() radius = 4 diff --git a/photutils/utils/tests/test_depths.py b/photutils/utils/tests/test_depths.py index 4df31f89cb..1431df925d 100644 --- a/photutils/utils/tests/test_depths.py +++ b/photutils/utils/tests/test_depths.py @@ -31,7 +31,7 @@ def setup_class(self): n_pixels = 10 threshold = 3.2 - finder = SourceFinder(n_pixels=n_pixels, progress_bar=False) + finder = SourceFinder(n_pixels=n_pixels) segment_map = finder(convolved_data, threshold) self.data = data self.mask = segment_map.make_source_mask()