Coverage for python / lsst / summit / extras / fastStarTrackerAnalysis.py: 0%
276 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:36 +0000
1# This file is part of summit_extras.
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/>.
22import glob
23import os
24import typing
25from dataclasses import dataclass
27import galsim
28import matplotlib
29import matplotlib.pyplot as plt
30import numpy as np
31from matplotlib.collections import PatchCollection
32from mpl_toolkits.axes_grid1 import make_axes_locatable
34import lsst.afw.image as afwImage
35import lsst.afw.math as afwMath
36import lsst.geom as geom
37from lsst.summit.utils.starTracker import (
38 dayObsSeqNumFrameNumFromFilename,
39 fastCam,
40 getRawDataDirForDayObs,
41 isStreamingModeFile,
42 openFile,
43)
44from lsst.summit.utils.utils import bboxToMatplotlibRectanle, detectObjectsInExp, getBboxAround, getSite
45from lsst.utils.iteration import ensure_iterable
47__all__ = (
48 "getStreamingSequences",
49 "getFlux",
50 "getBackgroundLevel",
51 "countOverThresholdPixels",
52 "sortSourcesByFlux",
53 "findFastStarTrackerImageSources",
54 "checkResultConsistency",
55 "plotSourceMovement",
56 "plotSource",
57 "plotSourcesOnImage",
58 "Source",
59 "NanSource",
60)
63@dataclass(slots=True)
64class Source:
65 """A dataclass for FastStarTracker analysis results."""
67 dayObs: int # mandatory attribute - the dayObs
68 seqNum: int # mandatory attribute - the seqNum
69 frameNum: int # mandatory attribute - the sub-sequence number, the position in the sequence
71 # raw numbers
72 centroidX: float = np.nan # in image coordinates
73 centroidY: float = np.nan # in image coordinates
74 rawFlux: float = np.nan
75 nPix: int | float = np.nan
76 bbox: geom.Box2I | None = None
77 cutout: np.ndarray | None = None
78 localCentroidX: float = np.nan # in cutout coordinates
79 localCentroidY: float = np.nan # in cutout coordinates
81 # numbers from the hsm moments fit
82 hsmFittedFlux: float = np.nan
83 hsmCentroidX: float = np.nan
84 hsmCentroidY: float = np.nan
85 moments: galsim.hsm.ShapeData | None = None # keep the full fit even though we pull some things out too
87 imageBackground: float = np.nan
88 imageStddev: float = np.nan
89 nSourcesInImage: int | float = np.nan
90 parentImageWidth: int | float = np.nan
91 parentImageHeight: int | float = np.nan
92 expTime: float = np.nan
94 def __repr__(self):
95 """Print everything except the full details of the moments."""
96 retStr = ""
97 for itemName in self.__slots__:
98 v = getattr(self, itemName)
99 if isinstance(v, int): # print ints as ints
100 retStr += f"{itemName} = {v}\n"
101 elif isinstance(v, float): # but round floats at 3dp
102 retStr += f"{itemName} = {v:.3f}\n"
103 elif itemName == "moments": # and don't spam the full moments
104 retStr += f"moments = {type(v)}\n"
105 elif itemName == "bbox": # and don't spam the full moments
106 retStr += f"bbox = lsst.geom.{repr(v)}\n"
107 elif itemName == "cutout": # and don't spam the full moments
108 if v is None:
109 retStr += "cutout = None\n"
110 else:
111 retStr += f"cutout = {type(v)}\n"
112 return retStr
115class NanSource:
116 def __getattribute__(self, name: str):
117 return np.nan
120def getStreamingSequences(dayObs: int) -> dict[int, list[str]]:
121 """Get the streaming sequences for a dayObs.
123 Note that this will need rewriting very soon once the way the data is
124 organised on disk is changed.
126 Parameters
127 ----------
128 dayObs : `int`
129 The dayObs.
131 Returns
132 -------
133 sequences : `dict` [`int`, `list`]
134 The streaming sequences in a dict, keyed by sequence number, with each
135 value being a list of the files in that sequence.
136 """
137 site = getSite()
138 if site in ["rubin-devl", "staff-rsp"]:
139 rootDataPath = "/sdf/data/rubin/offline/s3-backup/lfa/"
140 elif site == "summit":
141 rootDataPath = "/project"
142 else:
143 raise ValueError(f"Finding StarTracker data isn't supported at {site}")
145 dataDir = getRawDataDirForDayObs(rootDataPath, fastCam, dayObs)
146 files = glob.glob(os.path.join(dataDir, "*.fits"))
147 regularFiles = [f for f in files if not isStreamingModeFile(f)]
148 streamingFiles = [f for f in files if isStreamingModeFile(f)]
149 print(f"Found {len(regularFiles)} regular files on dayObs {dayObs}")
151 data = {}
152 if dayObs < 20240311:
153 # after this is when we changed the data layout on disk for streaming
154 # mode data in the GenericCamera
155 for filename in sorted(streamingFiles):
156 basename = os.path.basename(filename)
157 seqNum = int(basename.split("_")[3])
158 if seqNum not in data:
159 data[seqNum] = [filename]
160 else:
161 data[seqNum].append(filename)
162 else:
163 # dirNames here doesn't contain the full path, it's just the individual
164 # directory name and needs joining with dataDir for the full path
165 dirNames = sorted(d for d in os.listdir(dataDir) if os.path.isdir(os.path.join(dataDir, d)))
166 for d in dirNames:
167 files = sorted(glob.glob(os.path.join(dataDir, d, "*.fits")))
168 seqNum = int(d.split("_")[3])
169 data[seqNum] = files
171 print(f"Found {len(data)} streaming sequences on dayObs {dayObs}:")
172 for seqNum, files in data.items():
173 print(f"seqNum {seqNum} with {len(files)} frames")
175 return data
178def getFlux(cutout: np.ndarray[float], backgroundLevel: float = 0) -> float:
179 """Get the flux inside a cutout, subtracting the image-background.
181 Here the flux is simply summed, and if the image background level is
182 supplied, it is subtracted off, assuming it is constant over the cutout. A
183 more accurate(?) flux is obtained by the hsm model fit.
185 Parameters
186 ----------
187 cutout : `np.array`
188 The cutout as a raw array.
189 backgroundLevel : `float`, optional
190 If supplied, this is subtracted as a constant background level.
192 Returns
193 -------
194 flux : `float`
195 The flux of the source in the cutout.
196 """
197 rawFlux = np.sum(cutout)
198 if not backgroundLevel:
199 return rawFlux
201 return rawFlux - (cutout.size * backgroundLevel)
204def getBackgroundLevel(exp: afwImage.Exposure, nSigma: float = 3) -> tuple[float, float]:
205 """Calculate the clipped image mean and stddev of an exposure.
207 Testing shows on images like this, 2 rounds of sigma clipping is more than
208 enough so this is left fixed here.
210 Parameters
211 ----------
212 exp : `lsst.afw.image.Exposure`
213 The exposure.
214 nSigma : `float`, optional
215 The number of sigma to clip to for the background estimation.
217 Returns
218 -------
219 mean : `float`
220 The clipped mean, as an estimate of the background level
221 stddev : `float`
222 The clipped standard deviation, as an estimate of the background noise.
223 """
224 sctrl = afwMath.StatisticsControl()
225 sctrl.setNumSigmaClip(nSigma)
226 sctrl.setNumIter(2) # this is always plenty here
227 statTypes = afwMath.MEANCLIP | afwMath.STDEVCLIP
228 stats = afwMath.makeStatistics(exp.maskedImage, statTypes, sctrl)
229 std, _ = stats.getResult(afwMath.STDEVCLIP)
230 mean, _ = stats.getResult(afwMath.MEANCLIP)
231 return mean, std
234def countOverThresholdPixels(cutout: np.ndarray, bgMean: float, bgStd: float, nSigma: float = 15) -> int:
235 """Get the number of pixels in the cutout which are 'in the source'.
237 From the one image I've looked at so far, the drop-off is quite slow
238 probably due to some combination of focus, plate scale, star brightness,
239 pointing quality etc, so the default nSigma is 15 here as that looked about
240 right when I plotted it by eye.
242 Parameters
243 ----------
244 cutout : `np.array`
245 The cutout to measure.
246 bgMean : `float`
247 The background level.
248 bgStd : `float`
249 The clipped standard deviation in the image.
250 nSigma : `float`, optional
251 The number of sigma above background at which to count pixels as being
252 over threshold.
254 Returns
255 -------
256 nPix : `int`
257 The number of pixels above threshold.
258 """
259 inds = np.where(cutout > (bgMean + 0 * bgStd))
260 return len(inds[0])
263def sortSourcesByFlux(sources: list[Source], reverse: bool = False) -> list[Source]:
264 """Sort the sources by flux, returning the brightest first.
266 Parameters
267 ----------
268 sources : `list` of
269 `lsst.summit.extras.fastStarTrackerAnalysis.Source`
270 The list of sources to sort.
271 reverse : `bool`, optional
272 Return the brightest at the start of the list if ``reverse`` is
273 ``False``, or the brightest last if ``reverse`` is ``True``.
275 Returns
276 -------
277 sources : `list` of
278 `lsst.summit.extras.fastStarTrackerAnalysis.Source`
279 The sources, sorted by flux.
280 """
281 # invert reverse because we want brightest first by default, but want the
282 # reverse arg to still behave as one would expect
283 return sorted(sources, key=lambda s: s.rawFlux, reverse=not reverse)
286def findFastStarTrackerImageSources(
287 filename: str, boxSize: int, attachCutouts: bool = True
288) -> list[Source | NanSource]:
289 """Analyze a single FastStarTracker image.
291 Parameters
292 ----------
293 filename : `str`
294 The full name and path of the file.
295 boxSize : `int`
296 The size of the box to put around each source for measurement.
297 attachCutouts : `bool`, optional
298 Attach the cutouts to the ``Source`` objects? Useful for
299 debug/plotting but adds memory usage.
301 Returns
302 -------
303 sources : `list` of
304 `lsst.summit.extras.fastStarTrackerAnalysis.Source`
305 The sources in the image, sorted by rawFlux.
306 """
307 exp = openFile(filename)
308 # if the upstream exposure reading code hasn't set the
309 # visitInfo.exposureTime then this will return nan, as desired
310 expTime = exp.visitInfo.exposureTime
311 footprintSet = detectObjectsInExp(exp)
312 footprints = footprintSet.getFootprints()
313 bgMean, bgStd = getBackgroundLevel(exp)
315 dayObs, seqNum, frameNum = dayObsSeqNumFrameNumFromFilename(filename)
317 sources = []
318 if len(footprints) == 0:
319 sources = [NanSource()]
320 return sources
322 for footprint in footprints:
323 source = Source(dayObs=dayObs, seqNum=seqNum, frameNum=frameNum)
324 source.expTime = expTime
325 source.nSourcesInImage = len(footprints)
326 source.parentImageWidth, source.parentImageHeight = exp.getDimensions()
328 centroid = footprint.getCentroid()
329 bbox = getBboxAround(centroid, boxSize, exp)
330 source.bbox = bbox
331 cutout = exp.image[bbox].array
332 if attachCutouts:
333 source.cutout = cutout
334 source.centroidX = centroid[0]
335 source.centroidY = centroid[1]
336 source.rawFlux = getFlux(cutout, bgMean)
337 source.imageBackground = bgMean
338 source.imageStddev = bgStd
339 source.nPix = countOverThresholdPixels(cutout, bgMean, bgStd)
341 moments = galsim.hsm.FindAdaptiveMom(galsim.Image(cutout))
342 source.moments = moments
343 source.hsmFittedFlux = moments.moments_amp
344 source.hsmCentroidX = moments.moments_centroid.x + bbox.minX - 1
345 source.hsmCentroidY = moments.moments_centroid.y + bbox.minY - 1
346 source.localCentroidX = moments.moments_centroid.x - 1
347 source.localCentroidY = moments.moments_centroid.y - 1
348 sources.append(source)
349 return sortSourcesByFlux(sources)
352def checkResultConsistency(
353 results: dict[int, list[Source]],
354 maxAllowableShift: float = 5,
355 silent: bool = False,
356) -> bool:
357 """Check if a set of results are self-consistent.
359 Check the number of detected sources are the same in each image, that no
360 sources have been removed from each image's source list, and that all the
361 input images were the same size (because we read out sub frames, and
362 analyzing these with full frame data invalidates the centroid coordinates).
364 Also displays the maximum (x, y) movements between adjacent exposures, and
365 the mean and stddev of the main source's flux.
367 Parameters
368 ----------
369 results : `dict` of `list` of
370 `lsst.summit.extras.fastStarTrackerAnalysis.Source`
371 A dict, keyed by sequence number, with each value being a list of the
372 sources found in the image, e.g. as returned by
373 ``findFastStarTrackerImageSources()``.
374 maxAllowableShift : `float`
375 The biggest centroid shift between adjacent images allowable before
376 something is considered to have gone wrong.
377 silent : `bool`, optional
378 Print some useful checks and measurements if ``False``, otherwise just
379 return whether the results appear nominally OK silently (for use when
380 being called by other code rather than users).
382 Returns
383 -------
384 consistent : `bool`
385 Are the results nominally consistent?
386 """
387 if isinstance(results, typing.ValuesView): # in case we're passed a .values()
388 results = list(results)
390 sourceCounts = set([len(sourceSet) for sourceSet in results])
391 if sourceCounts == {0}: # none of the images contain any detections
392 if not silent:
393 print("No images contain any sources. Results are technically consistent, but also useless.")
394 # this is technically consistent, so return True, but any downstream
395 # code which tries to make plots with these will fail, of course.
396 return True
398 if 0 in ([len(sourceSet) for sourceSet in results]):
399 if not silent:
400 print(
401 "Some results contain no sources. Results are therefore fundamentally inconsistent"
402 " and other checks cannot be run"
403 )
404 return False
406 consistent = True
407 toPrint = []
408 nSources = set([sourceSet[0].nSourcesInImage for sourceSet in results])
409 if len(nSources) != 1:
410 toPrint.append(f"❌ Images contain a variable number of sources: {nSources}")
411 consistent = False
412 else:
413 n = nSources.pop()
414 toPrint.append(f"✅ All images contain the same nominal number of sources at detection stage: {n}")
416 nSourcesCounted = set([len(sourceSet) for sourceSet in results])
417 if len(nSourcesCounted) != 1:
418 toPrint.append(
419 f"❌ Number of actual sources in each sourceSet varies, got: {nSourcesCounted}."
420 " If some were manually removed you can ignore this"
421 )
422 consistent = False
423 else:
424 n = nSourcesCounted.pop()
425 toPrint.append(f"✅ All results contain the same number of actual sources per image: {n}")
427 widths = set([sourceSet[0].parentImageWidth for sourceSet in results])
428 heights = set([sourceSet[0].parentImageHeight for sourceSet in results])
429 if len(widths) != 1 or len(heights) != 1:
430 toPrint.append(f"❌ Input images were of variable dimenions! {widths=}, {heights=}")
431 consistent = False
432 else:
433 toPrint.append("✅ All input images were of the same dimensions")
435 if len(results) > 1: # can't np.diff an array of length 1 so these are not useful/defined
436 # now the basic checks have passed, do some sanity checks on the
437 # maximum deltas for the primary sources
438 sources = [sourceSet[0] for sourceSet in results]
439 dx = np.diff([s.centroidX for s in sources])
440 dy = np.diff([s.centroidY for s in sources])
441 maxMovementX = np.max(dx)
442 maxMovementY = np.max(dy)
443 happyOrSad = "✅"
444 if max(maxMovementX, maxMovementY) > maxAllowableShift:
445 consistent = False
446 happyOrSad = "❌"
448 toPrint.append(
449 f"{happyOrSad} Maximum centroid movement of brightest object between images in (x, y)"
450 f" = ({maxMovementX:.2f}, {maxMovementY:.2f}) pix"
451 )
453 fluxStd = np.nanstd([s.rawFlux for s in sources])
454 fluxMean = np.nanmean([s.rawFlux for s in sources])
455 toPrint.append(f"Mean and stddev of flux from brightest object = {fluxMean:.1f} ± {fluxStd:.1f} ADU")
457 if not silent:
458 for line in toPrint:
459 print(line)
461 return consistent
464def plotSourceMovement(
465 results: dict[int, list[Source]],
466 sourceIndex: int = 0,
467 allowInconsistent: bool = False,
468) -> list[matplotlib.figure.Figure]:
469 """Plot the centroid movements and fluxes etc for a set of results.
471 By default the brightest source in each image is plotted, but this can be
472 changed by setting ``sourceIndex`` to values greater than 0 to move through
473 the list of sources in each image.
475 Parameters
476 ----------
477 results : `dict` of `list` of
478 `lsst.summit.extras.fastStarTrackerAnalysis.Source`
479 A dict, keyed by sequence number, with each value being a list of the
480 sources found in the image, e.g. as returned by
481 ``findFastStarTrackerImageSources()``.
482 sourceIndex : `int`, optional
483 If there is more than one source in every image, which source number
484 should the plot be made for? Defaults to zero, which is the brightest
485 source by default.
486 allowInconsistent : `bool`, optional
487 Make the plots even if the input results appear to be inconsistent?
489 Returns
490 -------
491 figs : `list` of `matplotlib.figure.Figure`
492 The figures. The first is the source's flux and x, y movement over the
493 image sequence, and the second is a scatter plot of the x and y, with
494 the color showing the position in the sequence.
495 """
496 opts = {
497 "marker": "o",
498 "markersize": 6,
499 "linestyle": "-",
500 }
502 consistent = checkResultConsistency(results.values(), silent=True)
503 if not consistent and not allowInconsistent:
504 checkResultConsistency(results.values(), silent=False) # print the problem if we're raising
505 raise ValueError("The sources were found to be inconsistent and allowInconsistent=False")
507 sourceDict = {k: v[sourceIndex] for k, v in results.items()}
508 frameNums = [s.frameNum for s in sourceDict.values()]
509 sources = list(sourceDict.values())
511 allDayObs = set(s.dayObs for s in sources)
512 allSeqNums = set(s.seqNum for s in sources)
513 if len(allDayObs) > 1 or len(allSeqNums) > 1:
514 raise ValueError(
515 "The sources are from multiple days or sequences, found"
516 f" {allDayObs} dayObs and {allSeqNums} seqNum values."
517 )
518 dayObs = allDayObs.pop()
519 seqNum = allSeqNums.pop()
520 startFrame = min(frameNums)
521 endFrame = max(frameNums)
523 title = f"dayObs {dayObs}, seqNum {seqNum}, frames {startFrame}-{endFrame}"
525 axisLabelSize = 18
527 figs = []
528 fig = plt.figure(figsize=(10, 16))
529 ax1, ax2, ax3 = fig.subplots(3, sharex=True)
530 fig.subplots_adjust(hspace=0)
532 ax1.plot(frameNums, [s.rawFlux for s in sources], label="Raw Flux", **opts)
533 ax1.plot(frameNums, [s.hsmFittedFlux for s in sources], label="Fitted Flux", **opts)
534 ax1.set_ylabel("Flux (ADU)", size=axisLabelSize)
535 ax1.set_title(title)
536 ax1.legend()
538 ax2.plot(frameNums, [s.centroidX for s in sources], label="Raw centroid x", **opts)
539 ax2.plot(
540 frameNums,
541 [s.hsmCentroidX for s in sources],
542 label="Fitted centroid x",
543 **opts,
544 )
545 ax2.set_ylabel("x-centroid (pixels)", size=axisLabelSize)
546 ax2.legend()
548 ax3.plot(frameNums, [s.centroidY for s in sources], label="Raw centroid y", **opts)
549 ax3.plot(
550 frameNums,
551 [s.hsmCentroidY for s in sources],
552 label="Fitted centroid y",
553 **opts,
554 )
555 ax3.set_ylabel("y-centroid (pixels)", size=axisLabelSize)
556 ax3.set_xlabel("Frame number", size=axisLabelSize)
557 ax3.legend()
559 figs.append(fig)
561 fig = plt.figure(figsize=(10, 10))
562 ax4 = fig.subplots(1)
564 colors = np.arange(len(sources))
565 # gnuplot2 has a nice balance of nothing white, and having an intuitive
566 # progression of colours so the eye can pick out trends on the point cloud.
567 axRef = ax4.scatter(
568 [s.centroidX for s in sources],
569 [s.centroidY for s in sources],
570 c=colors,
571 cmap="gnuplot2",
572 )
573 ax4.set_xlabel("x-centroid (pixels)", size=axisLabelSize)
574 ax4.set_ylabel("y-centroid (pixels)", size=axisLabelSize)
575 ax4.set_aspect("equal", "box")
576 # move the colorbar
577 divider = make_axes_locatable(ax4)
578 cax = divider.append_axes("right", size="5%", pad=0.05)
579 cbar = plt.colorbar(axRef, cax=cax)
580 ax4.set_title(title)
581 cbar.set_label("Frame number in series", size=axisLabelSize * 0.75)
582 figs.append(fig)
584 return figs
587# -------------- plotting tools
590def plotSourcesOnImage(
591 parentFilename: str,
592 sources: Source | list[Source],
593) -> None:
594 """Plot one of more source on top of an image.
596 Parameters
597 ----------
598 parentFilename : `str`
599 The full path to the parent (.tif) file.
600 sources : `list` of
601 `lsst.summit.extras.fastStarTrackerAnalysis.Source` or
602 `lsst.summit.extras.fastStarTrackerAnalysis.Source`
603 The sources found in the image.
604 """
605 exp = openFile(parentFilename)
606 data = exp.image.array
608 fig = plt.figure(figsize=(16, 8))
609 ax = fig.subplots(1)
611 plt.imshow(data, interpolation="None", origin="lower")
613 sources = ensure_iterable(sources)
614 patches = []
615 for source in sources:
616 ax.scatter(source.centroidX, source.centroidY, color="red", marker="x") # mark the centroid
617 patch = bboxToMatplotlibRectanle(source.bbox)
618 patches.append(patch)
620 # move the colorbar
621 divider = make_axes_locatable(ax)
622 cax = divider.append_axes("right", size="5%", pad=0.05)
623 plt.colorbar(cax=cax)
625 # plot the bboxes on top
626 pc = PatchCollection(patches, edgecolor="r", facecolor="none")
627 ax.add_collection(pc)
629 plt.tight_layout()
632def plotSource(source: Source) -> None:
633 """Plot a single source.
635 Parameters
636 ----------
637 source : `lsst.summit.extras.fastStarTrackerAnalysis.Source`
638 The source to plot.
639 """
640 if source.cutout is None:
641 raise RuntimeError(
642 "Can only plot sources with attached cutouts. Either set attachCutouts=True "
643 "in findFastStarTrackerImageSources() or try using plotSourcesOnImage() instead"
644 )
646 fig = plt.figure(figsize=(16, 8))
647 ax = fig.subplots(1)
649 plt.imshow(source.cutout, interpolation="None", origin="lower") # plot the image
650 ax.scatter(source.localCentroidX, source.localCentroidY, color="red", marker="x", s=200) # mark centroid
652 # move the colorbar
653 divider = make_axes_locatable(ax)
654 cax = divider.append_axes("right", size="5%", pad=0.05)
655 plt.colorbar(cax=cax)
657 plt.tight_layout()