Add quantization noise + remove sky_image quantization

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

euclidlike_imsim/noise.py

@@ -14,14 +14,15 @@
 dark_current = 0.
 
 cfg_noise_key = [
-    "ignore_noise",
+    "use_noise",
     "stray_light",
     "thermal_background",
     "reciprocity_failure",
     "dark_current",
     "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."

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