Coverage for python / lsst / pipe / tasks / diffractionSpikeMask.py: 26%
155 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:11 +0000
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/>.
23__all__ = ["DiffractionSpikeMaskConfig", "DiffractionSpikeMaskTask"]
26import math
27import numpy as np
28import astropy.units as u
30import lsst.afw.geom as afwGeom
31from lsst.afw.image import abMagErrFromFluxErr
32import lsst.afw.table as afwTable
33import lsst.geom
34from lsst.pex.config import Config, ConfigField, Field
35from lsst.pipe.base import Struct, Task
36from lsst.meas.algorithms import getRefFluxField, LoadReferenceObjectsConfig
37from lsst.utils.timer import timeMethod
38import lsst.sphgeom
39from .colorterms import ColortermLibrary
42class DiffractionSpikeMaskConfig(Config):
43 """Config for BrightStarMaskTask.
44 """
45 refObjLoader = ConfigField(dtype=LoadReferenceObjectsConfig,
46 doc="Configuration of reference object loader")
47 applyColorTerms = Field(
48 dtype=bool,
49 default=False,
50 doc=("Apply photometric color terms to reference stars?\n"
51 "`True`: attempt to apply color terms; fail if color term data is "
52 "not available for the specified reference catalog and filter.\n"
53 "`False`: do not apply color terms."),
54 optional=True,
55 )
56 colorterms = ConfigField(
57 dtype=ColortermLibrary,
58 doc="Library of photometric reference catalog name: color term dict"
59 " (see also applyColorTerms).",
60 )
61 photoCatName = Field(
62 dtype=str,
63 optional=True,
64 doc=("Name of photometric reference catalog; used to select a color"
65 " term dict in colorterms. See also applyColorTerms."),
66 )
67 raKey = Field(
68 dtype=str,
69 default="coord_ra",
70 doc="RA column name in the reference catalog.",
71 )
72 decKey = Field(
73 dtype=str,
74 default="coord_dec",
75 doc="Declination column name in the reference catalog.",
76 )
77 angleMargin = Field(
78 dtype=float,
79 default=60.,
80 doc="Margin outside the exposure bounding box to include bright "
81 "sources. In arcseconds.",
82 )
83 magnitudeThreshold = Field(
84 dtype=float,
85 default=15,
86 doc="Threshold magnitude for treating a star from the reference catalog"
87 " as bright.",
88 )
89 diffractAngle = Field(
90 dtype=float,
91 default=45,
92 doc="Angle in degrees of the location of diffraction spikes with "
93 "respect to camera at 0 rotation angle.",
94 )
95 spikeAspectRatio = Field(
96 dtype=float,
97 default=10,
98 doc="Ratio of the length of a diffraction spike to it's width in the"
99 " core of the star.",
100 )
101 magSlope = Field(
102 dtype=float,
103 default=-0.12,
104 doc="Slope of the fit for the log(spike length) as a function of"
105 " magnitude.",
106 )
107 magOffset = Field(
108 dtype=float,
109 default=3.8,
110 doc="Intercept of the fit for the log(spike length) as a function of"
111 " magnitude.",
112 )
113 fallbackMagnitude = Field(
114 dtype=float,
115 default=12.,
116 doc="Default magnitude to use for sources in the reference catalog"
117 " with missing magnitudes that land in regions of saturated pixels.",
118 )
119 fallbackFluxField = Field(
120 dtype=str,
121 default="phot_g_mean",
122 doc="Fallback flux field in the reference catalog to use for sources"
123 " that don't have measurements in the science image's band.",
124 )
125 spikeMask = Field(
126 dtype=str,
127 default="SPIKE",
128 doc="Name of the mask plane indicating likely contamination from"
129 " a diffraction spike.",
130 )
131 saturatedMaskPlane = Field(
132 dtype=str,
133 default="SAT",
134 doc="Name of the mask plane indicating a saturated pixel.",
135 )
138class DiffractionSpikeMaskTask(Task):
139 """Load a reference catalog to identify bright stars that are likely to be
140 saturated and have visible diffraction spikes that need to be masked.
142 Attributes
143 ----------
144 angles : `numpy.ndarray`
145 Expected angles of diffraction spikes for bright sources.
146 refObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`
147 An instance of a reference object loader.
148 """
149 ConfigClass = DiffractionSpikeMaskConfig
150 _DefaultName = "diffractionSpikeMask"
152 def __init__(self, refObjLoader=None, **kwargs):
153 Task.__init__(self, **kwargs)
154 self.refObjLoader = refObjLoader
156 def setRefObjLoader(self, refObjLoader):
157 """Set the reference object loader for the task.
159 Parameters
160 ----------
161 refObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`
162 An instance of a reference object loader.
163 """
164 self.refObjLoader = refObjLoader
166 @timeMethod
167 def run(self, exposure):
168 """Load reference objects and mask bright stars on an exposure.
170 Parameters
171 ----------
172 exposure : `lsst.afw.image.Exposure`
173 Science exposure to set the SPIKE plane. Will be modified in place.
175 Returns
176 -------
177 spikeCat : `lsst.afw.table.SimpleCatalog`
178 The entries from the reference catalog selected as stars with
179 diffraction spikes.
180 """
181 if self.refObjLoader is None:
182 raise RuntimeError("Running diffraction spike mask task with no refObjLoader set in"
183 " __init__ or setRefObjLoader")
184 # First set the angles of the diffraction spikes from the exposure
185 # metadata.
186 self.set_diffraction_angle(exposure)
187 filterLabel = exposure.getFilter()
188 region = getRegion(exposure, lsst.sphgeom.Angle.fromDegrees(self.config.angleMargin/3600.))
189 refCat = self.refObjLoader.loadRegion(region, filterLabel.bandLabel).refCat
191 # Mask any sources with known or estimated magnitudes for the current
192 # filter, including sources off of the image which may have diffraction
193 # spikes that extend onto the image.
194 magnitudes = self.extractMagnitudes(refCat, filterLabel).refMag
195 radii = self.calculateReferenceRadius(magnitudes)
196 bright = (magnitudes < self.config.magnitudeThreshold)
198 nBright = np.count_nonzero(bright)
200 mask = exposure.mask
201 mask.addMaskPlane(self.config.spikeMask)
202 if nBright > 0:
203 xvals, yvals = exposure.wcs.skyToPixelArray(refCat[bright][self.config.raKey],
204 refCat[bright][self.config.decKey])
205 spikeCandidates = self.selectSources(xvals, yvals, radii, mask)
206 nSpike = len(spikeCandidates)
207 else:
208 nSpike = 0
209 if nSpike > 0:
210 self.log.info("Calculating mask for %d stars brighter than magnitude %f", nSpike,
211 self.config.magnitudeThreshold)
212 self.maskSources(xvals[spikeCandidates],
213 yvals[spikeCandidates],
214 radii[bright][spikeCandidates],
215 mask)
216 else:
217 self.log.info("No bright stars found in the reference catalog; not masking diffraction spikes.")
218 return afwTable.SimpleCatalog(refCat.schema)
220 return refCat[bright][spikeCandidates].copy(deep=True)
222 def selectSources(self, xvals, yvals, spikeRadii, mask):
223 """Select saturated sources, and bright sources that are off the image.
225 Parameters
226 ----------
227 xvals, yvals : `numpy.ndarray`
228 Array of x- and y-values of bright sources to mask.
229 mask : `lsst.afw.image.Mask`
230 The mask plane of the image to set the SPIKE mask plane.
231 spikeRadii : `numpy.ndarray`
232 Predicted lengths in pixels of the diffraction spikes for each
233 bright source.
235 Returns
236 -------
237 candidates : `numpy.ndarray`
238 Array of boolean flags indicating whether the given coordinates
239 should have a diffraction spike mask calculated.
240 """
241 candidates = np.zeros(len(xvals), dtype=bool)
242 saturatedBitMask = mask.getPlaneBitMask(self.config.saturatedMaskPlane)
243 bbox = mask.getBBox()
244 projection = np.max([abs(np.cos(np.deg2rad(self.angles))),
245 abs(np.sin(np.deg2rad(self.angles)))])
247 for i, (x, y, r) in enumerate(zip(xvals, yvals, spikeRadii)):
248 srcBox = lsst.geom.Box2I.makeCenteredBox(lsst.geom.Point2D(x, y), lsst.geom.Extent2I(3, 3))
249 if bbox.contains(srcBox):
250 candidates[i] = np.all((mask[srcBox].array & saturatedBitMask) > 0)
251 else:
252 # If the candidate bright source is off the image, only make a
253 # model for it if the diffraction spikes are long enough.
254 candidates[i] = boxSeparation(bbox, x, y) < r*projection
255 return candidates
257 def maskSources(self, xvals, yvals, spikeRadii, mask):
258 """Apply the SPIKE mask for a given set of coordinates. The mask plane
259 will be modified in place.
261 Parameters
262 ----------
263 xvals, yvals : `numpy.ndarray`
264 Array of x- and y-values of bright sources to mask.
265 spikeRadii : `numpy.ndarray`
266 Array of radius values for each bright source.
267 mask : `lsst.afw.image.Mask`
268 The mask plane of the image to set the SPIKE mask plane.
269 """
270 bbox = mask.getBBox()
271 for x, y, r in zip(xvals, yvals, spikeRadii):
272 maskSingle = self.makeSingleMask(x, y, r)
273 singleBBox = maskSingle.getBBox()
274 if bbox.overlaps(singleBBox):
275 singleBBox = singleBBox.clippedTo(bbox)
276 mask[singleBBox] |= maskSingle[singleBBox]
278 def makeSingleMask(self, x, y, r):
279 """Create a mask plane centered on a single source with the BRIGHT mask
280 set. This mask does not have to be fully contained in the bounding box
281 of the mask of the science image.
283 Parameters
284 ----------
285 x, y : `float`
286 Coordinates of the source to be masked.
287 r : `float`
288 Expected length of a diffraction spike for a source with a
289 magnitude
291 Returns
292 -------
293 mask : `lsst.afw.image.Mask`
294 A mask plane centered on the single source being modeled.
295 """
296 polygons = []
297 for angle in self.angles:
298 # Model the diffraction spike as an isosceles triangle with a long
299 # side along the spike and a short base. For efficiency, model the
300 # diffraction spike in the opposite direction at the same time,
301 # so the overall shape is a narrow diamond with equal length sides.
302 xLong = np.cos(np.deg2rad(angle))*r
303 yLong = np.sin(np.deg2rad(angle))*r
304 xShort = -np.sin(np.deg2rad(angle))*r/self.config.spikeAspectRatio
305 yShort = np.cos(np.deg2rad(angle))*r/self.config.spikeAspectRatio
307 corners = [lsst.geom.Point2D(x + xLong, y + yLong),
308 lsst.geom.Point2D(x + xShort, y + yShort),
309 lsst.geom.Point2D(x - xLong, y - yLong),
310 lsst.geom.Point2D(x - xShort, y - yShort)]
311 polygons.append(afwGeom.Polygon(corners))
312 # Combine all of the polygons into a single region
313 polygon = polygons[0]
314 for poly in polygons[1:]:
315 polygon = polygon.unionSingle(poly)
316 bbox = lsst.geom.Box2I(polygon.getBBox())
317 polyImage = polygon.createImage(bbox).array
318 mask = lsst.afw.image.Mask(bbox)
319 mask.array[polyImage > 0] = mask.getPlaneBitMask(self.config.spikeMask)
320 return mask
322 def set_diffraction_angle(self, exposure):
323 """Calculate the angle of diffration spikes on the image given the
324 camera rotation.
326 Parameters
327 ----------
328 exposure : `lsst.afw.image.Exposure`
329 The exposure to calculate the expected angle of diffraction spikes.
330 """
331 angle = (exposure.visitInfo.boresightParAngle.asDegrees()
332 - exposure.visitInfo.boresightRotAngle.asDegrees()
333 + self.config.diffractAngle)
334 self.angles = np.array([angle, angle + 90]) % 360
336 def calculateReferenceRadius(self, magnitudes):
337 """Calculate the size of the region to mask for each bright star.
339 Parameters
340 ----------
341 magnitudes : `numpy.ndarray`
342 Magnitudes of the bright stars.
344 Returns
345 -------
346 radii : `numpy.ndarray`
347 Array of radius values for the given magnitudes.
348 """
349 # In pixels
350 radii = [10**(self.config.magSlope*magnitude + self.config.magOffset) for magnitude in magnitudes]
351 return np.array(radii)
353 def extractMagnitudes(self, refCat, filterLabel):
354 """Extract magnitude and magnitude error arrays from the given catalog.
356 Parameters
357 ----------
358 refCat : `lsst.afw.table.SimpleCatalog`
359 The input reference catalog.
360 filterLabel : `str`
361 Label of filter being calibrated.
363 Returns
364 -------
365 result : `lsst.pipe.base.Struct`
366 Results as a struct with attributes:
368 ``refMag``
369 Reference magnitude (`np.array`).
370 ``refMagErr``
371 Reference magnitude error (`np.array`).
372 ``refFluxFieldList``
373 A list of field names of the reference catalog used for fluxes (1 or 2 strings) (`list`).
374 """
375 refSchema = refCat.schema
377 if self.config.applyColorTerms:
378 self.log.info("Applying color terms for filter=%r, config.photoCatName=%s",
379 filterLabel.physicalLabel, self.config.photoCatName)
380 colorterm = self.config.colorterms.getColorterm(filterLabel.physicalLabel,
381 self.config.photoCatName,
382 doRaise=True)
384 refMagArr, refMagErrArr = colorterm.getCorrectedMagnitudes(refCat)
385 fluxFieldList = [getRefFluxField(refSchema, filt) for filt in (colorterm.primary,
386 colorterm.secondary)]
387 else:
388 self.log.info("Not applying color terms.")
389 colorterm = None
391 fluxFieldList = [getRefFluxField(refSchema, filterLabel.bandLabel)]
392 fluxField = getRefFluxField(refSchema, filterLabel.bandLabel)
393 fluxKey = refSchema.find(fluxField).key
394 refFluxArr = np.array(refCat.get(fluxKey))
396 try:
397 fluxErrKey = refSchema.find(fluxField + "Err").key
398 refFluxErrArr = np.array(refCat.get(fluxErrKey))
399 except KeyError:
400 # Reference catalogue may not have flux uncertainties; HACK DM-2308
401 self.log.warning("Reference catalog does not have flux uncertainties for %s;"
402 " using sqrt(flux).", fluxField)
403 refFluxErrArr = np.sqrt(refFluxArr)
405 refMagArr = u.Quantity(refFluxArr, u.nJy).to_value(u.ABmag)
406 # HACK convert to Jy until we have a replacement for this (DM-16903)
407 refMagErrArr = abMagErrFromFluxErr(refFluxErrArr*1e-9, refFluxArr*1e-9)
408 # Not all sources in the reference catalog will have flux measurements
409 # in the current band. Since we are only masking sources use a fallback
410 # flux measurement in these cases to get an approximate magnitude to
411 # check for bright sources and to set the size of the mask.
412 badMagnitudes = np.isnan(refMagArr)
413 if np.count_nonzero(badMagnitudes):
414 fluxName = f"{self.config.fallbackFluxField}_flux"
415 # If the configured fallback flux field is missing, skip it.
416 try:
417 fallbackFluxKey = refSchema.find(fluxName).key
418 except KeyError:
419 self.log.warning(f"Field {fluxName} not found in photometric reference catalog."
420 " Not masking bright stars with missing magnitudes.")
421 else:
422 fallbackFluxArr = np.array(refCat.get(fallbackFluxKey))
423 fallbackMagArr = u.Quantity(fallbackFluxArr[badMagnitudes], u.nJy).to_value(u.ABmag)
424 refMagArr[badMagnitudes] = fallbackMagArr
426 return Struct(
427 refMag=refMagArr,
428 refMagErr=refMagErrArr,
429 refFluxFieldList=fluxFieldList,
430 )
433def getRegion(exposure, margin=None):
434 """Calculate an enveloping region for an exposure.
436 Parameters
437 ----------
438 exposure : `lsst.afw.image.Exposure`
439 Exposure object with calibrated WCS.
440 margin : `lsst.sphgeom.Angle`, optional
441 Grow the surrounding region of the exposure by this amount, in order to
442 mask the diffraction spikes of stars off the image.
444 Returns
445 -------
446 region : `lsst.sphgeom.Region`
447 Region enveloping an exposure.
448 """
449 # Bounding box needs to be a `Box2D` not a `Box2I` for `wcs.pixelToSky()`
450 bbox = lsst.geom.Box2D(exposure.getBBox())
451 wcs = exposure.getWcs()
453 region = lsst.sphgeom.ConvexPolygon([pp.getVector() for pp in wcs.pixelToSky(bbox.getCorners())])
454 if margin is not None:
455 # This is an ad-hoc, approximate implementation. It should be good
456 # enough for catalog loading, but is not a general-purpose solution.
457 center = lsst.geom.SpherePoint(region.getCentroid())
458 corners = [lsst.geom.SpherePoint(c) for c in region.getVertices()]
459 # Approximate the region as a Euclidian square
460 # geom.Angle(sphgeom.Angle) converter not pybind-wrapped???
461 diagonal_margin = lsst.geom.Angle(margin.asRadians() * math.sqrt(2.0))
462 padded = [c.offset(center.bearingTo(c), diagonal_margin) for c in corners]
463 return lsst.sphgeom.ConvexPolygon.convexHull([c.getVector() for c in padded])
465 return region
468def boxSeparation(bbox, x, y):
469 """Return the minimum horizontal or vertical distance from a point to the
470 outside edge of a bounding box.
472 Parameters
473 ----------
474 bbox : `lsst.geom.Box2I`
475 The bounding box to check.
476 x, y : `float`
477 Coordinates of the point.
479 Returns
480 -------
481 distance : `float`
482 The distance in pixels by which the point is outside the box, or 0 if
483 it is inside.
484 """
486 x0, y0 = bbox.getBegin()
487 x1, y1 = bbox.getEnd()
488 dx = 0 if x0 <= x <= x1 else min(abs(x0 - x), abs(x - x1))
489 dy = 0 if y0 <= y <= y1 else min(abs(y0 - y), abs(y - y1))
490 return max(dx, dy)