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