diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 2d7ac497a1..89f674251d 100644 --- a/CHANGELOG.rst +++ b/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 diff --git a/galsim/config/__init__.py b/galsim/config/__init__.py index 38c01e337b..c8a9e659ff 100644 --- a/galsim/config/__init__.py +++ b/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 diff --git a/galsim/config/image.py b/galsim/config/image.py index adbb4d25de..13aa4a2f00 100644 --- a/galsim/config/image.py +++ b/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): diff --git a/galsim/config/image_scattered.py b/galsim/config/image_scattered.py index a5d1d2f6c0..acc731ef4b 100644 --- a/galsim/config/image_scattered.py +++ b/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 diff --git a/galsim/config/image_tiled.py b/galsim/config/image_tiled.py index 28505bb191..5829ba0394 100644 --- a/galsim/config/image_tiled.py +++ b/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 diff --git a/galsim/config/input_image.py b/galsim/config/input_image.py new file mode 100644 index 0000000000..af0c3740db --- /dev/null +++ b/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()) diff --git a/galsim/config/stamp.py b/galsim/config/stamp.py index bb80226cc9..cae78f67fc 100644 --- a/galsim/config/stamp.py +++ b/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(): diff --git a/tests/config_input/sn.yaml b/tests/config_input/sn.yaml new file mode 100644 index 0000000000..ceb2e4ada9 --- /dev/null +++ b/tests/config_input/sn.yaml @@ -0,0 +1,83 @@ +# 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. +# + + +psf: + type: Moffat + beta: 3 + fwhm: 0.6 + ellip: '$galsim.Shear(e1=0.03, e2=-0.04)' + +image: + type: Single + size: 128 + dtype: float + random_seed: + - 12345 + pixel_scale: 0.2 + # No noise, which is unrealistic, but makes the test easier. + +--- + +gal: + type: Exponential + half_light_radius: 2.3 + flux: 3.e6 + ellip: + type: E1E2 + e1: 0.4 + e2: 0.7 + +output: + dir: output + file_name: ref.fits + +--- + +eval_variables: + # Model the time of the exposure as a linear function of the image_num with a 8 day cadence. + # More realistic versions of this could use datetime with observation times from + # a file or something similarly interesting. + ftime: '$(image_num - 2.3) * 7' # in days + fmax_flux: 1.e5 + +input: + initial_image: + dir: output + file_name: ref.fits + +gal: + # The "gal" here is the SN, which is a point source. + type: DeltaFunction + flux: '$0 if time < 0 else np.exp(-time/50) * max_flux' + +stamp: + image_pos: + type: XY + x: 80 + y: 60 + +output: + nfiles: 12 + dir: output + file_name: + type: NumberedFile + root: 'SN_image' + num: '$file_num' + digits: 2 + diff --git a/tests/test_config_image.py b/tests/test_config_image.py index d5cebbf62f..bc03da949e 100644 --- a/tests/test_config_image.py +++ b/tests/test_config_image.py @@ -214,7 +214,7 @@ def test_positions(): #logger.setLevel(logging.DEBUG) gal = galsim.Gaussian(sigma=1.7, flux=100) - im1= gal.drawImage(nx=21, ny=21, scale=1) + im1 = gal.drawImage(nx=21, ny=21, scale=1) im1.setCenter(39,43) im2 = galsim.config.BuildImage(config, logger=logger) @@ -279,6 +279,29 @@ def test_positions(): config['stamp']['world_pos'] = { 'type' : 'Random' } with assert_raises(galsim.GalSimConfigError): galsim.config.BuildImage(config) + del config['stamp']['world_pos'] + + # If the size is set in image, then image_pos will put the object offset in this image. + config['image'] = { + 'type' : 'Single', + 'size' : 64, + } + im11 = galsim.config.BuildImage(config) + assert im11.bounds == galsim.BoundsI(1,64,1,64) + assert im11[im2.bounds] == im2 + + # If the offset is larger, then only the overlap is included + config['image']['size'] = 50 + im12 = galsim.config.BuildImage(config) + assert im12.bounds == galsim.BoundsI(1,50,1,50) + b = im12.bounds & im2.bounds + assert im12[b] == im2[b] + + # If the offset is large enough, none of the stamp is included. + config['image']['size'] = 32 + im13 = galsim.config.BuildImage(config) + assert im13.bounds == galsim.BoundsI(1,32,1,32) + np.testing.assert_array_equal(im13.array, 0.) @timer @@ -3354,6 +3377,41 @@ def test_sensor(): with assert_raises(galsim.GalSimConfigError): galsim.config.BuildSensor(config, 'sensor', config) +@timer +def test_initial_image(): + # This test simulates a time series of a supernova going off near a big galaxy. + # It first makes a reference image of the galaxy alone. + # Then it makes 12 observations of the same scene with a SN. + # The SN appears between the 3rd and 4th images, and then decays exponentially. + + # Note: this can be run from the command line as just `galsim sn.yaml`, which does + # all the parts in one execution. + configs = galsim.config.ReadConfig('config_input/sn.yaml') + + # In Python, we have to run each config dict separately. + # First make the reference image. + galsim.config.Process(configs[0]) + ref_image = galsim.fits.read('output/ref.fits') + + # Now the supernovae. + galsim.config.ProcessInput(configs[1]) + sn_images = galsim.config.BuildImages(12, configs[1]) + + for i, sn_image in enumerate(sn_images): + print(i, sn_image.array.max(), sn_image.array.sum()) + t = (i - 2.3) * 7 # Time in day + if t > 0: + flux = np.exp(-t/50) * 1.e5 + else: + flux = 0. + + # The diff image is basically perfect here, since we didn't add noise. + # The only imprecision is that Moffats have flux off the edge of the stamp. + # It's pretty noticable with beta=2. But here with beta=3, it's pretty slight. + diff_image = sn_image - ref_image + print(i, t, flux, diff_image.array.sum()) + np.testing.assert_allclose(diff_image.array.sum(), flux, rtol=1.e-6) + if __name__ == "__main__": testfns = [v for k, v in vars().items() if k[:5] == 'test_' and callable(v)]