LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
subtractBrightStars.py
Go to the documentation of this file.
1# This file is part of pipe_tasks.
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/>.
21
22"""Retrieve extended PSF model and subtract bright stars at visit level."""
23
24__all__ = ["SubtractBrightStarsConnections", "SubtractBrightStarsConfig", "SubtractBrightStarsTask"]
25
26import logging
27from functools import reduce
28from operator import ior
29
30import numpy as np
31from lsst.afw.geom import SpanSet, Stencil
32from lsst.afw.image import Exposure, ExposureF, MaskedImageF
33from lsst.afw.math import (
34 StatisticsControl,
35 WarpingControl,
36 makeStatistics,
37 rotateImageBy90,
38 stringToStatisticsProperty,
39 warpImage,
40)
41from lsst.geom import Box2I, Point2D, Point2I
42from lsst.meas.algorithms import LoadReferenceObjectsConfig, ReferenceObjectLoader
43from lsst.meas.algorithms.brightStarStamps import BrightStarStamp, BrightStarStamps
44from lsst.pex.config import ChoiceField, ConfigField, Field, ListField
45from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct
46from lsst.pipe.base.connectionTypes import Input, Output, PrerequisiteInput
47from lsst.pipe.tasks.processBrightStars import ProcessBrightStarsTask
48
49logger = logging.getLogger(__name__)
50
51
53 PipelineTaskConnections,
54 dimensions=("instrument", "visit", "detector"),
55 defaultTemplates={
56 "outputExposureName": "brightStar_subtracted",
57 "outputBackgroundName": "brightStars",
58 "badStampsName": "brightStars",
59 },
60):
61 inputExposure = Input(
62 doc="Input exposure from which to subtract bright star stamps.",
63 name="calexp",
64 storageClass="ExposureF",
65 dimensions=(
66 "visit",
67 "detector",
68 ),
69 )
70 inputBrightStarStamps = Input(
71 doc="Set of preprocessed postage stamps, each centered on a single bright star.",
72 name="brightStarStamps",
73 storageClass="BrightStarStamps",
74 dimensions=(
75 "visit",
76 "detector",
77 ),
78 )
79 inputExtendedPsf = Input(
80 doc="Extended PSF model.",
81 name="extended_psf",
82 storageClass="ExtendedPsf",
83 dimensions=("band",),
84 )
85 skyCorr = Input(
86 doc="Input Sky Correction to be subtracted from the calexp if ``doApplySkyCorr``=True.",
87 name="skyCorr",
88 storageClass="Background",
89 dimensions=(
90 "instrument",
91 "visit",
92 "detector",
93 ),
94 )
95 refCat = PrerequisiteInput(
96 doc="Reference catalog that contains bright star positions",
97 name="gaia_dr3_20230707",
98 storageClass="SimpleCatalog",
99 dimensions=("skypix",),
100 multiple=True,
101 deferLoad=True,
102 )
103 outputExposure = Output(
104 doc="Exposure with bright stars subtracted.",
105 name="{outputExposureName}_calexp",
106 storageClass="ExposureF",
107 dimensions=(
108 "visit",
109 "detector",
110 ),
111 )
112 outputBackgroundExposure = Output(
113 doc="Exposure containing only the modelled bright stars.",
114 name="{outputBackgroundName}_calexp_background",
115 storageClass="ExposureF",
116 dimensions=(
117 "visit",
118 "detector",
119 ),
120 )
121 outputBadStamps = Output(
122 doc="The stamps that are not normalized and consequently not subtracted from the exposure.",
123 name="{badStampsName}_unsubtracted_stamps",
124 storageClass="BrightStarStamps",
125 dimensions=(
126 "visit",
127 "detector",
128 ),
129 )
130
131 def __init__(self, *, config=None):
132 super().__init__(config=config)
133 if not config.doApplySkyCorr:
134 self.inputs.remove("skyCorr")
135
136
137class SubtractBrightStarsConfig(PipelineTaskConfig, pipelineConnections=SubtractBrightStarsConnections):
138 """Configuration parameters for SubtractBrightStarsTask"""
139
140 doWriteSubtractor = Field[bool](
141 doc="Should an exposure containing all bright star models be written to disk?",
142 default=True,
143 )
144 doWriteSubtractedExposure = Field[bool](
145 doc="Should an exposure with bright stars subtracted be written to disk?",
146 default=True,
147 )
148 magLimit = Field[float](
149 doc="Magnitude limit, in Gaia G; all stars brighter than this value will be subtracted",
150 default=18,
151 )
152 minValidAnnulusFraction = Field[float](
153 doc="Minimum number of valid pixels that must fall within the annulus for the bright star to be "
154 "saved for subsequent generation of a PSF.",
155 default=0.0,
156 )
157 numSigmaClip = Field[float](
158 doc="Sigma for outlier rejection; ignored if annularFluxStatistic != 'MEANCLIP'.",
159 default=4,
160 )
161 numIter = Field[int](
162 doc="Number of iterations of outlier rejection; ignored if annularFluxStatistic != 'MEANCLIP'.",
163 default=3,
164 )
165 warpingKernelName = ChoiceField[str](
166 doc="Warping kernel",
167 default="lanczos5",
168 allowed={
169 "bilinear": "bilinear interpolation",
170 "lanczos3": "Lanczos kernel of order 3",
171 "lanczos4": "Lanczos kernel of order 4",
172 "lanczos5": "Lanczos kernel of order 5",
173 "lanczos6": "Lanczos kernel of order 6",
174 "lanczos7": "Lanczos kernel of order 7",
175 },
176 )
177 scalingType = ChoiceField[str](
178 doc="How the model should be scaled to each bright star; implemented options are "
179 "`annularFlux` to reuse the annular flux of each stamp, or `leastSquare` to perform "
180 "least square fitting on each pixel with no bad mask plane set.",
181 default="leastSquare",
182 allowed={
183 "annularFlux": "reuse BrightStarStamp annular flux measurement",
184 "leastSquare": "find least square scaling factor",
185 },
186 )
187 annularFluxStatistic = ChoiceField[str](
188 doc="Type of statistic to use to compute annular flux.",
189 default="MEANCLIP",
190 allowed={
191 "MEAN": "mean",
192 "MEDIAN": "median",
193 "MEANCLIP": "clipped mean",
194 },
195 )
196 badMaskPlanes = ListField[str](
197 doc="Mask planes that, if set, lead to associated pixels not being included in the computation of "
198 "the scaling factor (`BAD` should always be included). Ignored if scalingType is `annularFlux`, "
199 "as the stamps are expected to already be normalized.",
200 # Note that `BAD` should always be included, as secondary detected
201 # sources (i.e., detected sources other than the primary source of
202 # interest) also get set to `BAD`.
203 default=("BAD", "CR", "CROSSTALK", "EDGE", "NO_DATA", "SAT", "SUSPECT", "UNMASKEDNAN"),
204 )
205 subtractionBox = ListField[int](
206 doc="Size of the stamps to be extracted, in pixels.",
207 default=(250, 250),
208 )
209 subtractionBoxBuffer = Field[float](
210 doc=(
211 "'Buffer' (multiplicative) factor to be applied to determine the size of the stamp the "
212 "processed stars will be saved in. This is also the size of the extended PSF model. The buffer "
213 "region is masked and contain no data and subtractionBox determines the region where contains "
214 "the data."
215 ),
216 default=1.1,
217 )
218 doApplySkyCorr = Field[bool](
219 doc="Apply full focal plane sky correction before extracting stars?",
220 default=True,
221 )
222 refObjLoader = ConfigField[LoadReferenceObjectsConfig](
223 doc="Reference object loader for astrometric calibration.",
224 )
225
226
227class SubtractBrightStarsTask(PipelineTask):
228 """Use an extended PSF model to subtract bright stars from a calibrated
229 exposure (i.e. at single-visit level).
230
231 This task uses both a set of bright star stamps produced by
232 `~lsst.pipe.tasks.processBrightStars.ProcessBrightStarsTask`
233 and an extended PSF model produced by
234 `~lsst.pipe.tasks.extended_psf.MeasureExtendedPsfTask`.
235 """
236
237 ConfigClass = SubtractBrightStarsConfig
238 _DefaultName = "subtractBrightStars"
239
240 def __init__(self, *args, **kwargs):
241 super().__init__(*args, **kwargs)
242 # Placeholders to set up Statistics if scalingType is leastSquare.
243 self.statsControl, self.statsFlag = None, None
244 # Warping control; only contains shiftingALg provided in config.
245 self.warpControl = WarpingControl(self.config.warpingKernelName)
246
247 def runQuantum(self, butlerQC, inputRefs, outputRefs):
248 # Docstring inherited.
249 inputs = butlerQC.get(inputRefs)
250 if inputs["inputExtendedPsf"].default_extended_psf is None:
251 if not self._detectorInRegions(inputs["inputExposure"], inputs["inputExtendedPsf"]):
252 self.log.warn(
253 "Extended PSF model is not available for detector %i. Skipping withouth processing this "
254 "exposure.",
255 inputs["inputExposure"].detector.getId(),
256 )
257 return None
258 dataId = butlerQC.quantum.dataId
259 refObjLoader = ReferenceObjectLoader(
260 dataIds=[ref.datasetRef.dataId for ref in inputRefs.refCat],
261 refCats=inputs.pop("refCat"),
262 name=self.config.connections.refCat,
263 config=self.config.refObjLoader,
264 )
265 subtractor, _, badStamps = self.run(**inputs, dataId=dataId, refObjLoader=refObjLoader)
266 if self.config.doWriteSubtractedExposure:
267 outputExposure = inputs["inputExposure"].clone()
268 outputExposure.image -= subtractor.image
269 else:
270 outputExposure = None
271 outputBackgroundExposure = subtractor if self.config.doWriteSubtractor else None
272 # In its current state, the code produces outputBadStamps which are the
273 # stamps of stars that have not been subtracted from the image for any
274 # reason. If all the stars are subtracted from the calexp, the output
275 # is an empty fits file.
276 output = Struct(
277 outputExposure=outputExposure,
278 outputBackgroundExposure=outputBackgroundExposure,
279 outputBadStamps=badStamps,
280 )
281 butlerQC.put(output, outputRefs)
282
283 def run(
284 self, inputExposure, inputBrightStarStamps, inputExtendedPsf, dataId, skyCorr=None, refObjLoader=None
285 ):
286 """Iterate over all bright stars in an exposure to scale the extended
287 PSF model before subtracting bright stars.
288
289 Parameters
290 ----------
291 inputExposure : `~lsst.afw.image.ExposureF`
292 The image from which bright stars should be subtracted.
293 inputBrightStarStamps :
294 `~lsst.meas.algorithms.brightStarStamps.BrightStarStamps`
295 Set of stamps centered on each bright star to be subtracted,
296 produced by running
297 `~lsst.pipe.tasks.processBrightStars.ProcessBrightStarsTask`.
298 inputExtendedPsf : `~lsst.pipe.tasks.extended_psf.ExtendedPsf`
299 Extended PSF model, produced by
300 `~lsst.pipe.tasks.extended_psf.MeasureExtendedPsfTask`.
301 dataId : `dict` or `~lsst.daf.butler.DataCoordinate`
302 The dataId of the exposure (and detector) bright stars should be
303 subtracted from.
304 skyCorr : `~lsst.afw.math.backgroundList.BackgroundList`, optional
305 Full focal plane sky correction, obtained by running
306 `~lsst.pipe.tasks.skyCorrection.SkyCorrectionTask`. If
307 `doApplySkyCorr` is set to `True`, `skyCorr` cannot be `None`.
308 refObjLoader : `~lsst.meas.algorithms.ReferenceObjectLoader`, optional
309 Loader to find objects within a reference catalog.
310
311 Returns
312 -------
313 subtractorExp : `~lsst.afw.image.ExposureF`
314 An Exposure containing a scaled bright star model fit to every
315 bright star profile; its image can then be subtracted from the
316 input exposure.
317 invImages : `list` [`~lsst.afw.image.MaskedImageF`]
318 A list of small images ("stamps") containing the model, each scaled
319 to its corresponding input bright star.
320 """
321 self.inputExpBBox = inputExposure.getBBox()
322 if self.config.doApplySkyCorr and (skyCorr is not None):
323 self.log.info(
324 "Applying sky correction to exposure %s (exposure will be modified in-place).", dataId
325 )
326 self.applySkyCorr(inputExposure, skyCorr)
327
328 # Create an empty image the size of the exposure.
329 # TODO: DM-31085 (set mask planes).
330 subtractorExp = ExposureF(bbox=inputExposure.getBBox())
331 subtractor = subtractorExp.maskedImage
332
333 # Make a copy of the input model.
334 self.model = inputExtendedPsf(dataId["detector"]).clone()
335 self.modelStampSize = self.model.getDimensions()
336 # Number of 90 deg. rotations to reverse each stamp's rotation.
337 self.inv90Rots = 4 - inputBrightStarStamps.nb90Rots % 4
338 self.model = rotateImageBy90(self.model, self.inv90Rots)
339
340 brightStarList = self.makeBrightStarList(inputBrightStarStamps, inputExposure, refObjLoader)
341 invImages = []
342 subtractor, invImages = self.buildSubtractor(
343 inputBrightStarStamps, subtractor, invImages, multipleAnnuli=False
344 )
345 if brightStarList:
346 self.setMissedStarsStatsControl()
347 # This may change when multiple star bins are used for PSF
348 # creation.
349 innerRadius = inputBrightStarStamps._innerRadius
350 outerRadius = inputBrightStarStamps._outerRadius
351 brightStarStamps, badStamps = BrightStarStamps.initAndNormalize(
352 brightStarList,
353 innerRadius=innerRadius,
354 outerRadius=outerRadius,
355 nb90Rots=self.warpOutputs.nb90Rots,
356 imCenter=self.warper.modelCenter,
357 use_archive=True,
358 statsControl=self.missedStatsControl,
359 statsFlag=self.missedStatsFlag,
360 badMaskPlanes=self.warper.config.badMaskPlanes,
361 discardNanFluxObjects=False,
362 forceFindFlux=True,
363 )
364
365 self.psf_annular_fluxes = self.findPsfAnnularFluxes(brightStarStamps)
366 subtractor, invImages = self.buildSubtractor(
367 brightStarStamps, subtractor, invImages, multipleAnnuli=True
368 )
369 else:
370 badStamps = []
371 badStamps = BrightStarStamps(badStamps)
372
373 return subtractorExp, invImages, badStamps
374
375 def _detectorInRegions(self, inputExposure, inputExtendedPsf):
376 """Determine whether the input exposure's detector is in the region(s)
377 where the extended PSF model(s) is(are) available.
378
379 Parameters
380 ----------
381 inputExposure : `lsst.afw.image.ExposureF`
382 The image from which bright stars should be subtracted. The ID of
383 the detector will be used to determine whether the detector is in
384 the region(s) where the extended PSF model(s) is(are) available.
385 inputExtendedPsf: `~lsst.pipe.tasks.extended_psf.ExtendedPsf`
386 Extended PSF model(s), produced by
387 `~lsst.pipe.tasks.extended_psf.MeasureExtendedPsfTask`. The ID's of
388 the detectors in the region(s) where the extended PSF model(s)
389 is(are) available will be used to cross match with the ID of the
390 input exposure's detector.
391
392 Returns
393 -------
394 `bool`
395 True if the detector is in the region(s) where the extended PSF
396 model(s) is(are) available, False otherwise.
397 """
398 availableDetectors = [
399 detector
400 for detectorList in inputExtendedPsf.detectors_focal_plane_regions.values()
401 for detector in detectorList.detectors
402 ]
403 if inputExposure.detector.getId() in availableDetectors:
404 return True
405 else:
406 return False
407
408 def _setUpStatistics(self, exampleMask):
409 """Configure statistics control and flag, for use if ``scalingType`` is
410 `leastSquare`.
411 """
412 if self.config.scalingType == "leastSquare":
413 # Set the mask planes which will be ignored.
414 andMask = reduce(
415 ior,
416 (exampleMask.getPlaneBitMask(bm) for bm in self.config.badMaskPlanes),
417 )
418 self.statsControl = StatisticsControl(
419 andMask=andMask,
420 )
421 self.statsFlag = stringToStatisticsProperty("SUM")
422
423 def applySkyCorr(self, calexp, skyCorr):
424 """Apply correction to the sky background level.
425 Sky corrections can be generated via the SkyCorrectionTask within the
426 pipe_tools module. Because the sky model used by that code extends over
427 the entire focal plane, this can produce better sky subtraction.
428 The calexp is updated in-place.
429
430 Parameters
431 ----------
432 calexp : `~lsst.afw.image.Exposure` or `~lsst.afw.image.MaskedImage`
433 Calibrated exposure.
434 skyCorr : `~lsst.afw.math.backgroundList.BackgroundList`
435 Full focal plane sky correction, obtained by running
436 `~lsst.pipe.tasks.skyCorrection.SkyCorrectionTask`.
437 """
438 if isinstance(calexp, Exposure):
439 calexp = calexp.getMaskedImage()
440 calexp -= skyCorr.getImage()
441
442 def scaleModel(self, model, star, inPlace=True, nb90Rots=0, psf_annular_flux=1.0):
443 """Compute scaling factor to be applied to the extended PSF so that its
444 amplitude matches that of an individual star.
445
446 Parameters
447 ----------
448 model : `~lsst.afw.image.MaskedImageF`
449 The extended PSF model, shifted (and potentially warped) to match
450 the bright star position.
451 star : `~lsst.meas.algorithms.brightStarStamps.BrightStarStamp`
452 A stamp centered on the bright star to be subtracted.
453 inPlace : `bool`
454 Whether the model should be scaled in place. Default is `True`.
455 nb90Rots : `int`
456 The number of 90-degrees rotations to apply to the star stamp.
457 psf_annular_flux: `float`, optional
458 The annular flux of the PSF model at the radius where the flux of
459 the given star is determined. This is 1 for stars present in
460 inputBrightStarStamps, but can be different for stars that are
461 missing from inputBrightStarStamps.
462
463 Returns
464 -------
465 scalingFactor : `float`
466 The factor by which the model image should be multiplied for it
467 to be scaled to the input bright star.
468 """
469 if self.config.scalingType == "annularFlux":
470 scalingFactor = star.annularFlux * psf_annular_flux
471 elif self.config.scalingType == "leastSquare":
472 if self.statsControl is None:
473 self._setUpStatistics(star.stamp_im.mask)
474 starIm = star.stamp_im.clone()
475 # Rotate the star postage stamp.
476 starIm = rotateImageBy90(starIm, nb90Rots)
477 # Reverse the prior star flux normalization ("unnormalize").
478 starIm *= star.annularFlux
479 # The estimator of the scalingFactor (f) that minimizes (Y-fX)^2
480 # is E[XY]/E[XX].
481 xy = starIm.clone()
482 xy.image.array *= model.image.array
483 xx = starIm.clone()
484 xx.image.array = model.image.array**2
485 # Compute the least squares scaling factor.
486 xySum = makeStatistics(xy, self.statsFlag, self.statsControl).getValue()
487 xxSum = makeStatistics(xx, self.statsFlag, self.statsControl).getValue()
488 scalingFactor = xySum / xxSum if xxSum else 1
489 if inPlace:
490 model.image *= scalingFactor
491 return scalingFactor
492
493 def _overrideWarperConfig(self):
494 """Override the warper config with the config of this task.
495
496 This override is necessary for stars that are missing from the
497 inputBrightStarStamps object but still need to be subtracted.
498 """
499 # TODO: Replace these copied values with a warperConfig.
500 self.warper.config.minValidAnnulusFraction = self.config.minValidAnnulusFraction
501 self.warper.config.numSigmaClip = self.config.numSigmaClip
502 self.warper.config.numIter = self.config.numIter
503 self.warper.config.annularFluxStatistic = self.config.annularFluxStatistic
504 self.warper.config.badMaskPlanes = self.config.badMaskPlanes
505 self.warper.config.stampSize = self.config.subtractionBox
506 self.warper.modelStampBuffer = self.config.subtractionBoxBuffer
507 self.warper.config.magLimit = self.config.magLimit
508 self.warper.setModelStamp()
509
510 def setMissedStarsStatsControl(self):
511 """Configure statistics control for processing missing stars from
512 inputBrightStarStamps.
513 """
514 self.missedStatsControl = StatisticsControl(
515 numSigmaClip=self.warper.config.numSigmaClip,
516 numIter=self.warper.config.numIter,
517 )
518 self.missedStatsFlag = stringToStatisticsProperty(self.warper.config.annularFluxStatistic)
519
520 def setWarpTask(self):
521 """Create an instance of ProcessBrightStarsTask that will be used to
522 produce stamps of stars to be subtracted.
523 """
524 self.warper = ProcessBrightStarsTask()
525 self._overrideWarperConfig()
526 self.warper.modelCenter = self.modelStampSize[0] // 2, self.modelStampSize[1] // 2
527
528 def makeBrightStarList(self, inputBrightStarStamps, inputExposure, refObjLoader):
529 """Make a list of bright stars that are missing from
530 inputBrightStarStamps to be subtracted.
531
532 Parameters
533 ----------
534 inputBrightStarStamps :
535 `~lsst.meas.algorithms.brightStarStamps.BrightStarStamps`
536 Set of stamps centered on each bright star to be subtracted,
537 produced by running
538 `~lsst.pipe.tasks.processBrightStars.ProcessBrightStarsTask`.
539 inputExposure : `~lsst.afw.image.ExposureF`
540 The image from which bright stars should be subtracted.
541 refObjLoader : `~lsst.meas.algorithms.ReferenceObjectLoader`, optional
542 Loader to find objects within a reference catalog.
543
544 Returns
545 -------
546 brightStarList:
547 A list containing
548 `lsst.meas.algorithms.brightStarStamps.BrightStarStamp` of stars to
549 be subtracted.
550 """
551 self.setWarpTask()
552 missedStars = self.warper.extractStamps(
553 inputExposure, refObjLoader=refObjLoader, inputBrightStarStamps=inputBrightStarStamps
554 )
555 if missedStars.starStamps:
556 self.warpOutputs = self.warper.warpStamps(missedStars.starStamps, missedStars.pixCenters)
557 brightStarList = [
559 stamp_im=warp,
560 archive_element=transform,
561 position=self.warpOutputs.xy0s[j],
562 gaiaGMag=missedStars.gMags[j],
563 gaiaId=missedStars.gaiaIds[j],
564 minValidAnnulusFraction=self.warper.config.minValidAnnulusFraction,
565 )
566 for j, (warp, transform) in enumerate(
567 zip(self.warpOutputs.warpedStars, self.warpOutputs.warpTransforms)
568 )
569 ]
570 else:
571 brightStarList = []
572 return brightStarList
573
574 def initAnnulusImage(self):
575 """Initialize an annulus image of the given star.
576
577 Returns
578 -------
579 annulusImage : `~lsst.afw.image.MaskedImageF`
580 The initialized annulus image.
581 """
582 maskPlaneDict = self.model.mask.getMaskPlaneDict()
583 annulusImage = MaskedImageF(self.modelStampSize, planeDict=maskPlaneDict)
584 annulusImage.mask.array[:] = 2 ** maskPlaneDict["NO_DATA"]
585 return annulusImage
586
587 def createAnnulus(self, brightStarStamp):
588 """Create a circular annulus around the given star.
589
590 The circular annulus is set based on the inner and outer optimal radii.
591 These radii describe the annulus where the flux of the star is found.
592 The aim is to create the same annulus for the PSF model, eventually
593 measuring the model flux around that annulus.
594 An optimal radius usually differs from the radius where the PSF model
595 is normalized.
596
597 Parameters
598 ----------
599 brightStarStamp :
600 `~lsst.meas.algorithms.brightStarStamps.BrightStarStamp`
601 A stamp of a bright star to be subtracted.
602
603 Returns
604 -------
605 annulus : `~lsst.afw.image.MaskedImageF`
606 An annulus of the given star.
607 """
608 # Create SpanSet of annulus.
609 outerCircle = SpanSet.fromShape(
610 brightStarStamp.optimalOuterRadius, Stencil.CIRCLE, offset=self.warper.modelCenter
611 )
612 innerCircle = SpanSet.fromShape(
613 brightStarStamp.optimalInnerRadius, Stencil.CIRCLE, offset=self.warper.modelCenter
614 )
615 annulus = outerCircle.intersectNot(innerCircle)
616 return annulus
617
618 def applyStatsControl(self, annulusImage):
619 """Apply statistics control to the PSF annulus image.
620
621 Parameters
622 ----------
623 annulusImage : `~lsst.afw.image.MaskedImageF`
624 An image containing an annulus of the given model.
625
626 Returns
627 -------
628 annularFlux: float
629 The annular flux of the PSF model at the radius where the flux of
630 the given star is determined.
631 """
632 andMask = reduce(
633 ior, (annulusImage.mask.getPlaneBitMask(bm) for bm in self.warper.config.badMaskPlanes)
634 )
635 self.missedStatsControl.setAndMask(andMask)
636 annulusStat = makeStatistics(annulusImage, self.missedStatsFlag, self.missedStatsControl)
637 return annulusStat.getValue()
638
639 def findPsfAnnularFlux(self, brightStarStamp, maskedModel):
640 """Find the annular flux of the PSF model within a specified annulus.
641
642 This flux will be used for re-scaling the PSF to the level of stars
643 with bad stamps. Stars with bad stamps are those without a flux within
644 the normalization annulus.
645
646 Parameters
647 ----------
648 brightStarStamp :
649 `~lsst.meas.algorithms.brightStarStamps.BrightStarStamp`
650 A stamp of a bright star to be subtracted.
651 maskedModel : `~lsst.afw.image.MaskedImageF`
652 A masked image of the PSF model.
653
654 Returns
655 -------
656 annularFlux: float (between 0 and 1)
657 The annular flux of the PSF model at the radius where the flux of
658 the given star is determined.
659 """
660 annulusImage = self.initAnnulusImage()
661 annulus = self.createAnnulus(brightStarStamp)
662 annulus.copyMaskedImage(maskedModel, annulusImage)
663 annularFlux = self.applyStatsControl(annulusImage)
664 return annularFlux
665
666 def findPsfAnnularFluxes(self, brightStarStamps):
667 """Find the annular fluxes of the given PSF model.
668
669 Parameters
670 ----------
671 brightStarStamps :
672 `~lsst.meas.algorithms.brightStarStamps.BrightStarStamps`
673 The stamps of stars that will be subtracted from the exposure.
674
675 Returns
676 -------
677 PsfAnnularFluxes: numpy.array
678 A two column numpy.array containing annular fluxes of the PSF at
679 radii where the flux for stars exist (could be found).
680
681 Notes
682 -----
683 While the PSF model is normalized at a certain radius, the annular flux
684 of a star around that radius might be impossible to find. Therefore, we
685 have to scale the PSF model considering a radius where the star has an
686 identified flux. To do that, the flux of the model should be found and
687 used to adjust the scaling step.
688 """
689 outerRadii = []
690 annularFluxes = []
691 maskedModel = MaskedImageF(self.model.image)
692 # The model has wrong bbox values. Should be fixed in extended_psf.py?
693 maskedModel.setXY0(0, 0)
694 for star in brightStarStamps:
695 if star.optimalOuterRadius not in outerRadii:
696 annularFlux = self.findPsfAnnularFlux(star, maskedModel)
697 outerRadii.append(star.optimalOuterRadius)
698 annularFluxes.append(annularFlux)
699 return np.array([outerRadii, annularFluxes]).T
700
701 def preparePlaneModelStamp(self, brightStarStamp):
702 """Prepare the PSF plane model stamp.
703
704 It is called PlaneModel because, while it is a PSF model stamp that is
705 warped and rotated to the same orientation of a chosen star, it is not
706 yet scaled to the brightness level of the star.
707
708 Parameters
709 ----------
710 brightStarStamp :
711 `~lsst.meas.algorithms.brightStarStamps.BrightStarStamp`
712 The stamp of the star to which the PSF model will be scaled.
713
714 Returns
715 -------
716 bbox: `~lsst.geom.Box2I`
717 Contains the corner coordination and the dimensions of the model
718 stamp.
719
720 invImage: `~lsst.afw.image.MaskedImageF`
721 The extended PSF model, shifted (and potentially warped and
722 rotated) to match the bright star position.
723
724 Raises
725 ------
726 RuntimeError
727 Raised if warping of the model failed.
728
729 Notes
730 -----
731 Since detectors have different orientations, the PSF model should be
732 rotated to match the orientation of the detectors in some cases. To do
733 that, the code uses the inverse of the transform that is applied to the
734 bright star stamp to match the orientation of the detector.
735 """
736 # Set the origin.
737 self.model.setXY0(brightStarStamp.position)
738 # Create an empty destination image.
739 invTransform = brightStarStamp.archive_element.inverted()
740 invOrigin = Point2I(invTransform.applyForward(Point2D(brightStarStamp.position)))
741 bbox = Box2I(corner=invOrigin, dimensions=self.modelStampSize)
742 invImage = MaskedImageF(bbox)
743 # Apply inverse transform.
744 goodPix = warpImage(invImage, self.model, invTransform, self.warpControl)
745 if not goodPix:
746 # Do we want to find another way or just subtract the non-warped
747 # scaled model?
748 # Currently the code just leaves the failed ones un-subtracted.
749 raise RuntimeError(
750 f"Warping of a model failed for star {brightStarStamp.gaiaId}: no good pixel in output."
751 )
752 return bbox, invImage
753
754 def addScaledModel(self, subtractor, brightStarStamp, multipleAnnuli=False):
755 """Add the scaled model of the given star to the subtractor plane.
756
757 Parameters
758 ----------
759 subtractor : `~lsst.afw.image.MaskedImageF`
760 The full image containing the scaled model of bright stars to be
761 subtracted from the input exposure.
762 brightStarStamp :
763 `~lsst.meas.algorithms.brightStarStamps.BrightStarStamp`
764 The stamp of the star of which the PSF model will be scaled and
765 added to the subtractor.
766 multipleAnnuli : bool, optional
767 If true, the model should be scaled based on a flux at a radius
768 other than its normalization radius.
769
770 Returns
771 -------
772 subtractor : `~lsst.afw.image.MaskedImageF`
773 The input subtractor full image with the added scaled model at the
774 given star's location in the exposure.
775 invImage: `~lsst.afw.image.MaskedImageF`
776 The extended PSF model, shifted (and potentially warped) to match
777 the bright star position.
778 """
779 bbox, invImage = self.preparePlaneModelStamp(brightStarStamp)
780 bbox.clip(self.inputExpBBox)
781 if bbox.getArea() > 0:
782 if multipleAnnuli:
783 cond = self.psf_annular_fluxes[:, 0] == brightStarStamp.optimalOuterRadius
784 psf_annular_flux = self.psf_annular_fluxes[cond, 1][0]
785 self.scaleModel(
786 invImage,
787 brightStarStamp,
788 inPlace=True,
789 nb90Rots=self.inv90Rots,
790 psf_annular_flux=psf_annular_flux,
791 )
792 else:
793 self.scaleModel(invImage, brightStarStamp, inPlace=True, nb90Rots=self.inv90Rots)
794 # Replace NaNs before subtraction (all NaNs have the NO_DATA flag).
795 invImage.image.array[np.isnan(invImage.image.array)] = 0
796 subtractor[bbox] += invImage[bbox]
797 return subtractor, invImage
798
799 def buildSubtractor(self, brightStarStamps, subtractor, invImages, multipleAnnuli=False):
800 """Build an image containing potentially multiple scaled PSF models,
801 each at the location of a given bright star.
802
803 Parameters
804 ----------
805 brightStarStamps :
806 `~lsst.meas.algorithms.brightStarStamps.BrightStarStamps`
807 Set of stamps centered on each bright star to be subtracted,
808 produced by running
809 `~lsst.pipe.tasks.processBrightStars.ProcessBrightStarsTask`.
810 subtractor : `~lsst.afw.image.MaskedImageF`
811 The Exposure that will contain the scaled model of bright stars to
812 be subtracted from the exposure.
813 invImages : `list`
814 A list containing extended PSF models, shifted (and potentially
815 warped) to match the bright stars positions.
816 multipleAnnuli : bool, optional
817 This will be passed to addScaledModel method, by default False.
818
819 Returns
820 -------
821 subtractor : `~lsst.afw.image.MaskedImageF`
822 An Exposure containing a scaled bright star model fit to every
823 bright star profile; its image can then be subtracted from the
824 input exposure.
825 invImages: list
826 A list containing the extended PSF models, shifted (and potentially
827 warped) to match bright stars' positions.
828 """
829 for star in brightStarStamps:
830 if star.gaiaGMag < self.config.magLimit:
831 try:
832 # Add the scaled model at the star location to subtractor.
833 subtractor, invImage = self.addScaledModel(subtractor, star, multipleAnnuli)
834 invImages.append(invImage)
835 except RuntimeError as err:
836 logger.error(err)
837 return subtractor, invImages
Pass parameters to a Statistics object.
Definition Statistics.h:83
Parameters to control convolution.
An integer coordinate rectangle.
Definition Box.h:55