Coverage for python/lsst/source/injection/inject_engine.py: 66%
423 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-13 10:13 +0000
1# This file is part of source_injection.
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 __future__ import annotations
24__all__ = ["generate_galsim_objects", "inject_galsim_objects_into_exposure"]
26import numbers
27import os
28from collections import Counter
29from collections.abc import Generator, Iterable
30from typing import Any
32import galsim
33import galsim.errors as galsim_errors
34import numpy as np
35import numpy.ma as ma
36from astropy.io import fits
37from astropy.table import Table
38from galsim import GalSimFFTSizeError
39from pydantic import ConfigDict
40from pydantic.dataclasses import dataclass
42import lsst.afw.detection as afwDetection
43from lsst.afw.geom import SkyWcs
44from lsst.afw.image import ExposureF, PhotoCalib
45from lsst.geom import Box2I, Point2D, Point2I, SpherePoint, arcseconds, degrees
46from lsst.images import BoundsError, MaskPlane
47from lsst.images._mask import _guess_legacy_plane_map
48from lsst.images.cells import CellCoadd, CellIJ
49from lsst.pex.exceptions import InvalidParameterError, LogicError
51# This is off by default but we want it on since we catch the error
52galsim_errors.raise_fft_size_error = True
55def get_object_data(source_data: dict[str, Any], object_class: galsim.GSObject) -> dict[str, Any]:
56 """Assemble a dictionary of allowed keyword arguments and their
57 corresponding values to use when constructing a GSObject.
59 Parameters
60 ----------
61 source_data : `dict` [`str`, `Any`]
62 Dictionary of source data.
63 object_class : `galsim.gsobject.GSObject`
64 Class of GSObject to match against.
66 Returns
67 -------
68 object_data : `dict` [`str`, `Any`]
69 Dictionary of source data to pass to the GSObject constructor.
70 """
71 req = getattr(object_class, "_req_params", {})
72 opt = getattr(object_class, "_opt_params", {})
73 single = getattr(object_class, "_single_params", {})
75 object_data = {}
77 # Check required args.
78 for key in req:
79 if key not in source_data: 79 ↛ 80line 79 didn't jump to line 80 because the condition on line 79 was never true
80 raise ValueError(f"Required parameter {key} not found in input catalog.")
81 object_data[key] = source_data[key]
83 # Optional args.
84 for key in opt:
85 if key in source_data: 85 ↛ 86line 85 didn't jump to line 86 because the condition on line 85 was never true
86 object_data[key] = source_data[key]
88 # Single args.
89 for s in single:
90 count = 0
91 for key in s:
92 if key in source_data:
93 count += 1
94 if count > 1: 94 ↛ 95line 94 didn't jump to line 95 because the condition on line 94 was never true
95 raise ValueError(f"Only one of {s.keys()} allowed for type {object_class}.")
96 object_data[key] = source_data[key]
97 if count == 0: 97 ↛ 98line 97 didn't jump to line 98 because the condition on line 97 was never true
98 raise ValueError(f"One of the args {s.keys()} is required for type {object_class}.")
100 return object_data
103def get_shear_data(
104 source_data: dict[str, Any],
105 shear_attributes: list[str] = [
106 "g1",
107 "g2",
108 "g",
109 "e1",
110 "e2",
111 "e",
112 "eta1",
113 "eta2",
114 "eta",
115 "q",
116 "beta",
117 "shear",
118 ],
119) -> dict[str, Any]:
120 """Assemble a dictionary of allowed keyword arguments and their
121 corresponding values to use when constructing a Shear.
123 Parameters
124 ----------
125 source_data : `dict` [`str`, `Any`]
126 Dictionary of source data.
127 shear_attributes : `list` [`str`], optional
128 List of allowed shear attributes.
130 Returns
131 -------
132 shear_data : `dict` [`str`, `Any`]
133 Dictionary of source data to pass to the Shear constructor.
134 """
135 shear_params = set(shear_attributes) & set(source_data.keys())
136 shear_data = {shear_param: source_data[shear_param] for shear_param in shear_params}
137 # Return an empty dict if all of the values are masked
138 # This should only be the case for DeltaFunction, but validation of the
139 # return is left for the caller of this function.
140 # While None might perhaps have been a better choice here, changing the
141 # return type is unnecessary given that galsim.Shear can be initialized
142 # with no args and correctly returns a trivial shear.
143 if all(np.ma.is_masked(x) for x in shear_data.values()):
144 return {}
145 if (beta := shear_data.get("beta")) is not None: 145 ↛ 147line 145 didn't jump to line 147 because the condition on line 145 was always true
146 shear_data["beta"] = beta * galsim.degrees
147 return shear_data
150def generate_galsim_objects(
151 injection_catalog: Table,
152 photo_calib: PhotoCalib,
153 wcs: SkyWcs,
154 fits_alignment: str,
155 stamp_prefix: str = "",
156 logger: Any | None = None,
157) -> Generator[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None]:
158 """Generate GalSim objects from an injection catalog.
160 Parameters
161 ----------
162 injection_catalog : `astropy.table.Table`
163 Table of sources to be injected.
164 photo_calib : `lsst.afw.image.PhotoCalib`
165 Photometric calibration used to calibrate injected sources.
166 wcs : `lsst.afw.geom.SkyWcs`
167 WCS used to calibrate injected sources.
168 fits_alignment : `str`
169 Alignment of the FITS image to the WCS. Allowed values: "wcs", "pixel".
170 stamp_prefix : `str`
171 Prefix to add to the stamp name.
172 logger : `lsst.utils.logging.LsstLogAdapter`, optional
173 Logger to use for logging messages.
175 Yields
176 ------
177 sky_coords : `lsst.geom.SpherePoint`
178 RA/Dec coordinates of the source.
179 pixel_coords : `lsst.geom.Point2D`
180 Pixel coordinates of the source.
181 draw_size : `int`
182 Size of the stamp to draw.
183 object : `galsim.gsobject.GSObject`
184 A fully specified and transformed GalSim object.
185 """
186 if logger: 186 ↛ 187line 186 didn't jump to line 187 because the condition on line 186 was never true
187 source_types = Counter(injection_catalog["source_type"]) # type: ignore
188 grammar0 = "source" if len(injection_catalog) == 1 else "sources"
189 grammar1 = "type" if len(source_types) == 1 else "types"
190 logger.info(
191 "Generating %d injection %s consisting of %d unique %s: %s.",
192 len(injection_catalog),
193 grammar0,
194 len(source_types),
195 grammar1,
196 ", ".join(f"{k}({v})" for k, v in source_types.items()),
197 )
198 for source_data_full in injection_catalog:
199 items = dict(source_data_full).items() # type: ignore
200 source_data = {k: v for (k, v) in items if v is not ma.masked}
201 try:
202 sky_coords = SpherePoint(float(source_data["ra"]), float(source_data["dec"]), degrees)
203 except KeyError:
204 sky_coords = wcs.pixelToSky(float(source_data["x"]), float(source_data["y"]))
205 try:
206 pixel_coords = Point2D(source_data["x"], source_data["y"])
207 except KeyError:
208 pixel_coords = wcs.skyToPixel(sky_coords)
209 try:
210 inst_flux = photo_calib.magnitudeToInstFlux(source_data["mag"], pixel_coords)
211 except LogicError:
212 continue
213 try:
214 draw_size = int(source_data["draw_size"])
215 except KeyError:
216 draw_size = 0
218 if source_data["source_type"] == "Stamp": 218 ↛ 219line 218 didn't jump to line 219 because the condition on line 218 was never true
219 stamp_file = stamp_prefix + source_data["stamp"]
220 object = make_galsim_stamp(stamp_file, fits_alignment, wcs, sky_coords, inst_flux)
221 elif source_data["source_type"] == "Trail": 221 ↛ 222line 221 didn't jump to line 222 because the condition on line 221 was never true
222 object = make_galsim_trail(source_data, wcs, sky_coords, inst_flux)
223 else:
224 object = make_galsim_object(source_data, source_data["source_type"], inst_flux, logger)
226 yield sky_coords, pixel_coords, draw_size, object
229def make_galsim_object(
230 source_data: dict[str, Any],
231 source_type: str,
232 inst_flux: float,
233 logger: Any | None = None,
234) -> galsim.gsobject.GSObject:
235 """Make a generic GalSim object from a collection of source data.
237 Parameters
238 ----------
239 source_data : `dict` [`str`, `Any`]
240 Dictionary of source data.
241 source_type : `str`
242 Type of the source, corresponding to a GalSim class.
243 inst_flux : `float`
244 Instrumental flux of the source.
245 logger : `~lsst.utils.logging.LsstLogAdapter`, optional
247 Returns
248 -------
249 object : `galsim.gsobject.GSObject`
250 A fully specified and transformed GalSim object.
251 """
252 # Populate the non-shear and non-flux parameters.
253 object_class = getattr(galsim, source_type)
254 object_data = get_object_data(source_data, object_class)
255 object = object_class(**object_data)
256 # Create a version of the object with an area-preserving shear applied.
257 shear_data = get_shear_data(source_data)
258 if shear_data:
259 try:
260 object = object.shear(**shear_data)
261 except TypeError as err:
262 if logger:
263 logger.warning("Cannot apply shear to GalSim object: %s", err)
264 pass
265 # Apply the instrumental flux and return.
266 object = object.withFlux(inst_flux)
267 return object
270def make_galsim_trail(
271 source_data: dict[str, Any],
272 wcs: SkyWcs,
273 sky_coords: SpherePoint,
274 inst_flux: float,
275 trail_thickness: float = 1e-6,
276) -> galsim.gsobject.GSObject:
277 """Make a trail with GalSim from a collection of source data.
279 Parameters
280 ----------
281 source_data : `dict` [`str`, `Any`]
282 Dictionary of source data.
283 wcs : `lsst.afw.geom.SkyWcs`
284 World coordinate system.
285 sky_coords : `lsst.geom.SpherePoint`
286 Sky coordinates of the source.
287 inst_flux : `float`
288 Instrumental flux of the source.
289 trail_thickness : `float`
290 Thickness of the trail in pixels.
292 Returns
293 -------
294 object : `galsim.gsobject.GSObject`
295 A fully specified and transformed GalSim object.
296 """
297 # Make a 'thin' box to mimic a line surface brightness profile of default
298 # thickness = 1e-6 (i.e., much thinner than a pixel)
299 object = galsim.Box(source_data["trail_length"], trail_thickness)
300 try:
301 object = object.rotate(source_data["beta"] * galsim.degrees)
302 except KeyError:
303 pass
304 object = object.withFlux(inst_flux * source_data["trail_length"]) # type: ignore
305 # GalSim objects are assumed to be in sky-coords. As we want the trail to
306 # appear as defined above in image-coords, we must transform it here.
307 linear_wcs = wcs.linearizePixelToSky(sky_coords, arcseconds)
308 mat = linear_wcs.getMatrix()
309 object = object.transform(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1])
310 # Rescale flux; transformation preserves surface brightness, not flux.
311 object *= 1.0 / np.abs(mat[0, 0] * mat[1, 1] - mat[0, 1] * mat[1, 0])
312 return object # type: ignore
315def make_galsim_stamp(
316 stamp_file: str,
317 fits_alignment: str,
318 wcs: SkyWcs,
319 sky_coords: SpherePoint,
320 inst_flux: float,
321) -> galsim.gsobject.GSObject:
322 """Make a postage stamp with GalSim from a FITS file.
324 Parameters
325 ----------
326 stamp_file : `str`
327 Path to the FITS file containing the postage stamp.
328 fits_alignment : `str`
329 Alignment of the FITS image to the WCS. Allowed values: "wcs", "pixel".
330 wcs : `lsst.afw.geom.SkyWcs`
331 World coordinate system.
332 sky_coords : `lsst.geom.SpherePoint`
333 Sky coordinates of the source.
334 inst_flux : `float`
335 Instrumental flux of the source.
337 Returns
338 -------
339 object: `galsim.gsobject.GSObject`
340 A fully specified and transformed GalSim object.
341 """
342 stamp_file = stamp_file.strip()
343 if os.path.exists(stamp_file):
344 with fits.open(stamp_file) as hdul:
345 hdu_images = [hdu.is_image and hdu.size > 0 for hdu in hdul] # type: ignore
346 if any(hdu_images):
347 stamp_data = galsim.fits.read(stamp_file, read_header=True, hdu=np.where(hdu_images)[0][0])
348 else:
349 raise RuntimeError(f"Cannot find image in input FITS file {stamp_file}.")
350 match fits_alignment:
351 case "wcs":
352 # galsim.fits.read will always attach a WCS to its output.
353 # If it can't find a WCS in the FITS header, then it
354 # defaults to scale = 1.0 arcsec / pix. If that's the scale
355 # then we need to check if it was explicitly set or if it's
356 # just the default. If it's just the default then we should
357 # raise an exception.
358 if is_wcs_galsim_default(stamp_data.wcs, stamp_data.header): # type: ignore
359 raise RuntimeError(f"Cannot find WCS in input FITS file {stamp_file}")
360 case "pixel":
361 # We need to set stamp_data.wcs to the local target WCS.
362 linear_wcs = wcs.linearizePixelToSky(sky_coords, arcseconds)
363 mat = linear_wcs.getMatrix()
364 stamp_data.wcs = galsim.JacobianWCS( # type: ignore
365 mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]
366 )
367 object = galsim.InterpolatedImage(stamp_data, calculate_stepk=False)
368 object = object.withFlux(inst_flux)
369 return object # type: ignore
370 else:
371 raise RuntimeError(f"Cannot locate input FITS postage stamp {stamp_file}.")
374def is_wcs_galsim_default(
375 wcs: galsim.fitswcs.GSFitsWCS,
376 hdr: galsim.fits.FitsHeader,
377) -> bool:
378 """Decide if wcs = galsim.PixelScale(1.0) is explicitly present in header,
379 or if it's just the GalSim default.
381 Parameters
382 ----------
383 wcs : galsim.fitswcs.GSFitsWCS
384 Potentially default WCS.
385 hdr : galsim.fits.FitsHeader
386 Header as read in by GalSim.
388 Returns
389 -------
390 is_default : bool
391 True if default, False if explicitly set in header.
392 """
393 if wcs != galsim.PixelScale(1.0):
394 return False
395 if hdr.get("GS_WCS") is not None:
396 return False
397 if hdr.get("CTYPE1", "LINEAR") == "LINEAR":
398 return not any(k in hdr for k in ["CD1_1", "CDELT1"])
399 for wcs_type in galsim.fitswcs.fits_wcs_types:
400 # If one of these succeeds, then assume result is explicit.
401 try:
402 wcs_type._readHeader(hdr)
403 return False
404 except Exception:
405 pass
406 else:
407 return not any(k in hdr for k in ["CD1_1", "CDELT1"])
410def infer_gain_from_image(
411 exposure: ExposureF,
412 bad_mask_names: list[str] | None = None,
413 logger: Any | None = None,
414) -> float:
415 """Fit a straight line to image flux vs. estimated variance in each pixel
416 of an exposure, and from that infer an estimated average gain.
418 Parameters
419 ----------
420 exposure : `lsst.afw.image.ExposureF`
421 The exposure to inject synthetic sources into.
422 bad_mask_names : `list` [`str`], optional
423 Names of any mask plane bits that could indicate a pathological flux
424 vs. variance relationship.
425 logger : `lsst.utils.logging.LsstLogAdapter`, optional
426 Logger to use for logging messages.
428 Returns
429 -------
430 gain : `float`
431 The estimated average gain across all pixels.
432 """
433 image_arr = exposure.image.array
434 var_arr = exposure.variance.array
435 good = np.isfinite(image_arr) & np.isfinite(var_arr)
437 mask = exposure.mask
438 # Only look at mask names in the mask_plane_dict.
439 # Otherwise, getPlaneBitMask throws an exception.
440 mask_plane_dict = mask.getMaskPlaneDict()
441 if bad_mask_names is None:
442 bad_mask_names = []
443 missing = set(bad_mask_names) - set(mask_plane_dict)
444 if missing and logger:
445 logger.warning("Mask planes not found in exposure (skipping): %s", missing)
446 bad_mask_names = [name for name in bad_mask_names if name in mask_plane_dict]
447 if len(bad_mask_names) > 0:
448 bad_mask_bit_mask = mask.getPlaneBitMask(bad_mask_names)
449 good &= (mask.array.astype(int) & bad_mask_bit_mask) == 0
451 # Check for sufficient points to fit.
452 n_good = int(np.count_nonzero(good))
453 if n_good < 2:
454 if logger:
455 logger.warning("Too few valid pixels (%d) to infer gain from image; returning NaN.", n_good)
456 return float("nan")
458 fit = np.polyfit(image_arr[good], var_arr[good], deg=1)
459 slope = fit[0]
460 # The gain is 1 / slope, so a non-positive or non-finite slope gives a
461 # non-physical gain. Flag it as NaN rather than propagate a negative,
462 # zero, or infinite gain into the noise model.
463 if not np.isfinite(slope) or slope <= 0: 463 ↛ 464line 463 didn't jump to line 464 because the condition on line 463 was never true
464 if logger:
465 logger.warning("Inferred a non-physical variance-vs-flux slope (%g); returning NaN gain.", slope)
466 return float("nan")
467 return 1.0 / slope
470def get_gain_map(
471 exposure: ExposureF,
472 bad_mask_names: list[str] | None = None,
473 logger: Any | None = None,
474) -> galsim.Image:
475 """Retrieve gains for each pixel.
477 Parameters
478 ----------
479 exposure : `lsst.afw.image.ExposureF`
480 The exposure to inject synthetic sources into.
481 bad_mask_names : `list` [`str`], optional
482 Names of any mask plane bits that could indicate a pathological flux
483 vs. variance relationship.
484 logger : `lsst.utils.logging.LsstLogAdapter`, optional
485 Logger to use for logging messages.
487 Returns
488 -------
489 gain_map : `galsim.Image`
490 The gain to use in each pixel.
491 """
492 full_bbox = exposure.getBBox()
493 try:
494 amplifiers = exposure.getDetector().getAmplifiers()
495 bboxes = [amplifier.getBBox() for amplifier in amplifiers] if amplifiers else None
496 except AttributeError:
497 bboxes = None
498 if not bboxes: 498 ↛ 501line 498 didn't jump to line 501 because the condition on line 498 was always true
499 bboxes = [full_bbox]
501 gains = np.array(
502 [infer_gain_from_image(exposure.subset(bbox), bad_mask_names, logger) for bbox in bboxes],
503 dtype=float,
504 )
506 # Catch unphysical values and replace with fallback values
507 usable = np.isfinite(gains) & (gains > 0)
508 if not np.all(usable):
509 if np.any(usable): 509 ↛ 512line 509 didn't jump to line 512 because the condition on line 509 was never true
510 # Use the median of the regions that did fit (e.g. the other
511 # amplifiers of the same detector).
512 fallback = float(np.median(gains[usable]))
513 else:
514 # Nothing fit. Try a single fit over the whole exposure (which can
515 # succeed even when individual sub-regions do not), then fall back
516 # to unit gain as a last resort.
517 fallback = 1.0
518 if len(bboxes) > 1: 518 ↛ 519line 518 didn't jump to line 519 because the condition on line 518 was never true
519 whole = infer_gain_from_image(exposure, bad_mask_names, logger)
520 if np.isfinite(whole) and whole > 0:
521 fallback = whole
522 if logger: 522 ↛ 523line 522 didn't jump to line 523 because the condition on line 522 was never true
523 logger.warning(
524 "Could not infer gain for %d of %d region(s); using fallback gain %g.",
525 int(np.count_nonzero(~usable)),
526 len(gains),
527 fallback,
528 )
529 gains[~usable] = fallback
531 full_bounds = galsim.BoundsI(full_bbox.minX, full_bbox.maxX, full_bbox.minY, full_bbox.maxY)
532 gain_map = galsim.Image(full_bounds, dtype=float)
533 for bbox, gain in zip(bboxes, gains):
534 bounds = galsim.BoundsI(bbox.minX, bbox.maxX, bbox.minY, bbox.maxY)
535 gain_map[bounds].array[:] = 1.0 / gain
537 return gain_map
540def add_noise_to_galsim_image(
541 image_template: galsim.Image,
542 var_template: galsim.Image,
543 is_coadd: bool,
544 gain_map: galsim.Image,
545 object_common_bounds: galsim.BoundsI,
546 noise_seed: int | None = None,
547 logger: Any | None = None,
548) -> tuple[galsim.Image, galsim.Image]:
549 """Add noise to a supplied image, and adjust the estimated variance
550 accordingly. If the image represents a coadd, add Gaussian noise.
551 Otherwise add Poisson shot noise.
553 Parameters
554 ----------
555 image_template : `galsim.Image`
556 The image to add noise to.
557 var_template : `galsim.Image`
558 The variance to use for the noise generating process in each pixel.
559 is_coadd : `bool`
560 Whether the exposure is a coadd. If True, Gaussian noise is used;
561 if False, Poisson shot noise is used.
562 gain_map : `galsim.Image`
563 The gain to use in each pixel of the full exposure.
564 object_common_bounds : `galsim.BoundsI`
565 Specific region of the full exposure to which this image belongs.
566 noise_seed : `int`, optional
567 Random number generator seed to use for the noise realization.
568 logger : `lsst.utils.logging.LsstLogAdapter`, optional
569 Logger to use for logging messages.
571 Returns
572 -------
573 image_template : `galsim.Image`
574 The image with noise added.
575 var_template : `galsim.Image`
576 The estimated variance, adjusted in response to the added noise.
577 """
578 noise_template = var_template.copy()
579 negative_variance = var_template.array < 0
580 nonfinite_variance = ~np.isfinite(var_template.array)
581 if np.any(negative_variance):
582 if logger:
583 logger.debug("Setting negative-variance pixels to 0 variance for noise generation.")
584 noise_template.array[negative_variance] = 0
585 if np.any(nonfinite_variance):
586 if logger:
587 logger.debug("Setting non-finite-variance pixels to 0 variance for noise generation.")
588 noise_template.array[nonfinite_variance] = 0
590 rng = galsim.BaseDeviate(noise_seed)
592 if not is_coadd:
593 # Treat noise_template, which is the expected amount of
594 # injected flux divided by the gain, as the Poisson mean of
595 # the number of photons collected by each pixel.
596 noise = galsim.PoissonNoise(rng)
597 noise_template.addNoise(noise)
598 # Scale the photon count by the gain to get the amount of
599 # image flux to inject.
600 image_template = noise_template.copy()
601 image_template /= gain_map[object_common_bounds]
602 # Restore the original variance values
603 # of the bad-variance pixels.
604 bad_var = negative_variance | nonfinite_variance
605 noise_template.array[bad_var] = var_template.array[bad_var]
606 var_template = noise_template
607 else:
608 noise = galsim.VariableGaussianNoise(rng, noise_template)
609 image_template.addNoise(noise)
610 var_template = image_template.copy()
611 var_template *= gain_map[object_common_bounds]
613 return image_template, var_template
616def _get_galsim_psf(
617 *,
618 psf: afwDetection.Psf,
619 pixel_coords: Point2D,
620 bbox: Box2I,
621 calib_flux_radius: float,
622 galsim_wcs: galsim.JacobianWCS,
623 sky_coords: SpherePoint | None = None,
624 apply_aperture_correction: bool = True,
625):
626 """Return the GalSim PSF at given coordinates within a bbox.
628 Parameters
629 ----------
630 psf
631 The coadd PSF object to evaluate the PSF with.
632 pixel_coords
633 The pixel coordinates to get the PSF for.
634 bbox
635 A bounding box to search for an alternative PSF for. If pixel_coords
636 lies outside the box and PSF evaluation fails, it will retry using the
637 nearest point to pixel_coords within the box.
638 calib_flux_radius
639 The radius to use to compute the aperture flux.
640 galsim_wcs
641 The GalSim WCS to pass to the `galsim.InterpolatedImage` constructor.
642 sky_coords
643 The sky coordinates of pixel_coords. Only used for
644 apply_aperture_correction
645 If True, the PSF values will be divided by the aperture flux at the
646 average position of psf.
648 Returns
649 -------
650 galsim_psf
651 The GalSim PSF as a `galsim.InterpolatedImage`.
652 """
653 error_args = []
654 try:
655 psf_array = psf.computeKernelImage(pixel_coords).array
656 except InvalidParameterError:
657 # Try mapping to nearest point contained in bbox.
658 contained_point = Point2D(
659 np.clip(pixel_coords.x, bbox.minX, bbox.maxX), np.clip(pixel_coords.y, bbox.minY, bbox.maxY)
660 )
661 if pixel_coords == contained_point: # no difference, so skip immediately
662 error_args = ["Cannot compute PSF for object at %s; flagging and skipping.", sky_coords]
663 # Otherwise, try again with new point.
664 try:
665 psf_array = psf.computeKernelImage(contained_point).array
666 except InvalidParameterError:
667 error_args = ["Cannot compute PSF for object at %s; flagging and skipping.", sky_coords]
668 return None, error_args
670 # Compute the aperture corrected PSF interpolated image.
671 # Normalize first
672 # TODO: Should also deal with negative pixel values here as they are
673 # not valid for a PSF.
674 psf_array /= np.sum(psf_array)
675 if apply_aperture_correction: 675 ↛ 679line 675 didn't jump to line 679 because the condition on line 675 was always true
676 # TODO: Review whether we want to at least make this optional
677 aperture_correction = psf.computeApertureFlux(calib_flux_radius, psf.getAveragePosition())
678 psf_array /= aperture_correction
679 galsim_psf = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs)
681 return galsim_psf, error_args
684@dataclass(config=ConfigDict(arbitrary_types_allowed=True), kw_only=True)
685class CellBoundsInfo:
686 aperture_correction: float = 1.0
687 cell_ij: CellIJ
688 galsim_psf: galsim.InterpolatedImage
691def _inject_galsim_object_into_bounds(
692 *,
693 convolved_object,
694 position_d,
695 position_i,
696 object_common_bounds,
697 full_bounds,
698 galsim_image: galsim.Image,
699 galsim_variance: galsim.Image,
700 galsim_wcs,
701 gain_map,
702 exposure,
703 inject_variance: bool,
704 add_noise: bool,
705 noise_seed: int,
706 is_coadd: bool,
707 injection_core_size,
708 mask_plane_name,
709 mask_plane_core_name,
710 logger,
711):
712 common_image = galsim_image[object_common_bounds]
713 common_variance = galsim_variance[object_common_bounds]
715 offset = position_d - object_common_bounds.true_center
716 # Note, for preliminary_visit_image injection, pixel is already
717 # part of the PSF and for coadd injection, it's incorrect to
718 # include the output pixel.
719 # So for both cases, we draw using method='no_pixel'.
720 image_template = common_image.copy()
721 image_template = convolved_object.drawImage(
722 image_template, add_to_image=False, offset=offset, wcs=galsim_wcs, method="no_pixel"
723 )
725 # The amount of additional variance to inject,
726 # if inject_variance is true.
727 var_template = image_template.copy()
728 var_template *= gain_map[object_common_bounds]
730 if add_noise: 730 ↛ 731line 730 didn't jump to line 731 because the condition on line 730 was never true
731 image_template, var_template = add_noise_to_galsim_image(
732 image_template,
733 var_template,
734 is_coadd,
735 gain_map,
736 object_common_bounds,
737 noise_seed,
738 logger,
739 )
741 common_image += image_template
742 if inject_variance: 742 ↛ 745line 742 didn't jump to line 745 because the condition on line 742 was always true
743 common_variance += var_template
745 common_box = Box2I(
746 Point2I(object_common_bounds.xmin, object_common_bounds.ymin),
747 Point2I(object_common_bounds.xmax, object_common_bounds.ymax),
748 )
749 bitvalue = exposure.mask.getPlaneBitMask(mask_plane_name)
750 exposure[common_box].mask.array |= bitvalue
751 # Add a configurable pixel mask centered on the object. The mask must be
752 # large enough to always identify the core/peak of the injected
753 # source, but small enough that it rarely overlaps real sources.
754 sub_bounds_core = galsim.BoundsI(position_i).withBorder(injection_core_size // 2)
755 object_common_bounds_core = full_bounds & sub_bounds_core
756 if object_common_bounds_core.area() > 0:
757 common_box_core = Box2I(
758 Point2I(object_common_bounds_core.xmin, object_common_bounds_core.ymin),
759 Point2I(object_common_bounds_core.xmax, object_common_bounds_core.ymax),
760 )
761 bitvalue_core = exposure.mask.getPlaneBitMask(mask_plane_core_name)
762 exposure[common_box_core].mask.array |= bitvalue_core
765def inject_galsim_objects_into_exposure(
766 exposure: ExposureF | CellCoadd,
767 objects: Iterable[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None],
768 mask_plane_name: str = "INJECTED",
769 calib_flux_radius: float = 12.0,
770 draw_size_max: int = 1000,
771 inject_variance: bool = True,
772 add_noise: bool = True,
773 noise_seed: int = 0,
774 bad_mask_names: list[str] | None = None,
775 injection_core_size: int = 3,
776 logger: Any | None = None,
777) -> tuple[list[int], list[galsim.BoundsI], list[bool], list[bool]]:
778 """Inject sources into given exposure using GalSim.
780 Parameters
781 ----------
782 exposure : `lsst.afw.image.ExposureF` or `lsst.images.cells.CellCoadd`
783 The exposure to inject synthetic sources into.
784 objects : `Iterable` [`tuple`, None, None]
785 An iterator of tuples that contains (or generates) locations and object
786 surface brightness profiles to inject. The tuples should contain the
787 following elements: `lsst.geom.SpherePoint`, `lsst.geom.Point2D`,
788 `int`, `galsim.gsobject.GSObject`.
789 mask_plane_name : `str`
790 Name of the mask plane to use for the injected sources.
791 calib_flux_radius : `float`
792 Radius in pixels to use for the flux calibration. This is used to
793 produce the correct instrumental fluxes within the radius. The value
794 should match that of the field defined in slot_CalibFlux_instFlux.
795 draw_size_max : `int`
796 Maximum allowed size of the drawn object. If the object is larger than
797 this, the draw size will be clipped to this size.
798 inject_variance : `bool`
799 Whether, when injecting flux into the image plane, to inject a
800 corresponding amount of variance into the variance plane.
801 add_noise : `bool`
802 Whether to randomly vary the amount of injected image flux in each
803 pixel by an amount consistent with the amount of injected variance.
804 noise_seed : `int`
805 Seed for generating the noise of the first injected object. The seed
806 actually used increases by 1 for each subsequent object, to ensure
807 independent noise realizations.
808 injection_core_size : `int`
809 Size of the box around each injected source to mark in the injected
810 core mask. The value must be a positive odd integer.
811 bad_mask_names : `list[str]`, optional
812 List of mask plane names indicating pixels to ignore when fitting flux
813 vs variance in preparation for variance plane modification. If None,
814 then the all pixels are used.
815 logger : `lsst.utils.logging.LsstLogAdapter`, optional
816 Logger to use for logging messages.
818 Returns
819 -------
820 draw_sizes : `list` [`int`]
821 Draw sizes of the injected sources.
822 common_bounds : `list` [`galsim.BoundsI`]
823 Common bounds of the drawn objects.
824 fft_size_errors : `list` [`bool`]
825 Boolean flags indicating whether a GalSimFFTSizeError was raised.
826 psf_compute_errors : `list` [`bool`]
827 Boolean flags indicating whether a PSF computation error was raised.
828 """
829 is_cell = isinstance(exposure, CellCoadd)
830 if not is_cell and not isinstance(exposure, ExposureF): 830 ↛ 831line 830 didn't jump to line 831 because the condition on line 830 was never true
831 raise TypeError(f"Unsupported {type(exposure)=}")
832 if (
833 not isinstance(injection_core_size, numbers.Number)
834 or not (injection_core_size >= 1)
835 or not (injection_core_size % 2 == 1)
836 ):
837 raise ValueError(f"{injection_core_size=} must be a positive odd integer")
839 if is_cell:
840 cell_coadd = exposure
841 cell_psf = cell_coadd.psf
842 exposure = cell_coadd.to_legacy()
843 is_coadd = True
844 else:
845 is_coadd = exposure.info.getCoaddInputs() is not None
847 exposure.mask.addMaskPlane(mask_plane_name)
848 mask_plane_core_name = mask_plane_name + "_CORE"
849 exposure.mask.addMaskPlane(mask_plane_core_name)
850 wcs = exposure.getWcs()
851 bbox = exposure.getBBox()
852 psf = exposure.getPsf()
853 if logger: 853 ↛ 854line 853 didn't jump to line 854 because the condition on line 853 was never true
854 logger.info(
855 "Adding %s and %s mask planes to the exposure.",
856 mask_plane_name,
857 mask_plane_core_name,
858 )
859 full_bounds = galsim.BoundsI(bbox.minX, bbox.maxX, bbox.minY, bbox.maxY)
860 galsim_image = galsim.Image(exposure.image.array, bounds=full_bounds)
861 galsim_variance = galsim.Image(exposure.variance.array, bounds=full_bounds)
862 pixel_scale = wcs.getPixelScale(bbox.getCenter()).asArcseconds()
864 if is_cell:
865 all_bounds: list[tuple[galsim.BoundsI, CellBoundsInfo | None]] = []
866 cell_grid = cell_coadd.grid
867 contributions = cell_coadd.provenance.contributions
868 n_i, n_j = cell_grid.grid_size.as_tuple()
870 fallback_psf = None
871 cell_inputs: dict[CellIJ, str] = {}
873 # Iterate over all cells and find the ones with good PSFs
874 for cell_i in range(n_i):
875 for cell_j in range(n_j):
876 cell_ij = CellIJ(cell_i, cell_j)
877 bbox_cell = cell_grid.bbox_of(cell_ij)
878 bbox_cell_legacy = bbox_cell.to_legacy()
879 cen_bbox = bbox_cell_legacy.getCenter()
880 cen_x, cen_y = cen_bbox
881 try:
882 if cell_ij in cell_coadd.bounds.missing:
883 continue
884 is_cell_input = (contributions["cell_i"] == cell_i) & (contributions["cell_j"] == cell_j)
885 # Return a string of sorted visits to efficiently count
886 # unique combinations
887 cell_inputs[cell_ij] = ",".join(
888 str(visit) for visit in sorted(set(contributions["visit"][is_cell_input]))
889 )
891 # Use the native cell coadd PSF methods rather than
892 # _get_galsim_psf which uses the to_legacy PSF
893 galsim_psf = cell_psf.compute_kernel_image(x=cen_x, y=cen_y).array
894 mat = wcs.linearizePixelToSky(wcs.pixelToSky(x=cen_x, y=cen_y), arcseconds).getMatrix()
895 galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1])
896 aperture_correction = psf.computeApertureFlux(calib_flux_radius, cen_bbox)
897 galsim_psf = galsim.InterpolatedImage(
898 galsim.Image(galsim_psf / (np.sum(galsim_psf) * aperture_correction)),
899 wcs=galsim_wcs,
900 )
901 bounds_cell = galsim.BoundsI(
902 bbox_cell.min.x,
903 bbox_cell.max.x,
904 bbox_cell.min.y,
905 bbox_cell.max.y,
906 )
907 bounds_info = CellBoundsInfo(
908 aperture_correction=aperture_correction,
909 cell_ij=cell_ij,
910 galsim_psf=galsim_psf,
911 )
912 all_bounds.append((bounds_cell, bounds_info))
913 if fallback_psf is None:
914 fallback_psf = galsim_psf
915 except BoundsError:
916 # TODO: Should this ever happen? If it does, it probably
917 # means that cell_coadd.bounds.missing is wrong.
918 pass
919 cell_visits, cell_visits_index = np.unique(list(cell_inputs.values()), return_inverse=True)
920 # Store the index of the cell_visits to efficiently check if
921 # they're identical
922 cell_inputs_indices = dict(zip(cell_inputs.keys(), cell_visits_index, strict=True))
923 # Keep a dict for PSF/bounds retrieval prior to iteration over cells
924 cell_bounds = {bounds_info.cell_ij: bounds_info for _, bounds_info in all_bounds}
925 else:
926 all_bounds = [(full_bounds, None)]
928 # TODO: Refactor to not take exposure?
929 gain_map = get_gain_map(exposure, bad_mask_names, logger)
931 draw_sizes: list[int] = []
932 common_bounds: list[galsim.BoundsI] = []
933 fft_size_errors: list[bool] = []
934 psf_compute_errors: list[bool] = []
935 for i, (sky_coords, pixel_coords, draw_size, object) in enumerate(objects):
936 # Instantiate default returns in case of early exit from this loop.
937 draw_sizes.append(0)
938 common_bounds.append(galsim.BoundsI())
939 fft_size_errors.append(False)
940 psf_compute_errors.append(False)
942 # Get spatial coordinates and WCS.
943 posd = galsim.PositionD(pixel_coords.x, pixel_coords.y)
944 posi = galsim.PositionI(pixel_coords.x // 1, pixel_coords.y // 1)
945 if logger: 945 ↛ 946line 945 didn't jump to line 946 because the condition on line 945 was never true
946 logger.debug(f"Injecting synthetic source at {pixel_coords}.")
947 mat = wcs.linearizePixelToSky(sky_coords, arcseconds).getMatrix()
948 galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1])
950 # This check is here because sometimes the WCS is multivalued and
951 # objects that should not be included were being included.
952 galsim_pixel_scale = np.sqrt(galsim_wcs.pixelArea())
953 if galsim_pixel_scale < pixel_scale / 2 or galsim_pixel_scale > pixel_scale * 2: 953 ↛ 954line 953 didn't jump to line 954 because the condition on line 953 was never true
954 continue
956 # Get the fallback PSF (if available) to compute the injection box.
957 # Even if the PSF at the centroid of an object is not valid, nearby
958 # cells might be fine. For example, cells with saturated stars can
959 # end up with no visits below the masked pixel threshold and no PSF.
960 #
961 # Note that _get_galsim_psf already implements a fallback of getting
962 # the PSF at the nearest edge of the bbox for objects whose centroids
963 # land outside the exposure's bbox. Non-cell coadds and visits could
964 # also do something similar by searching for the nearest valid PSF,
965 # but that's not as trivial as picking an "average" cell.
966 if is_cell:
967 cell_ij_cen = cell_grid.index_of(x=pixel_coords.x, y=pixel_coords.y)
968 galsim_psf_centroid = cell_bounds.get(cell_ij_cen)
969 if galsim_psf_centroid is None:
970 galsim_psf_centroid, psf_error_args = fallback_psf, []
971 cell_psf_index_cen = None
972 else:
973 galsim_psf_centroid = galsim_psf_centroid.galsim_psf
974 cell_psf_index_cen = cell_inputs_indices[cell_ij_cen]
975 else:
976 # Get the PSF at the centroid of the object
977 galsim_psf_centroid, psf_error_args = _get_galsim_psf(
978 psf=psf,
979 pixel_coords=pixel_coords,
980 bbox=bbox,
981 calib_flux_radius=calib_flux_radius,
982 galsim_wcs=galsim_wcs,
983 sky_coords=sky_coords,
984 )
985 if galsim_psf_centroid is None: 985 ↛ 986line 985 didn't jump to line 986 because the condition on line 985 was never true
986 psf_compute_errors[i] = True
987 if logger:
988 logger.debug(*psf_error_args)
989 continue
991 # Convolve the object with the PSF and generate draw size.
992 conv = galsim.Convolve(object, galsim_psf_centroid)
993 if draw_size == 0: 993 ↛ 995line 993 didn't jump to line 995 because the condition on line 993 was always true
994 draw_size = conv.getGoodImageSize(galsim_wcs.minLinearScale()) # type: ignore
995 injection_draw_size = int(draw_size)
996 injection_core_size_obj = injection_core_size
997 if draw_size_max > 0 and injection_draw_size > draw_size_max: 997 ↛ 998line 997 didn't jump to line 998 because the condition on line 997 was never true
998 if logger:
999 logger.warning(
1000 "Clipping draw size for object at %s from %d to %d pixels.",
1001 sky_coords,
1002 injection_draw_size,
1003 draw_size_max,
1004 )
1005 injection_draw_size = draw_size_max
1006 draw_sizes[i] = injection_draw_size
1007 if injection_core_size_obj > injection_draw_size: 1007 ↛ 1008line 1007 didn't jump to line 1008 because the condition on line 1007 was never true
1008 if logger:
1009 logger.debug(
1010 "Clipping core size for object at %s from %d to %d pixels.",
1011 sky_coords,
1012 injection_core_size_obj,
1013 injection_draw_size,
1014 )
1015 injection_core_size_obj = injection_draw_size
1016 sub_bounds = galsim.BoundsI(posi).withBorder(injection_draw_size // 2)
1018 # These bounds may not be the same as what would be derived from
1019 # per-cell PSFs for large objects, but re-calculating them for every
1020 # cell is not worthwhile
1021 if is_cell:
1022 object_common_bounds = full_bounds & sub_bounds
1023 if object_common_bounds.area() > 0: 1023 ↛ 1026line 1023 didn't jump to line 1026 because the condition on line 1023 was always true
1024 common_bounds[i] = object_common_bounds # type: ignore
1026 area_to_inject = 0
1027 area_injected = 0
1028 bounds_fft_size_errors = {}
1029 for idx_bound, (injection_bounds, bounds_info) in enumerate(all_bounds):
1030 object_common_bounds = injection_bounds & sub_bounds
1032 # Inject the source if there is any overlap.
1033 if (common_area := object_common_bounds.area()) > 0:
1034 area_to_inject += common_area
1035 if is_cell:
1036 # Use the per-cell PSF
1037 galsim_psf = bounds_info.galsim_psf
1038 assert galsim_psf is not None
1039 if bounds_info.cell_ij == cell_ij_cen:
1040 # These have to be one and the same if the cell is not
1041 # missing, which it can't be if it's in all_bounds
1042 assert galsim_psf is galsim_psf_centroid
1043 elif cell_psf_index_cen == cell_inputs_indices[bounds_info.cell_ij]:
1044 # Use the PSF from the cell containing the centroid
1045 assert galsim_psf is not galsim_psf_centroid
1046 galsim_psf = galsim_psf_centroid
1047 conv = galsim.Convolve(object, galsim_psf)
1048 else:
1049 common_bounds[i] = object_common_bounds # type: ignore
1050 try:
1051 if is_cell:
1052 # only the cell containing the centroid should have
1053 # the core mask plane set
1054 injection_core_bounds = full_bounds & object_common_bounds.withBorder(
1055 injection_core_size // 2
1056 )
1057 else:
1058 injection_core_bounds = full_bounds
1059 _inject_galsim_object_into_bounds(
1060 convolved_object=conv,
1061 position_d=posd,
1062 position_i=posi,
1063 object_common_bounds=object_common_bounds,
1064 full_bounds=injection_core_bounds,
1065 galsim_image=galsim_image,
1066 galsim_variance=galsim_variance,
1067 galsim_wcs=galsim_wcs,
1068 gain_map=gain_map,
1069 exposure=exposure,
1070 inject_variance=inject_variance,
1071 add_noise=add_noise,
1072 noise_seed=noise_seed,
1073 is_coadd=is_coadd,
1074 injection_core_size=injection_core_size_obj,
1075 mask_plane_name=mask_plane_name,
1076 mask_plane_core_name=mask_plane_core_name,
1077 logger=logger,
1078 )
1079 area_injected += common_area
1080 except GalSimFFTSizeError as err:
1081 bounds_fft_size_errors[idx_bound] = err.size
1082 if logger and (area_to_inject == 0): 1082 ↛ 1083line 1082 didn't jump to line 1083 because the condition on line 1082 was never true
1083 logger.debug("No area overlap for object at %s; flagging and skipping.", sky_coords)
1084 if bounds_fft_size_errors: 1084 ↛ 1085line 1084 didn't jump to line 1085 because the condition on line 1084 was never true
1085 fft_size_errors[i] = True
1086 if logger:
1087 logger.debug(
1088 "GalSimFFTSizeError raised for object at index=%d and coords %s;"
1089 " flagging and skipping.\nbounds_fft_size_errors=%s",
1090 i,
1091 sky_coords,
1092 bounds_fft_size_errors,
1093 )
1095 # lsst.images and afw masks are incompatible due to dtypes alone, and so
1096 # they have to be converted explicitly (also to add plane descriptions)
1097 if is_cell:
1098 plane_map = _guess_legacy_plane_map(exposure.mask.getMaskPlaneDict())
1099 plane_map.update(
1100 {
1101 mask_plane_name: MaskPlane(
1102 name=mask_plane_name,
1103 description="Pixel had at least one source injected over it",
1104 ),
1105 mask_plane_core_name: MaskPlane(
1106 name=mask_plane_core_name,
1107 description=f"Pixel center is within the core region"
1108 f" ({injection_core_size}x{injection_core_size} box) of an injected source",
1109 ),
1110 }
1111 )
1112 converted = cell_coadd.mask.from_legacy(
1113 exposure.mask,
1114 plane_map=plane_map,
1115 sky_projection=cell_coadd.sky_projection,
1116 )
1117 cell_coadd.mask = converted
1118 return draw_sizes, common_bounds, fft_size_errors, psf_compute_errors