Coverage for python/lsst/summit/utils/guiders/detection.py: 20%
294 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-17 10:13 +0000
1# This file is part of summit_utils.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21from __future__ import annotations
23import logging
24from dataclasses import asdict, dataclass, field
26__all__ = [
27 "runSourceDetection",
28 "buildReferenceCatalog",
29 "trackStarAcrossStamp",
30 "makeBlankCatalog",
31 "runGalSim",
32]
34import galsim
35import numpy as np
36import pandas as pd
37from astropy.nddata import Cutout2D
38from astropy.stats import sigma_clipped_stats
40import lsst.afw.detection as afwDetect
41from lsst.afw.image import ExposureF, ImageF, MaskedImageF
42from lsst.afw.math import STDEVCLIP, makeStatistics
44from .reading import GuiderData
46log = logging.getLogger(__name__)
48_DEFAULT_COLUMNS: str = (
49 "trackid detector expid elapsed_time dalt daz dtheta dx dy "
50 "fwhm xroi yroi xccd yccd xroi_ref yroi_ref xccd_ref yccd_ref "
51 "dxfp dyfp xfp yfp alt az xfp_ref yfp_ref alt_ref az_ref "
52 "xerr yerr theta theta_err theta_ref flux flux_err magoffset snr "
53 "ixx iyy ixy e1 e2 e1_altaz e2_altaz "
54 "ampname timestamp stamp detid filter exptime "
55)
56DEFAULT_COLUMNS: tuple[str, ...] = tuple(_DEFAULT_COLUMNS.split())
59def makeBlankCatalog() -> pd.DataFrame:
60 """
61 Create a blank DataFrame with the default columns for a star catalog.
63 Returns
64 -------
65 catalog : `pd.DataFrame`
66 Empty catalog with the default schema.
67 """
68 return pd.DataFrame(columns=DEFAULT_COLUMNS)
71@dataclass(frozen=True, slots=True)
72class GuiderStarTrackerConfig:
73 """Configuration for the GuiderStarTracker.
75 Parameters
76 ----------
77 minSnr : `float`
78 Minimum signal-to-noise ratio for star detection.
79 minValidStampFraction : `float`
80 Minimum fraction of stamps of valid detection per detector.
81 If provided, this is used instead of `minStampDetections`.
82 edgeMargin : `int`
83 Margin in pixels to avoid edge effects in the image.
84 maxEllipticity : `float`
85 Maximum allowed ellipticity for a star to be considered valid.
86 cutOutSize : `int`
87 Size of the cutout around the star for tracking.
88 aperSizeArcsec : `float`
89 Aperture size in arcseconds for star detection.
90 gain : `float`
91 Gain factor for the guider data, used in flux calculations.
92 nFallbackStamps : `int`
93 Number of individual stamps to sample when coadd detection fails.
94 peakSnrThreshold : `float`
95 Minimum peak pixel SNR to consider a stamp non-empty during
96 single-stamp fallback detection.
97 """
99 minSnr: float = 10.0
100 minValidStampFraction: float = 0.5
101 edgeMargin: int = 5
102 maxEllipticity: float = 0.7
103 cutOutSize: int = 50
104 aperSizeArcsec: float = 5.0
105 gain: float = 1.0
106 nFallbackStamps: int = 10
107 peakSnrThreshold: float = 5.0
110def trackStarAcrossStamp(
111 refCenter: tuple[float, float],
112 guiderData: GuiderData,
113 guiderName: str,
114 config: GuiderStarTrackerConfig = GuiderStarTrackerConfig(),
115) -> pd.DataFrame:
116 """
117 Track a star across all guider stamps and compute centroid, shape, and
118 flux.
120 GalSim is used for centroid and shape measurements. Flux is measured with
121 aperture photometry.
123 Parameters
124 ----------
125 refCenter : `tuple[float, float]`
126 Reference position (x, y) in pixel coordinates for the star.
127 guiderData : `GuiderData`
128 Guider data containing image stamps and metadata.
129 guiderName : `str`
130 Name of the guider to process.
131 config : `GuiderStarTrackerConfig`
132 Configuration parameters for the star tracker.
134 Returns
135 -------
136 stars : `pd.DataFrame`
137 DataFrame containing the tracked star measurements across all stamps.
138 """
139 gd = guiderData
140 expid = gd.expid
141 wcs = gd.getWcs(guiderName)
142 pixelScale = wcs.getPixelScale().asArcseconds()
144 # Initialize parameters from config
145 apertureRadius = config.aperSizeArcsec / pixelScale
146 cutOutSize = config.cutOutSize
147 gain = config.gain
149 # check if the ref center is within the image bounds
150 stampShape = gd[guiderName, 0].shape
151 if not (0 <= refCenter[0] < stampShape[1]) or not (0 <= refCenter[1] < stampShape[0]):
152 return makeBlankCatalog()
154 # loop over stamps
155 results = []
156 for i in range(len(gd)):
157 data = gd[guiderName, i]
158 star = measureStarOnStamp(data, refCenter, cutOutSize, apertureRadius, gain=gain).toDataFrame()
160 # Add stamp index
161 if not star.empty:
162 star["stamp"] = i
163 results.append(star)
165 # 3) Concatenate
166 if not results:
167 return makeBlankCatalog()
168 stars = pd.concat(results, ignore_index=True)
170 # 4) Add metadata
171 stars["detector"] = guiderName
172 stars["expid"] = expid
173 stars["ampname"] = gd.getGuiderAmpName(guiderName)
174 stars["detid"] = gd.getGuiderDetNum(guiderName)
175 stars["filter"] = gd.header.get("filter", "UNKNOWN")
176 stars["exptime"] = gd.guiderDurationSec
177 return stars
180def annulusBackgroundSubtraction(data: np.ndarray, annulus: tuple[float, float]) -> tuple[np.ndarray, float]:
181 """
182 Subtract background from the data using an annulus.
184 Parameters
185 ----------
186 data : `np.ndarray`
187 Image cutout data.
188 annulus : `tuple[float, float]`
189 Inner and outer radii (pixels) defining the background annulus.
191 Returns
192 -------
193 dataBkgSub : `np.ndarray`
194 Background-subtracted data.
195 bkgStd : `float`
196 Standard deviation of the background estimation.
197 """
198 rin, rout = annulus
199 x0, y0 = data.shape[1] // 2, data.shape[0] // 2
200 x, y = np.indices(data.shape)
201 annMask = ((x - x0) ** 2 + (y - y0) ** 2 >= rin**2) & ((x - x0) ** 2 + (y - y0) ** 2 <= rout**2)
202 annMask &= np.isfinite(data)
203 _, bkgSub, bkgStd = sigma_clipped_stats(data[annMask], sigma=3.0)
204 dataBkgSub = data - bkgSub
205 return dataBkgSub, bkgStd
208@dataclass
209class StarMeasurement:
210 xroi: float = field(default=np.nan)
211 yroi: float = field(default=np.nan)
212 xerr: float = field(default=0.0)
213 yerr: float = field(default=0.0)
214 e1: float = field(default=np.nan)
215 e2: float = field(default=np.nan)
216 ixx: float = field(default=np.nan)
217 iyy: float = field(default=np.nan)
218 ixy: float = field(default=np.nan)
219 fwhm: float = field(default=np.nan)
220 flux: float = field(default=np.nan)
221 flux_err: float = field(default=0.0)
222 snr: float = field(default=0.0)
224 def toDataFrame(self) -> pd.DataFrame:
225 """
226 Convert this measurement to a single-row DataFrame.
228 Returns
229 -------
230 row : `pd.DataFrame`
231 Single-row DataFrame with measurement fields, or empty if invalid.
232 """
233 d = asdict(self)
234 # Only drop the column if xroi is NaN (i.e., measurement failed)
235 if not np.isfinite(d.get("xroi", np.nan)):
236 # Return an empty DataFrame with all the keys as columns,
237 return pd.DataFrame(columns=list(d.keys()))
238 # Otherwise, return all columns, even if some are NaN
239 return pd.DataFrame([d])
241 def runAperturePhotometry(
242 self, cutout: np.ndarray, radius: float, bkgStd: float = 1.0, gain: float = 1.0
243 ) -> None:
244 """
245 Perform aperture photometry on a cutout image.
247 Updates the flux, flux_err, and snr attributes of the StarMeasurement.
249 Parameters
250 ----------
251 cutout : `np.ndarray`
252 2D cutout image (background-subtracted).
253 radius : `float`
254 Aperture radius in pixels.
255 bkgStd : `float`
256 Background RMS per pixel.
257 gain : `float`
258 Detector gain (e-/ADU).
259 """
260 x0, y0 = self.xroi, self.yroi
261 if np.isfinite(x0) and np.isfinite(y0):
262 ny, nx = cutout.shape
263 y, x = np.indices((ny, nx))
264 x0, y0 = self.xroi, self.yroi
266 # Background mask
267 aperMask = (x - x0) ** 2 + (y - y0) ** 2 <= radius**2
269 # Aperture sum
270 fluxNet = np.nansum(cutout[aperMask])
271 fluxNet = np.clip(fluxNet, 0, None) # Ensure non-negative flux
272 npix = aperMask.sum()
274 # Flux error
275 fluxErr = np.sqrt(fluxNet / gain + npix * bkgStd**2)
276 snr = fluxNet / (fluxErr + 1e-9) if fluxErr > 0 else 0.0
278 # Update the measurement
279 self.flux = fluxNet
280 self.flux_err = fluxErr
281 self.snr = snr
284def runSourceDetection(
285 image: np.ndarray,
286 threshold: float = 10,
287 cutOutSize: int = 25,
288 apertureRadius: int = 5,
289 gain: float = 1.0,
290 nPixMin: int = 10,
291) -> pd.DataFrame:
292 """
293 Detect sources in an image and measure their properties.
295 Parameters
296 ----------
297 image : `np.ndarray`
298 2D image array.
299 threshold : `float`
300 Detection threshold in sigma units.
301 cutOutSize : `int`
302 Size of the cutout around each detected source (pixels).
303 apertureRadius : `int`
304 Aperture radius in pixels for photometry.
305 gain : `float`
306 Detector gain (e-/ADU).
307 nPixMin : `int`
308 Minimum number of pixels in a footprint for detection.
310 Returns
311 -------
312 sources : `pd.DataFrame`
313 DataFrame with detected source properties.
314 """
315 # Step 1: Convert numpy image to MaskedImage and Exposure
316 exposure = ExposureF(MaskedImageF(ImageF(image)))
318 # Step 2: Detect sources using STDEVCLIP for the background noise.
319 # The input coadd images are dithered (see GuiderData.getStampArrayCoadd)
320 # to prevent integer quantization from collapsing the pixel distribution
321 # and causing STDEVCLIP to return 0. (See DM-54263.)
322 footprints = None
323 if not isBlankImage(image):
324 median = np.nanmedian(image)
325 exposure.image -= median
326 imageStd = float(makeStatistics(exposure.getMaskedImage(), STDEVCLIP).getValue(STDEVCLIP))
327 if imageStd <= 0:
328 # Fallback: sigma68 is robust to quantization and bright stars.
329 p16, p84 = np.nanpercentile(image, [16, 84])
330 imageStd = (p84 - p16) / 2.0
331 if imageStd <= 0:
332 exposure.image += median
333 return pd.DataFrame(columns=DEFAULT_COLUMNS)
334 absThreshold = threshold * imageStd
335 thresh = afwDetect.Threshold(absThreshold, afwDetect.Threshold.VALUE)
336 footprints = afwDetect.FootprintSet(exposure.getMaskedImage(), thresh, "DETECTED", nPixMin)
337 exposure.image += median
339 if not footprints:
340 return pd.DataFrame(columns=DEFAULT_COLUMNS)
342 nFootprints = len(footprints.getFootprints())
343 results = []
344 for fp in footprints.getFootprints():
345 # Create a cutout of the image around the footprint
346 refCenter = tuple(fp.getCentroid())
347 star = measureStarOnStamp(image, refCenter, cutOutSize, apertureRadius, gain).toDataFrame()
348 if not star.empty:
349 results.append(star)
350 if not results:
351 if nFootprints > 0:
352 log.warning(
353 f"FootprintSet found {nFootprints} sources but GalSim "
354 f"failed on all of them (cutOutSize={cutOutSize})."
355 )
356 return pd.DataFrame(columns=DEFAULT_COLUMNS)
357 df = pd.concat([sf for sf in results], ignore_index=True)
358 return df
361def measureStarOnStamp(
362 stamp: np.ndarray,
363 refCenter: tuple[float, float],
364 cutOutSize: int,
365 apertureRadius: int,
366 gain: float = 1.0,
367) -> StarMeasurement:
368 """
369 Measure a star on a single stamp: background subtraction, shape, centroid,
370 photometry.
372 Parameters
373 ----------
374 stamp : `np.ndarray`
375 Full stamp array.
376 refCenter : `tuple[float, float]`
377 Reference (x, y) pixel position for the cutout center.
378 cutOutSize : `int`
379 Size of the cutout in pixels.
380 apertureRadius : `int`
381 Aperture radius in pixels for photometry.
382 gain : `float`
383 Detector gain (e-/ADU).
385 Returns
386 -------
387 measurement : `StarMeasurement`
388 StarMeasurement object with populated fields (may be empty on failure).
389 """
390 cutout = getCutouts(stamp, refCenter, cutoutSize=cutOutSize)
391 data = cutout.data.copy()
393 if np.all(data == 0):
394 return StarMeasurement()
396 # Replace NaN (from border-clipped cutouts) with the median so
397 # GalSim can still measure stars near the stamp edge.
398 nan_mask = ~np.isfinite(data)
399 if nan_mask.any():
400 data[nan_mask] = np.nanmedian(data)
402 # 1) Subtract the background
403 annulus = (apertureRadius * 1.0, apertureRadius * 2)
404 dataBkgSub, bkgStd = annulusBackgroundSubtraction(data, annulus)
406 # 2) Track the star across all stamps for this guider
407 star = runGalSim(dataBkgSub, gain=gain, bkgStd=bkgStd)
409 # 3) Make aperture photometry measurements
410 # Galsim flux is the normalization of the Gaussian, not w/ fixed aper.
411 star.runAperturePhotometry(dataBkgSub, apertureRadius, gain=gain, bkgStd=bkgStd)
413 # 4) Add centroid and shape in amplifier roi coordinates
414 star.xroi += cutout.xmin_original
415 star.yroi += cutout.ymin_original
416 return star
419def runGalSim(
420 imageArray: np.ndarray,
421 gain: float = 1.0,
422 bkgStd: float = 0.0,
423) -> StarMeasurement:
424 """
425 Measure star properties with GalSim adaptive moments.
427 Parameters
428 ----------
429 imageArray : `np.ndarray`
430 Background-subtracted image cutout.
431 gain : `float`
432 Detector gain (e-/ADU).
433 bkgStd : `float`
434 Background RMS per pixel.
436 Returns
437 -------
438 result : `StarMeasurement`
439 Resulting measurement (empty if measurement failed).
440 """
441 gsImg = galsim.Image(imageArray)
442 hsmRes = galsim.hsm.FindAdaptiveMom(gsImg, strict=False)
443 success = hsmRes.error_message == ""
445 if not success:
446 result = StarMeasurement()
447 else:
448 xCentroid = hsmRes.moments_centroid.x
449 yCentroid = hsmRes.moments_centroid.y
450 flux = hsmRes.moments_amp
451 sigma = hsmRes.moments_sigma
452 e1 = hsmRes.observed_shape.e1
453 e2 = hsmRes.observed_shape.e2
454 fwhm = 2.355 * sigma
456 # Calculate errors using GalSim's error estimation
457 xErr, yErr = calcGalsimError(imageArray, hsmRes, gain=gain, bkgStd=bkgStd, correctForGain=True)
459 # Calculate SNR and flux error
460 ellipticity = np.sqrt(e1**2 + e2**2)
461 nEff = 2 * np.pi * sigma**2 * np.sqrt(1 - ellipticity**2)
462 shotNoise = np.sqrt(nEff * bkgStd**2)
463 fluxErr = np.sqrt(max(0, flux / gain) + shotNoise**2)
464 snr = flux / (shotNoise + 1e-9) if shotNoise > 0 else 0.0
466 # Calculate second moments
467 ixx = sigma**2 * (1 + e1)
468 iyy = sigma**2 * (1 - e1)
469 ixy = sigma**2 * e2
471 result = StarMeasurement(
472 xroi=xCentroid,
473 yroi=yCentroid,
474 xerr=xErr,
475 yerr=yErr,
476 e1=e1,
477 e2=e2,
478 ixx=ixx,
479 iyy=iyy,
480 ixy=ixy,
481 fwhm=fwhm,
482 flux=flux,
483 flux_err=fluxErr,
484 snr=snr,
485 )
486 return result
489def calcGalsimError(
490 imageArray: np.ndarray,
491 shape: galsim.hsm.ShapeData,
492 gain: float = 1.0,
493 bkgStd: float = 0.0,
494 correctForGain: bool = False,
495) -> tuple[float, float]:
496 """
497 Estimate centroid errors from GalSim HSMShapeData.
499 Parameters
500 ----------
501 imageArray : `np.ndarray`
502 Image cutout used for measurement.
503 shape : `galsim.hsm`
504 GalSim HSM shape data result object.
505 gain : `float`
506 Detector gain (e-/ADU), ignored if `correctForGain` is False.
507 bkgStd : `float`
508 Background RMS per pixel.
509 correctForGain : `bool`
510 Whether to include gain-dependent weighting.
512 Returns
513 -------
514 xerr : `float`
515 Estimated x centroid uncertainty (pixels).
516 yerr : `float`
517 Estimated y centroid uncertainty (pixels).
518 """
519 if not shape or shape.error_message != "":
520 return 0.0, 0.0
522 x0 = shape.moments_centroid.x
523 y0 = shape.moments_centroid.y
524 sigma = shape.moments_sigma
525 e1 = shape.observed_shape.e1
526 e2 = shape.observed_shape.e2
527 flux = shape.moments_amp
529 kernel = makeEllipticalGaussianStar(
530 shape=(imageArray.shape[0], imageArray.shape[1]),
531 e1=e1,
532 e2=e2,
533 flux=1,
534 sigma=sigma,
535 center=(x0, y0),
536 )
538 weight = np.ones_like(imageArray) / (bkgStd**2 + 1e-9)
539 if correctForGain:
540 weight = np.ones_like(imageArray) / (bkgStd**2 + np.abs(flux * kernel / gain))
542 mask = weight == 0.0
543 data = imageArray.copy()
544 if np.any(mask):
545 kernelMasked = kernel.copy()
546 data[mask] = kernelMasked[mask] * np.sum(data[~mask]) / np.sum(kernelMasked[~mask])
548 u, v = np.meshgrid(np.arange(imageArray.shape[1]) - x0, np.arange(imageArray.shape[0]) - y0)
549 usq = u**2
550 vsq = v**2
551 WI = kernel * data
552 M00 = np.nansum(WI)
553 WV = (kernel**2).astype(float)
554 WV[~mask] /= weight[~mask]
555 WV[mask] /= np.median(weight[~mask])
556 WV = WV / float(M00**2)
558 varM10 = 4 * np.sum(WV * usq)
559 varM01 = 4 * np.sum(WV * vsq)
560 xerr = np.sqrt(varM10)
561 yerr = np.sqrt(varM01)
562 return xerr, yerr
565def makeEllipticalGaussianStar(
566 shape: tuple[int, int],
567 flux: float,
568 sigma: float,
569 e1: float,
570 e2: float,
571 center: tuple[float, float],
572) -> np.ndarray:
573 """
574 Create an elliptical 2D Gaussian star with specified parameters.
576 Parameters
577 ----------
578 shape : `tuple[int, int]`
579 (ny, nx) output array shape.
580 flux : `float`
581 Total flux (normalization).
582 sigma : `float`
583 Gaussian sigma (pixels).
584 e1 : `float`
585 Ellipticity component e1.
586 e2 : `float`
587 Ellipticity component e2.
588 center : `tuple[float, float]`
589 (x0, y0) centroid position in pixels.
591 Returns
592 -------
593 image : `np.ndarray`
594 Generated model image.
595 """
596 y, x = np.indices(shape)
597 x0, y0 = center
598 u = x - x0
599 v = y - y0
601 # Second-moment matrix elements
602 ixx = sigma**2 * (1 + e1)
603 iyy = sigma**2 * (1 - e1)
604 ixy = sigma**2 * e2
606 # Inverse covariance matrix
607 det = ixx * iyy - ixy**2
608 invIxx = iyy / det
609 invIyy = ixx / det
610 invIxy = -ixy / det
612 # Quadratic form: u^2 * invIxx + v^2 * invIyy + 2uv * invIxy
613 r2 = invIxx * u**2 + invIyy * v**2 + 2 * invIxy * u * v
615 e = np.sqrt(e1**2 + e2**2)
616 norm = flux / (2 * np.pi * sigma**2 * np.sqrt(1 - e**2))
617 image = norm * np.exp(-0.5 * r2)
618 return image
621def _detectOnSingleStamps(
622 guiderData: GuiderData,
623 guiderName: str,
624 config: GuiderStarTrackerConfig,
625 apertureRadius: int,
626 log: logging.Logger,
627) -> pd.DataFrame:
628 """Try source detection on individual stamps when coadd detection fails.
630 Samples up to ``config.nFallbackStamps`` stamps evenly spread across
631 the sequence. Returns the detection with the highest SNR, or an empty
632 DataFrame if no star is found on any stamp.
634 A stamp is considered truly empty if its peak pixel SNR
635 (max - median) / std is below ``config.peakSnrThreshold``.
636 """
637 nStamps = len(guiderData)
638 if nStamps == 0:
639 return pd.DataFrame(columns=DEFAULT_COLUMNS)
641 nSample = min(config.nFallbackStamps, nStamps)
642 indices = np.linspace(0, nStamps - 1, nSample, dtype=int)
644 bestSources = pd.DataFrame(columns=DEFAULT_COLUMNS)
645 bestSnr = 0.0
646 nTrulyEmpty = 0
648 for idx in indices:
649 stamp = guiderData[guiderName, idx].astype(np.float32)
650 arr = stamp - np.nanmin(stamp)
652 # Quick peak-SNR check before running full detection
653 med = np.nanmedian(arr)
654 std = np.nanstd(arr)
655 if std > 0:
656 peakSnr = (np.nanmax(arr) - med) / std
657 else:
658 peakSnr = 0.0
660 if peakSnr < config.peakSnrThreshold:
661 nTrulyEmpty += 1
662 continue
664 sources = runSourceDetection(
665 arr,
666 threshold=config.minSnr,
667 apertureRadius=apertureRadius,
668 cutOutSize=config.cutOutSize,
669 gain=config.gain,
670 )
671 if not sources.empty and sources["snr"].max() > bestSnr:
672 bestSnr = sources["snr"].max()
673 bestSources = sources
675 if not bestSources.empty:
676 log.info(
677 f"Single-stamp fallback recovered source on {guiderName} "
678 f"(SNR={bestSnr:.1f}, {nTrulyEmpty}/{nSample} stamps empty)"
679 )
681 return bestSources
684def buildReferenceCatalog(
685 guiderData: GuiderData,
686 log: logging.Logger,
687 config: GuiderStarTrackerConfig = GuiderStarTrackerConfig(),
688) -> pd.DataFrame:
689 """
690 Build a reference star catalog from each guider's coadded stamp.
692 Parameters
693 ----------
694 guiderData : `GuiderData`
695 Guider dataset containing stamps and metadata.
696 log : `logging.Logger`
697 Logger for warnings and diagnostics.
698 config : `GuiderStarTrackerConfig`
699 Star tracker configuration.
701 Returns
702 -------
703 refCatalog : `pd.DataFrame`
704 Concatenated reference catalog of brightest stars per guider.
705 """
706 expId = guiderData.expid
707 minSnr = config.minSnr
708 gain = config.gain
709 cutOutSize = config.cutOutSize
711 tableList = []
712 for guiderName in guiderData.guiderNames:
713 pixelScale = guiderData.getWcs(guiderName).getPixelScale().asArcseconds()
714 apertureRadius = int(config.aperSizeArcsec / pixelScale)
716 array = guiderData.getStampArrayCoadd(guiderName)
717 array = array - np.nanmin(array) # Ensure no negative values
718 sources = runSourceDetection(
719 array,
720 threshold=minSnr,
721 apertureRadius=apertureRadius,
722 cutOutSize=cutOutSize,
723 gain=gain,
724 )
725 detectionMethod = "coadd"
726 if sources.empty:
727 # Fallback: try detection on individual stamps. The coadd can
728 # wash out stars that drift between frames or have artifacts.
729 sources = _detectOnSingleStamps(guiderData, guiderName, config, apertureRadius, log)
730 detectionMethod = "single_stamp"
731 if sources.empty:
732 log.warning(f"No sources detected in `buildReferenceCatalog`for {guiderName} in {expId}. ")
733 continue
735 sources["detection_method"] = detectionMethod
737 sources.sort_values(by=["snr"], ascending=False, inplace=True)
738 sources.reset_index(drop=True, inplace=True)
740 detNum = guiderData.getGuiderDetNum(guiderName)
741 sources["detector"] = guiderName
742 sources["detid"] = detNum
743 sources["starid"] = detNum * 100
744 tableList.append(sources)
746 if len(tableList) == 0:
747 log.warning(f"buildReferenceCatalog failed - no stars detected in any guider for {expId}.")
748 return makeBlankCatalog()
750 refCatalog = pd.concat(tableList, ignore_index=True)
751 return refCatalog
754def getCutouts(imageArray: np.ndarray, refCenter: tuple[float, float], cutoutSize: int = 25) -> Cutout2D:
755 """
756 Get a cutout at the reference position from an image array.
758 Parameters
759 ----------
760 imageArray : `np.ndarray`
761 Full image array.
762 refCenter : `tuple[float, float]`
763 (x, y) center for the cutout in pixels.
764 cutoutSize : `int`
765 Size (pixels) of the square cutout.
767 Returns
768 -------
769 cutout : `Cutout2D`
770 Astropy Cutout2D object.
771 """
772 refX, refY = refCenter
773 return Cutout2D(imageArray, (refX, refY), size=cutoutSize, mode="partial", fill_value=np.nan)
776def isBlankImage(image: np.ndarray, fluxMin: float = 300, peakSnrMin: float = 5.0) -> bool:
777 """
778 Returns True if the image has no significant source (e.g., no star).
780 An image is considered non-blank if the peak flux above the median
781 exceeds ``fluxMin`` OR the peak pixel SNR (using MAD-based robust
782 std) exceeds ``peakSnrMin``.
784 Parameters
785 ----------
786 image : `np.ndarray`
787 2D image data.
788 fluxMin : `float`
789 Minimum peak flux above median (ADU) to consider non-blank.
790 peakSnrMin : `float`
791 Minimum peak pixel SNR to consider non-blank.
793 Returns
794 -------
795 bool
796 True if the image is blank, False otherwise.
797 """
798 med = np.nanmedian(image)
799 peakFlux = np.nanmax(image) - med
801 # Absolute flux check
802 if peakFlux > fluxMin:
803 return False
805 # SNR check with MAD-based robust std
806 mad = np.nanmedian(np.abs(image - med))
807 std = 1.4826 * mad
808 if std <= 0:
809 return True
810 peakSnr = peakFlux / std
811 return peakSnr < peakSnrMin