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