From b41dbd18d89115193ae13e8abc5e0eb6c582a91e Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Fri, 11 Jul 2025 14:25:05 -0400 Subject: [PATCH 1/4] Add function to split search hit regions --- proseco/acq.py | 59 +++++++++++++++++++++++++++++++++- proseco/characteristics_acq.py | 5 +++ 2 files changed, 63 insertions(+), 1 deletion(-) diff --git a/proseco/acq.py b/proseco/acq.py index 3d5aa31e..ca74bc83 100644 --- a/proseco/acq.py +++ b/proseco/acq.py @@ -1273,6 +1273,53 @@ def get_spoiler_stars(stars, acq, box_size): return spoilers +def split_labeled_regions(dc_labeled, max_size=8): + """ + Split labeled regions into smaller subregions. + + This is used to ensure that the labeled regions corresponding to "search hits" are + not too large, which can happen if the dark current is mostly above the threshold. + The algorithm splits each labeled region into smaller subregions of size at most + `max_size` x `max_size` pixels. This ensures that a large region of high background + does not turn into just one isolated search hit. + + Parameters + ---------- + dc_labeled : np.ndarray + Labeled array where each contiguous region has a unique label. + max_size : int, optional + Maximum size of subregions to create, by default 8. + + Returns + ------- + dc_labeled : np.ndarray + New labeled array with subregions of size at most `max_size`. + n_regions : int + Number of labeled regions found. + """ + # Find slices for each labeled region + slices = ndimage.find_objects(dc_labeled) + new_label = 1 + dc_labeled_new = np.zeros_like(dc_labeled) + for i, slc in enumerate(slices): + region = dc_labeled[slc] == (i + 1) + region_shape = region.shape + + # Split region into 8x8 subregions if needed + for r_start in range(0, region_shape[0], max_size): + for c_start in range(0, region_shape[1], max_size): + r_end = min(r_start + max_size, region_shape[0]) + c_end = min(c_start + max_size, region_shape[1]) + subregion = region[r_start:r_end, c_start:c_end] + if np.any(subregion): + mask = np.zeros_like(region, dtype=bool) + mask[r_start:r_end, c_start:c_end] = subregion + dc_labeled_new[slc][mask] = new_label + new_label += 1 + + return dc_labeled_new, new_label - 1 + + def get_imposter_stars( dark, star_row, @@ -1303,6 +1350,10 @@ def get_imposter_stars( :returns: numpy structured array of imposter stars """ + print( + f"get_imposter_stars: star_row={star_row}, star_col={star_col}, " + f"{thresh=} {box_size=} {maxmag=} {bgd=} {mag_limit=}" + ) # Convert row/col to array index coords unless testing. rc_off = 0 if test else 512 acq_row = int(star_row + rc_off) @@ -1346,8 +1397,14 @@ def get_imposter_stars( # contiguous regions above ``thresh`` labeled with a unique index. # This is a one-line way of doing the PEA merging process, roughly. dc_labeled, n_hits = ndimage.label(dc2x2 > thresh) + if ACQ.search_hit_max_size is not None: + dc_labeled, n_hits = split_labeled_regions( + dc_labeled, max_size=ACQ.search_hit_max_size + ) if test: - print(dc_labeled) + # Print the entire dc_labeled array without clipping + with np.printoptions(threshold=100000, linewidth=200): + print(dc_labeled) # If no hits just return empty list if n_hits == 0: diff --git a/proseco/characteristics_acq.py b/proseco/characteristics_acq.py index e034d9c4..2b57a94d 100644 --- a/proseco/characteristics_acq.py +++ b/proseco/characteristics_acq.py @@ -100,3 +100,8 @@ def _get_fid_acq_stages(): # Index template file name index_template_file = "index_template_acq.html" + +# Search hit max size in pixels. This is used to limit the size of search hits +# that are used in the acquisition algorithm. See `split_labeled_regions` for +# details. Set to None to disable the size limit. +search_hit_max_size = 8 From 8791b0c1ebb86e065acab8ec20bfef9a7f10f77d Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Sat, 12 Jul 2025 07:09:43 -0400 Subject: [PATCH 2/4] Add test of split_labeled_regions --- proseco/tests/test_acq.py | 63 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/proseco/tests/test_acq.py b/proseco/tests/test_acq.py index 489695f6..4871fa03 100644 --- a/proseco/tests/test_acq.py +++ b/proseco/tests/test_acq.py @@ -23,6 +23,7 @@ get_image_props, get_imposter_stars, get_p_man_err, + split_labeled_regions, ) from ..catalog import get_aca_catalog from ..core import ACABox, StarsTable @@ -1375,3 +1376,65 @@ def test_acq_include_ids_all_halfws_full_catalog_with_spoiler(): assert np.all(aca.acqs["id"] == ids) assert np.allclose(aca.acqs["mag"], mags) assert np.all(aca.acqs["halfw"] == halfws) + + +img_input = np.array( + [ + [4, 4, 0, 1, 1, 0], + [0, 0, 0, 1, 1, 0], + [1, 1, 1, 1, 1, 0], + [0, 2, 2, 3, 1, 0], + [0, 2, 2, 0, 1, 0], + [0, 0, 0, 0, 0, 0], + ] +) + +test_split_labeled_regions_cases = [ + { + "img_exp": [ + [4, 4, 0, 1, 1, 0], + [0, 0, 0, 1, 1, 0], + [1, 1, 1, 1, 1, 0], + [0, 2, 2, 3, 1, 0], + [0, 2, 2, 0, 1, 0], + [0, 0, 0, 0, 0, 0], + ], + "max_size": 5, + "n_region": 4, + }, + { + "img_exp": [ + [6, 6, 0, 1, 2, 0], + [0, 0, 0, 1, 2, 0], + [1, 1, 1, 1, 2, 0], + [0, 4, 4, 5, 2, 0], + [0, 4, 4, 0, 3, 0], + [0, 0, 0, 0, 0, 0], + ], + "max_size": 4, + "n_region": 6, + }, + { + "img_exp": [ + [9, 9, 0, 1, 2, 0], + [0, 0, 0, 1, 2, 0], + [3, 3, 4, 4, 5, 0], + [0, 7, 7, 8, 5, 0], + [0, 7, 7, 0, 6, 0], + [0, 0, 0, 0, 0, 0], + ], + "max_size": 2, + "n_region": 9, + }, +] + + +@pytest.mark.parametrize("test_case", test_split_labeled_regions_cases) +def test_split_labeled_regions(test_case): + img_exp = test_case["img_exp"] + max_size = test_case["max_size"] + n_region = test_case["n_region"] + + img_split, n_regions = split_labeled_regions(img_input, max_size=max_size) + assert n_regions == n_region + assert np.all(img_split == img_exp) From 21fe38b8248d1fd924914fe363f19802de7e237b Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Sat, 12 Jul 2025 07:30:43 -0400 Subject: [PATCH 3/4] Remove debug print statement --- proseco/acq.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/proseco/acq.py b/proseco/acq.py index ca74bc83..1d13984f 100644 --- a/proseco/acq.py +++ b/proseco/acq.py @@ -1350,10 +1350,6 @@ def get_imposter_stars( :returns: numpy structured array of imposter stars """ - print( - f"get_imposter_stars: star_row={star_row}, star_col={star_col}, " - f"{thresh=} {box_size=} {maxmag=} {bgd=} {mag_limit=}" - ) # Convert row/col to array index coords unless testing. rc_off = 0 if test else 512 acq_row = int(star_row + rc_off) From 8a84dbc16ff25de4352ad8ea4960c157e26e7538 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Thu, 11 Sep 2025 12:46:42 -0400 Subject: [PATCH 4/4] More changes --- proseco/acq.py | 15 ++++++++++----- ruff.toml | 4 +--- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/proseco/acq.py b/proseco/acq.py index 1d13984f..6a706b99 100644 --- a/proseco/acq.py +++ b/proseco/acq.py @@ -1343,13 +1343,15 @@ def get_imposter_stars( :param star_col: col of acq star (float) :param thresh: PEA search hit threshold for a 2x2 block (e-/sec) :param maxmag: Max mag (alternate way to specify search hit ``thresh``) - :param box_size: box size (arcsec) + :param box_size: float or ACABox box size (arcsec) :param bgd: assumed flat background (float, e-/sec) :param mag_limit: Max mag for imposter (using 6x6 readout) :param test: hook for convenience in algorithm testing :returns: numpy structured array of imposter stars """ + box_size = ACABox(box_size) + # Convert row/col to array index coords unless testing. rc_off = 0 if test else 512 acq_row = int(star_row + rc_off) @@ -1393,6 +1395,11 @@ def get_imposter_stars( # contiguous regions above ``thresh`` labeled with a unique index. # This is a one-line way of doing the PEA merging process, roughly. dc_labeled, n_hits = ndimage.label(dc2x2 > thresh) + + # If no hits just return empty list + if n_hits == 0: + return [] + if ACQ.search_hit_max_size is not None: dc_labeled, n_hits = split_labeled_regions( dc_labeled, max_size=ACQ.search_hit_max_size @@ -1402,10 +1409,6 @@ def get_imposter_stars( with np.printoptions(threshold=100000, linewidth=200): print(dc_labeled) - # If no hits just return empty list - if n_hits == 0: - return [] - outs = [] for idx in range(n_hits): # Get row and col index vals for each merged region of search hits @@ -1453,6 +1456,7 @@ def get_imposter_stars( yang, zang = pixels_to_yagzag(row, col, allow_bad=True) out = ( + len(vals), row, col, row - star_row, @@ -1470,6 +1474,7 @@ def get_imposter_stars( if len(outs) > 0: dtype = [ + ("n_pix", "