diff --git a/config/was.yaml b/config/was.yaml index a046d29..e165ade 100644 --- a/config/was.yaml +++ b/config/was.yaml @@ -85,15 +85,17 @@ image: # is even more pronounced in this case. draw_method: 'auto' - # These are all by default turned on, but you can turn any of them off if desired: - ignore_noise: False - stray_light: False - thermal_background: False - reciprocity_failure: False - dark_current: False - nonlinearity: False - ipc: False + # If "use_noise" is False, none of the noise sources will be added to the image. + # If "use_noise" is True, it will add sky background and associated noise. + use_noise: True + # stray_light: False # Not yet implemented + # thermal_background: False # Not yet implemented + # reciprocity_failure: False # Not yet implemented + # dark_current: False # Not yet implemented + # nonlinearity: False # Not yet implemented + # ipc: False # Not yet implemented read_noise: True + quantization_noise: True sky_subtract: False # nobjects: 1 diff --git a/euclidlike_imsim/noise.py b/euclidlike_imsim/noise.py index b6c5218..d057141 100644 --- a/euclidlike_imsim/noise.py +++ b/euclidlike_imsim/noise.py @@ -14,7 +14,7 @@ dark_current = 0. cfg_noise_key = [ - "ignore_noise", + "use_noise", "stray_light", "thermal_background", "reciprocity_failure", @@ -22,6 +22,7 @@ "nonlinearity", "ipc", "read_noise", + "quantization_noise", "sky_subtract", ] @@ -39,6 +40,19 @@ def parse_noise_config(params): def get_noise(cfg_noise, cfg_image, base, logger): + """ + Generate and apply noise to an image based on the provided configuration. + + Parameters: + cfg_noise (dict): Configuration dictionary for noise parameters. + + NOTE on the quantization: This prevents problems due to the rounding. + For example, given the very low amplitude of the sky background, we can have + spatial variations of only ~1 ADU. This will be impossible to properly pick + up by tools like SExtractor. Adding this noise prevents this problem and does not + change the signal in the image. This is discussed in [Cuillandre et al. 2025](https://arxiv.org/abs/2405.13496) + Sect. 4.2.7. + """ noise_img = base["current_image"].copy() noise_img.fill(0) @@ -97,8 +111,16 @@ def get_noise(cfg_noise, cfg_image, base, logger): # Make integer ADU now. noise_img.quantize() + # Add quantization noise (this avoid issues related to rounding e.g. + # background estimation) + if cfg_noise["quantization_noise"]: + quantization_noise = noise_img.copy() + quantization_noise.fill(0) + quantization_noise.addNoise(galsim.DeviateNoise(galsim.UniformDeviate(rng))) + noise_img += quantization_noise + noise_img -= 0.5 + sky_image /= gain - sky_image.quantize() base["noise_image"] = noise_img base["sky_image"] = sky_image @@ -142,7 +164,7 @@ def initialize(self, data, scratch, config, base, logger): self.final_data = None def _check_input(self): - if self.cfg_noise["ignore_noise"]: + if not self.cfg_noise["use_noise"]: raise GalSimConfigError( "You cannot ignore the noise and request the noise image at the same time." " Either active the noise or remove the output noise image." @@ -181,7 +203,7 @@ def processImage(self, index, obj_nums, config, base, logger): class SkyImageBuilder(NoiseImageBuilder): def _check_input(self): - if self.cfg_noise["ignore_noise"]: + if not self.cfg_noise["use_noise"]: raise GalSimConfigError( "You cannot ignore the noise and request the sky image at the same time." " Either activate the noise or remove the output sky image." @@ -213,7 +235,7 @@ def processImage(self, index, obj_nums, config, base, logger): class WeightImageBuilder(NoiseImageBuilder): def _check_input(self): - if self.cfg_noise["ignore_noise"]: + if not self.cfg_noise["use_noise"]: raise GalSimConfigError( "You cannot ignore the noise and request the weight image at the same time." " Either activate the noise or remove the output sky image." diff --git a/examples/end_to_end_demo.py b/examples/end_to_end_demo.py index 64aea0e..7f82161 100644 --- a/examples/end_to_end_demo.py +++ b/examples/end_to_end_demo.py @@ -395,10 +395,21 @@ def main(argv): # Finally, the analog-to-digital converter reads in an integer value. full_image.quantize() - sky_image.quantize() + # sky_image.quantize() # Quantizing the background will actually lead to some issues. # Note that the image type after this step is still a float. If we want to actually # get integer values, we can do new_img = galsim.Image(full_image, dtype=int) + # Here we add quantization noise to the image. This prevent to have problems due + # to the rounding. For example, given the very low amplitude of the sky background, + # we can have spatial variations of only ~1 ADU. This will be impossible to + # properly pick up by tools like SExtractor. Adding this noise prevents this problem + # and does not change the signal in the image. + quantization_noise = full_image.copy() + quantization_noise.fill(0) + quantization_noise.addNoise(galsim.DeviateNoise(galsim.UniformDeviate(image_rng))) + full_image += quantization_noise + full_image -= 0.5 + # Since many people are used to viewing background-subtracted images, we provide a # version with the background subtracted (also rounding that to an int). full_image -= sky_image