Add initial_image input type

CHANGELOG.rst

@@ -18,6 +18,8 @@ Config Updates
 
 - Fixed a bug in Scattered type, where it didn't respect world_pos being specified in the
   stamp field, as it should.  (#1190)
+- Added a new ``initial_image`` input type that lets you read in an existing image file
+  and draw onto it. (#1237)
 
 
 New Features

galsim/config/__init__.py

@@ -45,6 +45,7 @@
 from . import input_cosmos
 from . import input_nfw
 from . import input_powerspectrum
+from . import input_image
 
 from . import extra_psf
 from . import extra_weight

galsim/config/image.py

@@ -26,7 +26,7 @@
 from .wcs import BuildWCS
 from .sensor import BuildSensor
 from .bandpass import BuildBandpass
-from .stamp import BuildStamp, MakeStampTasks
+from .stamp import BuildStamp, MakeStampTasks, ParseDType
 from ..errors import GalSimConfigError, GalSimConfigValueError
 from ..position import PositionI, PositionD
 from ..bounds import BoundsI
@@ -181,6 +181,7 @@ def SetupConfigImageSize(config, xsize, ysize, logger=None):
     - Set config['wcs'] to be the built wcs
     - If wcs.isPixelScale(), also set config['pixel_scale'] for convenience.
     - Set config['world_center'] to either a given value or based on wcs and image_center
+    - Create a blank image if possible and store as config['current_image']
 
     Parameters:
         config:     The configuration dict.
@@ -206,7 +207,8 @@ def SetupConfigImageSize(config, xsize, ysize, logger=None):
 
     config['image_origin'] = PositionI(origin,origin)
     config['image_center'] = PositionD( origin + (xsize-1.)/2., origin + (ysize-1.)/2. )
-    config['image_bounds'] = BoundsI(origin, origin+xsize-1, origin, origin+ysize-1)
+    bounds = BoundsI(origin, origin+xsize-1, origin, origin+ysize-1)
+    config['image_bounds'] = bounds
 
     # Build the wcs
     wcs = BuildWCS(image, 'wcs', config, logger)
@@ -225,6 +227,12 @@ def SetupConfigImageSize(config, xsize, ysize, logger=None):
     else:
         config['world_center'] = wcs.toWorld(config['image_center'])
 
+    if bounds.isDefined():
+        dtype = ParseDType(image, config)
+        config['current_image'] = Image(bounds=bounds, dtype=dtype, wcs=wcs, init_value=0)
+    else:
+        config['current_image'] = None
+
 
 # Ignore these when parsing the parameters for specific Image types:
 from .stamp import stamp_image_keys
@@ -510,14 +518,18 @@ def buildImage(self, config, base, image_num, obj_num, logger):
         ysize = base['image_ysize']
         logger.debug('image %d: Single Image: size = %s, %s',image_num,xsize,ysize)
 
-        # In case there was one set from before, we don't want to confuse the stamp builder
-        # thinking that this is the full image onto which we are drawing this object.
-        base['current_image'] = None
-
         image, current_var = BuildStamp(
                 base, obj_num=obj_num, xsize=xsize, ysize=ysize, do_noise=True, logger=logger)
         if image is not None:
             image.wcs = base['wcs']   # in case stamp has a local jacobian.
+
+        current_image = base['current_image']
+        if current_image is not None and image is not None:
+            b = current_image.bounds & image.bounds
+            if b.isDefined():
+                current_image[b] += image[b]
+            image = current_image
+
         return image, current_var
 
     def makeTasks(self, config, base, jobs, logger):

galsim/config/image_scattered.py

@@ -21,7 +21,7 @@
 
 from .image import ImageBuilder, FlattenNoiseVariance, RegisterImageType
 from .value import ParseValue, GetAllParams
-from .stamp import BuildStamps, _ParseDType
+from .stamp import BuildStamps
 from .noise import AddSky, AddNoise
 from .input import ProcessInputNObjects
 from ..errors import GalSimConfigError, GalSimConfigValueError
@@ -92,16 +92,7 @@ def buildImage(self, config, base, image_num, obj_num, logger):
         Returns:
             the final image and the current noise variance in the image as a tuple
         """
-        full_xsize = base['image_xsize']
-        full_ysize = base['image_ysize']
-        wcs = base['wcs']
-
-        dtype = _ParseDType(config, base)
-        full_image = Image(full_xsize, full_ysize, dtype=dtype)
-        full_image.setOrigin(base['image_origin'])
-        full_image.wcs = wcs
-        full_image.setZero()
-        base['current_image'] = full_image
+        full_image = base['current_image']
 
         if 'image_pos' in config and 'world_pos' in config:
             raise GalSimConfigValueError(
@@ -111,6 +102,8 @@ def buildImage(self, config, base, image_num, obj_num, logger):
         if ('image_pos' not in config and 'world_pos' not in config and
                 not ('stamp' in base and
                     ('image_pos' in base['stamp'] or 'world_pos' in base['stamp']))):
+            full_xsize = base['image_xsize']
+            full_ysize = base['image_ysize']
             xmin = base['image_origin'].x
             xmax = xmin + full_xsize-1
             ymin = base['image_origin'].y

galsim/config/image_tiled.py

@@ -22,7 +22,7 @@
 from .image import ImageBuilder, FlattenNoiseVariance, RegisterImageType
 from .util import GetRNG
 from .value import ParseValue, GetAllParams
-from .stamp import BuildStamps, _ParseDType
+from .stamp import BuildStamps
 from .noise import AddSky, AddNoise
 from ..errors import GalSimConfigError, GalSimConfigValueError
 from ..image import Image
@@ -115,16 +115,7 @@ def buildImage(self, config, base, image_num, obj_num, logger):
         Returns:
             the final image and the current noise variance in the image as a tuple
         """
-        full_xsize = base['image_xsize']
-        full_ysize = base['image_ysize']
-        wcs = base['wcs']
-
-        dtype = _ParseDType(config, base)
-        full_image = Image(full_xsize, full_ysize, dtype=dtype)
-        full_image.setOrigin(base['image_origin'])
-        full_image.wcs = wcs
-        full_image.setZero()
-        base['current_image'] = full_image
+        full_image = base['current_image']
 
         nobjects = self.nx_tiles * self.ny_tiles
 

galsim/config/input_image.py

@@ -0,0 +1,63 @@
+# Copyright (c) 2012-2022 by the GalSim developers team on GitHub
+# https://github.com/GalSim-developers
+#
+# This file is part of GalSim: The modular galaxy image simulation toolkit.
+# https://github.com/GalSim-developers/GalSim
+#
+# GalSim is free software: redistribution and use in source and binary forms,
+# with or without modification, are permitted provided that the following
+# conditions are met:
+#
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions, and the disclaimer given in the accompanying LICENSE
+#    file.
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions, and the disclaimer given in the documentation
+#    and/or other materials provided with the distribution.
+#
+
+
+from .input import InputLoader, RegisterInputType
+from .value import GetAllParams
+from ..fits import read
+
+# This file adds input type initial_image.
+
+class InitialImageLoader(InputLoader):
+
+    def __init__(self):
+        self.init_func = read
+        self.has_nobj = False
+        self.file_scope = False
+        self.takes_logger = False
+        self.use_proxy = False
+
+    def getKwargs(self, config, base, logger):
+        """Parse the config dict and return the kwargs needed to build the PowerSpectrum object.
+
+        Parameters:
+            config:     The configuration dict for 'power_spectrum'
+            base:       The base configuration dict
+            logger:     If given, a logger object to log progress.
+
+        Returns:
+            kwargs, safe
+        """
+        req = { 'file_name': str }
+        opt = { 'dir': str, 'read_header': bool }
+        return GetAllParams(config, base, req=req, opt=opt)
+
+    def setupImage(self, input_obj, config, base, logger=None):
+        """Set up the PowerSpectrum input object's gridded values based on the
+        size of the image and the grid spacing.
+
+        Parameters:
+            input_obj:  The PowerSpectrum object to use
+            config:     The configuration dict for 'power_spectrum'
+            base:       The base configuration dict.
+            logger:     If given, a logger object to log progress.
+        """
+        base['current_image'] = input_obj.copy()
+
+# Register this as a valid input type
+RegisterInputType('initial_image', InitialImageLoader())

galsim/config/stamp.py

@@ -525,7 +525,7 @@ def DrawBasic(prof, image, method, offset, config, base, logger, **kwargs):
         kwargs['sensor'] = sensor
 
     if image is None:
-        kwargs['dtype'] = _ParseDType(config, base)
+        kwargs['dtype'] = ParseDType(config, base)
 
     if logger.isEnabledFor(logging.DEBUG):
         # Don't output the full image array.  Use str(image) for that kwarg.  And Bandpass.
@@ -584,7 +584,7 @@ def ParseWorldPos(config, param_name, base, logger):
     else:
         return ParseValue(config, param_name, base, PositionD)[0]
 
-def _ParseDType(config, base):
+def ParseDType(config, base):
     dtype = config.get('dtype', None)
     if isinstance(dtype, str):
         try:
@@ -637,25 +637,25 @@ def setup(self, config, base, xsize, ysize, ignore, logger):
 
         # Update the size if necessary
         image = base['image']
-        if not xsize:
-            if 'xsize' in config:
-                xsize = ParseValue(config,'xsize',base,int)[0]
-            elif 'size' in config:
-                xsize = ParseValue(config,'size',base,int)[0]
-            elif 'stamp_xsize' in image:
-                xsize = ParseValue(image,'stamp_xsize',base,int)[0]
-            elif 'stamp_size' in image:
-                xsize = ParseValue(image,'stamp_size',base,int)[0]
-
-        if not ysize:
-            if 'ysize' in config:
-                ysize = ParseValue(config,'ysize',base,int)[0]
-            elif 'size' in config:
-                ysize = ParseValue(config,'size',base,int)[0]
-            elif 'stamp_ysize' in image:
-                ysize = ParseValue(image,'stamp_ysize',base,int)[0]
-            elif 'stamp_size' in image:
-                ysize = ParseValue(image,'stamp_size',base,int)[0]
+        if 'xsize' in config:
+            xsize = ParseValue(config,'xsize',base,int)[0]
+        elif 'size' in config:
+            xsize = ParseValue(config,'size',base,int)[0]
+        elif 'stamp_xsize' in image:
+            xsize = ParseValue(image,'stamp_xsize',base,int)[0]
+        elif 'stamp_size' in image:
+            xsize = ParseValue(image,'stamp_size',base,int)[0]
+        # else use the input xsize
+
+        if 'ysize' in config:
+            ysize = ParseValue(config,'ysize',base,int)[0]
+        elif 'size' in config:
+            ysize = ParseValue(config,'size',base,int)[0]
+        elif 'stamp_ysize' in image:
+            ysize = ParseValue(image,'stamp_ysize',base,int)[0]
+        elif 'stamp_size' in image:
+            ysize = ParseValue(image,'stamp_size',base,int)[0]
+        # else use the input ysize
 
         # Determine where this object is going to go:
         if 'image_pos' in config:
@@ -891,9 +891,18 @@ def makeStamp(self, config, base, xsize, ysize, logger):
             the image
         """
         if xsize and ysize:
-            dtype = _ParseDType(config, base)
-            im = Image(xsize, ysize, dtype=dtype)
-            im.setZero()
+            dtype = ParseDType(config, base)
+            bounds = _BoundsI(1,xsize,1,ysize)
+
+            # Set the origin appropriately
+            stamp_center = base['stamp_center']
+            if stamp_center:
+                bounds = bounds.shift(stamp_center - bounds.center)
+            else:
+                bounds = bounds.shift(base.get('image_origin',PositionI(1,1)) - PositionI(1,1))
+
+            im = Image(bounds=bounds, dtype=dtype, init_value=0)
+
             return im
         else:
             return None
@@ -958,16 +967,15 @@ def updateSkip(self, prof, image, method, offset, config, base, logger):
                     prof = Convolution(prof, Pixel(1.))
                 N = prof.getGoodImageSize(1.)
                 bounds = _BoundsI(1,N,1,N)
-            else:
-                bounds = image.bounds
 
-            # Set the origin appropriately
-            stamp_center = base['stamp_center']
-            if stamp_center:
-                bounds = bounds.shift(stamp_center - bounds.center)
+                # Set the origin appropriately
+                stamp_center = base['stamp_center']
+                if stamp_center:
+                    bounds = bounds.shift(stamp_center - bounds.center)
+                else:
+                    bounds = bounds.shift(base.get('image_origin',PositionI(1,1)) - PositionI(1,1))
             else:
-                bounds = bounds.shift(base.get('image_origin',PositionI(1,1)) -
-                                      PositionI(bounds.xmin, bounds.ymin))
+                bounds = image.bounds
 
             overlap = bounds & base['current_image'].bounds
             if not overlap.isDefined():