Coverage for python/lsst/summit/extras/focusAnalysis.py: 14%
284 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-10 10:35 +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/>.
22from dataclasses import dataclass
24import matplotlib
25import matplotlib.pyplot as plt
26import numpy as np
27from matplotlib import gridspec
28from matplotlib.colors import LogNorm
29from matplotlib.patches import Arrow, Circle, Rectangle
30from scipy.linalg import norm
31from scipy.optimize import curve_fit
33import lsst.afw.image as afwImage
34import lsst.geom as geom
35from lsst.atmospec.utils import isDispersedExp
36from lsst.pipe.tasks.quickFrameMeasurement import QuickFrameMeasurementTask, QuickFrameMeasurementTaskConfig
37from lsst.summit.utils import ImageExaminer
39# TODO: change these back to local .imports
40from lsst.summit.utils.bestEffort import BestEffortIsr
41from lsst.summit.utils.butlerUtils import getExpRecordFromDataId, makeDefaultLatissButler
42from lsst.summit.utils.utils import FWHMTOSIGMA, SIGMATOFWHM
44__all__ = ["SpectralFocusAnalyzer", "NonSpectralFocusAnalyzer"]
47@dataclass
48class FitResult:
49 """A single 1D Gaussian fit result.
51 Attributes
52 ----------
53 amp : `float`
54 The fitted amplitude.
55 mean : `float`
56 The fitted mean (peak position).
57 sigma : `float`
58 The fitted standard deviation.
59 """
61 amp: float
62 mean: float
63 sigma: float
66def getFocusFromExposure(exp: afwImage.Exposure) -> float:
67 """Get the focus value from an exposure.
69 This was previously accessed via raw metadata but now lives inside the
70 visitInfo.
72 Parameters
73 ----------
74 exp : `lsst.afw.image.Exposure`
75 The exposure.
77 Returns
78 -------
79 focus : `float`
80 The focus value.
82 """
83 return float(exp.visitInfo.focusZ)
86class SpectralFocusAnalyzer:
87 """Analyze a focus sweep taken on dispersed (spectral) LATISS data.
89 Takes slices across the spectrum for each image, fits a Gaussian to
90 each slice to measure its width, and then fits a parabola to the
91 widths as a function of hexapod focus position to determine the
92 best-focus hexapod position. The number of slices and their
93 distances from the main star can be customized via
94 `setSpectrumBoxOffsets`.
96 Parameters
97 ----------
98 embargo : `bool`, optional
99 If `True`, use the embargo butler repo rather than the default
100 LATISS repo.
102 Examples
103 --------
104 >>> %matplotlib inline
105 >>> dayObs = 20210101
106 >>> seqNums = [100, 101, 102, 103, 104]
107 >>> focusAnalyzer = SpectralFocusAnalyzer()
108 >>> focusAnalyzer.setSpectrumBoxOffsets([500, 750, 1000, 1250])
109 >>> focusAnalyzer.getFocusData(dayObs, seqNums, doDisplay=True)
110 >>> focusAnalyzer.fitDataAndPlot()
112 `run` performs the two steps above in one call.
113 """
115 def __init__(self, embargo: bool = False):
116 self.butler = makeDefaultLatissButler(embargo=embargo)
117 self._bestEffort = BestEffortIsr(embargo=embargo)
118 qfmTaskConfig = QuickFrameMeasurementTaskConfig()
119 self._quickMeasure = QuickFrameMeasurementTask(config=qfmTaskConfig)
121 self.spectrumHalfWidth = 100
122 self.spectrumBoxLength = 20
123 self._spectrumBoxOffsets = [882, 1170, 1467]
124 self._setColors(len(self._spectrumBoxOffsets))
126 def setSpectrumBoxOffsets(self, offsets: list[int]) -> None:
127 """Set the spectrum slice offsets used by the analyzer.
129 Parameters
130 ----------
131 offsets : `list` [`int`]
132 The y-distances from the main star's centroid at which to
133 slice the spectrum, in pixels.
134 """
135 self._spectrumBoxOffsets = offsets
136 self._setColors(len(offsets))
138 def getSpectrumBoxOffsets(self) -> list[int]:
139 """Return the current spectrum slice offsets.
141 Returns
142 -------
143 offsets : `list` [`int`]
144 The y-distances from the main star's centroid at which the
145 spectrum is sliced, in pixels.
146 """
147 return self._spectrumBoxOffsets
149 def _setColors(self, nPoints: int) -> None:
150 """Populate ``self.COLORS`` with ``nPoints`` rainbow colors.
152 Parameters
153 ----------
154 nPoints : `int`
155 Number of colors to sample from the ``rainbow`` colormap.
156 """
157 self.COLORS = matplotlib.colormaps["rainbow"](np.linspace(0, 1, nPoints))
159 def _getBboxes(self, centroid: list[float]) -> list[geom.Box2I]:
160 """Build the spectrum-slice bounding boxes for this centroid.
162 Parameters
163 ----------
164 centroid : `list` [`float`] or `tuple` [`float`, `float`]
165 The ``(x, y)`` centroid of the main star.
167 Returns
168 -------
169 bboxes : `list` [`lsst.geom.Box2I`]
170 One bounding box per configured spectrum offset.
171 """
172 x, y = centroid
173 bboxes = []
175 for offset in self._spectrumBoxOffsets:
176 bbox = geom.Box2I(
177 geom.Point2I(x - self.spectrumHalfWidth, y + offset),
178 geom.Point2I(x + self.spectrumHalfWidth, y + offset + self.spectrumBoxLength),
179 )
180 bboxes.append(bbox)
181 return bboxes
183 def _bboxToMplRectangle(self, bbox: geom.Box2I, colorNum: int) -> matplotlib.patches.Rectangle:
184 """Convert a Box2I to a matplotlib `Rectangle` patch.
186 Parameters
187 ----------
188 bbox : `lsst.geom.Box2I`
189 The bounding box to draw.
190 colorNum : `int`
191 Index into ``self.COLORS`` to pick the edge color.
193 Returns
194 -------
195 rectangle : `matplotlib.patches.Rectangle`
196 A non-filled rectangle patch matching ``bbox``.
197 """
198 xmin = bbox.getBeginX()
199 ymin = bbox.getBeginY()
200 xsize = bbox.getWidth()
201 ysize = bbox.getHeight()
202 rectangle = Rectangle(
203 (xmin, ymin), xsize, ysize, alpha=1, facecolor="none", lw=2, edgecolor=self.COLORS[colorNum]
204 )
205 return rectangle
207 @staticmethod
208 def gauss(x: float | np.ndarray, *pars: float) -> float | np.ndarray:
209 """Evaluate a 1D Gaussian.
211 Parameters
212 ----------
213 x : `float` or `numpy.ndarray`
214 The abscissa at which to evaluate the Gaussian.
215 *pars : `float`
216 The Gaussian parameters in the order
217 ``(amplitude, mean, sigma)``, as required by
218 `scipy.optimize.curve_fit`.
220 Returns
221 -------
222 y : `float` or `numpy.ndarray`
223 The Gaussian evaluated at ``x``.
224 """
225 amp, mean, sigma = pars
226 return amp * np.exp(-((x - mean) ** 2) / (2.0 * sigma**2))
228 def run(
229 self,
230 dayObs: int,
231 seqNums: list[int],
232 doDisplay: bool = False,
233 hideFit: bool = False,
234 hexapodZeroPoint: float = 0,
235 ) -> list[float]:
236 """Perform a focus sweep analysis for spectral data.
238 For each seqNum for the specified dayObs, take a slice through the
239 spectrum at y-offsets as specified by the offsets
240 (see get/setSpectrumBoxOffsets() for setting these) and fit a Gaussian
241 to the spectrum slice to measure its width.
243 For each offset distance, fit a parabola to the fitted spectral widths
244 and return the hexapod position at which the best focus was achieved
245 for each.
247 Parameters
248 ----------
249 dayObs : `int`
250 The dayObs to use.
251 seqNums : `list` [`int`]
252 The seqNums for the focus sweep to analyze.
253 doDisplay : `bool`, optional
254 If `True`, render the per-image diagnostic plots. Designed
255 to be used in a notebook with ``%matplotlib inline``.
256 hideFit : `bool`, optional
257 If `True`, do not overlay the parabolic fit on the width
258 plot (data points are still shown).
259 hexapodZeroPoint : `float`, optional
260 Zero point subtracted from the hexapod focus axis before
261 plotting and fitting.
263 Returns
264 -------
265 bestFits : `list` [`float`]
266 Best-focus hexapod position for each spectral slice.
267 """
268 self.getFocusData(dayObs, seqNums, doDisplay=doDisplay)
269 bestFits = self.fitDataAndPlot(hideFit=hideFit, hexapodZeroPoint=hexapodZeroPoint)
270 return bestFits
272 def getFocusData(self, dayObs: int, seqNums: list[int], doDisplay: bool = False) -> None:
273 """Collect the per-image fit data for a spectral focus sweep.
275 For each seqNum of the specified dayObs, take a slice through
276 the spectrum at the y-offsets configured via
277 `setSpectrumBoxOffsets` and fit a Gaussian to each slice to
278 measure its width. Results are stored in ``self.fitData`` for
279 later use by `fitDataAndPlot`.
281 Parameters
282 ----------
283 dayObs : `int`
284 The dayObs to use.
285 seqNums : `list` [`int`]
286 The seqNums for the focus sweep to analyze.
287 doDisplay : `bool`, optional
288 If `True`, render per-image diagnostic plots showing the
289 exposure with the slice boxes overlaid and the Gaussian
290 fits to each slice. Designed to be used in a notebook with
291 ``%matplotlib inline``.
293 Notes
294 -----
295 Performs the focus analysis per-image, holding the data in the
296 class. Call `fitDataAndPlot` after this to fit a parabola to
297 the focus data itself.
298 """
299 fitData: dict[int, dict] = {}
300 filters: set[str] = set()
301 objects: set[str] = set()
303 axes = None
304 for seqNum in seqNums:
305 fitData[seqNum] = {}
306 dataId = {"day_obs": dayObs, "seq_num": seqNum, "detector": 0}
307 exp = self._bestEffort.getExposure(dataId)
309 # sanity checking
310 filt = exp.filter.physicalLabel
311 expRecord = getExpRecordFromDataId(self.butler, dataId)
312 obj = expRecord.target_name
313 objects.add(obj)
314 filters.add(filt)
315 assert isDispersedExp(exp), f"Image is not dispersed! (filter = {filt})"
316 assert len(filters) == 1, "You accidentally mixed filters!"
317 assert len(objects) == 1, "You accidentally mixed objects!"
319 quickMeasResult = self._quickMeasure.run(exp)
320 centroid = quickMeasResult.brightestObjCentroid
321 spectrumSliceBboxes = self._getBboxes(centroid) # inside the loop due to centroid shifts
323 if doDisplay:
324 fig, axes = plt.subplots(1, 2, figsize=(18, 9))
325 exp.image.array[exp.image.array <= 0] = 0.001
326 axes[0].imshow(exp.image.array, norm=LogNorm(), origin="lower", cmap="gray_r")
327 plt.tight_layout()
328 arrowX = centroid[0] - 400
329 arrowY = centroid[1]
330 arrowDx, arrowDy = 300, 0
331 arrow = Arrow(arrowX, arrowY, arrowDx, arrowDy, width=200.0, color="red")
332 circle = Circle(centroid, radius=25, facecolor="none", edgecolor="red")
333 axes[0].add_patch(arrow)
334 axes[0].add_patch(circle)
335 for i, bbox in enumerate(spectrumSliceBboxes):
336 rect = self._bboxToMplRectangle(bbox, i)
337 axes[0].add_patch(rect)
339 for i, bbox in enumerate(spectrumSliceBboxes):
340 data1d = np.mean(exp[bbox].image.array, axis=0) # flatten
341 data1d -= np.median(data1d)
342 xs = np.arange(len(data1d))
344 # get rough estimates for fit
345 # can't use sigma from quickMeasResult due to SDSS shape
346 # failing on saturated starts, and fp.getShape() is weird
347 amp = np.max(data1d)
348 mean = np.argmax(data1d)
349 sigma = 20
350 p0 = amp, mean, sigma
352 try:
353 coeffs, var_matrix = curve_fit(self.gauss, xs, data1d, p0=p0)
354 except RuntimeError:
355 coeffs = (np.nan, np.nan, np.nan)
357 fitData[seqNum][i] = FitResult(amp=abs(coeffs[0]), mean=coeffs[1], sigma=abs(coeffs[2]))
358 if doDisplay:
359 assert axes is not None
360 axes[1].plot(xs, data1d, "x", c=self.COLORS[i])
361 highResX = np.linspace(0, len(data1d), 1000)
362 if not np.isnan(coeffs[0]):
363 axes[1].plot(highResX, self.gauss(highResX, *coeffs), "k-")
365 if doDisplay: # show all color boxes together
366 plt.title(f"Fits to seqNum {seqNum}")
367 plt.show()
369 focuserPosition = getFocusFromExposure(exp)
370 fitData[seqNum]["focus"] = focuserPosition
372 self.fitData = fitData
373 self.filter = filters.pop()
374 self.object = objects.pop()
376 return
378 def fitDataAndPlot(self, hideFit: bool = False, hexapodZeroPoint: float = 0) -> list[float]:
379 """Fit a parabola to each spectral slice and plot the result.
381 For each configured slice offset, fit a parabola to the fitted
382 Gaussian widths vs. hexapod focus and return the hexapod
383 position at which each parabola is minimized.
385 Parameters
386 ----------
387 hideFit : `bool`, optional
388 If `True`, do not overlay the parabolic fit on the width
389 plot (data points are still shown).
390 hexapodZeroPoint : `float`, optional
391 Zero point subtracted from the hexapod focus axis before
392 plotting and fitting.
394 Returns
395 -------
396 bestFits : `list` [`float`]
397 Best-focus hexapod position for each spectral slice.
398 """
399 data = self.fitData
400 filt = self.filter
401 obj = self.object
403 bestFits = []
405 titleFontSize = 18
406 legendFontSize = 12
407 labelFontSize = 14
409 pixelsPerArcsec = 10
410 sigmaToFwhm = 2.355
412 f, axes = plt.subplots(2, 1, figsize=[10, 12])
413 focusPositions = [data[k]["focus"] - hexapodZeroPoint for k in sorted(data.keys())]
414 fineXs = np.linspace(np.min(focusPositions), np.max(focusPositions), 101)
415 seqNums = sorted(data.keys())
417 nSpectrumSlices = len(data[list(data.keys())[0]]) - 1
418 pointsForLegend = [0.0 for offset in range(nSpectrumSlices)]
419 for spectrumSlice in range(nSpectrumSlices): # the blue/green/red slices through the spectrum
420 # for scatter plots, the color needs to be a single-row 2d array
421 thisColor = np.array([self.COLORS[spectrumSlice]])
423 amps = [data[seqNum][spectrumSlice].amp for seqNum in seqNums]
424 widths = [data[seqNum][spectrumSlice].sigma / pixelsPerArcsec * sigmaToFwhm for seqNum in seqNums]
426 pointsForLegend[spectrumSlice] = axes[0].scatter(focusPositions, amps, c=thisColor)
427 axes[0].set_xlabel("Focus position (mm)", fontsize=labelFontSize)
428 axes[0].set_ylabel("Height (ADU)", fontsize=labelFontSize)
430 axes[1].scatter(focusPositions, widths, c=thisColor)
431 axes[1].set_xlabel("Focus position (mm)", fontsize=labelFontSize)
432 axes[1].set_ylabel("FWHM (arcsec)", fontsize=labelFontSize)
434 quadFitPars = np.polyfit(focusPositions, widths, 2)
435 fitMin = -quadFitPars[1] / (2.0 * quadFitPars[0])
436 bestFits.append(fitMin)
437 if not hideFit:
438 axes[1].plot(fineXs, np.poly1d(quadFitPars)(fineXs), c=self.COLORS[spectrumSlice])
439 axes[1].axvline(fitMin, color=self.COLORS[spectrumSlice])
440 msg = f"Best focus offset = {np.round(fitMin, 2)}"
441 axes[1].text(
442 fitMin,
443 np.mean(widths),
444 msg,
445 horizontalalignment="right",
446 verticalalignment="center",
447 rotation=90,
448 color=self.COLORS[spectrumSlice],
449 fontsize=legendFontSize,
450 )
452 titleText = f"Focus curve for {obj} w/ {filt}"
453 plt.suptitle(titleText, fontsize=titleFontSize)
454 legendText = self._generateLegendText(nSpectrumSlices)
455 axes[0].legend(pointsForLegend, legendText, fontsize=legendFontSize)
456 axes[1].legend(pointsForLegend, legendText, fontsize=legendFontSize)
457 f.tight_layout(rect=(0, 0.03, 1, 0.95))
459 for i, bestFit in enumerate(bestFits):
460 print(f"Best fit for spectrum slice {i} = {bestFit:.4f}mm")
461 return bestFits
463 def _generateLegendText(self, nSpectrumSlices: int) -> list[str]:
464 """Build legend entries labelling each spectrum slice.
466 Parameters
467 ----------
468 nSpectrumSlices : `int`
469 The number of slices to label.
471 Returns
472 -------
473 legendText : `list` [`str`]
474 One legend label per slice, running from the blue end to
475 the red end of the ``m=+1`` order.
476 """
477 if nSpectrumSlices == 1:
478 return ["m=+1 spectrum slice"]
479 if nSpectrumSlices == 2:
480 return ["m=+1 blue end", "m=+1 red end"]
482 legendText = []
483 legendText.append("m=+1 blue end")
484 for i in range(nSpectrumSlices - 2):
485 legendText.append("m=+1 redder...")
486 legendText.append("m=+1 red end")
487 return legendText
490class NonSpectralFocusAnalyzer:
491 """Analyze a focus sweep taken on direct (non-dispersed) imaging data.
493 For each image, measure the FWHM of the main star and the 50/80/90%
494 encircled-energy radii using
495 `lsst.summit.utils.ImageExaminer`, and fit a parabola to each
496 metric as a function of hexapod focus to determine the best-focus
497 position.
499 Parameters
500 ----------
501 embargo : `bool`, optional
502 If `True`, use the embargo butler repo rather than the default
503 LATISS repo.
505 Examples
506 --------
507 >>> %matplotlib inline
508 >>> dayObs = 20210101
509 >>> seqNums = [100, 101, 102, 103, 104]
510 >>> focusAnalyzer = NonSpectralFocusAnalyzer()
511 >>> focusAnalyzer.getFocusData(dayObs, seqNums, doDisplay=True)
512 >>> focusAnalyzer.fitDataAndPlot()
514 `run` performs the two steps above in one call.
515 """
517 def __init__(self, embargo: bool = False):
518 self.butler = makeDefaultLatissButler(embargo=embargo)
519 self._bestEffort = BestEffortIsr(embargo=embargo)
521 @staticmethod
522 def gauss(x: float | np.ndarray, *pars: float) -> float | np.ndarray:
523 """Evaluate a 1D Gaussian.
525 Parameters
526 ----------
527 x : `float` or `numpy.ndarray`
528 The abscissa at which to evaluate the Gaussian.
529 *pars : `float`
530 The Gaussian parameters in the order
531 ``(amplitude, mean, sigma)``.
533 Returns
534 -------
535 y : `float` or `numpy.ndarray`
536 The Gaussian evaluated at ``x``.
537 """
538 amp, mean, sigma = pars
539 return amp * np.exp(-((x - mean) ** 2) / (2.0 * sigma**2))
541 def run(
542 self,
543 dayObs: int,
544 seqNums: list[int],
545 *,
546 manualCentroid: tuple[float, float] | None = None,
547 doCheckDispersed: bool = True,
548 doDisplay: bool = False,
549 doForceCoM: bool = False,
550 ) -> dict:
551 """Perform a focus sweep analysis for direct imaging data.
553 For each seqNum for the specified dayObs, run the image through imExam
554 and collect the widths from the Gaussian fit and the 50/80/90%
555 encircled energy metrics, saving the data in the class for fitting.
557 For each of the [Gaussian fit, 50%, 80%, 90% encircled energy] metrics,
558 fit a parabola and return the focus value at which the minimum is
559 found.
561 Parameters
562 ----------
563 dayObs : `int`
564 The dayObs to use.
565 seqNums : `list` [`int`]
566 The seqNums for the focus sweep to analyze.
567 manualCentroid : `tuple` [`float`, `float`], optional
568 If provided, use this ``(x, y)`` centroid for every image
569 instead of fitting one per image.
570 doCheckDispersed : `bool`, optional
571 If `True`, assert that none of the exposures are dispersed.
572 doDisplay : `bool`, optional
573 If `True`, render per-image diagnostic plots. Designed to
574 be used in a notebook with ``%matplotlib inline``.
575 doForceCoM : `bool`, optional
576 If `True`, force the image examiner to use a centre-of-mass
577 centroid rather than its default.
579 Returns
580 -------
581 result : `dict` [`str`, `float`]
582 The parabolic-fit minima keyed by the metric name:
583 ``fwhmFitMin``, ``ee50FitMin``, ``ee80FitMin``,
584 ``ee90FitMin``.
585 """
586 self.getFocusData(
587 dayObs,
588 seqNums,
589 manualCentroid=manualCentroid,
590 doCheckDispersed=doCheckDispersed,
591 doDisplay=doDisplay,
592 doForceCoM=doForceCoM,
593 )
594 bestFit = self.fitDataAndPlot()
595 return bestFit
597 def getFocusData(
598 self,
599 dayObs: int,
600 seqNums: list[int],
601 *,
602 manualCentroid: tuple[float, float] | None = None,
603 doCheckDispersed: bool = True,
604 doDisplay: bool = False,
605 doForceCoM: bool = False,
606 ) -> None:
607 """Perform a focus sweep analysis for direct imaging data.
609 For each seqNum for the specified dayObs, run the image through imExam
610 and collect the widths from the Gaussian fit and the 50/80/90%
611 encircled energy metrics, saving the data in the class for fitting.
613 Parameters
614 ----------
615 dayObs : `int`
616 The dayObs to use.
617 seqNums : `list` [`int`]
618 The seqNums for the focus sweep to analyze.
619 manualCentroid : `tuple` [`float`, `float`], optional
620 If provided, use this ``(x, y)`` centroid for every image
621 instead of fitting one per image.
622 doCheckDispersed : `bool`, optional
623 If `True`, assert that none of the exposures are dispersed.
624 doDisplay : `bool`, optional
625 If `True`, render the per-image ImageExaminer diagnostic
626 plot. Designed to be used in a notebook with
627 ``%matplotlib inline``.
628 doForceCoM : `bool`, optional
629 If `True`, force the image examiner to use a centre-of-mass
630 centroid rather than its default.
632 Notes
633 -----
634 Performs the focus analysis per-image, holding the data in the
635 class. Call `fitDataAndPlot` after this to fit a parabola to
636 the focus data itself.
637 """
638 fitData: dict[int, dict] = {}
639 filters = set()
640 objects = set()
642 maxDistance = 200
643 firstCentroid = None
645 for seqNum in seqNums:
646 fitData[seqNum] = {}
647 dataId = {"day_obs": dayObs, "seq_num": seqNum, "detector": 0}
648 exp = self._bestEffort.getExposure(dataId)
650 # sanity/consistency checking
651 filt = exp.filter.physicalLabel
652 expRecord = getExpRecordFromDataId(self.butler, dataId)
653 obj = expRecord.target_name
654 objects.add(obj)
655 filters.add(filt)
656 if doCheckDispersed:
657 assert not isDispersedExp(exp), f"Image is dispersed! (filter = {filt})"
658 assert len(filters) == 1, "You accidentally mixed filters!"
659 assert len(objects) == 1, "You accidentally mixed objects!"
661 imExam = ImageExaminer(
662 exp, centroid=manualCentroid, doTweakCentroid=True, boxHalfSize=105, doForceCoM=doForceCoM
663 )
664 if doDisplay:
665 imExam.plot()
667 fwhm = imExam.imStats.fitFwhm
668 amp = imExam.imStats.fitAmp
669 gausMean = imExam.imStats.fitGausMean
670 centroid = imExam.centroid
672 if seqNum == seqNums[0]:
673 firstCentroid = centroid
675 dist = norm(np.array(centroid) - np.array(firstCentroid))
676 if dist > maxDistance:
677 print(f"Skipping {seqNum} because distance {dist}> maxDistance {maxDistance}")
678 del fitData[seqNum]
679 continue
681 fitData[seqNum]["fitResult"] = FitResult(amp=amp, mean=gausMean, sigma=fwhm * FWHMTOSIGMA)
682 fitData[seqNum]["eeRadius50"] = imExam.imStats.eeRadius50
683 fitData[seqNum]["eeRadius80"] = imExam.imStats.eeRadius80
684 fitData[seqNum]["eeRadius90"] = imExam.imStats.eeRadius90
686 focuserPosition = getFocusFromExposure(exp)
687 fitData[seqNum]["focus"] = focuserPosition
689 self.fitData = fitData
690 self.filter = filters.pop()
691 self.object = objects.pop()
693 return
695 def fitDataAndPlot(self) -> dict:
696 """Fit a parabola to each width metric and plot the result.
698 Fits a parabola to the Gaussian FWHM and the 50%, 80% and 90%
699 encircled-energy radii as functions of the hexapod focus
700 position. Returns the focus position at which each parabola is
701 minimized.
703 Returns
704 -------
705 result : `dict` [`str`, `float`]
706 The parabolic-fit minima keyed by the metric name:
707 ``fwhmFitMin``, ``ee50FitMin``, ``ee80FitMin``,
708 ``ee90FitMin``.
709 """
710 fitData = self.fitData
712 labelFontSize = 14
714 pixelsPerArcsec = 10
716 fig = plt.figure(figsize=(10, 10)) # noqa
717 gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1])
719 seqNums = sorted(fitData.keys())
720 widths = [fitData[seqNum]["fitResult"].sigma * SIGMATOFWHM / pixelsPerArcsec for seqNum in seqNums]
721 focusPositions = [fitData[seqNum]["focus"] for seqNum in seqNums]
722 fineXs = np.linspace(np.min(focusPositions), np.max(focusPositions), 101)
724 fwhmFitPars = np.polyfit(focusPositions, widths, 2)
725 fwhmFitMin = -fwhmFitPars[1] / (2.0 * fwhmFitPars[0])
727 ax0 = plt.subplot(gs[0])
728 ax0.scatter(focusPositions, widths, c="k")
729 ax0.set_ylabel("FWHM (arcsec)", fontsize=labelFontSize)
730 ax0.plot(fineXs, np.poly1d(fwhmFitPars)(fineXs), "b-")
731 ax0.axvline(fwhmFitMin, c="r", ls="--")
733 ee90s = [fitData[seqNum]["eeRadius90"] for seqNum in seqNums]
734 ee80s = [fitData[seqNum]["eeRadius80"] for seqNum in seqNums]
735 ee50s = [fitData[seqNum]["eeRadius50"] for seqNum in seqNums]
736 ax1 = plt.subplot(gs[1], sharex=ax0)
737 ax1.scatter(focusPositions, ee90s, c="r", label="Encircled energy 90%")
738 ax1.scatter(focusPositions, ee80s, c="g", label="Encircled energy 80%")
739 ax1.scatter(focusPositions, ee50s, c="b", label="Encircled energy 50%")
741 ee90FitPars = np.polyfit(focusPositions, ee90s, 2)
742 ee90FitMin = -ee90FitPars[1] / (2.0 * ee90FitPars[0])
743 ee80FitPars = np.polyfit(focusPositions, ee80s, 2)
744 ee80FitMin = -ee80FitPars[1] / (2.0 * ee80FitPars[0])
745 ee50FitPars = np.polyfit(focusPositions, ee50s, 2)
746 ee50FitMin = -ee50FitPars[1] / (2.0 * ee50FitPars[0])
748 ax1.plot(fineXs, np.poly1d(ee90FitPars)(fineXs), "r-")
749 ax1.plot(fineXs, np.poly1d(ee80FitPars)(fineXs), "g-")
750 ax1.plot(fineXs, np.poly1d(ee50FitPars)(fineXs), "b-")
752 ax1.axvline(ee90FitMin, c="r", ls="--")
753 ax1.axvline(ee80FitMin, c="g", ls="--")
754 ax1.axvline(ee50FitMin, c="b", ls="--")
756 ax1.set_xlabel("User-applied focus offset (mm)", fontsize=labelFontSize)
757 ax1.set_ylabel("Radius (pixels)", fontsize=labelFontSize)
759 ax1.legend()
761 plt.subplots_adjust(hspace=0.0)
763 results = {
764 "fwhmFitMin": fwhmFitMin,
765 "ee90FitMin": ee90FitMin,
766 "ee80FitMin": ee80FitMin,
767 "ee50FitMin": ee50FitMin,
768 }
770 return results
773if __name__ == "__main__": 773 ↛ 775line 773 didn't jump to line 775 because the condition on line 773 was never true
774 # TODO: DM-34239 Move this to be a butler-driven test
775 analyzer = SpectralFocusAnalyzer()
776 # dataId = {'dayObs': '2020-02-20', 'seqNum': 485} # direct image
777 dataId = {"day_obs": 20200312}
778 seqNums = [121, 122]
779 analyzer.getFocusData(dataId["day_obs"], seqNums, doDisplay=True)
780 analyzer.fitDataAndPlot()