Coverage for python / lsst / summit / utils / guiders / detection.py: 19%
233 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 09:17 +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
40from lsst.afw.image import ExposureF, ImageF, MaskedImageF
41from lsst.summit.utils.utils import detectObjectsInExp
43from .reading import GuiderData
45_DEFAULT_COLUMNS: str = (
46 "trackid detector expid elapsed_time dalt daz dtheta dx dy "
47 "fwhm xroi yroi xccd yccd xroi_ref yroi_ref xccd_ref yccd_ref "
48 "dxfp dyfp xfp yfp alt az xfp_ref yfp_ref alt_ref az_ref "
49 "xerr yerr theta theta_err theta_ref flux flux_err magoffset snr "
50 "ixx iyy ixy e1 e2 e1_altaz e2_altaz "
51 "ampname timestamp stamp detid filter exptime "
52)
53DEFAULT_COLUMNS: tuple[str, ...] = tuple(_DEFAULT_COLUMNS.split())
56def makeBlankCatalog() -> pd.DataFrame:
57 """
58 Create a blank DataFrame with the default columns for a star catalog.
60 Returns
61 -------
62 catalog : `pd.DataFrame`
63 Empty catalog with the default schema.
64 """
65 return pd.DataFrame(columns=DEFAULT_COLUMNS)
68@dataclass(frozen=True, slots=True)
69class GuiderStarTrackerConfig:
70 """Configuration for the GuiderStarTracker.
72 Parameters
73 ----------
74 minSnr : `float`
75 Minimum signal-to-noise ratio for star detection.
76 minValidStampFraction : `float`
77 Minimum fraction of stamps of valid detection per detector.
78 If provided, this is used instead of `minStampDetections`.
79 edgeMargin : `int`
80 Margin in pixels to avoid edge effects in the image.
81 maxEllipticity : `float`
82 Maximum allowed ellipticity for a star to be considered valid.
83 cutOutSize : `int`
84 Size of the cutout around the star for tracking.
85 aperSizeArcsec : `float`
86 Aperture size in arcseconds for star detection.
87 gain : `float`
88 Gain factor for the guider data, used in flux calculations.
89 """
91 minSnr: float = 10.0
92 minValidStampFraction: float = 0.5
93 edgeMargin: int = 5
94 maxEllipticity: float = 0.5
95 cutOutSize: int = 50
96 aperSizeArcsec: float = 3.0
97 gain: float = 1.0
100def trackStarAcrossStamp(
101 refCenter: tuple[float, float],
102 guiderData: GuiderData,
103 guiderName: str,
104 config: GuiderStarTrackerConfig = GuiderStarTrackerConfig(),
105) -> pd.DataFrame:
106 """
107 Track a star across all guider stamps and compute centroid, shape, and
108 flux.
110 GalSim is used for centroid and shape measurements. Flux is measured with
111 aperture photometry.
113 Parameters
114 ----------
115 refCenter : `tuple[float, float]`
116 Reference position (x, y) in pixel coordinates for the star.
117 guiderData : `GuiderData`
118 Guider data containing image stamps and metadata.
119 guiderName : `str`
120 Name of the guider to process.
121 config : `GuiderStarTrackerConfig`
122 Configuration parameters for the star tracker.
124 Returns
125 -------
126 stars : `pd.DataFrame`
127 DataFrame containing the tracked star measurements across all stamps.
128 """
129 gd = guiderData
130 expid = gd.expid
131 wcs = gd.getWcs(guiderName)
132 pixelScale = wcs.getPixelScale().asArcseconds()
134 # Initialize parameters from config
135 apertureRadius = config.aperSizeArcsec / pixelScale
136 cutOutSize = config.cutOutSize
137 gain = config.gain
139 # check if the ref center is within the image bounds
140 stampShape = gd[guiderName, 0].shape
141 if not (0 <= refCenter[0] < stampShape[1]) or not (0 <= refCenter[1] < stampShape[0]):
142 return makeBlankCatalog()
144 # loop over stamps
145 results = []
146 for i in range(len(gd)):
147 data = gd[guiderName, i]
148 star = measureStarOnStamp(data, refCenter, cutOutSize, apertureRadius, gain=gain).toDataFrame()
150 # Add stamp index
151 if not star.empty:
152 star["stamp"] = i
153 results.append(star)
155 # 3) Concatenate
156 if not results:
157 return makeBlankCatalog()
158 stars = pd.concat(results, ignore_index=True)
160 # 4) Add metadata
161 stars["detector"] = guiderName
162 stars["expid"] = expid
163 stars["ampname"] = gd.getGuiderAmpName(guiderName)
164 stars["detid"] = gd.getGuiderDetNum(guiderName)
165 stars["filter"] = gd.header.get("filter", "UNKNOWN")
166 stars["exptime"] = gd.guiderDurationSec
167 return stars
170def annulusBackgroundSubtraction(data: np.ndarray, annulus: tuple[float, float]) -> tuple[np.ndarray, float]:
171 """
172 Subtract background from the data using an annulus.
174 Parameters
175 ----------
176 data : `np.ndarray`
177 Image cutout data.
178 annulus : `tuple[float, float]`
179 Inner and outer radii (pixels) defining the background annulus.
181 Returns
182 -------
183 dataBkgSub : `np.ndarray`
184 Background-subtracted data.
185 bkgStd : `float`
186 Standard deviation of the background estimation.
187 """
188 rin, rout = annulus
189 x0, y0 = data.shape[1] // 2, data.shape[0] // 2
190 x, y = np.indices(data.shape)
191 annMask = ((x - x0) ** 2 + (y - y0) ** 2 >= rin**2) & ((x - x0) ** 2 + (y - y0) ** 2 <= rout**2)
192 annMask &= np.isfinite(data)
193 _, bkgSub, bkgStd = sigma_clipped_stats(data[annMask], sigma=3.0)
194 dataBkgSub = data - bkgSub
195 return dataBkgSub, bkgStd
198@dataclass
199class StarMeasurement:
200 xroi: float = field(default=np.nan)
201 yroi: float = field(default=np.nan)
202 xerr: float = field(default=0.0)
203 yerr: float = field(default=0.0)
204 e1: float = field(default=np.nan)
205 e2: float = field(default=np.nan)
206 ixx: float = field(default=np.nan)
207 iyy: float = field(default=np.nan)
208 ixy: float = field(default=np.nan)
209 fwhm: float = field(default=np.nan)
210 flux: float = field(default=np.nan)
211 flux_err: float = field(default=0.0)
212 snr: float = field(default=0.0)
214 def toDataFrame(self) -> pd.DataFrame:
215 """
216 Convert this measurement to a single-row DataFrame.
218 Returns
219 -------
220 row : `pd.DataFrame`
221 Single-row DataFrame with measurement fields, or empty if invalid.
222 """
223 d = asdict(self)
224 # Only drop the column if xroi is NaN (i.e., measurement failed)
225 if not np.isfinite(d.get("xroi", np.nan)):
226 # Return an empty DataFrame with all the keys as columns,
227 return pd.DataFrame(columns=list(d.keys()))
228 # Otherwise, return all columns, even if some are NaN
229 return pd.DataFrame([d])
231 def runAperturePhotometry(
232 self, cutout: np.ndarray, radius: float, bkgStd: float = 1.0, gain: float = 1.0
233 ) -> None:
234 """
235 Perform aperture photometry on a cutout image.
237 Updates the flux, flux_err, and snr attributes of the StarMeasurement.
239 Parameters
240 ----------
241 cutout : `np.ndarray`
242 2D cutout image (background-subtracted).
243 radius : `float`
244 Aperture radius in pixels.
245 bkgStd : `float`
246 Background RMS per pixel.
247 gain : `float`
248 Detector gain (e-/ADU).
249 """
250 x0, y0 = self.xroi, self.yroi
251 if np.isfinite(x0) and np.isfinite(y0):
252 ny, nx = cutout.shape
253 y, x = np.indices((ny, nx))
254 x0, y0 = self.xroi, self.yroi
256 # Background mask
257 aperMask = (x - x0) ** 2 + (y - y0) ** 2 <= radius**2
259 # Aperture sum
260 fluxNet = np.nansum(cutout[aperMask])
261 fluxNet = np.clip(fluxNet, 0, None) # Ensure non-negative flux
262 npix = aperMask.sum()
264 # Flux error
265 fluxErr = np.sqrt(fluxNet / gain + npix * bkgStd**2)
266 snr = fluxNet / (fluxErr + 1e-9) if fluxErr > 0 else 0.0
268 # Update the measurement
269 self.flux = fluxNet
270 self.flux_err = fluxErr
271 self.snr = snr
274def runSourceDetection(
275 image: np.ndarray,
276 threshold: float = 10,
277 cutOutSize: int = 25,
278 apertureRadius: int = 5,
279 gain: float = 1.0,
280) -> pd.DataFrame:
281 """
282 Detect sources in an image and measure their properties.
284 Parameters
285 ----------
286 image : `np.ndarray`
287 2D image array.
288 threshold : `float`
289 Detection threshold in sigma units.
290 cutOutSize : `int`
291 Size of the cutout around each detected source (pixels).
292 apertureRadius : `int`
293 Aperture radius in pixels for photometry.
294 gain : `float`
295 Detector gain (e-/ADU).
297 Returns
298 -------
299 sources : `pd.DataFrame`
300 DataFrame with detected source properties.
301 """
302 # Step 1: Convert numpy image to MaskedImage and Exposure
303 exposure = ExposureF(MaskedImageF(ImageF(image)))
305 # Step 2: Detect sources
306 # we assume that we have bright stars
307 # filter out stamps with no stars
308 if not isBlankImage(image):
309 footprints = detectObjectsInExp(exposure, nSigma=threshold)
310 else:
311 footprints = None
313 if not footprints:
314 return pd.DataFrame(columns=DEFAULT_COLUMNS)
316 results = []
317 for fp in footprints.getFootprints():
318 # Create a cutout of the image around the footprint
319 refCenter = tuple(fp.getCentroid())
320 star = measureStarOnStamp(image, refCenter, cutOutSize, apertureRadius, gain).toDataFrame()
321 if not star.empty:
322 results.append(star)
323 if not results:
324 return pd.DataFrame(columns=DEFAULT_COLUMNS)
325 df = pd.concat([sf for sf in results], ignore_index=True)
326 return df
329def measureStarOnStamp(
330 stamp: np.ndarray,
331 refCenter: tuple[float, float],
332 cutOutSize: int,
333 apertureRadius: int,
334 gain: float = 1.0,
335) -> StarMeasurement:
336 """
337 Measure a star on a single stamp: background subtraction, shape, centroid,
338 photometry.
340 Parameters
341 ----------
342 stamp : `np.ndarray`
343 Full stamp array.
344 refCenter : `tuple[float, float]`
345 Reference (x, y) pixel position for the cutout center.
346 cutOutSize : `int`
347 Size of the cutout in pixels.
348 apertureRadius : `int`
349 Aperture radius in pixels for photometry.
350 gain : `float`
351 Detector gain (e-/ADU).
353 Returns
354 -------
355 measurement : `StarMeasurement`
356 StarMeasurement object with populated fields (may be empty on failure).
357 """
358 cutout = getCutouts(stamp, refCenter, cutoutSize=cutOutSize)
359 data = cutout.data
361 if np.all(data == 0) | (not np.isfinite(data).all()):
362 return StarMeasurement()
364 # 1) Subtract the background
365 annulus = (apertureRadius * 1.0, apertureRadius * 2)
366 dataBkgSub, bkgStd = annulusBackgroundSubtraction(data, annulus)
368 # 2) Track the star across all stamps for this guider
369 star = runGalSim(dataBkgSub, gain=gain, bkgStd=bkgStd)
371 # 3) Make aperture photometry measurements
372 # Galsim flux is the normalization of the Gaussian, not w/ fixed aper.
373 star.runAperturePhotometry(dataBkgSub, apertureRadius, gain=gain, bkgStd=bkgStd)
375 # 4) Add centroid and shape in amplifier roi coordinates
376 star.xroi += cutout.xmin_original
377 star.yroi += cutout.ymin_original
378 return star
381def runGalSim(
382 imageArray: np.ndarray,
383 gain: float = 1.0,
384 bkgStd: float = 0.0,
385) -> StarMeasurement:
386 """
387 Measure star properties with GalSim adaptive moments.
389 Parameters
390 ----------
391 imageArray : `np.ndarray`
392 Background-subtracted image cutout.
393 gain : `float`
394 Detector gain (e-/ADU).
395 bkgStd : `float`
396 Background RMS per pixel.
398 Returns
399 -------
400 result : `StarMeasurement`
401 Resulting measurement (empty if measurement failed).
402 """
403 gsImg = galsim.Image(imageArray)
404 hsmRes = galsim.hsm.FindAdaptiveMom(gsImg, strict=False)
405 success = hsmRes.error_message == ""
407 if not success:
408 result = StarMeasurement()
409 else:
410 xCentroid = hsmRes.moments_centroid.x
411 yCentroid = hsmRes.moments_centroid.y
412 flux = hsmRes.moments_amp
413 sigma = hsmRes.moments_sigma
414 e1 = hsmRes.observed_shape.e1
415 e2 = hsmRes.observed_shape.e2
416 fwhm = 2.355 * sigma
418 # Calculate errors using GalSim's error estimation
419 xErr, yErr = calcGalsimError(imageArray, hsmRes, gain=gain, bkgStd=bkgStd, correctForGain=True)
421 # Calculate SNR and flux error
422 ellipticity = np.sqrt(e1**2 + e2**2)
423 nEff = 2 * np.pi * sigma**2 * np.sqrt(1 - ellipticity**2)
424 shotNoise = np.sqrt(nEff * bkgStd**2)
425 fluxErr = np.sqrt(flux / gain + shotNoise**2)
426 snr = flux / (shotNoise + 1e-9) if shotNoise > 0 else 0.0
428 # Calculate second moments
429 ixx = sigma**2 * (1 + e1)
430 iyy = sigma**2 * (1 - e1)
431 ixy = sigma**2 * e2
433 result = StarMeasurement(
434 xroi=xCentroid,
435 yroi=yCentroid,
436 xerr=xErr,
437 yerr=yErr,
438 e1=e1,
439 e2=e2,
440 ixx=ixx,
441 iyy=iyy,
442 ixy=ixy,
443 fwhm=fwhm,
444 flux=flux,
445 flux_err=fluxErr,
446 snr=snr,
447 )
448 return result
451def calcGalsimError(
452 imageArray: np.ndarray,
453 shape: galsim.hsm.ShapeData,
454 gain: float = 1.0,
455 bkgStd: float = 0.0,
456 correctForGain: bool = False,
457) -> tuple[float, float]:
458 """
459 Estimate centroid errors from GalSim HSMShapeData.
461 Parameters
462 ----------
463 imageArray : `np.ndarray`
464 Image cutout used for measurement.
465 shape : `galsim.hsm`
466 GalSim HSM shape data result object.
467 gain : `float`
468 Detector gain (e-/ADU), ignored if `correctForGain` is False.
469 bkgStd : `float`
470 Background RMS per pixel.
471 correctForGain : `bool`
472 Whether to include gain-dependent weighting.
474 Returns
475 -------
476 xerr : `float`
477 Estimated x centroid uncertainty (pixels).
478 yerr : `float`
479 Estimated y centroid uncertainty (pixels).
480 """
481 if not shape or shape.error_message != "":
482 return 0.0, 0.0
484 x0 = shape.moments_centroid.x
485 y0 = shape.moments_centroid.y
486 sigma = shape.moments_sigma
487 e1 = shape.observed_shape.e1
488 e2 = shape.observed_shape.e2
489 flux = shape.moments_amp
491 kernel = makeEllipticalGaussianStar(
492 shape=(imageArray.shape[0], imageArray.shape[1]),
493 e1=e1,
494 e2=e2,
495 flux=1,
496 sigma=sigma,
497 center=(x0, y0),
498 )
500 weight = np.ones_like(imageArray) / (bkgStd**2 + 1e-9)
501 if correctForGain:
502 weight = np.ones_like(imageArray) / (bkgStd**2 + np.abs(flux * kernel / gain))
504 mask = weight == 0.0
505 data = imageArray.copy()
506 if np.any(mask):
507 kernelMasked = kernel.copy()
508 data[mask] = kernelMasked[mask] * np.sum(data[~mask]) / np.sum(kernelMasked[~mask])
510 u, v = np.meshgrid(np.arange(imageArray.shape[1]) - x0, np.arange(imageArray.shape[0]) - y0)
511 usq = u**2
512 vsq = v**2
513 WI = kernel * data
514 M00 = np.nansum(WI)
515 WV = (kernel**2).astype(float)
516 WV[~mask] /= weight[~mask]
517 WV[mask] /= np.median(weight[~mask])
518 WV = WV / float(M00**2)
520 varM10 = 4 * np.sum(WV * usq)
521 varM01 = 4 * np.sum(WV * vsq)
522 xerr = np.sqrt(varM10)
523 yerr = np.sqrt(varM01)
524 return xerr, yerr
527def makeEllipticalGaussianStar(
528 shape: tuple[int, int],
529 flux: float,
530 sigma: float,
531 e1: float,
532 e2: float,
533 center: tuple[float, float],
534) -> np.ndarray:
535 """
536 Create an elliptical 2D Gaussian star with specified parameters.
538 Parameters
539 ----------
540 shape : `tuple[int, int]`
541 (ny, nx) output array shape.
542 flux : `float`
543 Total flux (normalization).
544 sigma : `float`
545 Gaussian sigma (pixels).
546 e1 : `float`
547 Ellipticity component e1.
548 e2 : `float`
549 Ellipticity component e2.
550 center : `tuple[float, float]`
551 (x0, y0) centroid position in pixels.
553 Returns
554 -------
555 image : `np.ndarray`
556 Generated model image.
557 """
558 y, x = np.indices(shape)
559 x0, y0 = center
560 u = x - x0
561 v = y - y0
563 # Second-moment matrix elements
564 ixx = sigma**2 * (1 + e1)
565 iyy = sigma**2 * (1 - e1)
566 ixy = sigma**2 * e2
568 # Inverse covariance matrix
569 det = ixx * iyy - ixy**2
570 invIxx = iyy / det
571 invIyy = ixx / det
572 invIxy = -ixy / det
574 # Quadratic form: u^2 * invIxx + v^2 * invIyy + 2uv * invIxy
575 r2 = invIxx * u**2 + invIyy * v**2 + 2 * invIxy * u * v
577 e = np.sqrt(e1**2 + e2**2)
578 norm = flux / (2 * np.pi * sigma**2 * np.sqrt(1 - e**2))
579 image = norm * np.exp(-0.5 * r2)
580 return image
583def buildReferenceCatalog(
584 guiderData: GuiderData,
585 log: logging.Logger,
586 config: GuiderStarTrackerConfig = GuiderStarTrackerConfig(),
587) -> pd.DataFrame:
588 """
589 Build a reference star catalog from each guider's coadded stamp.
591 Parameters
592 ----------
593 guiderData : `GuiderData`
594 Guider dataset containing stamps and metadata.
595 log : `logging.Logger`
596 Logger for warnings and diagnostics.
597 config : `GuiderStarTrackerConfig`
598 Star tracker configuration.
600 Returns
601 -------
602 refCatalog : `pd.DataFrame`
603 Concatenated reference catalog of brightest stars per guider.
604 """
605 expId = guiderData.expid
606 minSnr = config.minSnr
607 gain = config.gain
608 cutOutSize = config.cutOutSize
610 tableList = []
611 for guiderName in guiderData.guiderNames:
612 pixelScale = guiderData.getWcs(guiderName).getPixelScale().asArcseconds()
613 apertureRadius = int(config.aperSizeArcsec / pixelScale)
615 array = guiderData.getStampArrayCoadd(guiderName)
616 # array = np.where(array < 0, 0, array) # Ensure no negative values
617 array = array - np.nanmin(array) # Ensure no negative values
618 sources = runSourceDetection(
619 array,
620 threshold=minSnr,
621 apertureRadius=apertureRadius,
622 cutOutSize=cutOutSize,
623 gain=gain,
624 )
625 if sources.empty:
626 log.warning(f"No sources detected in `buildReferenceCatalog`" f"for {guiderName} in {expId}. ")
627 continue
629 sources.sort_values(by=["snr"], ascending=False, inplace=True)
630 sources.reset_index(drop=True, inplace=True)
632 detNum = guiderData.getGuiderDetNum(guiderName)
633 sources["detector"] = guiderName
634 sources["detid"] = detNum
635 sources["starid"] = detNum * 100
636 tableList.append(sources)
638 if len(tableList) == 0:
639 log.warning(f"buildReferenceCatalog failed - no stars detected in any guider for {expId}.")
640 return makeBlankCatalog()
642 refCatalog = pd.concat(tableList, ignore_index=True)
643 return refCatalog
646def getCutouts(imageArray: np.ndarray, refCenter: tuple[float, float], cutoutSize: int = 25) -> Cutout2D:
647 """
648 Get a cutout at the reference position from an image array.
650 Parameters
651 ----------
652 imageArray : `np.ndarray`
653 Full image array.
654 refCenter : `tuple[float, float]`
655 (x, y) center for the cutout in pixels.
656 cutoutSize : `int`
657 Size (pixels) of the square cutout.
659 Returns
660 -------
661 cutout : `Cutout2D`
662 Astropy Cutout2D object.
663 """
664 refX, refY = refCenter
665 return Cutout2D(imageArray, (refX, refY), size=cutoutSize, mode="partial", fill_value=np.nan)
668def isBlankImage(image: np.ndarray, fluxMin: int = 300) -> bool:
669 """
670 Returns True if the image has no significant source (e.g., no star).
672 Parameters
673 ----------
674 image : 2D array
675 Image data (float or int).
676 fluxMin : float
677 Minimum deviation from the background median to be considered a source.
679 Returns
680 -------
681 bool
682 True if the image is blank, False otherwise.
683 (no pixel deviates more than flux_min)
684 """
685 med = np.nanmedian(image)
686 diff = np.abs(image - med)
687 return not np.any(diff > fluxMin)