Coverage for python/lsst/source/injection/inject_engine.py: 50%
321 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 09:15 +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 os
27from collections import Counter
28from collections.abc import Generator
29from typing import Any
31import galsim
32import numpy as np
33import numpy.ma as ma
34from astropy.io import fits
35from astropy.table import Table
36from galsim import GalSimFFTSizeError
38from lsst.afw.geom import SkyWcs
39from lsst.afw.image import ExposureF, PhotoCalib
40from lsst.geom import Box2I, Point2D, Point2I, SpherePoint, arcseconds, degrees
41from lsst.pex.exceptions import InvalidParameterError, LogicError
44def get_object_data(source_data: dict[str, Any], object_class: galsim.GSObject) -> dict[str, Any]:
45 """Assemble a dictionary of allowed keyword arguments and their
46 corresponding values to use when constructing a GSObject.
48 Parameters
49 ----------
50 source_data : `dict` [`str`, `Any`]
51 Dictionary of source data.
52 object_class : `galsim.gsobject.GSObject`
53 Class of GSObject to match against.
55 Returns
56 -------
57 object_data : `dict` [`str`, `Any`]
58 Dictionary of source data to pass to the GSObject constructor.
59 """
60 req = getattr(object_class, "_req_params", {})
61 opt = getattr(object_class, "_opt_params", {})
62 single = getattr(object_class, "_single_params", {})
64 object_data = {}
66 # Check required args.
67 for key in req: 67 ↛ 68line 67 didn't jump to line 68 because the loop on line 67 never started
68 if key not in source_data:
69 raise ValueError(f"Required parameter {key} not found in input catalog.")
70 object_data[key] = source_data[key]
72 # Optional args.
73 for key in opt:
74 if key in source_data: 74 ↛ 75line 74 didn't jump to line 75 because the condition on line 74 was never true
75 object_data[key] = source_data[key]
77 # Single args.
78 for s in single: 78 ↛ 79line 78 didn't jump to line 79 because the loop on line 78 never started
79 count = 0
80 for key in s:
81 if key in source_data:
82 count += 1
83 if count > 1:
84 raise ValueError(f"Only one of {s.keys()} allowed for type {object_class}.")
85 object_data[key] = source_data[key]
86 if count == 0:
87 raise ValueError(f"One of the args {s.keys()} is required for type {object_class}.")
89 return object_data
92def get_shear_data(
93 source_data: dict[str, Any],
94 shear_attributes: list[str] = [
95 "g1",
96 "g2",
97 "g",
98 "e1",
99 "e2",
100 "e",
101 "eta1",
102 "eta2",
103 "eta",
104 "q",
105 "beta",
106 "shear",
107 ],
108) -> dict[str, Any]:
109 """Assemble a dictionary of allowed keyword arguments and their
110 corresponding values to use when constructing a Shear.
112 Parameters
113 ----------
114 source_data : `dict` [`str`, `Any`]
115 Dictionary of source data.
116 shear_attributes : `list` [`str`], optional
117 List of allowed shear attributes.
119 Returns
120 -------
121 shear_data : `dict` [`str`, `Any`]
122 Dictionary of source data to pass to the Shear constructor.
123 """
124 shear_params = set(shear_attributes) & set(source_data.keys())
125 shear_data = {}
126 for shear_param in shear_params: 126 ↛ 127line 126 didn't jump to line 127 because the loop on line 126 never started
127 if shear_param == "beta":
128 shear_data.update({shear_param: source_data[shear_param] * galsim.degrees})
129 else:
130 shear_data.update({shear_param: source_data[shear_param]})
131 return shear_data
134def generate_galsim_objects(
135 injection_catalog: Table,
136 photo_calib: PhotoCalib,
137 wcs: SkyWcs,
138 fits_alignment: str,
139 stamp_prefix: str = "",
140 logger: Any | None = None,
141) -> Generator[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None]:
142 """Generate GalSim objects from an injection catalog.
144 Parameters
145 ----------
146 injection_catalog : `astropy.table.Table`
147 Table of sources to be injected.
148 photo_calib : `lsst.afw.image.PhotoCalib`
149 Photometric calibration used to calibrate injected sources.
150 wcs : `lsst.afw.geom.SkyWcs`
151 WCS used to calibrate injected sources.
152 fits_alignment : `str`
153 Alignment of the FITS image to the WCS. Allowed values: "wcs", "pixel".
154 stamp_prefix : `str`
155 Prefix to add to the stamp name.
156 logger : `lsst.utils.logging.LsstLogAdapter`, optional
157 Logger to use for logging messages.
159 Yields
160 ------
161 sky_coords : `lsst.geom.SpherePoint`
162 RA/Dec coordinates of the source.
163 pixel_coords : `lsst.geom.Point2D`
164 Pixel coordinates of the source.
165 draw_size : `int`
166 Size of the stamp to draw.
167 object : `galsim.gsobject.GSObject`
168 A fully specified and transformed GalSim object.
169 """
170 if logger: 170 ↛ 171line 170 didn't jump to line 171 because the condition on line 170 was never true
171 source_types = Counter(injection_catalog["source_type"]) # type: ignore
172 grammar0 = "source" if len(injection_catalog) == 1 else "sources"
173 grammar1 = "type" if len(source_types) == 1 else "types"
174 logger.info(
175 "Generating %d injection %s consisting of %d unique %s: %s.",
176 len(injection_catalog),
177 grammar0,
178 len(source_types),
179 grammar1,
180 ", ".join(f"{k}({v})" for k, v in source_types.items()),
181 )
182 for source_data_full in injection_catalog:
183 items = dict(source_data_full).items() # type: ignore
184 source_data = {k: v for (k, v) in items if v is not ma.masked}
185 try:
186 sky_coords = SpherePoint(float(source_data["ra"]), float(source_data["dec"]), degrees)
187 except KeyError:
188 sky_coords = wcs.pixelToSky(float(source_data["x"]), float(source_data["y"]))
189 try:
190 pixel_coords = Point2D(source_data["x"], source_data["y"])
191 except KeyError:
192 pixel_coords = wcs.skyToPixel(sky_coords)
193 try:
194 inst_flux = photo_calib.magnitudeToInstFlux(source_data["mag"], pixel_coords)
195 except LogicError:
196 continue
197 try:
198 draw_size = int(source_data["draw_size"])
199 except KeyError:
200 draw_size = 0
202 if source_data["source_type"] == "Stamp": 202 ↛ 203line 202 didn't jump to line 203 because the condition on line 202 was never true
203 stamp_file = stamp_prefix + source_data["stamp"]
204 object = make_galsim_stamp(stamp_file, fits_alignment, wcs, sky_coords, inst_flux)
205 elif source_data["source_type"] == "Trail": 205 ↛ 206line 205 didn't jump to line 206 because the condition on line 205 was never true
206 object = make_galsim_trail(source_data, wcs, sky_coords, inst_flux)
207 else:
208 object = make_galsim_object(source_data, source_data["source_type"], inst_flux, logger)
210 yield sky_coords, pixel_coords, draw_size, object
213def make_galsim_object(
214 source_data: dict[str, Any],
215 source_type: str,
216 inst_flux: float,
217 logger: Any | None = None,
218) -> galsim.gsobject.GSObject:
219 """Make a generic GalSim object from a collection of source data.
221 Parameters
222 ----------
223 source_data : `dict` [`str`, `Any`]
224 Dictionary of source data.
225 source_type : `str`
226 Type of the source, corresponding to a GalSim class.
227 inst_flux : `float`
228 Instrumental flux of the source.
229 logger : `~lsst.utils.logging.LsstLogAdapter`, optional
231 Returns
232 -------
233 object : `galsim.gsobject.GSObject`
234 A fully specified and transformed GalSim object.
235 """
236 # Populate the non-shear and non-flux parameters.
237 object_class = getattr(galsim, source_type)
238 object_data = get_object_data(source_data, object_class)
239 object = object_class(**object_data)
240 # Create a version of the object with an area-preserving shear applied.
241 shear_data = get_shear_data(source_data)
242 if shear_data: 242 ↛ 243line 242 didn't jump to line 243 because the condition on line 242 was never true
243 try:
244 object = object.shear(**shear_data)
245 except TypeError as err:
246 if logger:
247 logger.warning("Cannot apply shear to GalSim object: %s", err)
248 pass
249 # Apply the instrumental flux and return.
250 object = object.withFlux(inst_flux)
251 return object
254def make_galsim_trail(
255 source_data: dict[str, Any],
256 wcs: SkyWcs,
257 sky_coords: SpherePoint,
258 inst_flux: float,
259 trail_thickness: float = 1e-6,
260) -> galsim.gsobject.GSObject:
261 """Make a trail with GalSim from a collection of source data.
263 Parameters
264 ----------
265 source_data : `dict` [`str`, `Any`]
266 Dictionary of source data.
267 wcs : `lsst.afw.geom.SkyWcs`
268 World coordinate system.
269 sky_coords : `lsst.geom.SpherePoint`
270 Sky coordinates of the source.
271 inst_flux : `float`
272 Instrumental flux of the source.
273 trail_thickness : `float`
274 Thickness of the trail in pixels.
276 Returns
277 -------
278 object : `galsim.gsobject.GSObject`
279 A fully specified and transformed GalSim object.
280 """
281 # Make a 'thin' box to mimic a line surface brightness profile of default
282 # thickness = 1e-6 (i.e., much thinner than a pixel)
283 object = galsim.Box(source_data["trail_length"], trail_thickness)
284 try:
285 object = object.rotate(source_data["beta"] * galsim.degrees)
286 except KeyError:
287 pass
288 object = object.withFlux(inst_flux * source_data["trail_length"]) # type: ignore
289 # GalSim objects are assumed to be in sky-coords. As we want the trail to
290 # appear as defined above in image-coords, we must transform it here.
291 linear_wcs = wcs.linearizePixelToSky(sky_coords, arcseconds)
292 mat = linear_wcs.getMatrix()
293 object = object.transform(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1])
294 # Rescale flux; transformation preserves surface brightness, not flux.
295 object *= 1.0 / np.abs(mat[0, 0] * mat[1, 1] - mat[0, 1] * mat[1, 0])
296 return object # type: ignore
299def make_galsim_stamp(
300 stamp_file: str,
301 fits_alignment: str,
302 wcs: SkyWcs,
303 sky_coords: SpherePoint,
304 inst_flux: float,
305) -> galsim.gsobject.GSObject:
306 """Make a postage stamp with GalSim from a FITS file.
308 Parameters
309 ----------
310 stamp_file : `str`
311 Path to the FITS file containing the postage stamp.
312 fits_alignment : `str`
313 Alignment of the FITS image to the WCS. Allowed values: "wcs", "pixel".
314 wcs : `lsst.afw.geom.SkyWcs`
315 World coordinate system.
316 sky_coords : `lsst.geom.SpherePoint`
317 Sky coordinates of the source.
318 inst_flux : `float`
319 Instrumental flux of the source.
321 Returns
322 -------
323 object: `galsim.gsobject.GSObject`
324 A fully specified and transformed GalSim object.
325 """
326 stamp_file = stamp_file.strip()
327 if os.path.exists(stamp_file):
328 with fits.open(stamp_file) as hdul:
329 hdu_images = [hdu.is_image and hdu.size > 0 for hdu in hdul] # type: ignore
330 if any(hdu_images):
331 stamp_data = galsim.fits.read(stamp_file, read_header=True, hdu=np.where(hdu_images)[0][0])
332 else:
333 raise RuntimeError(f"Cannot find image in input FITS file {stamp_file}.")
334 match fits_alignment:
335 case "wcs":
336 # galsim.fits.read will always attach a WCS to its output.
337 # If it can't find a WCS in the FITS header, then it
338 # defaults to scale = 1.0 arcsec / pix. If that's the scale
339 # then we need to check if it was explicitly set or if it's
340 # just the default. If it's just the default then we should
341 # raise an exception.
342 if is_wcs_galsim_default(stamp_data.wcs, stamp_data.header): # type: ignore
343 raise RuntimeError(f"Cannot find WCS in input FITS file {stamp_file}")
344 case "pixel":
345 # We need to set stamp_data.wcs to the local target WCS.
346 linear_wcs = wcs.linearizePixelToSky(sky_coords, arcseconds)
347 mat = linear_wcs.getMatrix()
348 stamp_data.wcs = galsim.JacobianWCS( # type: ignore
349 mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1]
350 )
351 object = galsim.InterpolatedImage(stamp_data, calculate_stepk=False)
352 object = object.withFlux(inst_flux)
353 return object # type: ignore
354 else:
355 raise RuntimeError(f"Cannot locate input FITS postage stamp {stamp_file}.")
358def is_wcs_galsim_default(
359 wcs: galsim.fitswcs.GSFitsWCS,
360 hdr: galsim.fits.FitsHeader,
361) -> bool:
362 """Decide if wcs = galsim.PixelScale(1.0) is explicitly present in header,
363 or if it's just the GalSim default.
365 Parameters
366 ----------
367 wcs : galsim.fitswcs.GSFitsWCS
368 Potentially default WCS.
369 hdr : galsim.fits.FitsHeader
370 Header as read in by GalSim.
372 Returns
373 -------
374 is_default : bool
375 True if default, False if explicitly set in header.
376 """
377 if wcs != galsim.PixelScale(1.0):
378 return False
379 if hdr.get("GS_WCS") is not None:
380 return False
381 if hdr.get("CTYPE1", "LINEAR") == "LINEAR":
382 return not any(k in hdr for k in ["CD1_1", "CDELT1"])
383 for wcs_type in galsim.fitswcs.fits_wcs_types:
384 # If one of these succeeds, then assume result is explicit.
385 try:
386 wcs_type._readHeader(hdr)
387 return False
388 except Exception:
389 pass
390 else:
391 return not any(k in hdr for k in ["CD1_1", "CDELT1"])
394def infer_gain_from_image(
395 exposure: ExposureF,
396 bad_mask_names: list[str] | None = None,
397 logger: Any | None = None,
398) -> float:
399 """Fit a straight line to image flux vs. estimated variance in each pixel
400 of an exposure, and from that infer an estimated average gain.
402 Parameters
403 ----------
404 exposure : `lsst.afw.image.ExposureF`
405 The exposure to inject synthetic sources into.
406 bad_mask_names : `list` [`str`], optional
407 Names of any mask plane bits that could indicate a pathological flux
408 vs. variance relationship.
409 logger : `lsst.utils.logging.LsstLogAdapter`, optional
410 Logger to use for logging messages.
412 Returns
413 -------
414 gain : `float`
415 The estimated average gain across all pixels.
416 """
417 image_arr = exposure.image.array
418 var_arr = exposure.variance.array
419 good = np.isfinite(image_arr) & np.isfinite(var_arr)
421 mask = exposure.mask
422 # Only look at mask names in the mask_plane_dict.
423 # Otherwise, getPlaneBitMask throws an exception.
424 mask_plane_dict = mask.getMaskPlaneDict()
425 if bad_mask_names is None:
426 bad_mask_names = []
427 missing = set(bad_mask_names) - set(mask_plane_dict)
428 if missing and logger:
429 logger.warning("Mask planes not found in exposure (skipping): %s", missing)
430 bad_mask_names = [name for name in bad_mask_names if name in mask_plane_dict]
431 if len(bad_mask_names) > 0:
432 bad_mask_bit_mask = mask.getPlaneBitMask(bad_mask_names)
433 good &= (mask.array.astype(int) & bad_mask_bit_mask) == 0
435 # Check for sufficient points to fit.
436 n_good = int(np.count_nonzero(good))
437 if n_good < 2:
438 if logger:
439 logger.warning("Too few valid pixels (%d) to infer gain from image; returning NaN.", n_good)
440 return float("nan")
442 fit = np.polyfit(image_arr[good], var_arr[good], deg=1)
443 slope = fit[0]
444 # The gain is 1 / slope, so a non-positive or non-finite slope gives a
445 # non-physical gain. Flag it as NaN rather than propagate a negative,
446 # zero, or infinite gain into the noise model.
447 if not np.isfinite(slope) or slope <= 0: 447 ↛ 448line 447 didn't jump to line 448 because the condition on line 447 was never true
448 if logger:
449 logger.warning("Inferred a non-physical variance-vs-flux slope (%g); returning NaN gain.", slope)
450 return float("nan")
451 return 1.0 / slope
454def get_gain_map(
455 exposure: ExposureF,
456 bad_mask_names: list[str] | None = None,
457 logger: Any | None = None,
458) -> galsim.Image:
459 """Retrieve gains for each pixel.
461 Parameters
462 ----------
463 exposure : `lsst.afw.image.ExposureF`
464 The exposure to inject synthetic sources into.
465 bad_mask_names : `list` [`str`], optional
466 Names of any mask plane bits that could indicate a pathological flux
467 vs. variance relationship.
468 logger : `lsst.utils.logging.LsstLogAdapter`, optional
469 Logger to use for logging messages.
471 Returns
472 -------
473 gain_map : `galsim.Image`
474 The gain to use in each pixel.
475 """
476 full_bbox = exposure.getBBox()
477 try:
478 amplifiers = exposure.getDetector().getAmplifiers()
479 bboxes = [amplifier.getBBox() for amplifier in amplifiers] if amplifiers else None
480 except AttributeError:
481 bboxes = None
482 if not bboxes: 482 ↛ 485line 482 didn't jump to line 485 because the condition on line 482 was always true
483 bboxes = [full_bbox]
485 gains = np.array(
486 [infer_gain_from_image(exposure.subset(bbox), bad_mask_names, logger) for bbox in bboxes],
487 dtype=float,
488 )
490 # Catch unphysical values and replace with fallback values
491 usable = np.isfinite(gains) & (gains > 0)
492 if not np.all(usable):
493 if np.any(usable): 493 ↛ 496line 493 didn't jump to line 496 because the condition on line 493 was never true
494 # Use the median of the regions that did fit (e.g. the other
495 # amplifiers of the same detector).
496 fallback = float(np.median(gains[usable]))
497 else:
498 # Nothing fit. Try a single fit over the whole exposure (which can
499 # succeed even when individual sub-regions do not), then fall back
500 # to unit gain as a last resort.
501 fallback = 1.0
502 if len(bboxes) > 1: 502 ↛ 503line 502 didn't jump to line 503 because the condition on line 502 was never true
503 whole = infer_gain_from_image(exposure, bad_mask_names, logger)
504 if np.isfinite(whole) and whole > 0:
505 fallback = whole
506 if logger: 506 ↛ 507line 506 didn't jump to line 507 because the condition on line 506 was never true
507 logger.warning(
508 "Could not infer gain for %d of %d region(s); using fallback gain %g.",
509 int(np.count_nonzero(~usable)),
510 len(gains),
511 fallback,
512 )
513 gains[~usable] = fallback
515 full_bounds = galsim.BoundsI(full_bbox.minX, full_bbox.maxX, full_bbox.minY, full_bbox.maxY)
516 gain_map = galsim.Image(full_bounds, dtype=float)
517 for bbox, gain in zip(bboxes, gains):
518 bounds = galsim.BoundsI(bbox.minX, bbox.maxX, bbox.minY, bbox.maxY)
519 gain_map[bounds].array[:] = 1.0 / gain
521 return gain_map
524def add_noise_to_galsim_image(
525 image_template: galsim.Image,
526 var_template: galsim.Image,
527 is_coadd: bool,
528 gain_map: galsim.Image,
529 object_common_bounds: galsim.BoundsI,
530 noise_seed: int | None = None,
531 logger: Any | None = None,
532) -> tuple[galsim.Image, galsim.Image]:
533 """Add noise to a supplied image, and adjust the estimated variance
534 accordingly. If the image represents a coadd, add Gaussian noise.
535 Otherwise add Poisson shot noise.
537 Parameters
538 ----------
539 image_template : `galsim.Image`
540 The image to add noise to.
541 var_template : `galsim.Image`
542 The variance to use for the noise generating process in each pixel.
543 is_coadd : `bool`
544 Whether the exposure is a coadd. If True, Gaussian noise is used;
545 if False, Poisson shot noise is used.
546 gain_map : `galsim.Image`
547 The gain to use in each pixel of the full exposure.
548 object_common_bounds : `galsim.BoundsI`
549 Specific region of the full exposure to which this image belongs.
550 noise_seed : `int`, optional
551 Random number generator seed to use for the noise realization.
552 logger : `lsst.utils.logging.LsstLogAdapter`, optional
553 Logger to use for logging messages.
555 Returns
556 -------
557 image_template : `galsim.Image`
558 The image with noise added.
559 var_template : `galsim.Image`
560 The estimated variance, adjusted in response to the added noise.
561 """
562 noise_template = var_template.copy()
563 negative_variance = var_template.array < 0
564 nonfinite_variance = ~np.isfinite(var_template.array)
565 if np.any(negative_variance):
566 if logger:
567 logger.debug("Setting negative-variance pixels to 0 variance for noise generation.")
568 noise_template.array[negative_variance] = 0
569 if np.any(nonfinite_variance):
570 if logger:
571 logger.debug("Setting non-finite-variance pixels to 0 variance for noise generation.")
572 noise_template.array[nonfinite_variance] = 0
574 rng = galsim.BaseDeviate(noise_seed)
576 if not is_coadd:
577 # Treat noise_template, which is the expected amount of
578 # injected flux divided by the gain, as the Poisson mean of
579 # the number of photons collected by each pixel.
580 noise = galsim.PoissonNoise(rng)
581 noise_template.addNoise(noise)
582 # Scale the photon count by the gain to get the amount of
583 # image flux to inject.
584 image_template = noise_template.copy()
585 image_template /= gain_map[object_common_bounds]
586 # Restore the original variance values
587 # of the bad-variance pixels.
588 bad_var = negative_variance | nonfinite_variance
589 noise_template.array[bad_var] = var_template.array[bad_var]
590 var_template = noise_template
591 else:
592 noise = galsim.VariableGaussianNoise(rng, noise_template)
593 image_template.addNoise(noise)
594 var_template = image_template.copy()
595 var_template *= gain_map[object_common_bounds]
597 return image_template, var_template
600def inject_galsim_objects_into_exposure(
601 exposure: ExposureF,
602 objects: Generator[tuple[SpherePoint, Point2D, int, galsim.gsobject.GSObject], None, None],
603 mask_plane_name: str = "INJECTED",
604 calib_flux_radius: float = 12.0,
605 draw_size_max: int = 1000,
606 inject_variance: bool = True,
607 add_noise: bool = True,
608 noise_seed: int = 0,
609 bad_mask_names: list[str] | None = None,
610 logger: Any | None = None,
611) -> tuple[list[int], list[galsim.BoundsI], list[bool], list[bool]]:
612 """Inject sources into given exposure using GalSim.
614 Parameters
615 ----------
616 exposure : `lsst.afw.image.ExposureF`
617 The exposure to inject synthetic sources into.
618 objects : `Generator` [`tuple`, None, None]
619 An iterator of tuples that contains (or generates) locations and object
620 surface brightness profiles to inject. The tuples should contain the
621 following elements: `lsst.geom.SpherePoint`, `lsst.geom.Point2D`,
622 `int`, `galsim.gsobject.GSObject`.
623 mask_plane_name : `str`
624 Name of the mask plane to use for the injected sources.
625 calib_flux_radius : `float`
626 Radius in pixels to use for the flux calibration. This is used to
627 produce the correct instrumental fluxes within the radius. The value
628 should match that of the field defined in slot_CalibFlux_instFlux.
629 draw_size_max : `int`
630 Maximum allowed size of the drawn object. If the object is larger than
631 this, the draw size will be clipped to this size.
632 inject_variance : `bool`
633 Whether, when injecting flux into the image plane, to inject a
634 corresponding amount of variance into the variance plane.
635 add_noise : `bool`
636 Whether to randomly vary the amount of injected image flux in each
637 pixel by an amount consistent with the amount of injected variance.
638 noise_seed : `int`
639 Seed for generating the noise of the first injected object. The seed
640 actually used increases by 1 for each subsequent object, to ensure
641 independent noise realizations.
642 bad_mask_names : `list[str]`, optional
643 List of mask plane names indicating pixels to ignore when fitting flux
644 vs variance in preparation for variance plane modification. If None,
645 then the all pixels are used.
646 logger : `lsst.utils.logging.LsstLogAdapter`, optional
647 Logger to use for logging messages.
649 Returns
650 -------
651 draw_sizes : `list` [`int`]
652 Draw sizes of the injected sources.
653 common_bounds : `list` [`galsim.BoundsI`]
654 Common bounds of the drawn objects.
655 fft_size_errors : `list` [`bool`]
656 Boolean flags indicating whether a GalSimFFTSizeError was raised.
657 psf_compute_errors : `list` [`bool`]
658 Boolean flags indicating whether a PSF computation error was raised.
659 """
660 exposure.mask.addMaskPlane(mask_plane_name)
661 mask_plane_core_name = mask_plane_name + "_CORE"
662 exposure.mask.addMaskPlane(mask_plane_core_name)
663 if logger: 663 ↛ 664line 663 didn't jump to line 664 because the condition on line 663 was never true
664 logger.info(
665 "Adding %s and %s mask planes to the exposure.",
666 mask_plane_name,
667 mask_plane_core_name,
668 )
669 psf = exposure.getPsf()
670 wcs = exposure.getWcs()
671 bbox = exposure.getBBox()
672 full_bounds = galsim.BoundsI(bbox.minX, bbox.maxX, bbox.minY, bbox.maxY)
673 galsim_image = galsim.Image(exposure.image.array, bounds=full_bounds)
674 galsim_variance = galsim.Image(exposure.variance.array, bounds=full_bounds)
675 pixel_scale = wcs.getPixelScale(bbox.getCenter()).asArcseconds()
677 gain_map = get_gain_map(exposure, bad_mask_names, logger)
678 is_coadd = exposure.info.getCoaddInputs() is not None
680 draw_sizes: list[int] = []
681 common_bounds: list[galsim.BoundsI] = []
682 fft_size_errors: list[bool] = []
683 psf_compute_errors: list[bool] = []
684 for i, (sky_coords, pixel_coords, draw_size, object) in enumerate(objects):
685 # Instantiate default returns in case of early exit from this loop.
686 draw_sizes.append(0)
687 common_bounds.append(galsim.BoundsI())
688 fft_size_errors.append(False)
689 psf_compute_errors.append(False)
691 # Get spatial coordinates and WCS.
692 posd = galsim.PositionD(pixel_coords.x, pixel_coords.y)
693 posi = galsim.PositionI(pixel_coords.x // 1, pixel_coords.y // 1)
694 if logger: 694 ↛ 695line 694 didn't jump to line 695 because the condition on line 694 was never true
695 logger.debug(f"Injecting synthetic source at {pixel_coords}.")
696 mat = wcs.linearizePixelToSky(sky_coords, arcseconds).getMatrix()
697 galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1])
699 # This check is here because sometimes the WCS is multivalued and
700 # objects that should not be included were being included.
701 galsim_pixel_scale = np.sqrt(galsim_wcs.pixelArea())
702 if galsim_pixel_scale < pixel_scale / 2 or galsim_pixel_scale > pixel_scale * 2: 702 ↛ 703line 702 didn't jump to line 703 because the condition on line 702 was never true
703 continue
705 # Compute the PSF at the object location.
706 try:
707 psf_array = psf.computeKernelImage(pixel_coords).array
708 except InvalidParameterError:
709 # Try mapping to nearest point contained in bbox.
710 contained_point = Point2D(
711 np.clip(pixel_coords.x, bbox.minX, bbox.maxX), np.clip(pixel_coords.y, bbox.minY, bbox.maxY)
712 )
713 if pixel_coords == contained_point: # no difference, so skip immediately
714 psf_compute_errors[i] = True
715 if logger:
716 logger.debug("Cannot compute PSF for object at %s; flagging and skipping.", sky_coords)
717 continue
718 # Otherwise, try again with new point.
719 try:
720 psf_array = psf.computeKernelImage(contained_point).array
721 except InvalidParameterError:
722 psf_compute_errors[i] = True
723 if logger:
724 logger.debug("Cannot compute PSF for object at %s; flagging and skipping.", sky_coords)
725 continue
727 # Compute the aperture corrected PSF interpolated image.
728 aperture_correction = psf.computeApertureFlux(calib_flux_radius, psf.getAveragePosition())
729 psf_array /= aperture_correction
730 galsim_psf = galsim.InterpolatedImage(galsim.Image(psf_array), wcs=galsim_wcs)
732 # Convolve the object with the PSF and generate draw size.
733 conv = galsim.Convolve(object, galsim_psf)
734 if draw_size == 0: 734 ↛ 736line 734 didn't jump to line 736 because the condition on line 734 was always true
735 draw_size = conv.getGoodImageSize(galsim_wcs.minLinearScale()) # type: ignore
736 injection_draw_size = int(draw_size)
737 injection_core_size = 3
738 if draw_size_max > 0 and injection_draw_size > draw_size_max: 738 ↛ 739line 738 didn't jump to line 739 because the condition on line 738 was never true
739 if logger:
740 logger.warning(
741 "Clipping draw size for object at %s from %d to %d pixels.",
742 sky_coords,
743 injection_draw_size,
744 draw_size_max,
745 )
746 injection_draw_size = draw_size_max
747 draw_sizes[i] = injection_draw_size
748 if injection_core_size > injection_draw_size: 748 ↛ 749line 748 didn't jump to line 749 because the condition on line 748 was never true
749 if logger:
750 logger.debug(
751 "Clipping core size for object at %s from %d to %d pixels.",
752 sky_coords,
753 injection_core_size,
754 injection_draw_size,
755 )
756 injection_core_size = injection_draw_size
757 sub_bounds = galsim.BoundsI(posi).withBorder(injection_draw_size // 2)
758 object_common_bounds = full_bounds & sub_bounds
759 common_bounds[i] = object_common_bounds # type: ignore
761 # Inject the source if there is any overlap.
762 if object_common_bounds.area() > 0: 762 ↛ 829line 762 didn't jump to line 829 because the condition on line 762 was always true
763 common_image = galsim_image[object_common_bounds]
764 common_variance = galsim_variance[object_common_bounds]
766 offset = posd - object_common_bounds.true_center
767 # Note, for preliminary_visit_image injection, pixel is already
768 # part of the PSF and for coadd injection, it's incorrect to
769 # include the output pixel.
770 # So for both cases, we draw using method='no_pixel'.
771 image_template = common_image.copy()
772 try:
773 image_template = conv.drawImage(
774 image_template, add_to_image=False, offset=offset, wcs=galsim_wcs, method="no_pixel"
775 )
776 except GalSimFFTSizeError as err:
777 fft_size_errors[i] = True
778 if logger:
779 logger.debug(
780 "GalSimFFTSizeError raised for object at %s; flagging and skipping.\n%s",
781 sky_coords,
782 err,
783 )
784 continue
786 # The amount of additional variance to inject,
787 # if inject_variance is true.
788 var_template = image_template.copy()
789 var_template *= gain_map[object_common_bounds]
791 if add_noise: 791 ↛ 792line 791 didn't jump to line 792 because the condition on line 791 was never true
792 image_template, var_template = add_noise_to_galsim_image(
793 image_template,
794 var_template,
795 is_coadd,
796 gain_map,
797 object_common_bounds,
798 noise_seed,
799 logger,
800 )
802 common_image += image_template
803 if inject_variance: 803 ↛ 808line 803 didn't jump to line 808 because the condition on line 803 was always true
804 common_variance += var_template
806 # Increment the seed so different noise is generated for different
807 # objects.
808 noise_seed += 1
810 common_box = Box2I(
811 Point2I(object_common_bounds.xmin, object_common_bounds.ymin),
812 Point2I(object_common_bounds.xmax, object_common_bounds.ymax),
813 )
814 bitvalue = exposure.mask.getPlaneBitMask(mask_plane_name)
815 exposure[common_box].mask.array |= bitvalue
816 # Add a 3 x 3 pixel mask centered on the object. The mask must be
817 # large enough to always identify the core/peak of the injected
818 # source, but small enough that it rarely overlaps real sources.
819 sub_bounds_core = galsim.BoundsI(posi).withBorder(injection_core_size // 2)
820 object_common_bounds_core = full_bounds & sub_bounds_core
821 if object_common_bounds_core.area() > 0: 821 ↛ 684line 821 didn't jump to line 684 because the condition on line 821 was always true
822 common_box_core = Box2I(
823 Point2I(object_common_bounds_core.xmin, object_common_bounds_core.ymin),
824 Point2I(object_common_bounds_core.xmax, object_common_bounds_core.ymax),
825 )
826 bitvalue_core = exposure.mask.getPlaneBitMask(mask_plane_core_name)
827 exposure[common_box_core].mask.array |= bitvalue_core
828 else:
829 if logger:
830 logger.debug("No area overlap for object at %s; flagging and skipping.", sky_coords)
832 return draw_sizes, common_bounds, fft_size_errors, psf_compute_errors