Coverage for python/lsst/summit/utils/guiders/detection.py: 19%
229 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-11 02:56 -0700
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 minStampDetections : `int`
77 Minimum number of detections across all stamps for a star to be
78 considered valid.
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 minStampDetections: int = 30
93 edgeMargin: int = 25
94 maxEllipticity: float = 0.3
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 footprints = detectObjectsInExp(exposure, nSigma=threshold)
308 if not footprints:
309 return pd.DataFrame(columns=DEFAULT_COLUMNS)
311 results = []
312 for fp in footprints.getFootprints():
313 # only single peaked stars
314 if len(fp.getPeaks()) > 1:
315 continue
317 # Create a cutout of the image around the footprint
318 refCenter = tuple(fp.getCentroid())
319 star = measureStarOnStamp(image, refCenter, cutOutSize, apertureRadius, gain).toDataFrame()
320 if not star.empty:
321 results.append(star)
322 if not results:
323 return pd.DataFrame(columns=DEFAULT_COLUMNS)
324 df = pd.concat([sf for sf in results], ignore_index=True)
325 return df
328def measureStarOnStamp(
329 stamp: np.ndarray,
330 refCenter: tuple[float, float],
331 cutOutSize: int,
332 apertureRadius: int,
333 gain: float = 1.0,
334) -> StarMeasurement:
335 """
336 Measure a star on a single stamp: background subtraction, shape, centroid,
337 photometry.
339 Parameters
340 ----------
341 stamp : `np.ndarray`
342 Full stamp array.
343 refCenter : `tuple[float, float]`
344 Reference (x, y) pixel position for the cutout center.
345 cutOutSize : `int`
346 Size of the cutout in pixels.
347 apertureRadius : `int`
348 Aperture radius in pixels for photometry.
349 gain : `float`
350 Detector gain (e-/ADU).
352 Returns
353 -------
354 measurement : `StarMeasurement`
355 StarMeasurement object with populated fields (may be empty on failure).
356 """
357 cutout = getCutouts(stamp, refCenter, cutoutSize=cutOutSize)
358 data = cutout.data
360 if np.all(data == 0) | (not np.isfinite(data).all()):
361 return StarMeasurement()
363 # 1) Subtract the background
364 annulus = (apertureRadius * 1.0, apertureRadius * 2)
365 dataBkgSub, bkgStd = annulusBackgroundSubtraction(data, annulus)
367 # 2) Track the star across all stamps for this guider
368 star = runGalSim(dataBkgSub, gain=gain, bkgStd=bkgStd)
370 # 3) Make aperture photometry measurements
371 # Galsim flux is the normalization of the Gaussian, not w/ fixed aper.
372 star.runAperturePhotometry(dataBkgSub, apertureRadius, gain=gain, bkgStd=bkgStd)
374 # 4) Add centroid and shape in amplifier roi coordinates
375 star.xroi += cutout.xmin_original
376 star.yroi += cutout.ymin_original
377 return star
380def runGalSim(
381 imageArray: np.ndarray,
382 gain: float = 1.0,
383 bkgStd: float = 0.0,
384) -> StarMeasurement:
385 """
386 Measure star properties with GalSim adaptive moments.
388 Parameters
389 ----------
390 imageArray : `np.ndarray`
391 Background-subtracted image cutout.
392 gain : `float`
393 Detector gain (e-/ADU).
394 bkgStd : `float`
395 Background RMS per pixel.
397 Returns
398 -------
399 result : `StarMeasurement`
400 Resulting measurement (empty if measurement failed).
401 """
402 gsImg = galsim.Image(imageArray)
403 hsmRes = galsim.hsm.FindAdaptiveMom(gsImg, strict=False)
404 success = hsmRes.error_message == ""
406 if not success:
407 result = StarMeasurement()
408 else:
409 xCentroid = hsmRes.moments_centroid.x
410 yCentroid = hsmRes.moments_centroid.y
411 flux = hsmRes.moments_amp
412 sigma = hsmRes.moments_sigma
413 e1 = hsmRes.observed_shape.e1
414 e2 = hsmRes.observed_shape.e2
415 fwhm = 2.355 * sigma
417 # Calculate errors using GalSim's error estimation
418 xErr, yErr = calcGalsimError(imageArray, hsmRes, gain=gain, bkgStd=bkgStd, correctForGain=True)
420 # Calculate SNR and flux error
421 ellipticity = np.sqrt(e1**2 + e2**2)
422 nEff = 2 * np.pi * sigma**2 * np.sqrt(1 - ellipticity**2)
423 shotNoise = np.sqrt(nEff * bkgStd**2)
424 fluxErr = np.sqrt(flux / gain + shotNoise**2)
425 snr = flux / (shotNoise + 1e-9) if shotNoise > 0 else 0.0
427 # Calculate second moments
428 ixx = sigma**2 * (1 + e1)
429 iyy = sigma**2 * (1 - e1)
430 ixy = sigma**2 * e2
432 result = StarMeasurement(
433 xroi=xCentroid,
434 yroi=yCentroid,
435 xerr=xErr,
436 yerr=yErr,
437 e1=e1,
438 e2=e2,
439 ixx=ixx,
440 iyy=iyy,
441 ixy=ixy,
442 fwhm=fwhm,
443 flux=flux,
444 flux_err=fluxErr,
445 snr=snr,
446 )
447 return result
450def calcGalsimError(
451 imageArray: np.ndarray,
452 shape: galsim.hsm.ShapeData,
453 gain: float = 1.0,
454 bkgStd: float = 0.0,
455 correctForGain: bool = False,
456) -> tuple[float, float]:
457 """
458 Estimate centroid errors from GalSim HSMShapeData.
460 Parameters
461 ----------
462 imageArray : `np.ndarray`
463 Image cutout used for measurement.
464 shape : `galsim.hsm`
465 GalSim HSM shape data result object.
466 gain : `float`
467 Detector gain (e-/ADU), ignored if `correctForGain` is False.
468 bkgStd : `float`
469 Background RMS per pixel.
470 correctForGain : `bool`
471 Whether to include gain-dependent weighting.
473 Returns
474 -------
475 xerr : `float`
476 Estimated x centroid uncertainty (pixels).
477 yerr : `float`
478 Estimated y centroid uncertainty (pixels).
479 """
480 if not shape or shape.error_message != "":
481 return 0.0, 0.0
483 x0 = shape.moments_centroid.x
484 y0 = shape.moments_centroid.y
485 sigma = shape.moments_sigma
486 e1 = shape.observed_shape.e1
487 e2 = shape.observed_shape.e2
488 flux = shape.moments_amp
490 kernel = makeEllipticalGaussianStar(
491 shape=(imageArray.shape[0], imageArray.shape[1]),
492 e1=e1,
493 e2=e2,
494 flux=1,
495 sigma=sigma,
496 center=(x0, y0),
497 )
499 weight = np.ones_like(imageArray) / (bkgStd**2 + 1e-9)
500 if correctForGain:
501 weight = np.ones_like(imageArray) / (bkgStd**2 + np.abs(flux * kernel / gain))
503 mask = weight == 0.0
504 data = imageArray.copy()
505 if np.any(mask):
506 kernelMasked = kernel.copy()
507 data[mask] = kernelMasked[mask] * np.sum(data[~mask]) / np.sum(kernelMasked[~mask])
509 u, v = np.meshgrid(np.arange(imageArray.shape[1]) - x0, np.arange(imageArray.shape[0]) - y0)
510 usq = u**2
511 vsq = v**2
512 WI = kernel * data
513 M00 = np.nansum(WI)
514 WV = (kernel**2).astype(float)
515 WV[~mask] /= weight[~mask]
516 WV[mask] /= np.median(weight[~mask])
517 WV = WV / float(M00**2)
519 varM10 = 4 * np.sum(WV * usq)
520 varM01 = 4 * np.sum(WV * vsq)
521 xerr = np.sqrt(varM10)
522 yerr = np.sqrt(varM01)
523 return xerr, yerr
526def makeEllipticalGaussianStar(
527 shape: tuple[int, int],
528 flux: float,
529 sigma: float,
530 e1: float,
531 e2: float,
532 center: tuple[float, float],
533) -> np.ndarray:
534 """
535 Create an elliptical 2D Gaussian star with specified parameters.
537 Parameters
538 ----------
539 shape : `tuple[int, int]`
540 (ny, nx) output array shape.
541 flux : `float`
542 Total flux (normalization).
543 sigma : `float`
544 Gaussian sigma (pixels).
545 e1 : `float`
546 Ellipticity component e1.
547 e2 : `float`
548 Ellipticity component e2.
549 center : `tuple[float, float]`
550 (x0, y0) centroid position in pixels.
552 Returns
553 -------
554 image : `np.ndarray`
555 Generated model image.
556 """
557 y, x = np.indices(shape)
558 x0, y0 = center
559 u = x - x0
560 v = y - y0
562 # Second-moment matrix elements
563 ixx = sigma**2 * (1 + e1)
564 iyy = sigma**2 * (1 - e1)
565 ixy = sigma**2 * e2
567 # Inverse covariance matrix
568 det = ixx * iyy - ixy**2
569 invIxx = iyy / det
570 invIyy = ixx / det
571 invIxy = -ixy / det
573 # Quadratic form: u^2 * invIxx + v^2 * invIyy + 2uv * invIxy
574 r2 = invIxx * u**2 + invIyy * v**2 + 2 * invIxy * u * v
576 e = np.sqrt(e1**2 + e2**2)
577 norm = flux / (2 * np.pi * sigma**2 * np.sqrt(1 - e**2))
578 image = norm * np.exp(-0.5 * r2)
579 return image
582def buildReferenceCatalog(
583 guiderData: GuiderData,
584 log: logging.Logger,
585 config: GuiderStarTrackerConfig = GuiderStarTrackerConfig(),
586) -> pd.DataFrame:
587 """
588 Build a reference star catalog from each guider's coadded stamp.
590 Parameters
591 ----------
592 guiderData : `GuiderData`
593 Guider dataset containing stamps and metadata.
594 log : `logging.Logger`
595 Logger for warnings and diagnostics.
596 config : `GuiderStarTrackerConfig`
597 Star tracker configuration.
599 Returns
600 -------
601 refCatalog : `pd.DataFrame`
602 Concatenated reference catalog of brightest stars per guider.
603 """
604 expId = guiderData.expid
605 minSnr = config.minSnr
606 gain = config.gain
607 cutOutSize = config.cutOutSize
609 tableList = []
610 for guiderName in guiderData.guiderNames:
611 pixelScale = guiderData.getWcs(guiderName).getPixelScale().asArcseconds()
612 apertureRadius = int(config.aperSizeArcsec / pixelScale)
614 array = guiderData.getStampArrayCoadd(guiderName)
615 # array = np.where(array < 0, 0, array) # Ensure no negative values
616 array = array - np.nanmin(array) # Ensure no negative values
617 sources = runSourceDetection(
618 array,
619 threshold=minSnr,
620 apertureRadius=apertureRadius,
621 cutOutSize=cutOutSize,
622 gain=gain,
623 )
624 if sources.empty:
625 log.warning(f"No sources detected in `buildReferenceCatalog`" f"for {guiderName} in {expId}. ")
626 continue
628 sources.sort_values(by=["snr"], ascending=False, inplace=True)
629 sources.reset_index(drop=True, inplace=True)
631 detNum = guiderData.getGuiderDetNum(guiderName)
632 sources["detector"] = guiderName
633 sources["detid"] = detNum
634 sources["starid"] = detNum * 100
635 tableList.append(sources)
637 if len(tableList) == 0:
638 log.warning(f"buildReferenceCatalog failed - no stars detected in any guider for {expId}.")
639 return makeBlankCatalog()
641 refCatalog = pd.concat(tableList, ignore_index=True)
642 return refCatalog
645def getCutouts(imageArray: np.ndarray, refCenter: tuple[float, float], cutoutSize: int = 25) -> Cutout2D:
646 """
647 Get a cutout at the reference position from an image array.
649 Parameters
650 ----------
651 imageArray : `np.ndarray`
652 Full image array.
653 refCenter : `tuple[float, float]`
654 (x, y) center for the cutout in pixels.
655 cutoutSize : `int`
656 Size (pixels) of the square cutout.
658 Returns
659 -------
660 cutout : `Cutout2D`
661 Astropy Cutout2D object.
662 """
663 refX, refY = refCenter
664 return Cutout2D(imageArray, (refX, refY), size=cutoutSize, mode="partial", fill_value=np.nan)