From 9818ff7337d1bf34f65654916536ad9bbf2b86d8 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Tue, 18 Feb 2025 10:09:06 +0000 Subject: [PATCH 01/36] Machinery to write photons outside image bounds as extra output to FITS. --- config/imsim-config-instcat.yaml | 2 +- config/imsim-config-photon-pooling.yaml | 13 +++++++- config/imsim-config-skycat.yaml | 2 +- imsim/__init__.py | 1 + imsim/photon_pooling.py | 13 +++++++- imsim/photon_scattering.py | 42 +++++++++++++++++++++++++ 6 files changed, 69 insertions(+), 4 deletions(-) create mode 100644 imsim/photon_scattering.py diff --git a/config/imsim-config-instcat.yaml b/config/imsim-config-instcat.yaml index 82e44fc6..9f4038e2 100644 --- a/config/imsim-config-instcat.yaml +++ b/config/imsim-config-instcat.yaml @@ -12,7 +12,7 @@ # Get most of the configuration from the base imSim config template. modules: - imsim -template: imsim-config +template: imsim-config-photon-pooling #################################################################### # The following entries are added to the base configuration above. diff --git a/config/imsim-config-photon-pooling.yaml b/config/imsim-config-photon-pooling.yaml index 541813e6..2be4b981 100644 --- a/config/imsim-config-photon-pooling.yaml +++ b/config/imsim-config-photon-pooling.yaml @@ -65,4 +65,15 @@ output.photon_pooling_truth: nominal_flux: "@nominal_flux" # The nominal "expectation value" flux. phot_flux: "@phot_flux" # The realized flux for photon shooting. fft_flux: "@fft_flux" # If FFT rendering, then the flux used, including vignetting. - incident_flux: "@incident_flux" # The sum of the fluxes of all photons from the object. \ No newline at end of file + incident_flux: "@incident_flux" # The sum of the fluxes of all photons from the object. +output.scattered_photons: + dir: output + file_name: + type: FormattedStr + format : scattered_%08d-%1d-%s-%s-det%03d.fits + items: + - { type: OpsimData, field: observationId } + - { type: OpsimData, field: snap } + - $band + - $det_name + - "@output.det_num" \ No newline at end of file diff --git a/config/imsim-config-skycat.yaml b/config/imsim-config-skycat.yaml index 57ba1ab3..ebcba212 100644 --- a/config/imsim-config-skycat.yaml +++ b/config/imsim-config-skycat.yaml @@ -13,7 +13,7 @@ # Get most of the configuration from the base imSim config template. modules: - imsim -template: imsim-config +template: imsim-config-photon-pooling #################################################################### # The following entires are added to the base configuration above. diff --git a/imsim/__init__.py b/imsim/__init__.py index b982ffde..33c6540f 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -40,3 +40,4 @@ from .process_info import * from .table_row import * from .photon_pooling import * +from .photon_scattering import * diff --git a/imsim/photon_pooling.py b/imsim/photon_pooling.py index 04001169..dc890e74 100644 --- a/imsim/photon_pooling.py +++ b/imsim/photon_pooling.py @@ -138,6 +138,8 @@ def buildImage(self, config, base, image_num, _obj_num, logger): local_wcs = base['wcs'].local(full_image.true_center) if sensor is None: sensor = Sensor() + if 'scattered_photons' in base['output']: + base['scattered_photons'] = [] # Initialize the scattered_photons list for batch_num, batch in enumerate(phot_batches, start=current_photon_batch_num): if not batch: continue @@ -157,8 +159,17 @@ def buildImage(self, config, base, image_num, _obj_num, logger): # Now accumulate the photons onto the sensor. Resume is true for all calls but the first. Recalculate the pixel # boundaries on the first subbatch of each full batch. self.accumulate_photons(photons, full_image, sensor, resume=(batch_num > current_photon_batch_num or subbatch_num > 0), recalc=(subbatch_num == 0)) - del photons # As with the stamps above, let the garbage collector know we don't need the photons anymore. + # Gather non-accumulated photons if we're going to be outputting them. + if 'scattered_photons' in base['output']: + scattered_indices = [i for i in range(len(photons)) if not imview.bounds.includes(photons.x[i], photons.y[i])] + if len(scattered_indices) > 0: + scattered_photons = galsim.PhotonArray(len(scattered_indices)) + scattered_photons.copyFrom(photons, target_indices=slice(len(scattered_indices)), source_indices=scattered_indices) + base['scattered_photons'].append(scattered_photons) + + del photons # As with the stamps above, let the garbage collector know we don't need the photons anymore. + # Note: in typical imsim usage, all current_vars will be 0. So this normally doesn't # add much to the checkpointing data. nz_var = np.nonzero(current_vars)[0] diff --git a/imsim/photon_scattering.py b/imsim/photon_scattering.py new file mode 100644 index 00000000..88af3951 --- /dev/null +++ b/imsim/photon_scattering.py @@ -0,0 +1,42 @@ +from galsim.config.extra import ExtraOutputBuilder, RegisterExtraOutput +from galsim import PhotonArray +from astropy.io.fits import BinTableHDU + +# An extra output type to enable two-pass photon scattering. +# This output should be used in the first pass to write scattered photons to file, +# while the second pass will read the scattered photons then go over all the images +# again and add the photons. +class ScatteredPhotonsBuilder(ExtraOutputBuilder): + """Build pickled photon arrays containing photons which were not accumulated + to the sensor during each image's construction. + """ + + def processImage(self, index, obj_nums, config, base, logger): + """After each image is drawn, concatenate the list of scattered photon arrays + and store in data ready to be pickled to file later on. + """ + # If scattered photons have been stored in the base config, concatenate these + # photon arrays to a single one and store in data, indexing by image number. + if 'scattered_photons' in base and len(base['scattered_photons']) > 1: + self.data[index] = PhotonArray.concatenate(base['scattered_photons']) + else: + self.data[index] = PhotonArray(N=0) + + def finalize(self, config, base, main_data, logger): + return self.data + + def writeFile(self, file_name, config, base, logger): + """Write the photon array to fits file. + Might it be faster to write to a pickle file? + """ + for photon_array in self.final_data: + photon_array.write(file_name) + + def writeHdu(self, config, base, logger): + """We don't want to write the scattered photons to FITS, so return an empty + BinTable in lieu of something else. + """ + return BinTableHDU(data=None) + + +RegisterExtraOutput('scattered_photons', ScatteredPhotonsBuilder()) From 9152ec540fc518f1d31e617c5c33454ffed828a1 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Mon, 24 Feb 2025 12:50:57 +0000 Subject: [PATCH 02/36] Put scattered photons in focal plane coords. --- imsim/photon_pooling.py | 2 +- imsim/photon_scattering.py | 15 +++++++++++---- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/imsim/photon_pooling.py b/imsim/photon_pooling.py index dc890e74..5ab74927 100644 --- a/imsim/photon_pooling.py +++ b/imsim/photon_pooling.py @@ -169,7 +169,7 @@ def buildImage(self, config, base, image_num, _obj_num, logger): base['scattered_photons'].append(scattered_photons) del photons # As with the stamps above, let the garbage collector know we don't need the photons anymore. - + # Note: in typical imsim usage, all current_vars will be 0. So this normally doesn't # add much to the checkpointing data. nz_var = np.nonzero(current_vars)[0] diff --git a/imsim/photon_scattering.py b/imsim/photon_scattering.py index 88af3951..8d4d6c85 100644 --- a/imsim/photon_scattering.py +++ b/imsim/photon_scattering.py @@ -2,6 +2,9 @@ from galsim import PhotonArray from astropy.io.fits import BinTableHDU +from .utils import pixel_to_focal +from .camera import get_camera + # An extra output type to enable two-pass photon scattering. # This output should be used in the first pass to write scattered photons to file, # while the second pass will read the scattered photons then go over all the images @@ -19,22 +22,26 @@ def processImage(self, index, obj_nums, config, base, logger): # photon arrays to a single one and store in data, indexing by image number. if 'scattered_photons' in base and len(base['scattered_photons']) > 1: self.data[index] = PhotonArray.concatenate(base['scattered_photons']) + # We need to store the pixels' using focal plane coordinates so they + # can be accumulated in the second pass on an arbitrary sensor. + detector = get_camera(base['output']['camera'])[base['det_name']] + self.data[index].x, self.data[index].y = pixel_to_focal(self.data[index].x, self.data[index].y, detector) else: self.data[index] = PhotonArray(N=0) def finalize(self, config, base, main_data, logger): return self.data - + def writeFile(self, file_name, config, base, logger): """Write the photon array to fits file. - Might it be faster to write to a pickle file? + Might it be faster/more space efficient to pickle? """ for photon_array in self.final_data: photon_array.write(file_name) def writeHdu(self, config, base, logger): - """We don't want to write the scattered photons to FITS, so return an empty - BinTable in lieu of something else. + """We don't want to write the scattered photons to FITS, so return an + empty BinTable in lieu of something else. """ return BinTableHDU(data=None) From 2937c80be2844198084ef5c8fed9fea7fbee1936 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Fri, 9 May 2025 12:46:12 +0100 Subject: [PATCH 03/36] First partial attempt at reading in the photons from file. --- imsim/photon_scattering.py | 107 +++++++++++++++++++++++++++++++++++-- 1 file changed, 102 insertions(+), 5 deletions(-) diff --git a/imsim/photon_scattering.py b/imsim/photon_scattering.py index 8d4d6c85..13dac7d8 100644 --- a/imsim/photon_scattering.py +++ b/imsim/photon_scattering.py @@ -1,8 +1,10 @@ +import galsim from galsim.config.extra import ExtraOutputBuilder, RegisterExtraOutput -from galsim import PhotonArray +from galsim.config import ImageBuilder, RegisterImageType, InputLoader, RegisterInputType from astropy.io.fits import BinTableHDU -from .utils import pixel_to_focal +# from .lsst_image import LSST_ImageBuilderBase +from .utils import pixel_to_focal, focal_to_pixel from .camera import get_camera # An extra output type to enable two-pass photon scattering. @@ -21,13 +23,13 @@ def processImage(self, index, obj_nums, config, base, logger): # If scattered photons have been stored in the base config, concatenate these # photon arrays to a single one and store in data, indexing by image number. if 'scattered_photons' in base and len(base['scattered_photons']) > 1: - self.data[index] = PhotonArray.concatenate(base['scattered_photons']) - # We need to store the pixels' using focal plane coordinates so they + self.data[index] = galsim.PhotonArray.concatenate(base['scattered_photons']) + # We need to store the photons using focal plane coordinates so they # can be accumulated in the second pass on an arbitrary sensor. detector = get_camera(base['output']['camera'])[base['det_name']] self.data[index].x, self.data[index].y = pixel_to_focal(self.data[index].x, self.data[index].y, detector) else: - self.data[index] = PhotonArray(N=0) + self.data[index] = galsim.PhotonArray(N=0) def finalize(self, config, base, main_data, logger): return self.data @@ -46,4 +48,99 @@ def writeHdu(self, config, base, logger): return BinTableHDU(data=None) +class ScatteredPhotonsInput(object): + """A class to read in the scattered photons that were written to file + during an earlier first pass. + """ + + def __init__(self, file_name, wcs, det, xsize=4096, ysize=4096, logger=None): + """ + Initialize the scattered photons input class. + + Parameters: + file_name: str + The name of the file to read. + wcs: galsim.WCS + The WCS object to use for the image. + det: lsst.afw.cameraGeom.Detector + The detector for the current sensor. + xsize: int + The x size in pixels of the CCD. (default: 4096) + ysize: int + The y size in pixels of the CCD. (default: 4096) + logger: logging.logger + A logger object. (default: None) + """ + self.file_name = file_name + self.wcs = wcs + self.det = det + self.xsize = xsize + self.ysize = ysize + self.logger = logger + self._photons = None + + def read_photons(self): + """Read the scattered photons from the file. + """ + # Read the scattered photons from the file then convert them from + # focal plane coordinates to pixel coordinates. + self.photons = galsim.PhotonArray.read(self.file_name) + self.photons.x, self.photons.y = focal_to_pixel(self.photons.x, self.photons.y, self.det) + + @property + def photons(self): + """Get the scattered photons. + """ + if self._photons is None: + self.read_photons() + return self._photons + + +class ScatteredPhotonsLoader(InputLoader): + """ + Class to load scattered photons from file. + """ + def getKwargs(self, config, base, logger): + req = {'file_name': str} + opt = {} + kwargs, safe = galsim.config.GetAllParams(config, base, req=req, + opt=opt) + wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) + kwargs['wcs'] = wcs + kwargs['xsize'] = base.get('det_xsize', 4096) + kwargs['ysize'] = base.get('det_ysize', 4096) + kwargs['logger'] = logger + + # For now assume the scattered photons on file can be used for all detectors. + safe = True + return kwargs, safe + + +# class LSST_ScatteredPhotonsImageBuilder(LSST_ImageBuilderBase): + +# def setup(self, config, base, logger): +# """Set up the scattered photons image type. +# """ +# # We need to set the pixel scale to be the same as the camera's pixel scale +# # so that we can convert between focal plane coordinates and pixel coordinates. +# self.pixel_scale = base['output']['pixel_scale'] +# self.camera = get_camera(base['output']['camera'])[base['det_name']] +# self.focal_plane = self.camera.get_focal_plane() +# self.focal_plane.set_pixel_scale(self.pixel_scale) + +# def draw(self, config, base, logger): +# """Draw the scattered photons to the image. +# """ +# # Get the scattered photons from the base config. +# if 'scattered_photons' in base and len(base['scattered_photons']) > 1: +# self.data = PhotonArray.concatenate(base['scattered_photons']) +# # We need to convert the focal plane coordinates to pixel coordinates +# # so that we can draw them on the image. +# self.data.x, self.data.y = self.focal_plane.focal_to_pixel(self.data.x, self.data.y) +# else: +# self.data = PhotonArray(N=0) + + RegisterExtraOutput('scattered_photons', ScatteredPhotonsBuilder()) +RegisterInputType('scattered_photons', ScatteredPhotonsLoader()) +# RegisterImageType('LSST_ScatteredPhotonsImage', LSST_ScatteredPhotonsImageBuilder) From eaf79d32cfa1683b328877fa1538bc5ad20fccd6 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Fri, 13 Jun 2025 17:13:24 +0100 Subject: [PATCH 04/36] Working scattered photon input --- examples/imsim-user-instcat-ffp-pass1.yaml | 38 +++++++++++++++++++ examples/imsim-user-instcat-ffp-pass2.yaml | 44 ++++++++++++++++++++++ imsim/photon_pooling.py | 7 +++- imsim/photon_scattering.py | 44 +++++++++------------- 4 files changed, 106 insertions(+), 27 deletions(-) create mode 100644 examples/imsim-user-instcat-ffp-pass1.yaml create mode 100644 examples/imsim-user-instcat-ffp-pass2.yaml diff --git a/examples/imsim-user-instcat-ffp-pass1.yaml b/examples/imsim-user-instcat-ffp-pass1.yaml new file mode 100644 index 00000000..2ea985f1 --- /dev/null +++ b/examples/imsim-user-instcat-ffp-pass1.yaml @@ -0,0 +1,38 @@ +# Use imSim custom modules +modules: + - imsim + +# Get most of the configuration from the imSim config-template +# for instance catalogs. +template: imsim-config-instcat + +################################################################ +# Make your changes below. +################################################################ + +# Put your own commands that override the defaults below here. For example +# input.instance_catalog.file_name: ./imsim_cat_197356.txt +# input.instance_catalog.sort_mag: False +# input.tree_rings.only_dets: [R22_S11] +# image.nobjects: 5 + +input.instance_catalog.file_name: $os.environ.get('IMSIM_HOME')+'/imSim/examples/example_instance_catalog_ffp.txt' +input.instance_catalog.sort_mag: False + +input.tree_rings.only_dets: [R22_S11] +input.atm_psf: "" +psf: + type: Convolve + items: + - + type: Gaussian + fwhm: 0.3 + +image.nobjects: 10 + +stamp.fft_sb_thresh: 1.e30 + +output.det_num.first: 94 +output.nfiles: 1 + +output.dir: output diff --git a/examples/imsim-user-instcat-ffp-pass2.yaml b/examples/imsim-user-instcat-ffp-pass2.yaml new file mode 100644 index 00000000..937424ae --- /dev/null +++ b/examples/imsim-user-instcat-ffp-pass2.yaml @@ -0,0 +1,44 @@ +# Use imSim custom modules +modules: + - imsim + +# Get most of the configuration from the imSim config-template +# for instance catalogs. +template: imsim-config-instcat + +################################################################ +# Make your changes below. +################################################################ + +# Put your own commands that override the defaults below here. For example +# input.instance_catalog.file_name: ./imsim_cat_197356.txt +# input.instance_catalog.sort_mag: False +# input.tree_rings.only_dets: [R22_S11] +# image.nobjects: 5 + +input.instance_catalog.file_name: $os.environ.get('IMSIM_HOME')+'/imSim/examples/example_instance_catalog_ffp.txt' +input.instance_catalog.sort_mag: False + +input.tree_rings.only_dets: [R22_S12] +input.atm_psf: "" + +input.scattered_photons: + file_name: "output/scattered_00398414-0-r-R22_S11-det094.fits" + camera: $camera + det_name: $det_name + +psf: + type: Convolve + items: + - + type: Gaussian + fwhm: 0.3 + +image.nobjects: 10 + +stamp.fft_sb_thresh: 1.e30 + +output.det_num.first: 95 +output.nfiles: 1 + +output.dir: output diff --git a/imsim/photon_pooling.py b/imsim/photon_pooling.py index 5ab74927..12e31777 100644 --- a/imsim/photon_pooling.py +++ b/imsim/photon_pooling.py @@ -165,7 +165,12 @@ def buildImage(self, config, base, image_num, _obj_num, logger): scattered_indices = [i for i in range(len(photons)) if not imview.bounds.includes(photons.x[i], photons.y[i])] if len(scattered_indices) > 0: scattered_photons = galsim.PhotonArray(len(scattered_indices)) - scattered_photons.copyFrom(photons, target_indices=slice(len(scattered_indices)), source_indices=scattered_indices) + scattered_photons.copyFrom(photons, + target_indices=slice(len(scattered_indices)), + source_indices=scattered_indices, + do_xy=True, + do_flux=True, + do_other=False) base['scattered_photons'].append(scattered_photons) del photons # As with the stamps above, let the garbage collector know we don't need the photons anymore. diff --git a/imsim/photon_scattering.py b/imsim/photon_scattering.py index 13dac7d8..2368eb38 100644 --- a/imsim/photon_scattering.py +++ b/imsim/photon_scattering.py @@ -48,22 +48,23 @@ def writeHdu(self, config, base, logger): return BinTableHDU(data=None) -class ScatteredPhotonsInput(object): - """A class to read in the scattered photons that were written to file - during an earlier first pass. +class ScatteredPhotons(object): + """A class to hold the photons which fell outside the sensor being drawn + by the task that createed them during the first pass. They were saved to file + then, and will now be read in this second pass to be accumulated on other sensors. """ - def __init__(self, file_name, wcs, det, xsize=4096, ysize=4096, logger=None): + def __init__(self, file_name, camera, det_name, xsize=4096, ysize=4096, logger=None): """ Initialize the scattered photons input class. Parameters: file_name: str The name of the file to read. - wcs: galsim.WCS - The WCS object to use for the image. - det: lsst.afw.cameraGeom.Detector - The detector for the current sensor. + camera: str + The name of the camera containing the detector. + det_name: str + The name of the detector to use for photon coordinate transformations. xsize: int The x size in pixels of the CCD. (default: 4096) ysize: int @@ -72,28 +73,21 @@ def __init__(self, file_name, wcs, det, xsize=4096, ysize=4096, logger=None): A logger object. (default: None) """ self.file_name = file_name - self.wcs = wcs - self.det = det + self.det = get_camera(camera)[det_name] self.xsize = xsize self.ysize = ysize self.logger = logger self._photons = None + self._photons = galsim.PhotonArray.read(self.file_name) + self._photons.x, self._photons.y = focal_to_pixel(self._photons.x, self._photons.y, self.det) def read_photons(self): """Read the scattered photons from the file. """ - # Read the scattered photons from the file then convert them from - # focal plane coordinates to pixel coordinates. - self.photons = galsim.PhotonArray.read(self.file_name) - self.photons.x, self.photons.y = focal_to_pixel(self.photons.x, self.photons.y, self.det) - - @property - def photons(self): - """Get the scattered photons. - """ - if self._photons is None: - self.read_photons() - return self._photons + # Read the scattered photons from the file then transform them from + # focal plane coordinates to this detector's pixel coordinates. + self._photons = galsim.PhotonArray.read(self.file_name) + self._photons.x, self._photons.y = focal_to_pixel(self._photons.x, self._photons.y, self.det) class ScatteredPhotonsLoader(InputLoader): @@ -101,12 +95,10 @@ class ScatteredPhotonsLoader(InputLoader): Class to load scattered photons from file. """ def getKwargs(self, config, base, logger): - req = {'file_name': str} + req = {'file_name': str, 'camera': str, 'det_name': str} opt = {} kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) - wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) - kwargs['wcs'] = wcs kwargs['xsize'] = base.get('det_xsize', 4096) kwargs['ysize'] = base.get('det_ysize', 4096) kwargs['logger'] = logger @@ -142,5 +134,5 @@ def getKwargs(self, config, base, logger): RegisterExtraOutput('scattered_photons', ScatteredPhotonsBuilder()) -RegisterInputType('scattered_photons', ScatteredPhotonsLoader()) +RegisterInputType('scattered_photons', ScatteredPhotonsLoader(ScatteredPhotons)) # RegisterImageType('LSST_ScatteredPhotonsImage', LSST_ScatteredPhotonsImageBuilder) From 43840a988b8a56efb78c121926120f27f4b7a02f Mon Sep 17 00:00:00 2001 From: William Lucas Date: Thu, 10 Jul 2025 09:38:06 +0100 Subject: [PATCH 05/36] Initial full working pipeline --- examples/imsim-user-instcat-ffp-pass1.yaml | 10 ++- examples/imsim-user-instcat-ffp-pass2.yaml | 16 ++++- imsim/photon_pooling.py | 10 +-- imsim/photon_scattering.py | 73 +++++++++++++--------- 4 files changed, 71 insertions(+), 38 deletions(-) diff --git a/examples/imsim-user-instcat-ffp-pass1.yaml b/examples/imsim-user-instcat-ffp-pass1.yaml index 2ea985f1..fa7aafb3 100644 --- a/examples/imsim-user-instcat-ffp-pass1.yaml +++ b/examples/imsim-user-instcat-ffp-pass1.yaml @@ -19,7 +19,12 @@ template: imsim-config-instcat input.instance_catalog.file_name: $os.environ.get('IMSIM_HOME')+'/imSim/examples/example_instance_catalog_ffp.txt' input.instance_catalog.sort_mag: False -input.tree_rings.only_dets: [R22_S11] +# For now, disable tree rings to make the diffraction spikes +# easier to see in the images. +input.tree_rings: "" +image.sensor.treering_center: "" +image.sensor.treering_func: "" + input.atm_psf: "" psf: type: Convolve @@ -33,6 +38,7 @@ image.nobjects: 10 stamp.fft_sb_thresh: 1.e30 output.det_num.first: 94 -output.nfiles: 1 +output.nfiles: 2 +output.nproc: 1 output.dir: output diff --git a/examples/imsim-user-instcat-ffp-pass2.yaml b/examples/imsim-user-instcat-ffp-pass2.yaml index 937424ae..077031e9 100644 --- a/examples/imsim-user-instcat-ffp-pass2.yaml +++ b/examples/imsim-user-instcat-ffp-pass2.yaml @@ -16,12 +16,19 @@ template: imsim-config-instcat # input.tree_rings.only_dets: [R22_S11] # image.nobjects: 5 +input.checkpoint: "" +# input.instance_catalog: "" input.instance_catalog.file_name: $os.environ.get('IMSIM_HOME')+'/imSim/examples/example_instance_catalog_ffp.txt' input.instance_catalog.sort_mag: False +# input.opsim_data: "" -input.tree_rings.only_dets: [R22_S12] +input.tree_rings: "" input.atm_psf: "" +input.initial_image: + file_name: "eimage_00398414-0-r-R22_S12-det095.fits" + dir: "output" + input.scattered_photons: file_name: "output/scattered_00398414-0-r-R22_S11-det094.fits" camera: $camera @@ -34,11 +41,14 @@ psf: type: Gaussian fwhm: 0.3 -image.nobjects: 10 +stamp: "" -stamp.fft_sb_thresh: 1.e30 +image.type: LSST_ScatteredPhotonsImage +image.sensor.treering_center: "" +image.sensor.treering_func: "" output.det_num.first: 95 output.nfiles: 1 +output.file_name.format: "eimage_final_%08d-%1d-%s-%s-det%03d.fits" output.dir: output diff --git a/imsim/photon_pooling.py b/imsim/photon_pooling.py index 12e31777..34862495 100644 --- a/imsim/photon_pooling.py +++ b/imsim/photon_pooling.py @@ -156,13 +156,9 @@ def buildImage(self, config, base, image_num, _obj_num, logger): for op in photon_ops: op.applyTo(photons, local_wcs, rng) - # Now accumulate the photons onto the sensor. Resume is true for all calls but the first. Recalculate the pixel - # boundaries on the first subbatch of each full batch. - self.accumulate_photons(photons, full_image, sensor, resume=(batch_num > current_photon_batch_num or subbatch_num > 0), recalc=(subbatch_num == 0)) - # Gather non-accumulated photons if we're going to be outputting them. if 'scattered_photons' in base['output']: - scattered_indices = [i for i in range(len(photons)) if not imview.bounds.includes(photons.x[i], photons.y[i])] + scattered_indices = [i for i in range(len(photons)) if not full_image.bounds.includes(photons.x[i], photons.y[i])] if len(scattered_indices) > 0: scattered_photons = galsim.PhotonArray(len(scattered_indices)) scattered_photons.copyFrom(photons, @@ -173,6 +169,10 @@ def buildImage(self, config, base, image_num, _obj_num, logger): do_other=False) base['scattered_photons'].append(scattered_photons) + # Now accumulate the photons onto the sensor. Resume is true for all calls but the first. Recalculate the pixel + # boundaries on the first subbatch of each full batch. + self.accumulate_photons(photons, full_image, sensor, resume=(batch_num > current_photon_batch_num or subbatch_num > 0), recalc=(subbatch_num == 0)) + del photons # As with the stamps above, let the garbage collector know we don't need the photons anymore. # Note: in typical imsim usage, all current_vars will be 0. So this normally doesn't diff --git a/imsim/photon_scattering.py b/imsim/photon_scattering.py index 2368eb38..13d86e70 100644 --- a/imsim/photon_scattering.py +++ b/imsim/photon_scattering.py @@ -1,9 +1,10 @@ import galsim from galsim.config.extra import ExtraOutputBuilder, RegisterExtraOutput from galsim.config import ImageBuilder, RegisterImageType, InputLoader, RegisterInputType +from galsim.sensor import Sensor, SiliconSensor from astropy.io.fits import BinTableHDU -# from .lsst_image import LSST_ImageBuilderBase +from .lsst_image import LSST_ImageBuilderBase from .utils import pixel_to_focal, focal_to_pixel from .camera import get_camera @@ -77,17 +78,17 @@ def __init__(self, file_name, camera, det_name, xsize=4096, ysize=4096, logger=N self.xsize = xsize self.ysize = ysize self.logger = logger - self._photons = None - self._photons = galsim.PhotonArray.read(self.file_name) - self._photons.x, self._photons.y = focal_to_pixel(self._photons.x, self._photons.y, self.det) + self.photons = None + self.photons = galsim.PhotonArray.read(self.file_name) + self.photons.x, self.photons.y = focal_to_pixel(self.photons.x, self.photons.y, self.det) def read_photons(self): """Read the scattered photons from the file. """ # Read the scattered photons from the file then transform them from # focal plane coordinates to this detector's pixel coordinates. - self._photons = galsim.PhotonArray.read(self.file_name) - self._photons.x, self._photons.y = focal_to_pixel(self._photons.x, self._photons.y, self.det) + self.photons = galsim.PhotonArray.read(self.file_name) + self.photons.x, self.photons.y = focal_to_pixel(self.photons.x, self.photons.y, self.det) class ScatteredPhotonsLoader(InputLoader): @@ -108,31 +109,47 @@ def getKwargs(self, config, base, logger): return kwargs, safe -# class LSST_ScatteredPhotonsImageBuilder(LSST_ImageBuilderBase): +class LSST_ScatteredPhotonsImageBuilder(LSST_ImageBuilderBase): -# def setup(self, config, base, logger): -# """Set up the scattered photons image type. -# """ -# # We need to set the pixel scale to be the same as the camera's pixel scale -# # so that we can convert between focal plane coordinates and pixel coordinates. -# self.pixel_scale = base['output']['pixel_scale'] -# self.camera = get_camera(base['output']['camera'])[base['det_name']] -# self.focal_plane = self.camera.get_focal_plane() -# self.focal_plane.set_pixel_scale(self.pixel_scale) + # def setup(self, config, base, image_num, obj_num, ignore, logger): + # """Set up the scattered photons image type. + # """ + # # We need to set the pixel scale to be the same as the camera's pixel scale + # # so that we can convert between focal plane coordinates and pixel coordinates. + # self.pixel_scale = base['output']['pixel_scale'] + # self.camera = get_camera(base['output']['camera'])[base['det_name']] + # self.focal_plane = self.camera.get_focal_plane() + # self.focal_plane.set_pixel_scale(self.pixel_scale) -# def draw(self, config, base, logger): -# """Draw the scattered photons to the image. -# """ -# # Get the scattered photons from the base config. -# if 'scattered_photons' in base and len(base['scattered_photons']) > 1: -# self.data = PhotonArray.concatenate(base['scattered_photons']) -# # We need to convert the focal plane coordinates to pixel coordinates -# # so that we can draw them on the image. -# self.data.x, self.data.y = self.focal_plane.focal_to_pixel(self.data.x, self.data.y) -# else: -# self.data = PhotonArray(N=0) + def buildImage(self, config, base, image_num, _obj_num, logger): + """Draw the scattered photons to the image. + """ + # Make sure we have an input image on which to draw. + image = base['current_image'] + + # Make sure we have a scattered_photons input. + if not isinstance(base['_input_objs']['scattered_photons'][0], ScatteredPhotons): + raise galsim.config.GalSimConfigError( + "When using LSST_ScatteredPhotonsImage, you must provide a scattered_photons input.",) + + # Get the scattered photons from the base config. + scattered_photons = base['_input_objs']['scattered_photons'][0].photons + + if len(scattered_photons) > 0: + # There are photons to be accumulated. + # They should already have been transformed to pixel coordinates when they + # were read in from file, so go ahead and accumulate them. + + # Just for testing! + scattered_photons.flux *= 1.e3 + + sensor = base.get('sensor', Sensor()) + + sensor.accumulate(scattered_photons, image) + + return image, [] RegisterExtraOutput('scattered_photons', ScatteredPhotonsBuilder()) RegisterInputType('scattered_photons', ScatteredPhotonsLoader(ScatteredPhotons)) -# RegisterImageType('LSST_ScatteredPhotonsImage', LSST_ScatteredPhotonsImageBuilder) +RegisterImageType('LSST_ScatteredPhotonsImage', LSST_ScatteredPhotonsImageBuilder()) From 12db575c79c4c3a973500302aca78aebd69441e8 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Tue, 15 Jul 2025 18:13:53 +0100 Subject: [PATCH 06/36] Full focal plane tests --- tests/test_photon_scattering.py | 126 ++++++++++++++++++++++++++++++++ 1 file changed, 126 insertions(+) create mode 100644 tests/test_photon_scattering.py diff --git a/tests/test_photon_scattering.py b/tests/test_photon_scattering.py new file mode 100644 index 00000000..aea40f1f --- /dev/null +++ b/tests/test_photon_scattering.py @@ -0,0 +1,126 @@ +import numpy as np +import galsim + +import os + +from imsim.camera import get_camera +from imsim.photon_scattering import gather_scattered_photons +from imsim.utils import focal_to_pixel + + +def test_gather_scattered_photons(): + """Make sure that photons inside and outside the image bounds are gathered correctly.""" + xymin = 0 + xymax = 100 + nphotons = 10 + bounds = galsim.BoundsI(xymin, xymax, xymin, xymax) + + rng = np.random.default_rng(42) + + # Create test photon array + photons = galsim.PhotonArray(nphotons) + photons.flux = np.ones(nphotons) + + # Test all photons with random positions inside the bounds -- none should be gathered. + photons.x = rng.uniform(xymin, xymax, size=nphotons) + photons.y = rng.uniform(xymin, xymax, size=nphotons) + gathered_photons = gather_scattered_photons(bounds, photons) + assert len(gathered_photons) == 0, "Expected 0 photons to be gathered when all are inside the bounds." + + # Test all photons with random positions outside the bounds -- all should be gathered. + photons.x = rng.uniform(xymax+100, xymax+200, size=nphotons) + photons.y = rng.uniform(xymin+100, xymax+200, size=nphotons) + gathered_photons = gather_scattered_photons(bounds, photons) + assert len(gathered_photons) == nphotons, "Expected all {} photons to be gathered when they are outside the bounds.".format(nphotons) + + # Test a mix of photons inside and outside the bounds. + n_inside = nphotons // 2 + n_outside = nphotons - n_inside + photons.x[:n_inside] = rng.uniform(xymin, xymax, size=n_inside) + photons.x[n_inside:] = rng.uniform(xymax+100, xymax+200, size=n_outside) + photons.y[:n_inside] = rng.uniform(xymin, xymax, size=n_inside) + photons.y[n_inside:] = rng.uniform(xymin+100, xymax+200, size=n_outside) + gathered_photons = gather_scattered_photons(bounds, photons) + assert len(gathered_photons) == n_outside, "Of {} photons total, expected the {} outside the bounds to be gathered.".format(nphotons, n_outside) + + +def process_scattered_photons_output(pa_list): + # Given a list of photon arrays (as would be generated by gather_scattered_photons across an image's photon pools), + # create a config and process it to write the photons to file. + base = { + "det_name": "R22_S11", + "output": { + "camera": "LsstCam", + "scattered_photons": { + "dir": "output", + "file_name": "scattered_photons.fits", + }, + }, + "file_num": 0, + "scattered_photons": pa_list, + } + + # This should now generate a concatenated photon array in output/scattered_photons.fits + # with positions transformed to focal plane coordinates. + galsim.config.SetupExtraOutput(base) + galsim.config.ProcessExtraOutputsForImage(base) + galsim.config.WriteExtraOutputs(base, []) + + # Check that the output file exists. + assert os.path.exists("output/scattered_photons.fits"), "Output file 'output/scattered_photons.fits' was not created." + + # Assert that the photons stored in output match. They should need to be transformed + # focal plane to detector coordinates first. + output_photons = galsim.PhotonArray.read("output/scattered_photons.fits") + if len(output_photons) > 0: + detector = get_camera(base['output']['camera'])[base['det_name']] + output_photons.x, output_photons.y = focal_to_pixel(output_photons.x, output_photons.y, detector) + focal_plane_photons = galsim.PhotonArray.concatenate(pa_list) + assert np.allclose(output_photons.x, focal_plane_photons.x), "X coordinates of output photons do not match original photons in focal plane coordinates." + assert np.allclose(output_photons.y, focal_plane_photons.y), "Y coordinates of output photons do not match original photons in focal plane coordinates." + assert np.allclose(output_photons.flux, focal_plane_photons.flux), "Flux of output photons does not match original photons." + + # Final cleanup. + os.remove("output/scattered_photons.fits") + + +def test_scattered_photons_output(): + # Create a list of photon arrays which are all of > 0 length, corresponding to the case + # for which all photon pools found and stored scattered photons. + nphotons_total = 100 + nbatch = 10 + nphot_per_batch = nphotons_total // nbatch + pa_list = [galsim.PhotonArray(nphot_per_batch, + x=np.arange(i*nphot_per_batch, (i+1)*nphot_per_batch), + y=np.arange(i*nphot_per_batch, (i+1)*nphot_per_batch), + flux=np.ones(nphot_per_batch)) + for i in range(nbatch)] + process_scattered_photons_output(pa_list) + + +def test_scattered_photons_output_with_none(): + # A list of empty photon arrays should produce output with an empty photon array. + # This is the case if no photons were found outside the image bounds for the entire image. + nbatch = 10 + pa_list = [galsim.PhotonArray(N=0) for _ in range(nbatch)] + process_scattered_photons_output(pa_list) + + +def test_scattered_photons_output_with_mix(): + # Create a list of photon arrays of different sizes, including some empty ones. + rng = np.random.default_rng(42) + nscattered = [0, 10, 20, 15, 0, 10, 0, 15, 0, 20, 10, 0] + xymin = 0.0 + xymax = 100.0 + pa_list = [galsim.PhotonArray(nphot, + x=rng.uniform(xymin, xymax, size=nphot), + y=rng.uniform(xymin, xymax, size=nphot), + flux=np.ones(nphot)) + for nphot in nscattered] + process_scattered_photons_output(pa_list) + + +if __name__ == "__main__": + testfns = [v for k, v in vars().items() if k.startswith("test_") and callable(v)] + for testfn in testfns: + testfn() From 06ca377f41a26735ebb7c894db3af7e722f4c84b Mon Sep 17 00:00:00 2001 From: William Lucas Date: Wed, 16 Jul 2025 14:14:14 +0100 Subject: [PATCH 07/36] Move gathering of scattered photons into new function in photon_scattering.py --- imsim/photon_pooling.py | 12 ++-------- imsim/photon_scattering.py | 47 +++++++++++++++++++++++--------------- 2 files changed, 30 insertions(+), 29 deletions(-) diff --git a/imsim/photon_pooling.py b/imsim/photon_pooling.py index 34862495..1fd56a0e 100644 --- a/imsim/photon_pooling.py +++ b/imsim/photon_pooling.py @@ -13,6 +13,7 @@ from .stamp import ProcessingMode, ObjectInfo, build_obj from .lsst_image import LSST_ImageBuilderBase +from .photon_scattering import gather_scattered_photons class LSST_PhotonPoolingImageBuilder(LSST_ImageBuilderBase): @@ -158,16 +159,7 @@ def buildImage(self, config, base, image_num, _obj_num, logger): # Gather non-accumulated photons if we're going to be outputting them. if 'scattered_photons' in base['output']: - scattered_indices = [i for i in range(len(photons)) if not full_image.bounds.includes(photons.x[i], photons.y[i])] - if len(scattered_indices) > 0: - scattered_photons = galsim.PhotonArray(len(scattered_indices)) - scattered_photons.copyFrom(photons, - target_indices=slice(len(scattered_indices)), - source_indices=scattered_indices, - do_xy=True, - do_flux=True, - do_other=False) - base['scattered_photons'].append(scattered_photons) + base['scattered_photons'].append(gather_scattered_photons(full_image.bounds, photons)) # Now accumulate the photons onto the sensor. Resume is true for all calls but the first. Recalculate the pixel # boundaries on the first subbatch of each full batch. diff --git a/imsim/photon_scattering.py b/imsim/photon_scattering.py index 13d86e70..02f6ba66 100644 --- a/imsim/photon_scattering.py +++ b/imsim/photon_scattering.py @@ -8,6 +8,33 @@ from .utils import pixel_to_focal, focal_to_pixel from .camera import get_camera +def gather_scattered_photons(image_bounds, photons): + """ Given image bounds and list of photons, gather the photons + that fall outside the bounds, copy them into a new PhotonArray, + and return them. + + Parameters: + image_bounds: galsim.Bounds + The bounds used to determine which photons fall within the image area. + + Returns: + scattered_photons: galsim.PhotonArray + A new PhotonArray containing the photons falling outside the image. + If no photons are found to fall outside, the PhotonArray will be empty. + """ + scattered_indices = [i for i in range(len(photons)) if not image_bounds.includes(photons.x[i], photons.y[i])] + if len(scattered_indices) > 0: + scattered_photons = galsim.PhotonArray(len(scattered_indices)) + scattered_photons.copyFrom(photons, + target_indices=slice(len(scattered_indices)), + source_indices=scattered_indices, + do_xy=True, + do_flux=True, + do_other=False) + else: + scattered_photons = galsim.PhotonArray(N=0) + return scattered_photons + # An extra output type to enable two-pass photon scattering. # This output should be used in the first pass to write scattered photons to file, # while the second pass will read the scattered photons then go over all the images @@ -23,7 +50,7 @@ def processImage(self, index, obj_nums, config, base, logger): """ # If scattered photons have been stored in the base config, concatenate these # photon arrays to a single one and store in data, indexing by image number. - if 'scattered_photons' in base and len(base['scattered_photons']) > 1: + if 'scattered_photons' in base and len(base['scattered_photons']) > 0: self.data[index] = galsim.PhotonArray.concatenate(base['scattered_photons']) # We need to store the photons using focal plane coordinates so they # can be accumulated in the second pass on an arbitrary sensor. @@ -82,14 +109,6 @@ def __init__(self, file_name, camera, det_name, xsize=4096, ysize=4096, logger=N self.photons = galsim.PhotonArray.read(self.file_name) self.photons.x, self.photons.y = focal_to_pixel(self.photons.x, self.photons.y, self.det) - def read_photons(self): - """Read the scattered photons from the file. - """ - # Read the scattered photons from the file then transform them from - # focal plane coordinates to this detector's pixel coordinates. - self.photons = galsim.PhotonArray.read(self.file_name) - self.photons.x, self.photons.y = focal_to_pixel(self.photons.x, self.photons.y, self.det) - class ScatteredPhotonsLoader(InputLoader): """ @@ -111,16 +130,6 @@ def getKwargs(self, config, base, logger): class LSST_ScatteredPhotonsImageBuilder(LSST_ImageBuilderBase): - # def setup(self, config, base, image_num, obj_num, ignore, logger): - # """Set up the scattered photons image type. - # """ - # # We need to set the pixel scale to be the same as the camera's pixel scale - # # so that we can convert between focal plane coordinates and pixel coordinates. - # self.pixel_scale = base['output']['pixel_scale'] - # self.camera = get_camera(base['output']['camera'])[base['det_name']] - # self.focal_plane = self.camera.get_focal_plane() - # self.focal_plane.set_pixel_scale(self.pixel_scale) - def buildImage(self, config, base, image_num, _obj_num, logger): """Draw the scattered photons to the image. """ From ed76d3ecb9b25dfb98ea4fc8f49eeee8fe6fc6e2 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Wed, 23 Jul 2025 15:40:11 +0100 Subject: [PATCH 08/36] Initial go at code to find a science detectors neighbours as part of the Camera class. --- imsim/camera.py | 38 ++++++++++++++++ tests/test_camera.py | 104 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 142 insertions(+) diff --git a/imsim/camera.py b/imsim/camera.py index d3cb45dc..baa82ff4 100644 --- a/imsim/camera.py +++ b/imsim/camera.py @@ -1,6 +1,8 @@ import os +import re import json from collections import defaultdict +import numpy as np import galsim import lsst.utils from .meta_data import data_dir @@ -214,3 +216,39 @@ def update(self, other): def __getattr__(self, attr): """Provide access to the attributes of the underlying lsst_camera.""" return getattr(self.lsst_camera, attr) + + def get_adjacent_detectors(self, det_name): + """ + Given a science detector's name, return a list of the names of science + detectors within the 3x3 grid around it, including itself. + """ + if det_name not in self: + raise galsim.GalSimValueError("Detector is not in camera", det_name) + match = re.match(r"R(\d{2})_S(\d{2})", det_name) + if not match: + # If the given detector is in the list but couldn't be matched, then + # it's a wavefront or guide detector. + raise galsim.GalSimValueError("Arg must be a science detector, not wavefront or guide", det_name) + r00 = (int(match.group(1)[0]), int(match.group(1)[1])) + s00 = (int(match.group(2)[0]), int(match.group(2)[1])) + + # The S indices aren't too bad. We'll always have a full set of the nine + # index pairs, but they need to be permuted to match the position within + # the raft. + s_rows = np.roll(np.array([0, 1, 2]), -s00[0] + 1) + s_cols = np.roll(np.array([0, 1, 2]), -s00[1] + 1) + + # The R indices are more complex. If det_name is positioned to one side + # of the raft, we'll need to decrement the lowest index or increment the + # highest (and never both) for either or both of the row or column + # indices. This integer arithmatic makes it work. + r_rows = np.array([r00[0]-s_rows[0]//2, r00[0], r00[0]+(2-s_rows[2])//2]) + r_cols = np.array([r00[1]-s_cols[0]//2, r00[1], r00[1]+(2-s_cols[2])//2]) + + # That out of the way, loop through all index combinations and add to the + # the list if they exist in the camera. + adjacent_detectors = [f"R{ri}{rj}_S{si}{sj}" + for ri, si in zip(r_rows, s_rows) + for rj, sj in zip(r_cols, s_cols) + if f"R{ri}{rj}_S{si}{sj}" in self] + return adjacent_detectors diff --git a/tests/test_camera.py b/tests/test_camera.py index 47a1a318..d5426581 100644 --- a/tests/test_camera.py +++ b/tests/test_camera.py @@ -6,6 +6,7 @@ import json from pathlib import Path import imsim +from galsim import GalSimValueError DATA_DIR = str(Path(__file__).parent.parent / 'data') @@ -43,6 +44,109 @@ def test_bias_levels(self): for amp_name, amp in ccd.items(): self.assertEqual(bias_level, amp.bias_level) + def check_adjacent_detectors_match(self, camera, det_name, expected_detectors): + camera = imsim.Camera(camera) + adjacent_detectors = camera.get_adjacent_detectors(det_name) + self.assertCountEqual(expected_detectors, adjacent_detectors, + f"Adjacent detectors for {det_name} do not match expected values.") + + def test_get_adjacent_detectors_central(self): + # Case of detectors adjacent to one central in a raft. + det_name = 'R22_S11' + expected_detectors = [ + 'R22_S00', 'R22_S01', 'R22_S02', + 'R22_S10', 'R22_S11', 'R22_S12', + 'R22_S20', 'R22_S21', 'R22_S22', + ] + self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + + def test_get_adjacent_detectors_offset(self): + # Case of detectors adjacent to one offset to a corner in a raft. + # Top left corner. + det_name = 'R22_S20' + expected_detectors = [ + 'R31_S02', 'R32_S00', 'R32_S01', + 'R21_S22', 'R22_S20', 'R22_S21', + 'R21_S12', 'R22_S10', 'R22_S11', + ] + self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + # Bottom right corner. + det_name = 'R22_S02' + expected_detectors = [ + 'R22_S11', 'R22_S12', 'R23_S10', + 'R22_S01', 'R22_S02', 'R23_S00', + 'R12_S21', 'R12_S22', 'R13_S20', + ] + self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + + def test_get_adjacent_detectors_edge(self): + # Case of detectors adjacent to one on the edge of the camera. + # Bottom edge. + det_name = 'R02_S01' + expected_detectors = [ + 'R02_S10', 'R02_S11', 'R02_S12', + 'R02_S00', 'R02_S01', 'R02_S02', + ] + self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + # Left edge. + det_name = 'R20_S10' + expected_detectors = [ + 'R20_S20', 'R20_S21', + 'R20_S10', 'R20_S11', + 'R20_S00', 'R20_S01', + ] + # Top edge. + det_name = 'R42_S21' + expected_detectors = [ + 'R42_S20', 'R42_S21', 'R42_S22', + 'R42_S10', 'R42_S11', 'R42_S12', + ] + # Right edge. + det_name = 'R24_S12' + expected_detectors = [ + 'R24_S21', 'R24_S22', + 'R24_S11', 'R24_S12', + 'R24_S01', 'R24_S02', + ] + self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + + def test_get_adjacent_detectors_corner(self): + # Case of detectors adjacent to one in a corner of the camera. + # We don't expect wavefront or guider detectors to be returned. + # On the side of a corner cutout. + det_name = 'R01_S10' + expected_detectors = [ + 'R01_S20', 'R01_S21', + 'R01_S10', 'R01_S11', + 'R01_S00', 'R01_S01', + ] + self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + # On the corner of a corner cutout. + det_name = 'R11_S00' + expected_detectors = [ + 'R10_S12', 'R11_S10', 'R11_S11', + 'R10_S02', 'R11_S00', 'R11_S01', + 'R01_S20', 'R01_S21', + ] + self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + + def test_get_adjacent_detectors_exceptions(self): + # Ensure that get_adjacent_detectors raises errors for invalid detecor + # names -- either ones which don't exist, or non science ones. + camera = imsim.Camera('LsstCam') + + # Non-existent detectors. + with self.assertRaises(GalSimValueError): + camera.get_adjacent_detectors('R99_S99') + with self.assertRaises(GalSimValueError): + camera.get_adjacent_detectors('R00_S00') + + # Non-science detectors. + with self.assertRaises(GalSimValueError): + camera.get_adjacent_detectors('R00_SG0') + with self.assertRaises(GalSimValueError): + camera.get_adjacent_detectors('R00_SW0') + if __name__ == '__main__': unittest.main() From cb004103644793a70f79782e1be8c214336beea7 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Wed, 23 Jul 2025 17:38:44 +0100 Subject: [PATCH 09/36] Make references to 'scattered' photons 'off-detector' photons instead -- more general, as we don't want to presume how they got off the detector. Be a bit more general. --- imsim/photon_pooling.py | 13 +- imsim/photon_scattering.py | 132 +++++++++--------- ...scattering.py => test_full_focal_plane.py} | 71 ++++++---- 3 files changed, 114 insertions(+), 102 deletions(-) rename tests/{test_photon_scattering.py => test_full_focal_plane.py} (65%) diff --git a/imsim/photon_pooling.py b/imsim/photon_pooling.py index 1fd56a0e..62d0d467 100644 --- a/imsim/photon_pooling.py +++ b/imsim/photon_pooling.py @@ -13,7 +13,7 @@ from .stamp import ProcessingMode, ObjectInfo, build_obj from .lsst_image import LSST_ImageBuilderBase -from .photon_scattering import gather_scattered_photons +from .full_focal_plane import gather_out_of_bounds_photons class LSST_PhotonPoolingImageBuilder(LSST_ImageBuilderBase): @@ -139,8 +139,8 @@ def buildImage(self, config, base, image_num, _obj_num, logger): local_wcs = base['wcs'].local(full_image.true_center) if sensor is None: sensor = Sensor() - if 'scattered_photons' in base['output']: - base['scattered_photons'] = [] # Initialize the scattered_photons list + if 'off_detector_photons' in base['output']: + base['off_detector_photons'] = [] # Initialize the off_detector_photons list for batch_num, batch in enumerate(phot_batches, start=current_photon_batch_num): if not batch: continue @@ -157,9 +157,10 @@ def buildImage(self, config, base, image_num, _obj_num, logger): for op in photon_ops: op.applyTo(photons, local_wcs, rng) - # Gather non-accumulated photons if we're going to be outputting them. - if 'scattered_photons' in base['output']: - base['scattered_photons'].append(gather_scattered_photons(full_image.bounds, photons)) + # Gather off-detector photons if we're going to be outputting them + # (likely to draw them on the other sensors in a second pass). + if 'off_detector_photons' in base['output']: + base['off_detector_photons'].append(gather_out_of_bounds_photons(full_image.bounds, photons)) # Now accumulate the photons onto the sensor. Resume is true for all calls but the first. Recalculate the pixel # boundaries on the first subbatch of each full batch. diff --git a/imsim/photon_scattering.py b/imsim/photon_scattering.py index 02f6ba66..5c5484c6 100644 --- a/imsim/photon_scattering.py +++ b/imsim/photon_scattering.py @@ -1,82 +1,85 @@ import galsim from galsim.config.extra import ExtraOutputBuilder, RegisterExtraOutput -from galsim.config import ImageBuilder, RegisterImageType, InputLoader, RegisterInputType -from galsim.sensor import Sensor, SiliconSensor +from galsim.config import RegisterImageType, InputLoader, RegisterInputType +from galsim.sensor import Sensor from astropy.io.fits import BinTableHDU from .lsst_image import LSST_ImageBuilderBase from .utils import pixel_to_focal, focal_to_pixel from .camera import get_camera -def gather_scattered_photons(image_bounds, photons): - """ Given image bounds and list of photons, gather the photons - that fall outside the bounds, copy them into a new PhotonArray, - and return them. +def gather_out_of_bounds_photons(image_bounds, photons): + """ Given image bounds and a list of photon arrays, gather the photons + that fall outside the bounds, copy them into a new PhotonArray, and return + them. Parameters: image_bounds: galsim.Bounds - The bounds used to determine which photons fall within the image area. + The bounds used to determine which photons fall outside. + photons: List(galsim.PhotonArray) + A list of PhotonArrays containing all the photons to check. Returns: - scattered_photons: galsim.PhotonArray - A new PhotonArray containing the photons falling outside the image. - If no photons are found to fall outside, the PhotonArray will be empty. + out_of_bounds_photons: galsim.PhotonArray + A new PhotonArray containing the photons falling outside the bounds. + If no such photons are found, the PhotonArray will be empty (N=0). """ - scattered_indices = [i for i in range(len(photons)) if not image_bounds.includes(photons.x[i], photons.y[i])] - if len(scattered_indices) > 0: - scattered_photons = galsim.PhotonArray(len(scattered_indices)) - scattered_photons.copyFrom(photons, - target_indices=slice(len(scattered_indices)), - source_indices=scattered_indices, - do_xy=True, - do_flux=True, - do_other=False) + out_of_bounds_indices = [i for i in range(len(photons)) + if not image_bounds.includes(photons.x[i], photons.y[i])] + if len(out_of_bounds_indices) > 0: + out_of_bounds_photons = galsim.PhotonArray(len(out_of_bounds_indices)) + out_of_bounds_photons.copyFrom(photons, + target_indices=slice(len(out_of_bounds_indices)), + source_indices=out_of_bounds_indices, + do_xy=True, + do_flux=True, + do_other=False) else: - scattered_photons = galsim.PhotonArray(N=0) - return scattered_photons - -# An extra output type to enable two-pass photon scattering. -# This output should be used in the first pass to write scattered photons to file, -# while the second pass will read the scattered photons then go over all the images -# again and add the photons. -class ScatteredPhotonsBuilder(ExtraOutputBuilder): - """Build pickled photon arrays containing photons which were not accumulated - to the sensor during each image's construction. + out_of_bounds_photons = galsim.PhotonArray(N=0) + return out_of_bounds_photons + +# An extra output type to enable two-pass drawing of off-detector photons. This +# output should be used in the first pass to write off-detector photons to file, +# while the second pass will read the off-detector photons then go through all +# the images again and draw the photons on top of the first pass image. +class OffDetectorPhotonsBuilder(ExtraOutputBuilder): + """Build photon arrays containing the off-detector photons found during an + image's construction and write them to file. """ def processImage(self, index, obj_nums, config, base, logger): - """After each image is drawn, concatenate the list of scattered photon arrays - and store in data ready to be pickled to file later on. + """After each image is drawn, concatenate the list of off-detector + photon arrays found from each batch/sub-batch and store in data ready + to be written to file later on. """ - # If scattered photons have been stored in the base config, concatenate these - # photon arrays to a single one and store in data, indexing by image number. - if 'scattered_photons' in base and len(base['scattered_photons']) > 0: - self.data[index] = galsim.PhotonArray.concatenate(base['scattered_photons']) + if 'off_detector_photons' in base and len(base['off_detector_photons']) > 0: + self.data[index] = galsim.PhotonArray.concatenate(base['off_detector_photons']) # We need to store the photons using focal plane coordinates so they # can be accumulated in the second pass on an arbitrary sensor. detector = get_camera(base['output']['camera'])[base['det_name']] self.data[index].x, self.data[index].y = pixel_to_focal(self.data[index].x, self.data[index].y, detector) else: + # If we've been told to write off-detector photons but found none, + # let's at least write an empty PhotonArray. self.data[index] = galsim.PhotonArray(N=0) def finalize(self, config, base, main_data, logger): return self.data def writeFile(self, file_name, config, base, logger): - """Write the photon array to fits file. - Might it be faster/more space efficient to pickle? + """Write the photon array to fits. """ for photon_array in self.final_data: photon_array.write(file_name) def writeHdu(self, config, base, logger): - """We don't want to write the scattered photons to FITS, so return an + """We don't want to write the off-detector photons to FITS, so return an empty BinTable in lieu of something else. """ return BinTableHDU(data=None) -class ScatteredPhotons(object): +class OffDetectorPhotons(object): """A class to hold the photons which fell outside the sensor being drawn by the task that createed them during the first pass. They were saved to file then, and will now be read in this second pass to be accumulated on other sensors. @@ -84,7 +87,7 @@ class ScatteredPhotons(object): def __init__(self, file_name, camera, det_name, xsize=4096, ysize=4096, logger=None): """ - Initialize the scattered photons input class. + Initialize the off-detector photons input class. Parameters: file_name: str @@ -105,14 +108,14 @@ def __init__(self, file_name, camera, det_name, xsize=4096, ysize=4096, logger=N self.xsize = xsize self.ysize = ysize self.logger = logger - self.photons = None + # self.photons = None self.photons = galsim.PhotonArray.read(self.file_name) self.photons.x, self.photons.y = focal_to_pixel(self.photons.x, self.photons.y, self.det) -class ScatteredPhotonsLoader(InputLoader): +class OffDetectorPhotonsLoader(InputLoader): """ - Class to load scattered photons from file. + Class to load off-detector photons from file. """ def getKwargs(self, config, base, logger): req = {'file_name': str, 'camera': str, 'det_name': str} @@ -123,42 +126,41 @@ def getKwargs(self, config, base, logger): kwargs['ysize'] = base.get('det_ysize', 4096) kwargs['logger'] = logger - # For now assume the scattered photons on file can be used for all detectors. + # For now assume the off-detector photons on file can be used for + # all detectors. safe = True return kwargs, safe -class LSST_ScatteredPhotonsImageBuilder(LSST_ImageBuilderBase): +# A class to build one of the images making up a full focal plane, used in a +# second pass which uses as input images from the first pass along with any +# photons which fell off-detector. +class LSST_FocalPlaneImageBuilder(LSST_ImageBuilderBase): def buildImage(self, config, base, image_num, _obj_num, logger): - """Draw the scattered photons to the image. + """Draw the off-detector photons to the image. """ - # Make sure we have an input image on which to draw. + # Make sure we have an input image and off-detector photons to draw. + # Without them, this type of image won't work. image = base['current_image'] - - # Make sure we have a scattered_photons input. - if not isinstance(base['_input_objs']['scattered_photons'][0], ScatteredPhotons): + if not isinstance(base['_input_objs']['off_detector_photons'][0], OffDetectorPhotons): raise galsim.config.GalSimConfigError( - "When using LSST_ScatteredPhotonsImage, you must provide a scattered_photons input.",) - - # Get the scattered photons from the base config. - scattered_photons = base['_input_objs']['scattered_photons'][0].photons - - if len(scattered_photons) > 0: - # There are photons to be accumulated. - # They should already have been transformed to pixel coordinates when they - # were read in from file, so go ahead and accumulate them. + "When using LSST_FocalPlaneImage, you must provide an off_detector_photons input.",) + + off_detector_photons = base['_input_objs']['off_detector_photons'][0].photons - # Just for testing! - scattered_photons.flux *= 1.e3 + if len(off_detector_photons) > 0: + # There are photons to be accumulated. They should already have been + # transformed to pixel coordinates when they were read from file, so + # go ahead and accumulate them. sensor = base.get('sensor', Sensor()) - sensor.accumulate(scattered_photons, image) + sensor.accumulate(off_detector_photons, image) return image, [] -RegisterExtraOutput('scattered_photons', ScatteredPhotonsBuilder()) -RegisterInputType('scattered_photons', ScatteredPhotonsLoader(ScatteredPhotons)) -RegisterImageType('LSST_ScatteredPhotonsImage', LSST_ScatteredPhotonsImageBuilder()) +RegisterExtraOutput('off_detector_photons', OffDetectorPhotonsBuilder()) +RegisterInputType('off_detector_photons', OffDetectorPhotonsLoader(OffDetectorPhotons)) +RegisterImageType('LSST_FocalPlaneImage', LSST_FocalPlaneImageBuilder()) diff --git a/tests/test_photon_scattering.py b/tests/test_full_focal_plane.py similarity index 65% rename from tests/test_photon_scattering.py rename to tests/test_full_focal_plane.py index aea40f1f..c48ea793 100644 --- a/tests/test_photon_scattering.py +++ b/tests/test_full_focal_plane.py @@ -4,11 +4,11 @@ import os from imsim.camera import get_camera -from imsim.photon_scattering import gather_scattered_photons +from imsim.full_focal_plane import gather_out_of_bounds_photons from imsim.utils import focal_to_pixel -def test_gather_scattered_photons(): +def test_gather_out_of_bounds_photons(): """Make sure that photons inside and outside the image bounds are gathered correctly.""" xymin = 0 xymax = 100 @@ -24,13 +24,13 @@ def test_gather_scattered_photons(): # Test all photons with random positions inside the bounds -- none should be gathered. photons.x = rng.uniform(xymin, xymax, size=nphotons) photons.y = rng.uniform(xymin, xymax, size=nphotons) - gathered_photons = gather_scattered_photons(bounds, photons) + gathered_photons = gather_out_of_bounds_photons(bounds, photons) assert len(gathered_photons) == 0, "Expected 0 photons to be gathered when all are inside the bounds." # Test all photons with random positions outside the bounds -- all should be gathered. photons.x = rng.uniform(xymax+100, xymax+200, size=nphotons) photons.y = rng.uniform(xymin+100, xymax+200, size=nphotons) - gathered_photons = gather_scattered_photons(bounds, photons) + gathered_photons = gather_out_of_bounds_photons(bounds, photons) assert len(gathered_photons) == nphotons, "Expected all {} photons to be gathered when they are outside the bounds.".format(nphotons) # Test a mix of photons inside and outside the bounds. @@ -40,38 +40,40 @@ def test_gather_scattered_photons(): photons.x[n_inside:] = rng.uniform(xymax+100, xymax+200, size=n_outside) photons.y[:n_inside] = rng.uniform(xymin, xymax, size=n_inside) photons.y[n_inside:] = rng.uniform(xymin+100, xymax+200, size=n_outside) - gathered_photons = gather_scattered_photons(bounds, photons) + gathered_photons = gather_out_of_bounds_photons(bounds, photons) assert len(gathered_photons) == n_outside, "Of {} photons total, expected the {} outside the bounds to be gathered.".format(nphotons, n_outside) -def process_scattered_photons_output(pa_list): - # Given a list of photon arrays (as would be generated by gather_scattered_photons across an image's photon pools), - # create a config and process it to write the photons to file. +def process_off_detector_photons_output(pa_list): + # Given a list of photon arrays (as would be generated by + # gather_out_of_bounds_photons across an image's photon pools), create a + # config and process it to write the photons to file. base = { "det_name": "R22_S11", "output": { "camera": "LsstCam", - "scattered_photons": { + "off_detector_photons": { "dir": "output", - "file_name": "scattered_photons.fits", + "file_name": "off_det_photons.fits", }, }, "file_num": 0, - "scattered_photons": pa_list, + "off_detector_photons": pa_list, } - # This should now generate a concatenated photon array in output/scattered_photons.fits - # with positions transformed to focal plane coordinates. + # This should now generate a concatenated photon array in + # output/off_det_photons.fits with positions transformed to focal plane + # coordinates. galsim.config.SetupExtraOutput(base) galsim.config.ProcessExtraOutputsForImage(base) galsim.config.WriteExtraOutputs(base, []) # Check that the output file exists. - assert os.path.exists("output/scattered_photons.fits"), "Output file 'output/scattered_photons.fits' was not created." + assert os.path.exists("output/off_det_photons.fits"), "Output file 'output/off_det_photons.fits' was not created." - # Assert that the photons stored in output match. They should need to be transformed - # focal plane to detector coordinates first. - output_photons = galsim.PhotonArray.read("output/scattered_photons.fits") + # Assert that the photons stored in the output match. They should need to be + # transformed from focal plane to detector coordinates first. + output_photons = galsim.PhotonArray.read("output/off_det_photons.fits") if len(output_photons) > 0: detector = get_camera(base['output']['camera'])[base['det_name']] output_photons.x, output_photons.y = focal_to_pixel(output_photons.x, output_photons.y, detector) @@ -81,12 +83,13 @@ def process_scattered_photons_output(pa_list): assert np.allclose(output_photons.flux, focal_plane_photons.flux), "Flux of output photons does not match original photons." # Final cleanup. - os.remove("output/scattered_photons.fits") + os.remove("output/off_det_photons.fits") -def test_scattered_photons_output(): - # Create a list of photon arrays which are all of > 0 length, corresponding to the case - # for which all photon pools found and stored scattered photons. +def test_off_detector_photons_output(): + # Create a list of photon arrays which are all of > 0 length, corresponding + # to the case for which all photon pools found and stored off-detector + # photons. nphotons_total = 100 nbatch = 10 nphot_per_batch = nphotons_total // nbatch @@ -95,29 +98,35 @@ def test_scattered_photons_output(): y=np.arange(i*nphot_per_batch, (i+1)*nphot_per_batch), flux=np.ones(nphot_per_batch)) for i in range(nbatch)] - process_scattered_photons_output(pa_list) + process_off_detector_photons_output(pa_list) -def test_scattered_photons_output_with_none(): - # A list of empty photon arrays should produce output with an empty photon array. - # This is the case if no photons were found outside the image bounds for the entire image. +def test_off_detector_photons_output_with_none(): + # A list of empty photon arrays should produce output with an empty photon + # array. This is the case if no photons were found outside the image bounds + # for the entire image. nbatch = 10 pa_list = [galsim.PhotonArray(N=0) for _ in range(nbatch)] - process_scattered_photons_output(pa_list) + process_off_detector_photons_output(pa_list) -def test_scattered_photons_output_with_mix(): - # Create a list of photon arrays of different sizes, including some empty ones. +def test_off_detector_photons_output_with_mix(): + # Create a list of photon arrays of different sizes, including some empty + # ones. This is the case when some photon pools contained off-detector + # photons and some didn't. rng = np.random.default_rng(42) - nscattered = [0, 10, 20, 15, 0, 10, 0, 15, 0, 20, 10, 0] + n_off_det = [0, 10, 20, 15, 0, 10, 0, 15, 0, 20, 10, 0] xymin = 0.0 xymax = 100.0 pa_list = [galsim.PhotonArray(nphot, x=rng.uniform(xymin, xymax, size=nphot), y=rng.uniform(xymin, xymax, size=nphot), flux=np.ones(nphot)) - for nphot in nscattered] - process_scattered_photons_output(pa_list) + for nphot in n_off_det] + process_off_detector_photons_output(pa_list) + +# TODO: Need some tests for the OffDetectorPhotonsLoader and +# LSST_FocalPlaneImageBuilder. if __name__ == "__main__": From 582c76d0848da43c6bd5bdc3bb2879785f9aa4c8 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Wed, 23 Jul 2025 17:39:31 +0100 Subject: [PATCH 10/36] Rename photon_scattering.py to full_focal_plane.py --- imsim/__init__.py | 2 +- imsim/camera.py | 27 ++++++++++--------- ...oton_scattering.py => full_focal_plane.py} | 0 3 files changed, 15 insertions(+), 14 deletions(-) rename imsim/{photon_scattering.py => full_focal_plane.py} (100%) diff --git a/imsim/__init__.py b/imsim/__init__.py index 33c6540f..69965189 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -40,4 +40,4 @@ from .process_info import * from .table_row import * from .photon_pooling import * -from .photon_scattering import * +from .full_focal_plane import * diff --git a/imsim/camera.py b/imsim/camera.py index baa82ff4..fa34353d 100644 --- a/imsim/camera.py +++ b/imsim/camera.py @@ -228,25 +228,26 @@ def get_adjacent_detectors(self, det_name): if not match: # If the given detector is in the list but couldn't be matched, then # it's a wavefront or guide detector. - raise galsim.GalSimValueError("Arg must be a science detector, not wavefront or guide", det_name) + raise galsim.GalSimValueError( + "Arg must be a science detector, not wavefront or guide", det_name) r00 = (int(match.group(1)[0]), int(match.group(1)[1])) s00 = (int(match.group(2)[0]), int(match.group(2)[1])) - # The S indices aren't too bad. We'll always have a full set of the nine - # index pairs, but they need to be permuted to match the position within - # the raft. - s_rows = np.roll(np.array([0, 1, 2]), -s00[0] + 1) - s_cols = np.roll(np.array([0, 1, 2]), -s00[1] + 1) + # The following method of manipulating the indices in the detector name + # is about 5000 times faster than using lsst.afw.cameraGeom to find the + # detector's center and then the names of those shifted away by one + # detector width/height. - # The R indices are more complex. If det_name is positioned to one side - # of the raft, we'll need to decrement the lowest index or increment the - # highest (and never both) for either or both of the row or column - # indices. This integer arithmatic makes it work. + # We'll always have a full set of the nine S index pairs, but they need + # to be permuted to match the position within the raft. + s_rows = np.roll(np.array([0, 1, 2]), 1 - s00[0]) + s_cols = np.roll(np.array([0, 1, 2]), 1 - s00[1]) + + # Use some integer arithmetic to increment or decrement the raft indices + # if the detector is positioned on an edge of the raft. r_rows = np.array([r00[0]-s_rows[0]//2, r00[0], r00[0]+(2-s_rows[2])//2]) r_cols = np.array([r00[1]-s_cols[0]//2, r00[1], r00[1]+(2-s_cols[2])//2]) - - # That out of the way, loop through all index combinations and add to the - # the list if they exist in the camera. + adjacent_detectors = [f"R{ri}{rj}_S{si}{sj}" for ri, si in zip(r_rows, s_rows) for rj, sj in zip(r_cols, s_cols) diff --git a/imsim/photon_scattering.py b/imsim/full_focal_plane.py similarity index 100% rename from imsim/photon_scattering.py rename to imsim/full_focal_plane.py From 32e57c17b9d50b17cf4b9c23eb673ca3c83b9303 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Fri, 8 Aug 2025 15:06:15 +0100 Subject: [PATCH 11/36] Changes to variable names in get_adjacent_detectors() --- imsim/camera.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/imsim/camera.py b/imsim/camera.py index fa34353d..0537ad21 100644 --- a/imsim/camera.py +++ b/imsim/camera.py @@ -230,8 +230,8 @@ def get_adjacent_detectors(self, det_name): # it's a wavefront or guide detector. raise galsim.GalSimValueError( "Arg must be a science detector, not wavefront or guide", det_name) - r00 = (int(match.group(1)[0]), int(match.group(1)[1])) - s00 = (int(match.group(2)[0]), int(match.group(2)[1])) + r_det = (int(match.group(1)[0]), int(match.group(1)[1])) + s_det = (int(match.group(2)[0]), int(match.group(2)[1])) # The following method of manipulating the indices in the detector name # is about 5000 times faster than using lsst.afw.cameraGeom to find the @@ -240,13 +240,13 @@ def get_adjacent_detectors(self, det_name): # We'll always have a full set of the nine S index pairs, but they need # to be permuted to match the position within the raft. - s_rows = np.roll(np.array([0, 1, 2]), 1 - s00[0]) - s_cols = np.roll(np.array([0, 1, 2]), 1 - s00[1]) + s_rows = np.roll(np.array([0, 1, 2]), 1 - s_det[0]) + s_cols = np.roll(np.array([0, 1, 2]), 1 - s_det[1]) - # Use some integer arithmetic to increment or decrement the raft indices - # if the detector is positioned on an edge of the raft. - r_rows = np.array([r00[0]-s_rows[0]//2, r00[0], r00[0]+(2-s_rows[2])//2]) - r_cols = np.array([r00[1]-s_cols[0]//2, r00[1], r00[1]+(2-s_cols[2])//2]) + # A little integer arithmetic is used to increment or decrement the raft + # indices if the detector is positioned on an edge of the raft. + r_rows = np.array([r_det[0]-s_rows[0]//2, r_det[0], r_det[0]+(2-s_rows[2])//2]) + r_cols = np.array([r_det[1]-s_cols[0]//2, r_det[1], r_det[1]+(2-s_cols[2])//2]) adjacent_detectors = [f"R{ri}{rj}_S{si}{sj}" for ri, si in zip(r_rows, s_rows) From b9817915d83b568556c5c16675aa9d0203e2c353 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Fri, 22 Aug 2025 10:15:36 +0100 Subject: [PATCH 12/36] Move off detector photons from photon pooling config to FFP pass 1 example --- config/imsim-config-photon-pooling.yaml | 11 ----------- config/imsim-config-skycat.yaml | 2 +- examples/imsim-user-instcat-ffp-pass1.yaml | 12 ++++++++++++ examples/imsim-user-instcat-ffp-pass2.yaml | 8 ++++---- 4 files changed, 17 insertions(+), 16 deletions(-) diff --git a/config/imsim-config-photon-pooling.yaml b/config/imsim-config-photon-pooling.yaml index 2be4b981..58a56a1c 100644 --- a/config/imsim-config-photon-pooling.yaml +++ b/config/imsim-config-photon-pooling.yaml @@ -66,14 +66,3 @@ output.photon_pooling_truth: phot_flux: "@phot_flux" # The realized flux for photon shooting. fft_flux: "@fft_flux" # If FFT rendering, then the flux used, including vignetting. incident_flux: "@incident_flux" # The sum of the fluxes of all photons from the object. -output.scattered_photons: - dir: output - file_name: - type: FormattedStr - format : scattered_%08d-%1d-%s-%s-det%03d.fits - items: - - { type: OpsimData, field: observationId } - - { type: OpsimData, field: snap } - - $band - - $det_name - - "@output.det_num" \ No newline at end of file diff --git a/config/imsim-config-skycat.yaml b/config/imsim-config-skycat.yaml index ebcba212..57ba1ab3 100644 --- a/config/imsim-config-skycat.yaml +++ b/config/imsim-config-skycat.yaml @@ -13,7 +13,7 @@ # Get most of the configuration from the base imSim config template. modules: - imsim -template: imsim-config-photon-pooling +template: imsim-config #################################################################### # The following entires are added to the base configuration above. diff --git a/examples/imsim-user-instcat-ffp-pass1.yaml b/examples/imsim-user-instcat-ffp-pass1.yaml index fa7aafb3..e8860e30 100644 --- a/examples/imsim-user-instcat-ffp-pass1.yaml +++ b/examples/imsim-user-instcat-ffp-pass1.yaml @@ -42,3 +42,15 @@ output.nfiles: 2 output.nproc: 1 output.dir: output + +output.off_detector_photons: + dir: output + file_name: + type: FormattedStr + format : off_det_%08d-%1d-%s-%s-det%03d.fits + items: + - { type: OpsimData, field: observationId } + - { type: OpsimData, field: snap } + - $band + - $det_name + - "@output.det_num" diff --git a/examples/imsim-user-instcat-ffp-pass2.yaml b/examples/imsim-user-instcat-ffp-pass2.yaml index 077031e9..7a0bb4a7 100644 --- a/examples/imsim-user-instcat-ffp-pass2.yaml +++ b/examples/imsim-user-instcat-ffp-pass2.yaml @@ -29,8 +29,8 @@ input.initial_image: file_name: "eimage_00398414-0-r-R22_S12-det095.fits" dir: "output" -input.scattered_photons: - file_name: "output/scattered_00398414-0-r-R22_S11-det094.fits" +input.off_detector_photons: + file_name: "output/off_det_00398414-0-r-R22_S11-det094.fits" camera: $camera det_name: $det_name @@ -43,12 +43,12 @@ psf: stamp: "" -image.type: LSST_ScatteredPhotonsImage +image.type: LSST_FocalPlaneImage image.sensor.treering_center: "" image.sensor.treering_func: "" output.det_num.first: 95 output.nfiles: 1 -output.file_name.format: "eimage_final_%08d-%1d-%s-%s-det%03d.fits" +output.file_name.format: "eimage_full_focal_%08d-%1d-%s-%s-det%03d.fits" output.dir: output From 9062a7564bb9f06f037cfe02c32278fa1249a487 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Fri, 22 Aug 2025 16:33:04 +0100 Subject: [PATCH 13/36] Doing some horrible initial experimentation in instcat.py. --- imsim/instcat.py | 91 +++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 86 insertions(+), 5 deletions(-) diff --git a/imsim/instcat.py b/imsim/instcat.py index ff6b006e..d66dc54a 100644 --- a/imsim/instcat.py +++ b/imsim/instcat.py @@ -13,7 +13,11 @@ from galsim import CelestialCoord import galsim import pickle -from .utils import RUBIN_AREA +from .camera import Camera +from .utils import RUBIN_AREA, focal_to_pixel +import lsst.afw.cameraGeom as cameraGeom + +from copy import deepcopy def clarify_radec_limits( @@ -168,8 +172,9 @@ class InstCatalog(object): # cached SEDs. fnu = (0 * u.ABmag).to(u.erg/u.s/u.cm**2/u.Hz) _flux_density = fnu.to_value(u.ph/u.nm/u.s/u.cm**2, u.spectral_density(500*u.nm)) - def __init__(self, file_name, wcs, xsize=4096, ysize=4096, sed_dir=None, - edge_pix=100, sort_mag=True, flip_g2=True, + def __init__(self, file_name, wcs, fiducial_wcs=None, xsize=4096, ysize=4096, sed_dir=None, + edge_pix=200, camera=None, det_name=None, + sort_mag=True, flip_g2=True, pupil_area=RUBIN_AREA, min_source=None, skip_invalid=True, logger=None): logger = galsim.config.LoggerWrapper(logger) @@ -178,6 +183,22 @@ def __init__(self, file_name, wcs, xsize=4096, ysize=4096, sed_dir=None, self.xsize = xsize self.ysize = ysize self.edge_pix = edge_pix + if camera is not None and det_name is not None and fiducial_wcs is not None: + self.multi_det = True + self.camera = Camera(camera_class=camera) + self.det_name = det_name + self.fiducial_wcs = fiducial_wcs + elif camera is None and det_name is None and fiducial_wcs is None: + # Fine, just set everything to None and move on. + self.multi_det = False + self.camera = None + self.det_name = None + self.fiducial_wcs = None + else: + # We're in a bad state where we have some but not all of the + # information needd for a multi-detector run. + # I don't think this should be possible, but just in case... + raise galsim.GalSimNotImplementedError("For multi-detector runs, the camera, detector name, and fiducial_wcs must all be provided to InstCatalog.") self.sort_mag = sort_mag self.flip_g2 = flip_g2 self.min_source = min_source @@ -227,6 +248,19 @@ def _load(self, logger=None): nuse = 0 ntot = 0 + if self.multi_det: + # Get adjacent detectors; then their centers in pixel coordinates, then + # store after converting to RA and dec. + adjacent_detectors = self.camera.get_adjacent_detectors(self.det_name) + detector_centers = {} + transform_mapping = self.camera.lsst_camera[self.det_name].getTransform(cameraGeom.FOCAL_PLANE, cameraGeom.PIXELS).getMapping() + + # Get the focal plane centers for each detector in the list. + for detector in adjacent_detectors: + focal_plane_center = self.camera.lsst_camera[detector].getCenter(cameraGeom.FOCAL_PLANE) + pixel_center = galsim.PositionD(transform_mapping.applyForward(focal_plane_center)) + detector_centers[detector] = self.wcs.toWorld(pixel_center) + with fopen(self.file_name, mode='rt') as _input: for line in _input: if ' inf ' in line: continue @@ -244,6 +278,21 @@ def _load(self, logger=None): and min_dec <= dec <= max_dec ): continue + + # # Next check is whether the object is closer to this detector's center + # # than any of the adjacent detectors. If not, then skip it. + # if self.multi_det: + # obj_dist = {} + # for detector in adjacent_detectors: + # det_center = detector_centers[detector] + # obj_dist[detector] = math.sqrt( + # (ra - det_center.ra)**2 + + # (dec - det_center.dec)**2 + # ) + # closest_detector = min(obj_dist, key=obj_dist.get) + # if closest_detector != self.det_name: + # continue + world_pos = galsim.CelestialCoord(ra, dec) #logger.debug('world_pos = %s',world_pos) try: @@ -642,10 +691,42 @@ def getKwargs(self, config, base, logger): 'skip_invalid' : bool, } kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) - wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) - kwargs['wcs'] = wcs + kwargs['xsize'] = base.get('det_xsize', 4096) kwargs['ysize'] = base.get('det_ysize', 4096) + kwargs['camera'] = base['output'].get('camera', None) + kwargs['det_name'] = base.get('det_name', None) + + if kwargs['camera'] is not None and kwargs['det_name'] is not None: + # A horrible hack to get the Batoid WCS for R22_S11. + # Need a nicer way to do this, and if possible only + # calculated once pre multiprocessing then inherited + # by all the processes. + # orig_det_name = base.get('det_name', None) + # base['det_name'] = 'R22_S11' + # base['image']['det_name'] = 'R22_S11' + if base['image']['wcs'].wcs_type == 'Batoid': + base['image']['wcs']['det_name'] = 'R22_S11' + central_wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) + # Now we have the TAN-SIP WCS for the central detector, + # as determined by batoid, pull out the CD matrix and use + # it to make a new TAN WCS (ie. without the SIP + # distortion terms). + central_affine = galsim.AffineTransform(central_wcs.cd[0,0], + central_wcs.cd[0,1], + central_wcs.cd[1,0], + central_wcs.cd[1,1], + origin=central_wcs.origin) + fiducial_wcs = galsim.TanWCS(central_affine, central_wcs.center) + kwargs['fiducial_wcs'] = fiducial_wcs + del base['image']['wcs']['current'] + # base['det_name'] = orig_det_name + # base['image']['det_name'] = orig_det_name + base['image']['wcs']['det_name'] = '$det_name' + + wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) + + kwargs['wcs'] = wcs kwargs['logger'] = logger return kwargs, False From c205bd9f1859b40685f6b85b6cbadca87b472126 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Tue, 26 Aug 2025 14:18:53 +0100 Subject: [PATCH 14/36] Make the fiducial WCS a component of the instcat in the YAML config. --- config/imsim-config-instcat.yaml | 27 ++++++++++++++++ imsim/instcat.py | 55 ++++++++++++-------------------- 2 files changed, 48 insertions(+), 34 deletions(-) diff --git a/config/imsim-config-instcat.yaml b/config/imsim-config-instcat.yaml index 9f4038e2..2f51f5c5 100644 --- a/config/imsim-config-instcat.yaml +++ b/config/imsim-config-instcat.yaml @@ -25,6 +25,33 @@ input.instance_catalog: file_name: default_catalog_file.txt sed_dir: $os.environ.get('SIMS_SED_LIBRARY_DIR') pupil_area: $pupil_area + # A single WCS calculated used when determining which objects are to be drawn + # by which detector. We could in principle use any WCS. Here we start from the + # same Batoid WCS as in the image. In the code, any SIP components of the + # fiducial WCS, once it's been calculated, are removed. + fiducial_wcs: + type: Batoid + + # These are required: + camera: $camera + boresight: $boresight + + obstime: + type: Eval + str: "astropy.time.Time(mjd_val, format='mjd', scale='tai')" + fmjd_val: { type: OpsimData, field: mjd } + + # The fiducial WCS will use this detector. We default here to the + # central detector, but it can be overrdden if needed. + det_name: "R22_S11" + wavelength: $(@image.bandpass).effective_wavelength + + # The rest can be omitted, since these are the default values, but shown here + # for reference. + temperature: 280 # Kelvin + pressure: 72.7 # kPa + H2O_pressure: 1.0 # kPa + order: 3 # Order of the SIP polynomial input.opsim_data: # Read the visit meta data. By default, we use the same file as the above diff --git a/imsim/instcat.py b/imsim/instcat.py index d66dc54a..049d2c5d 100644 --- a/imsim/instcat.py +++ b/imsim/instcat.py @@ -198,7 +198,7 @@ def __init__(self, file_name, wcs, fiducial_wcs=None, xsize=4096, ysize=4096, se # We're in a bad state where we have some but not all of the # information needd for a multi-detector run. # I don't think this should be possible, but just in case... - raise galsim.GalSimNotImplementedError("For multi-detector runs, the camera, detector name, and fiducial_wcs must all be provided to InstCatalog.") + raise galsim.GalSimNotImplementedError("For multi-detector runs, the camera, det_name, and fiducial_wcs must all be available to InstCatalog.") self.sort_mag = sort_mag self.flip_g2 = flip_g2 self.min_source = min_source @@ -690,43 +690,30 @@ def getKwargs(self, config, base, logger): 'min_source' : int, 'skip_invalid' : bool, } - kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) - + ignore = ['fiducial_wcs'] # Handled separately. + kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore) + wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) + kwargs['wcs'] = wcs + if config.get('fiducial_wcs', None) is not None: + fiducial_wcs = galsim.config.BuildWCS(config, 'fiducial_wcs', base, logger=logger) + if fiducial_wcs.wcs_type == "TAN-SIP": + # If it's a TAN-SIP WCS (such as from Batoid) then we use it to + # create a TAN WCS without the SIP distortions, since we'll + # potentially use this WCS a long way from the det where the + # polynomial was fitted. + fiducial_affine = galsim.AffineTransform(fiducial_wcs.cd[0,0], + fiducial_wcs.cd[0,1], + fiducial_wcs.cd[1,0], + fiducial_wcs.cd[1,1], + origin=fiducial_wcs.origin) + fiducial_wcs = galsim.TanWCS(fiducial_affine, fiducial_wcs.center) + kwargs['fiducial_wcs'] = fiducial_wcs + else: + kwargs['fiducial_wcs'] = None kwargs['xsize'] = base.get('det_xsize', 4096) kwargs['ysize'] = base.get('det_ysize', 4096) kwargs['camera'] = base['output'].get('camera', None) kwargs['det_name'] = base.get('det_name', None) - - if kwargs['camera'] is not None and kwargs['det_name'] is not None: - # A horrible hack to get the Batoid WCS for R22_S11. - # Need a nicer way to do this, and if possible only - # calculated once pre multiprocessing then inherited - # by all the processes. - # orig_det_name = base.get('det_name', None) - # base['det_name'] = 'R22_S11' - # base['image']['det_name'] = 'R22_S11' - if base['image']['wcs'].wcs_type == 'Batoid': - base['image']['wcs']['det_name'] = 'R22_S11' - central_wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) - # Now we have the TAN-SIP WCS for the central detector, - # as determined by batoid, pull out the CD matrix and use - # it to make a new TAN WCS (ie. without the SIP - # distortion terms). - central_affine = galsim.AffineTransform(central_wcs.cd[0,0], - central_wcs.cd[0,1], - central_wcs.cd[1,0], - central_wcs.cd[1,1], - origin=central_wcs.origin) - fiducial_wcs = galsim.TanWCS(central_affine, central_wcs.center) - kwargs['fiducial_wcs'] = fiducial_wcs - del base['image']['wcs']['current'] - # base['det_name'] = orig_det_name - # base['image']['det_name'] = orig_det_name - base['image']['wcs']['det_name'] = '$det_name' - - wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) - - kwargs['wcs'] = wcs kwargs['logger'] = logger return kwargs, False From 536f5aa0b98bbc3b1908352e561a2f0fb788aeae Mon Sep 17 00:00:00 2001 From: William Lucas Date: Fri, 29 Aug 2025 15:30:50 +0100 Subject: [PATCH 15/36] Big improvements to the instcat multi detector object handling. And it now works. --- imsim/instcat.py | 56 ++++++++++++++++-------------------------------- 1 file changed, 19 insertions(+), 37 deletions(-) diff --git a/imsim/instcat.py b/imsim/instcat.py index 049d2c5d..956aad9f 100644 --- a/imsim/instcat.py +++ b/imsim/instcat.py @@ -172,7 +172,8 @@ class InstCatalog(object): # cached SEDs. fnu = (0 * u.ABmag).to(u.erg/u.s/u.cm**2/u.Hz) _flux_density = fnu.to_value(u.ph/u.nm/u.s/u.cm**2, u.spectral_density(500*u.nm)) - def __init__(self, file_name, wcs, fiducial_wcs=None, xsize=4096, ysize=4096, sed_dir=None, + def __init__(self, file_name, wcs, fiducial_wcs=None, fiducial_det_name=None, + xsize=4096, ysize=4096, sed_dir=None, edge_pix=200, camera=None, det_name=None, sort_mag=True, flip_g2=True, pupil_area=RUBIN_AREA, min_source=None, skip_invalid=True, @@ -188,17 +189,10 @@ def __init__(self, file_name, wcs, fiducial_wcs=None, xsize=4096, ysize=4096, se self.camera = Camera(camera_class=camera) self.det_name = det_name self.fiducial_wcs = fiducial_wcs - elif camera is None and det_name is None and fiducial_wcs is None: - # Fine, just set everything to None and move on. - self.multi_det = False - self.camera = None - self.det_name = None - self.fiducial_wcs = None + self.fiducial_det_name = fiducial_det_name else: - # We're in a bad state where we have some but not all of the - # information needd for a multi-detector run. - # I don't think this should be possible, but just in case... - raise galsim.GalSimNotImplementedError("For multi-detector runs, the camera, det_name, and fiducial_wcs must all be available to InstCatalog.") + # Set multi_det to False and move on. + self.multi_det = False self.sort_mag = sort_mag self.flip_g2 = flip_g2 self.min_source = min_source @@ -253,13 +247,13 @@ def _load(self, logger=None): # store after converting to RA and dec. adjacent_detectors = self.camera.get_adjacent_detectors(self.det_name) detector_centers = {} - transform_mapping = self.camera.lsst_camera[self.det_name].getTransform(cameraGeom.FOCAL_PLANE, cameraGeom.PIXELS).getMapping() + transform_mapping = self.camera.lsst_camera[self.fiducial_det_name].getTransform(cameraGeom.FOCAL_PLANE, cameraGeom.PIXELS).getMapping() # Get the focal plane centers for each detector in the list. for detector in adjacent_detectors: focal_plane_center = self.camera.lsst_camera[detector].getCenter(cameraGeom.FOCAL_PLANE) pixel_center = galsim.PositionD(transform_mapping.applyForward(focal_plane_center)) - detector_centers[detector] = self.wcs.toWorld(pixel_center) + detector_centers[detector] = self.fiducial_wcs.toWorld(pixel_center) with fopen(self.file_name, mode='rt') as _input: for line in _input: @@ -279,22 +273,15 @@ def _load(self, logger=None): ): continue - # # Next check is whether the object is closer to this detector's center - # # than any of the adjacent detectors. If not, then skip it. - # if self.multi_det: - # obj_dist = {} - # for detector in adjacent_detectors: - # det_center = detector_centers[detector] - # obj_dist[detector] = math.sqrt( - # (ra - det_center.ra)**2 + - # (dec - det_center.dec)**2 - # ) - # closest_detector = min(obj_dist, key=obj_dist.get) - # if closest_detector != self.det_name: - # continue - world_pos = galsim.CelestialCoord(ra, dec) #logger.debug('world_pos = %s',world_pos) + if self.multi_det: + # Check whether the object is closer to this detector's center + # than any of the adjacent detectors. If not, then skip it. + obj_dist = {det: detector_centers[det].distanceTo(world_pos) for det in adjacent_detectors} + closest_detector = min(obj_dist, key=obj_dist.get) + if closest_detector != self.det_name: + continue try: image_pos = self.wcs.toImage(world_pos) except RuntimeError as e: @@ -697,17 +684,12 @@ def getKwargs(self, config, base, logger): if config.get('fiducial_wcs', None) is not None: fiducial_wcs = galsim.config.BuildWCS(config, 'fiducial_wcs', base, logger=logger) if fiducial_wcs.wcs_type == "TAN-SIP": - # If it's a TAN-SIP WCS (such as from Batoid) then we use it to - # create a TAN WCS without the SIP distortions, since we'll - # potentially use this WCS a long way from the det where the - # polynomial was fitted. - fiducial_affine = galsim.AffineTransform(fiducial_wcs.cd[0,0], - fiducial_wcs.cd[0,1], - fiducial_wcs.cd[1,0], - fiducial_wcs.cd[1,1], - origin=fiducial_wcs.origin) - fiducial_wcs = galsim.TanWCS(fiducial_affine, fiducial_wcs.center) + # If it's a TAN-SIP WCS (such as from Batoid) then disable the + # SIP distortions, since we'll potentially use this WCS a long + # way from the detector where the polynomial was fitted. + fiducial_wcs.ab = None kwargs['fiducial_wcs'] = fiducial_wcs + kwargs['fiducial_det_name'] = config['fiducial_wcs']['det_name'] else: kwargs['fiducial_wcs'] = None kwargs['xsize'] = base.get('det_xsize', 4096) From 95499b15624be7ca20c22fae9f6fbfe8ca19b619 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Wed, 10 Sep 2025 15:52:37 +0100 Subject: [PATCH 16/36] First go at getting skyCatalog interface to accept/reject objects by proximity to detector. --- imsim/skycat.py | 57 +++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 53 insertions(+), 4 deletions(-) diff --git a/imsim/skycat.py b/imsim/skycat.py index 054af0b5..ba14f2bf 100644 --- a/imsim/skycat.py +++ b/imsim/skycat.py @@ -8,13 +8,17 @@ RegisterObjectType from skycatalogs import skyCatalogs from skycatalogs.utils import PolygonalRegion +from .camera import Camera from .utils import RUBIN_AREA +import lsst.afw.cameraGeom as cameraGeom class SkyCatalogInterface: """Interface to skyCatalogs package.""" - def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, - obj_types=None, skycatalog_root=None, edge_pix=100, + def __init__(self, file_name, wcs, band, mjd, fiducial_wcs=None, + fiducial_det_name=None, xsize=4096, ysize=4096, + obj_types=None, skycatalog_root=None, edge_pix=200, + camera=None, det_name=None, pupil_area=RUBIN_AREA, max_flux=None, logger=None, apply_dc2_dilation=False, approx_nobjects=None): """ @@ -28,6 +32,11 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, LSST band, one of ugrizy mjd : float MJD of the midpoint of the exposure. + fiducial_wcs : galsim.WCS [None] + The fiducial WCS, used across all detectors to determine which + should handle which objects. + fiducial_det_name : str [None] + The name of the detector used to create the fiducial WCS. xsize : int Size in pixels of CCD in x-direction. ysize : int @@ -39,9 +48,13 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, Root directory containing skyCatalogs files. If None, then set skycatalog_root to the same directory containing the yaml config file. - edge_pix : float [100] + edge_pix : float [200] Size in pixels of the buffer region around nominal image to consider objects. + camera : str [None] + Name of the imSim camera to use for a multi-detector run. + det_name : str [None] + Name of the current detector. pupil_area : float [RUBIN_AREA] The area of the telescope pupil in cm**2. The default value uses R_outer=418 cm and R_inner=255 cm. @@ -73,6 +86,15 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, else: self.skycatalog_root = skycatalog_root self.edge_pix = edge_pix + if camera is not None and det_name is not None and fiducial_wcs is not None: + self.multi_det = True + self.camera = Camera(camera_class=camera) + self.det_name = det_name + self.fiducial_wcs = fiducial_wcs + self.fiducial_det_name = fiducial_det_name + else: + # Set multi_det to False and move on. + self.multi_det = False self.pupil_area = pupil_area self.max_flux = max_flux self.logger = galsim.config.LoggerWrapper(logger) @@ -87,6 +109,19 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, @property def objects(self): if self._objects is None: + if self.multi_det: + # Get adjacent detectors; then their centers in pixel coordinates, then + # store after converting to RA and dec. + adjacent_detectors = self.camera.get_adjacent_detectors(self.det_name) + detector_centers = {} + transform_mapping = self.camera.lsst_camera[self.fiducial_det_name].getTransform(cameraGeom.FOCAL_PLANE, cameraGeom.PIXELS).getMapping() + + # Get the focal plane centers for each detector in the list. + for detector in adjacent_detectors: + focal_plane_center = self.camera.lsst_camera[detector].getCenter(cameraGeom.FOCAL_PLANE) + pixel_center = galsim.PositionD(transform_mapping.applyForward(focal_plane_center)) + detector_centers[detector] = self.fiducial_wcs.toWorld(pixel_center) + # Select objects from polygonal region bounded by CCD edges corners = ((-self.edge_pix, -self.edge_pix), (self.xsize + self.edge_pix, -self.edge_pix), @@ -218,12 +253,26 @@ def getKwargs(self, config, base, logger): 'mjd': float, 'pupil_area': float, } + ignore = ['fiducial_wcs'] # Handled separately. kwargs, safe = galsim.config.GetAllParams(config, base, req=req, - opt=opt) + opt=opt, ignore=ignore) wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) kwargs['wcs'] = wcs + if config.get('fiducial_wcs', None) is not None: + fiducial_wcs = galsim.config.BuildWCS(config, 'fiducial_wcs', base, logger=logger) + if fiducial_wcs.wcs_type == "TAN-SIP": + # If it's a TAN-SIP WCS (such as from Batoid) then disable the + # SIP distortions, since we'll potentially use this WCS a long + # way from the detector where the polynomial was fitted. + fiducial_wcs.ab = None + kwargs['fiducial_wcs'] = fiducial_wcs + kwargs['fiducial_det_name'] = config['fiducial_wcs']['det_name'] + else: + kwargs['fiducial_wcs'] = None kwargs['xsize'] = base.get('det_xsize', 4096) kwargs['ysize'] = base.get('det_ysize', 4096) + kwargs['camera'] = base['output'].get('camera', None) + kwargs['det_name'] = base.get('det_name', None) kwargs['logger'] = logger # Sky catalog object lists are created per CCD, so they are From b1d9bf8a5041f87005271e3a0d7d9d58e3d289a1 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Thu, 11 Sep 2025 11:57:37 +0100 Subject: [PATCH 17/36] Add example_instance_catalog_ffp.txt --- examples/example_instance_catalog_ffp.txt | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 examples/example_instance_catalog_ffp.txt diff --git a/examples/example_instance_catalog_ffp.txt b/examples/example_instance_catalog_ffp.txt new file mode 100644 index 00000000..a2dadb84 --- /dev/null +++ b/examples/example_instance_catalog_ffp.txt @@ -0,0 +1,23 @@ +rightascension 60.49045502638662697 +declination -38.16437495898705379 +mjd 60143.42156961110595148 +altitude 53.16185928082865786 +azimuth 114.38574818197552929 +filter 2 +rotskypos 131.23506207988341998 +dist2moon 88.2886136 +moonalt -28.7317781 +moondec 24.6999775 +moonphase 0.5804871 +moonra 127.0562204 +nsnap 1 +obshistid 398414 +rottelpos 40.0386345 +seed 398414 +seeing 0.7501810 +sunalt -19.2243392 +vistime 30.0000000 +seqnum 0 +object 32166149840900 60.5808675 -38.07410111 13.8487927 starSED/phoSimMLT/lte035-4.5-1.0a+0.4.BT-Settl.spec.gz 0 0 0 0 0 0 point none CCM 0.0112685 3.1 +object 32166139536388 60.55213375 -38.04768278 13.2693929 starSED/phoSimMLT/lte035-4.5-1.0a+0.4.BT-Settl.spec.gz 0 0 0 0 0 0 point none CCM 0.01014165 3.1 +object 32166139536388 60.55989542 -38.038889 13.2693929 starSED/phoSimMLT/lte035-4.5-1.0a+0.4.BT-Settl.spec.gz 0 0 0 0 0 0 point none CCM 0.01014165 3.1 \ No newline at end of file From 63dbf79e121467f54ad28e6d6375581b19c07ba5 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Thu, 11 Sep 2025 11:59:41 +0100 Subject: [PATCH 18/36] Add the fiducial WCS to imsim-config-skycat.yaml. --- config/imsim-config-skycat.yaml | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/config/imsim-config-skycat.yaml b/config/imsim-config-skycat.yaml index 57ba1ab3..f13d8a52 100644 --- a/config/imsim-config-skycat.yaml +++ b/config/imsim-config-skycat.yaml @@ -26,6 +26,33 @@ input.sky_catalog: band: $band mjd: { type: OpsimData, field: mjd } pupil_area: $pupil_area + # A single WCS calculated used when determining which objects are to be drawn + # by which detector. We could in principle use any WCS. Here we start from the + # same Batoid WCS as in the image. In the code, any SIP components of the + # fiducial WCS, once it's been calculated, are removed. + fiducial_wcs: + type: Batoid + + # These are required: + camera: $camera + boresight: $boresight + + obstime: + type: Eval + str: "astropy.time.Time(mjd_val, format='mjd', scale='tai')" + fmjd_val: { type: OpsimData, field: mjd } + + # The fiducial WCS will use this detector. We default here to the + # central detector, but it can be overrdden if needed. + det_name: "R22_S11" + wavelength: $(@image.bandpass).effective_wavelength + + # The rest can be omitted, since these are the default values, but shown here + # for reference. + temperature: 280 # Kelvin + pressure: 72.7 # kPa + H2O_pressure: 1.0 # kPa + order: 3 # Order of the SIP polynomial input.opsim_data: file_name: default_opsim.db From 2678932511750a44d0e2dac989e827a5075d708e Mon Sep 17 00:00:00 2001 From: William Lucas Date: Thu, 11 Sep 2025 12:35:27 +0100 Subject: [PATCH 19/36] First attempt at logic to filter out objects closer to other dets in skycat.py --- imsim/skycat.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/imsim/skycat.py b/imsim/skycat.py index ba14f2bf..24c2370e 100644 --- a/imsim/skycat.py +++ b/imsim/skycat.py @@ -137,6 +137,16 @@ def objects(self): self.file_name, skycatalog_root=self.skycatalog_root) self._objects = sky_cat.get_objects_by_region( region, obj_type_set=self.obj_types, mjd=self.mjd) + if self.multi_det: + # Work out the which of the adjacent detectors are closest + # to each object found in the region. + closest_detector = {} + for obj in self._objects: + world_pos = galsim.CelestialCoord(obj.ra*galsim.degrees, + obj.dec*galsim.degrees) + obj_dist = {det: detector_centers[det].distanceTo(world_pos) for det in adjacent_detectors} + closest_detector[obj] = min(obj_dist, key=obj_dist.get) + self._objects = filter(lambda object: closest_detector[object] == self.det_name, self._objects) if not self._objects: self.logger.warning("No objects found on image.") return self._objects From b050370cb09d5461b1a526d8ff976df7f30b9a2d Mon Sep 17 00:00:00 2001 From: Jim Chiang Date: Mon, 15 Sep 2025 13:05:01 -0700 Subject: [PATCH 20/36] keep track of objects to skip instead of trying to use an iterator provided by the filter built-in --- imsim/skycat.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/imsim/skycat.py b/imsim/skycat.py index 24c2370e..d5f1b778 100644 --- a/imsim/skycat.py +++ b/imsim/skycat.py @@ -105,6 +105,7 @@ def __init__(self, file_name, wcs, band, mjd, fiducial_wcs=None, self.logger.info(f'Object types restricted to {obj_types}') self.ccd_center = wcs.toWorld(galsim.PositionD(xsize/2.0, ysize/2.0)) self._objects = None + self._objects_to_skip = set() @property def objects(self): @@ -139,14 +140,18 @@ def objects(self): region, obj_type_set=self.obj_types, mjd=self.mjd) if self.multi_det: # Work out the which of the adjacent detectors are closest - # to each object found in the region. - closest_detector = {} + # to each object found in the region. Add the object to the + # set of objects to skip if the current detector isn't the + # closest one. for obj in self._objects: world_pos = galsim.CelestialCoord(obj.ra*galsim.degrees, obj.dec*galsim.degrees) obj_dist = {det: detector_centers[det].distanceTo(world_pos) for det in adjacent_detectors} - closest_detector[obj] = min(obj_dist, key=obj_dist.get) - self._objects = filter(lambda object: closest_detector[object] == self.det_name, self._objects) + closest_detector = min(obj_dist, key=obj_dist.get) + if closest_detector != self.det_name: + self._objects_to_skip.add(obj.id) + self.logger.info("Skipping {len(self._objects_to_skip)} " + "objects for this detector.") if not self._objects: self.logger.warning("No objects found on image.") return self._objects @@ -210,6 +215,8 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): raise RuntimeError("Trying to get an object from an empty sky catalog") skycat_obj = self.objects[index] + if skycat_obj.id in self._objects_to_skip: + return None gsobjs = skycat_obj.get_gsobject_components(gsparams) if self.apply_dc2_dilation and skycat_obj.object_type == 'galaxy': From bfcb934c43b63d91df6b92a8b951ea06535f3201 Mon Sep 17 00:00:00 2001 From: Jim Chiang Date: Tue, 16 Sep 2025 09:53:49 -0700 Subject: [PATCH 21/36] include f-string prefix --- imsim/skycat.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imsim/skycat.py b/imsim/skycat.py index d5f1b778..2b0cec58 100644 --- a/imsim/skycat.py +++ b/imsim/skycat.py @@ -150,7 +150,7 @@ def objects(self): closest_detector = min(obj_dist, key=obj_dist.get) if closest_detector != self.det_name: self._objects_to_skip.add(obj.id) - self.logger.info("Skipping {len(self._objects_to_skip)} " + self.logger.info(f"Skipping {len(self._objects_to_skip)} " "objects for this detector.") if not self._objects: self.logger.warning("No objects found on image.") From 36605bd717b3c09556f0f924925af0c042ad710a Mon Sep 17 00:00:00 2001 From: Jim Chiang Date: Mon, 29 Sep 2025 15:23:33 -0700 Subject: [PATCH 22/36] don't run AddNoise for pass 2 rendering --- imsim/full_focal_plane.py | 8 ++++++++ imsim/lsst_image.py | 5 ++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/imsim/full_focal_plane.py b/imsim/full_focal_plane.py index 5c5484c6..d5f6bf83 100644 --- a/imsim/full_focal_plane.py +++ b/imsim/full_focal_plane.py @@ -137,6 +137,14 @@ def getKwargs(self, config, base, logger): # photons which fell off-detector. class LSST_FocalPlaneImageBuilder(LSST_ImageBuilderBase): + def setup(self, config, base, image_num, obj_num, ignore, logger): + xsize, ysize = super().setup(config, base, image_num, obj_num, ignore, logger) + # Disable application of the addNoise function in + # LSST_ImageBuilderBase.addNoise since it was already applied + # in pass1 and not needed for the added off-detector photons + self.add_noise = False + return xsize, ysize + def buildImage(self, config, base, image_num, _obj_num, logger): """Draw the off-detector photons to the image. """ diff --git a/imsim/lsst_image.py b/imsim/lsst_image.py index da85cf99..b936f266 100644 --- a/imsim/lsst_image.py +++ b/imsim/lsst_image.py @@ -88,6 +88,8 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): else: self.boresight = params.get('boresight') + self.add_noise = True + self.camera_name = params.get('camera', 'LsstCamSim') # If using a SiliconSensor and not overridden in config, determine sensor type to use. @@ -197,7 +199,8 @@ def addNoise(self, image, config, base, image_num, obj_num, current_var, logger) sky.array[:] *= fringing_map image += sky - AddNoise(base,image,current_var,logger) + if self.add_noise: + AddNoise(base,image,current_var,logger) @staticmethod def _create_full_image(config, base): From b6fe71c390bd1b356f8e6bff60b97f11ebd155ef Mon Sep 17 00:00:00 2001 From: William Lucas Date: Fri, 17 Oct 2025 15:20:50 +0100 Subject: [PATCH 23/36] Config now keys off-detector photons to det_num/name. --- examples/imsim-user-instcat-ffp-pass2.yaml | 34 ++++++++++++++++++++-- imsim/full_focal_plane.py | 14 +++++++-- 2 files changed, 42 insertions(+), 6 deletions(-) diff --git a/examples/imsim-user-instcat-ffp-pass2.yaml b/examples/imsim-user-instcat-ffp-pass2.yaml index 7a0bb4a7..82faf26a 100644 --- a/examples/imsim-user-instcat-ffp-pass2.yaml +++ b/examples/imsim-user-instcat-ffp-pass2.yaml @@ -26,11 +26,39 @@ input.tree_rings: "" input.atm_psf: "" input.initial_image: - file_name: "eimage_00398414-0-r-R22_S12-det095.fits" - dir: "output" + dir: output + file_name: + type: FormattedStr + format : eimage_%08d-%1d-%s-%s-det%03d.fits + items: + - { type: OpsimData, field: observationId } + - { type: OpsimData, field: snap } + - $band + - $det_name + - "@output.det_num" +# Simple off detector photon input reading one file. +# input.off_detector_photons: +# file_name: "output/off_det_00398414-0-r-R22_S11-det094.fits" +# camera: $camera +# det_name: $det_name + +# If we're doing the preprocessing step then we've already reorganised the +# photons so they live in one file per destination detector, so we only need to +# read that one. But, we might feasibly also wish to directly specify a list of +# off-detector photon files, which presumably haven't been preprocessed, so I +# might want to add the capability to work from a file list too. input.off_detector_photons: - file_name: "output/off_det_00398414-0-r-R22_S11-det094.fits" + dir: output + file_name: + type: FormattedStr + format : off_det_%08d-%1d-%s-%s-det%03d.fits + items: + - { type: OpsimData, field: observationId } + - { type: OpsimData, field: snap } + - $band + - $det_name + - "@output.det_num" camera: $camera det_name: $det_name diff --git a/imsim/full_focal_plane.py b/imsim/full_focal_plane.py index d5f6bf83..d109ffa8 100644 --- a/imsim/full_focal_plane.py +++ b/imsim/full_focal_plane.py @@ -1,3 +1,5 @@ +import os + import galsim from galsim.config.extra import ExtraOutputBuilder, RegisterExtraOutput from galsim.config import RegisterImageType, InputLoader, RegisterInputType @@ -85,7 +87,7 @@ class OffDetectorPhotons(object): then, and will now be read in this second pass to be accumulated on other sensors. """ - def __init__(self, file_name, camera, det_name, xsize=4096, ysize=4096, logger=None): + def __init__(self, file_name, camera, det_name, dir=None, xsize=4096, ysize=4096, logger=None): """ Initialize the off-detector photons input class. @@ -96,6 +98,8 @@ def __init__(self, file_name, camera, det_name, xsize=4096, ysize=4096, logger=N The name of the camera containing the detector. det_name: str The name of the detector to use for photon coordinate transformations. + dir: str + The directory where the file is. (default: None) xsize: int The x size in pixels of the CCD. (default: 4096) ysize: int @@ -104,6 +108,8 @@ def __init__(self, file_name, camera, det_name, xsize=4096, ysize=4096, logger=N A logger object. (default: None) """ self.file_name = file_name + if dir is not None: + self.file_name = os.path.join(dir, self.file_name) self.det = get_camera(camera)[det_name] self.xsize = xsize self.ysize = ysize @@ -119,7 +125,7 @@ class OffDetectorPhotonsLoader(InputLoader): """ def getKwargs(self, config, base, logger): req = {'file_name': str, 'camera': str, 'det_name': str} - opt = {} + opt = {'dir': str} kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) kwargs['xsize'] = base.get('det_xsize', 4096) @@ -128,7 +134,9 @@ def getKwargs(self, config, base, logger): # For now assume the off-detector photons on file can be used for # all detectors. - safe = True + # Change assumption, so now we consider that we have a per-detector + # input file that cannot be reused for other detectors. + safe = False return kwargs, safe From c1466fdbb1001f6b72f5c6885f101807fe8f0f73 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Tue, 17 Feb 2026 16:08:09 +0000 Subject: [PATCH 24/36] Fix test_trimmer.py to set edge_pix=100, the old default, to keep a test useful. (Default is now 200 pixels to cover the gap between detectors.) --- config/imsim-config-instcat.yaml | 2 +- tests/test_trimmer.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/config/imsim-config-instcat.yaml b/config/imsim-config-instcat.yaml index 2f51f5c5..3e28ff11 100644 --- a/config/imsim-config-instcat.yaml +++ b/config/imsim-config-instcat.yaml @@ -12,7 +12,7 @@ # Get most of the configuration from the base imSim config template. modules: - imsim -template: imsim-config-photon-pooling +template: imsim-config #################################################################### # The following entries are added to the base configuration above. diff --git a/tests/test_trimmer.py b/tests/test_trimmer.py index 698b9e23..13faf4e2 100644 --- a/tests/test_trimmer.py +++ b/tests/test_trimmer.py @@ -45,8 +45,8 @@ def test_InstCatTrimmer(self): instcat = imsim.InstCatalog(instcat_file, wcs, edge_pix=1000, skip_invalid=False) self.assertEqual(instcat.nobjects, 24) - # With the default edge_pix=100, only 17 make the cut. - instcat = imsim.InstCatalog(instcat_file, wcs, skip_invalid=False) + # With a smaller than default edge_pix=100, only 17 make the cut. + instcat = imsim.InstCatalog(instcat_file, wcs, edge_pix=100, skip_invalid=False) self.assertEqual(instcat.nobjects, 17) # Check the application of min_source. From 7a365bed2a80bcb9f1e76c3c71bc58b7148d21c6 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Fri, 27 Feb 2026 14:00:08 +0000 Subject: [PATCH 25/36] New input.instance_catalog.full_focal_plane enables and disables FFP. --- imsim/instcat.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/imsim/instcat.py b/imsim/instcat.py index 956aad9f..04623da5 100644 --- a/imsim/instcat.py +++ b/imsim/instcat.py @@ -177,6 +177,7 @@ def __init__(self, file_name, wcs, fiducial_wcs=None, fiducial_det_name=None, edge_pix=200, camera=None, det_name=None, sort_mag=True, flip_g2=True, pupil_area=RUBIN_AREA, min_source=None, skip_invalid=True, + full_focal_plane=False, logger=None): logger = galsim.config.LoggerWrapper(logger) self.file_name = file_name @@ -184,15 +185,16 @@ def __init__(self, file_name, wcs, fiducial_wcs=None, fiducial_det_name=None, self.xsize = xsize self.ysize = ysize self.edge_pix = edge_pix - if camera is not None and det_name is not None and fiducial_wcs is not None: - self.multi_det = True + if full_focal_plane: + if camera is None or det_name is None or fiducial_wcs is None: + raise galsim.GalSimConfigValueError("full_focal_plane is set to True: camera, det_name and fiducial_wcs must be provided.") + self.full_focal_plane = True self.camera = Camera(camera_class=camera) self.det_name = det_name self.fiducial_wcs = fiducial_wcs self.fiducial_det_name = fiducial_det_name else: - # Set multi_det to False and move on. - self.multi_det = False + self.full_focal_plane = False self.sort_mag = sort_mag self.flip_g2 = flip_g2 self.min_source = min_source @@ -242,7 +244,7 @@ def _load(self, logger=None): nuse = 0 ntot = 0 - if self.multi_det: + if self.full_focal_plane: # Get adjacent detectors; then their centers in pixel coordinates, then # store after converting to RA and dec. adjacent_detectors = self.camera.get_adjacent_detectors(self.det_name) @@ -275,7 +277,7 @@ def _load(self, logger=None): world_pos = galsim.CelestialCoord(ra, dec) #logger.debug('world_pos = %s',world_pos) - if self.multi_det: + if self.full_focal_plane: # Check whether the object is closer to this detector's center # than any of the adjacent detectors. If not, then skip it. obj_dist = {det: detector_centers[det].distanceTo(world_pos) for det in adjacent_detectors} @@ -676,6 +678,7 @@ def getKwargs(self, config, base, logger): 'pupil_area' : float, 'min_source' : int, 'skip_invalid' : bool, + 'full_focal_plane' : bool, } ignore = ['fiducial_wcs'] # Handled separately. kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore) From b78dfc64580745eb8cd30af39d985ceda6d57a9f Mon Sep 17 00:00:00 2001 From: William Lucas Date: Mon, 9 Mar 2026 14:35:04 +0000 Subject: [PATCH 26/36] Have test_image_nobjects set edge_pix=100 as the new default of 200 (to fill detector gaps) changes results. --- tests/test_lsst_image.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_lsst_image.py b/tests/test_lsst_image.py index 8f59a353..9fdf1c28 100644 --- a/tests/test_lsst_image.py +++ b/tests/test_lsst_image.py @@ -35,6 +35,7 @@ def test_image_nobjects(): config = {'modules': ['imsim'], 'template': template, 'input.instance_catalog.file_name': instcat_file, + 'input.instance_catalog.edge_pix': 100, 'input.opsim_data.file_name': instcat_file, 'input.tree_rings': '', 'input.atm_psf': '', From 75cd0622ac7f633aaee86b79238c7597e29e151f Mon Sep 17 00:00:00 2001 From: William Lucas Date: Mon, 9 Mar 2026 14:35:53 +0000 Subject: [PATCH 27/36] Calculate fiducial WCS before image WCS in order to correctly store _icrf_to_field in base config. --- imsim/instcat.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/imsim/instcat.py b/imsim/instcat.py index 04623da5..34e66236 100644 --- a/imsim/instcat.py +++ b/imsim/instcat.py @@ -682,8 +682,6 @@ def getKwargs(self, config, base, logger): } ignore = ['fiducial_wcs'] # Handled separately. kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore) - wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) - kwargs['wcs'] = wcs if config.get('fiducial_wcs', None) is not None: fiducial_wcs = galsim.config.BuildWCS(config, 'fiducial_wcs', base, logger=logger) if fiducial_wcs.wcs_type == "TAN-SIP": @@ -695,6 +693,10 @@ def getKwargs(self, config, base, logger): kwargs['fiducial_det_name'] = config['fiducial_wcs']['det_name'] else: kwargs['fiducial_wcs'] = None + # Must build the 'real' WCS after the fiducial WCS as BatoidWCSBuilder.buildWCS() + # stores base['_icrf_to_field'], and we want that to be the 'real' one. + wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) + kwargs['wcs'] = wcs kwargs['xsize'] = base.get('det_xsize', 4096) kwargs['ysize'] = base.get('det_ysize', 4096) kwargs['camera'] = base['output'].get('camera', None) From 68b36b54d507430552a34c01c60ec89c0e16de85 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Wed, 11 Mar 2026 14:15:53 +0000 Subject: [PATCH 28/36] Add full_focal_plane parameter to skycat and ensure correct _icrf_to_field is stored. --- imsim/skycat.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/imsim/skycat.py b/imsim/skycat.py index 2b0cec58..c8c511f4 100644 --- a/imsim/skycat.py +++ b/imsim/skycat.py @@ -20,7 +20,8 @@ def __init__(self, file_name, wcs, band, mjd, fiducial_wcs=None, obj_types=None, skycatalog_root=None, edge_pix=200, camera=None, det_name=None, pupil_area=RUBIN_AREA, max_flux=None, logger=None, - apply_dc2_dilation=False, approx_nobjects=None): + apply_dc2_dilation=False, approx_nobjects=None, + full_focal_plane=False): """ Parameters ---------- @@ -52,7 +53,7 @@ def __init__(self, file_name, wcs, band, mjd, fiducial_wcs=None, Size in pixels of the buffer region around nominal image to consider objects. camera : str [None] - Name of the imSim camera to use for a multi-detector run. + Name of the imSim camera to use for a full_focal_plane run. det_name : str [None] Name of the current detector. pupil_area : float [RUBIN_AREA] @@ -73,6 +74,8 @@ def __init__(self, file_name, wcs, band, mjd, fiducial_wcs=None, Approximate number of objects per CCD used by galsim to set up the image processing. If None, then the actual number of objects found by skyCatalogs, via .getNObjects, will be used. + full_focal_plane : bool [False] + Enables or disables two-pass full focal plane processing. """ self.file_name = file_name self.wcs = wcs @@ -86,15 +89,16 @@ def __init__(self, file_name, wcs, band, mjd, fiducial_wcs=None, else: self.skycatalog_root = skycatalog_root self.edge_pix = edge_pix - if camera is not None and det_name is not None and fiducial_wcs is not None: - self.multi_det = True + if full_focal_plane: + if camera is None or det_name is None or fiducial_wcs is None: + raise galsim.GalSimConfigValueError("full_focal_plane is set to True: camera, det_name and fiducial_wcs must be provided.") + self.full_focal_plane = True self.camera = Camera(camera_class=camera) self.det_name = det_name self.fiducial_wcs = fiducial_wcs self.fiducial_det_name = fiducial_det_name else: - # Set multi_det to False and move on. - self.multi_det = False + self.full_focal_plane = False self.pupil_area = pupil_area self.max_flux = max_flux self.logger = galsim.config.LoggerWrapper(logger) @@ -110,7 +114,7 @@ def __init__(self, file_name, wcs, band, mjd, fiducial_wcs=None, @property def objects(self): if self._objects is None: - if self.multi_det: + if self.full_focal_plane: # Get adjacent detectors; then their centers in pixel coordinates, then # store after converting to RA and dec. adjacent_detectors = self.camera.get_adjacent_detectors(self.det_name) @@ -138,7 +142,7 @@ def objects(self): self.file_name, skycatalog_root=self.skycatalog_root) self._objects = sky_cat.get_objects_by_region( region, obj_type_set=self.obj_types, mjd=self.mjd) - if self.multi_det: + if self.full_focal_plane: # Work out the which of the adjacent detectors are closest # to each object found in the region. Add the object to the # set of objects to skip if the current detector isn't the @@ -269,12 +273,11 @@ def getKwargs(self, config, base, logger): 'approx_nobjects': int, 'mjd': float, 'pupil_area': float, + 'full_focal_plane': bool, } ignore = ['fiducial_wcs'] # Handled separately. kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore) - wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) - kwargs['wcs'] = wcs if config.get('fiducial_wcs', None) is not None: fiducial_wcs = galsim.config.BuildWCS(config, 'fiducial_wcs', base, logger=logger) if fiducial_wcs.wcs_type == "TAN-SIP": @@ -286,6 +289,10 @@ def getKwargs(self, config, base, logger): kwargs['fiducial_det_name'] = config['fiducial_wcs']['det_name'] else: kwargs['fiducial_wcs'] = None + # Must build the 'real' WCS after the fiducial WCS as BatoidWCSBuilder.buildWCS() + # stores base['_icrf_to_field'], and we want that to be the 'real' one. + wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) + kwargs['wcs'] = wcs kwargs['xsize'] = base.get('det_xsize', 4096) kwargs['ysize'] = base.get('det_ysize', 4096) kwargs['camera'] = base['output'].get('camera', None) From 5f8d821b02f97cf2c8d65423a776d78b0427d22b Mon Sep 17 00:00:00 2001 From: William Lucas Date: Mon, 16 Mar 2026 15:18:11 +0000 Subject: [PATCH 29/36] Allow off_detector_photons to be read from multiple files. --- examples/imsim-user-instcat-ffp-pass2.yaml | 15 +++-- imsim/full_focal_plane.py | 72 +++++++++++++--------- 2 files changed, 55 insertions(+), 32 deletions(-) diff --git a/examples/imsim-user-instcat-ffp-pass2.yaml b/examples/imsim-user-instcat-ffp-pass2.yaml index 82faf26a..89d153c0 100644 --- a/examples/imsim-user-instcat-ffp-pass2.yaml +++ b/examples/imsim-user-instcat-ffp-pass2.yaml @@ -50,15 +50,22 @@ input.initial_image: # might want to add the capability to work from a file list too. input.off_detector_photons: dir: output - file_name: + # file_name: + # type: FormattedStr + # format: off_det_%08d-%1d-%s-%s-det%03d.fits + # items: + # - { type: OpsimData, field: observationId } + # - { type: OpsimData, field: snap } + # - $band + # - $det_name + # - "@output.det_num" + file_pattern: type: FormattedStr - format : off_det_%08d-%1d-%s-%s-det%03d.fits + format: off_det_%08d-%1d-%s-DETNAME-detDETNUM.fits items: - { type: OpsimData, field: observationId } - { type: OpsimData, field: snap } - $band - - $det_name - - "@output.det_num" camera: $camera det_name: $det_name diff --git a/imsim/full_focal_plane.py b/imsim/full_focal_plane.py index d109ffa8..ee2f46c9 100644 --- a/imsim/full_focal_plane.py +++ b/imsim/full_focal_plane.py @@ -87,55 +87,73 @@ class OffDetectorPhotons(object): then, and will now be read in this second pass to be accumulated on other sensors. """ - def __init__(self, file_name, camera, det_name, dir=None, xsize=4096, ysize=4096, logger=None): + def __init__(self, camera, det_name, file_name=None, file_pattern=None, dir=None, logger=None): """ Initialize the off-detector photons input class. Parameters: - file_name: str - The name of the file to read. camera: str The name of the camera containing the detector. det_name: str The name of the detector to use for photon coordinate transformations. + file_name: str + The name of the file to read. (default: None) + file_pattern: str + The file pattern to use if up to 189 are to be read. (default: None) dir: str The directory where the file is. (default: None) - xsize: int - The x size in pixels of the CCD. (default: 4096) - ysize: int - The y size in pixels of the CCD. (default: 4096) logger: logging.logger A logger object. (default: None) """ - self.file_name = file_name - if dir is not None: - self.file_name = os.path.join(dir, self.file_name) - self.det = get_camera(camera)[det_name] - self.xsize = xsize - self.ysize = ysize self.logger = logger - # self.photons = None - self.photons = galsim.PhotonArray.read(self.file_name) + camera = get_camera(camera) + self.det = camera[det_name] + # The Loader should have ensured we have one but not both of file_name or file_pattern. + if file_name is not None: + # Only have a single file to read. + file_names = [file_name] + elif file_pattern is not None: + # We want to read up to 188 files: potentially photons from every + # detector except this one (as those photons were drawn normally in + # the first pass). We may want to allow det_num_start and + # det_num_end to be set in the config in a similar way to output + # where it's controlled with det_num.first, nfiles and only_dets. + det_num_start = 0 + det_num_end = 189 + file_names = [file_pattern.replace("DETNAME", camera[i].getName()).replace("DETNUM", "{:03d}".format(i)) + for i in range(det_num_start, det_num_end) + if camera[i].getName() != det_name + ] + if dir is not None: + file_names = [os.path.join(dir, fname) for fname in file_names] + # Allow for the possibility that not all the 188 files exist. + self.file_names = [fname for fname in file_names if os.path.isfile(fname)] + if len(self.file_names) == 0: + raise galsim.GalSimConfigValueError(f"Detector {det_name} did not find any expected off detector photon files: {file_names}.") + logger.info(f"Detector {det_name} reading {len(self.file_names)} files for off-detector photons.") + # Would it be better to make photons a @property and read lazily when + # first accessed? + photon_storage = [] + for file_name in self.file_names: + photon_storage.append(galsim.PhotonArray.read(file_name)) + self.photons = galsim.PhotonArray.concatenate(photon_storage) self.photons.x, self.photons.y = focal_to_pixel(self.photons.x, self.photons.y, self.det) class OffDetectorPhotonsLoader(InputLoader): - """ - Class to load off-detector photons from file. + """Class to load off-detector photons from file. """ def getKwargs(self, config, base, logger): - req = {'file_name': str, 'camera': str, 'det_name': str} - opt = {'dir': str} + req = {'camera': str, 'det_name': str} + opt = {'file_name': str, 'file_pattern': str, 'dir': str} kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) - kwargs['xsize'] = base.get('det_xsize', 4096) - kwargs['ysize'] = base.get('det_ysize', 4096) + if ('file_name' in kwargs) == ('file_pattern' in kwargs): + raise galsim.GalSimConfigValueError('off_detector_photons must contain one of either file_name or file_pattern') kwargs['logger'] = logger - # For now assume the off-detector photons on file can be used for - # all detectors. - # Change assumption, so now we consider that we have a per-detector - # input file that cannot be reused for other detectors. + # Base assumption (for now) is that we have a per-detector set of inputs + # that cannot be reused for other detectors. safe = False return kwargs, safe @@ -162,16 +180,14 @@ def buildImage(self, config, base, image_num, _obj_num, logger): if not isinstance(base['_input_objs']['off_detector_photons'][0], OffDetectorPhotons): raise galsim.config.GalSimConfigError( "When using LSST_FocalPlaneImage, you must provide an off_detector_photons input.",) - + off_detector_photons = base['_input_objs']['off_detector_photons'][0].photons if len(off_detector_photons) > 0: # There are photons to be accumulated. They should already have been # transformed to pixel coordinates when they were read from file, so # go ahead and accumulate them. - sensor = base.get('sensor', Sensor()) - sensor.accumulate(off_detector_photons, image) return image, [] From 85272cff5253e8c0eb8dee6e129dfa338e166974 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Mon, 16 Mar 2026 17:04:59 +0000 Subject: [PATCH 30/36] Update tests to use LsstCamSim. --- tests/test_camera.py | 16 ++++++++-------- tests/test_full_focal_plane.py | 2 +- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/test_camera.py b/tests/test_camera.py index d5426581..ba129fc8 100644 --- a/tests/test_camera.py +++ b/tests/test_camera.py @@ -58,7 +58,7 @@ def test_get_adjacent_detectors_central(self): 'R22_S10', 'R22_S11', 'R22_S12', 'R22_S20', 'R22_S21', 'R22_S22', ] - self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) def test_get_adjacent_detectors_offset(self): # Case of detectors adjacent to one offset to a corner in a raft. @@ -69,7 +69,7 @@ def test_get_adjacent_detectors_offset(self): 'R21_S22', 'R22_S20', 'R22_S21', 'R21_S12', 'R22_S10', 'R22_S11', ] - self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) # Bottom right corner. det_name = 'R22_S02' expected_detectors = [ @@ -77,7 +77,7 @@ def test_get_adjacent_detectors_offset(self): 'R22_S01', 'R22_S02', 'R23_S00', 'R12_S21', 'R12_S22', 'R13_S20', ] - self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) def test_get_adjacent_detectors_edge(self): # Case of detectors adjacent to one on the edge of the camera. @@ -87,7 +87,7 @@ def test_get_adjacent_detectors_edge(self): 'R02_S10', 'R02_S11', 'R02_S12', 'R02_S00', 'R02_S01', 'R02_S02', ] - self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) # Left edge. det_name = 'R20_S10' expected_detectors = [ @@ -108,7 +108,7 @@ def test_get_adjacent_detectors_edge(self): 'R24_S11', 'R24_S12', 'R24_S01', 'R24_S02', ] - self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) def test_get_adjacent_detectors_corner(self): # Case of detectors adjacent to one in a corner of the camera. @@ -120,7 +120,7 @@ def test_get_adjacent_detectors_corner(self): 'R01_S10', 'R01_S11', 'R01_S00', 'R01_S01', ] - self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) # On the corner of a corner cutout. det_name = 'R11_S00' expected_detectors = [ @@ -128,12 +128,12 @@ def test_get_adjacent_detectors_corner(self): 'R10_S02', 'R11_S00', 'R11_S01', 'R01_S20', 'R01_S21', ] - self.check_adjacent_detectors_match('LsstCam', det_name, expected_detectors) + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) def test_get_adjacent_detectors_exceptions(self): # Ensure that get_adjacent_detectors raises errors for invalid detecor # names -- either ones which don't exist, or non science ones. - camera = imsim.Camera('LsstCam') + camera = imsim.Camera('LsstCamSim') # Non-existent detectors. with self.assertRaises(GalSimValueError): diff --git a/tests/test_full_focal_plane.py b/tests/test_full_focal_plane.py index c48ea793..220e4717 100644 --- a/tests/test_full_focal_plane.py +++ b/tests/test_full_focal_plane.py @@ -51,7 +51,7 @@ def process_off_detector_photons_output(pa_list): base = { "det_name": "R22_S11", "output": { - "camera": "LsstCam", + "camera": "LsstCamSim", "off_detector_photons": { "dir": "output", "file_name": "off_det_photons.fits", From 0292416e501675d27d6c39156bb8b2498a9397cb Mon Sep 17 00:00:00 2001 From: William Lucas Date: Tue, 17 Mar 2026 09:40:24 +0000 Subject: [PATCH 31/36] Fix edge_pix values used in test_stamp.py to old value of 100. --- tests/test_stamp.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_stamp.py b/tests/test_stamp.py index 9b691b5d..60a56236 100644 --- a/tests/test_stamp.py +++ b/tests/test_stamp.py @@ -204,6 +204,7 @@ def test_stamp_sizes(): input.sky_catalog.pupil_area: type: Eval str: "0.25 * np.pi * 649**2" + input.sky_catalog.edge_pix: 100 input.opsim_data.file_name: data/small_opsim_9683.db input.opsim_data.visit: 449053 input.tree_rings.only_dets: [R22_S11] @@ -426,6 +427,7 @@ def test_faint_high_redshift_stamp(): input.sky_catalog.file_name: data/sky_cat_9683.yaml input.sky_catalog.obj_types: [galaxy] input.sky_catalog.band: u + input.sky_catalog.edge_pix: 100 input.opsim_data.file_name: data/small_opsim_9683.db input.opsim_data.visit: 449053 input.tree_rings.only_dets: [R22_S11] From 5342e861a304ded89fa8af7ccfd15f5017810f5e Mon Sep 17 00:00:00 2001 From: William Lucas Date: Tue, 24 Mar 2026 09:51:08 +0000 Subject: [PATCH 32/36] Add tests for off-detector photon loader. --- imsim/full_focal_plane.py | 10 ++- tests/test_full_focal_plane.py | 152 +++++++++++++++++++++++++++++++-- 2 files changed, 153 insertions(+), 9 deletions(-) diff --git a/imsim/full_focal_plane.py b/imsim/full_focal_plane.py index ee2f46c9..e07d0fd8 100644 --- a/imsim/full_focal_plane.py +++ b/imsim/full_focal_plane.py @@ -46,7 +46,8 @@ def gather_out_of_bounds_photons(image_bounds, photons): # the images again and draw the photons on top of the first pass image. class OffDetectorPhotonsBuilder(ExtraOutputBuilder): """Build photon arrays containing the off-detector photons found during an - image's construction and write them to file. + image's construction and write them to file. This item will typically be + included in the config for the first pass of a full focal plane run. """ def processImage(self, index, obj_nums, config, base, logger): @@ -129,7 +130,8 @@ def __init__(self, camera, det_name, file_name=None, file_pattern=None, dir=None # Allow for the possibility that not all the 188 files exist. self.file_names = [fname for fname in file_names if os.path.isfile(fname)] if len(self.file_names) == 0: - raise galsim.GalSimConfigValueError(f"Detector {det_name} did not find any expected off detector photon files: {file_names}.") + raise galsim.GalSimConfigValueError( + f"Detector {det_name} did not find any expected off detector photon files: {file_names}.") logger.info(f"Detector {det_name} reading {len(self.file_names)} files for off-detector photons.") # Would it be better to make photons a @property and read lazily when # first accessed? @@ -142,6 +144,8 @@ def __init__(self, camera, det_name, file_name=None, file_pattern=None, dir=None class OffDetectorPhotonsLoader(InputLoader): """Class to load off-detector photons from file. + These options should be provided in the config for the second pass of a full + focal plane run. """ def getKwargs(self, config, base, logger): req = {'camera': str, 'det_name': str} @@ -171,7 +175,7 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): self.add_noise = False return xsize, ysize - def buildImage(self, config, base, image_num, _obj_num, logger): + def buildImage(self, config, base, image_num, obj_num, logger): """Draw the off-detector photons to the image. """ # Make sure we have an input image and off-detector photons to draw. diff --git a/tests/test_full_focal_plane.py b/tests/test_full_focal_plane.py index 220e4717..8d0b8e69 100644 --- a/tests/test_full_focal_plane.py +++ b/tests/test_full_focal_plane.py @@ -2,11 +2,15 @@ import galsim import os +from pathlib import Path +import pytest from imsim.camera import get_camera from imsim.full_focal_plane import gather_out_of_bounds_photons from imsim.utils import focal_to_pixel +DATA_DIR = Path(__file__).parent / 'data' +OUTPUT_DIR = Path(__file__).parent / 'output' def test_gather_out_of_bounds_photons(): """Make sure that photons inside and outside the image bounds are gathered correctly.""" @@ -53,7 +57,7 @@ def process_off_detector_photons_output(pa_list): "output": { "camera": "LsstCamSim", "off_detector_photons": { - "dir": "output", + "dir": OUTPUT_DIR, "file_name": "off_det_photons.fits", }, }, @@ -69,11 +73,11 @@ def process_off_detector_photons_output(pa_list): galsim.config.WriteExtraOutputs(base, []) # Check that the output file exists. - assert os.path.exists("output/off_det_photons.fits"), "Output file 'output/off_det_photons.fits' was not created." + assert os.path.exists(OUTPUT_DIR / "off_det_photons.fits"), "Output file 'output/off_det_photons.fits' was not created." # Assert that the photons stored in the output match. They should need to be # transformed from focal plane to detector coordinates first. - output_photons = galsim.PhotonArray.read("output/off_det_photons.fits") + output_photons = galsim.PhotonArray.read(OUTPUT_DIR / "off_det_photons.fits") if len(output_photons) > 0: detector = get_camera(base['output']['camera'])[base['det_name']] output_photons.x, output_photons.y = focal_to_pixel(output_photons.x, output_photons.y, detector) @@ -83,7 +87,7 @@ def process_off_detector_photons_output(pa_list): assert np.allclose(output_photons.flux, focal_plane_photons.flux), "Flux of output photons does not match original photons." # Final cleanup. - os.remove("output/off_det_photons.fits") + os.remove(OUTPUT_DIR / "off_det_photons.fits") def test_off_detector_photons_output(): @@ -125,8 +129,144 @@ def test_off_detector_photons_output_with_mix(): for nphot in n_off_det] process_off_detector_photons_output(pa_list) -# TODO: Need some tests for the OffDetectorPhotonsLoader and -# LSST_FocalPlaneImageBuilder. + +def create_base_off_detector_config(det_name="R22_S11"): + base = { + "input": { + "off_detector_photons": { + "camera": "LsstCamSim", + "det_name": det_name, + "dir": DATA_DIR, + } + }, + } + return base + + +def process_single_off_detector_load(det_name, file_name, expected_photons): + config = create_base_off_detector_config(det_name) + config["input"]["off_detector_photons"]["file_name"] = file_name + galsim.config.ProcessInput(config) + photons = config['_input_objs']['off_detector_photons'][0].photons + np.testing.assert_allclose(photons.x, expected_photons.x, rtol=1.e-10) + np.testing.assert_allclose(photons.y, expected_photons.y, rtol=1.e-10) + np.testing.assert_allclose(photons.flux, expected_photons.flux, rtol=1.e-10) + + +def test_load_photons_from_single_file(): + # Test that the OffDetectorPhotonsLoader can load photons using the 'file_name' + # config option, working as both R22_S11 and R22_S12 to ensure photons + # are transformed from focal plane to detector coordinates. + expected_photons = galsim.PhotonArray.fromArrays(x=np.arange(1000., 1010.), + y=np.arange(2000., 2010.), + flux=np.ones(10), + ) + process_single_off_detector_load("R22_S11", "off_det_00398414-0-r-R22_S11-det094.fits", expected_photons) + expected_photons = galsim.PhotonArray.fromArrays(x=np.arange(-3000., -3010., -1.), + y=np.arange(-4000., -4010., -1), + flux=np.ones(10), + ) + process_single_off_detector_load("R22_S12", "off_det_00398414-0-r-R22_S12-det095.fits", expected_photons) + + +def test_load_photons_from_multiple_files(): + # Test that the OffDetectorPhotonsLoader can load photons from multiple files + # given a file name pattern. + config = create_base_off_detector_config() + config["input"]["off_detector_photons"]["file_pattern"] = "off_det_multi_00398414-0-r-DETNAME-detDETNUM.fits" + # Processing this config should read the photons from the two files + # data/off_det_multi_00398414-0-r-R22_S10-det093.fits and + # data/off_det_multi_00398414-0-r-R22_S12-det095.fits, concatenate them, and + # transform them to the R22_S11 detector coordinates. + galsim.config.ProcessInput(config) + photons = config['_input_objs']['off_detector_photons'][0].photons + + # We expect the photons from the two files to have been transformed to the + # following in the R22_S11 detector coordinates. + x = np.concatenate((np.arange(1000., 1010.), np.arange(5000., 5010.))) + y = np.concatenate((np.arange(2000., 2010.), np.arange(10000., 10010.))) + expected_photons = galsim.PhotonArray.fromArrays(x=x, + y=y, + flux=np.ones(20), + ) + np.testing.assert_allclose(photons.x, expected_photons.x, rtol=1.e-10) + np.testing.assert_allclose(photons.y, expected_photons.y, rtol=1.e-10) + np.testing.assert_allclose(photons.flux, expected_photons.flux, rtol=1.e-10) + + +def test_off_detector_loader_file_naming(): + # Make sure that OffDetectorPhotonsLoader raises an error if neither a file + # name nor a file pattern is provided. + config = create_base_off_detector_config() + with pytest.raises(galsim.GalSimConfigError, + match='off_detector_photons must contain one of either file_name or file_pattern in config'): + galsim.config.ProcessInput(config) + + # Make sure the same error is raised if both are provided. + config = create_base_off_detector_config() + config["input"]["off_detector_photons"]["file_name"] = "some_file_name.fits" + config["input"]["off_detector_photons"]["file_pattern"] = "some_file_pattern.fits" + with pytest.raises(galsim.GalSimConfigError, + match='off_detector_photons must contain one of either file_name or file_pattern in config'): + galsim.config.ProcessInput(config) + + # Now ensure that, if a file pattern is provided, the config includes the + # DETNAME and DETNUM fields that are used to read the range of files. + config = create_base_off_detector_config() + config["input"]["off_detector_photons"]["file_pattern"] = "some_file_pattern.fits" + with pytest.raises(galsim.GalSimConfigError, + match='file_pattern must contain both "DETNAME" and "DETNUM" placeholders to be replaced ' + 'with the detector name and number when reading off-detector photon files.'): + galsim.config.ProcessInput(config) + + +def test_load_photons_raises_error_with_no_files(): + # Photons are read lazily, so check that an error is raised when we try to access them + # but the files they are to be read from don't exist. + + # Firstly check the case when a single non-existent file_name is provided. + with pytest.raises(galsim.GalSimConfigError): + off_detector_photons_input = OffDetectorPhotons("LsstCamSim", "R22_S11", file_name="non_existent_file.fits") + photons = off_detector_photons_input.photons + + # The same but with a file_pattern that doesn't match any existing files. + with pytest.raises(galsim.GalSimConfigError): + off_detector_photons_input = OffDetectorPhotons("LsstCamSim", "R22_S11", file_pattern="non_existent_file_pattern_DETNAME_DETNUM.fits") + photons = off_detector_photons_input.photons + + +def test_draw_off_detector_photons(): + # Test that the LSST_FocalPlaneImageBuilder can draw photons loaded by the + # OffDetectorPhotonsLoader on top of the image. + xsize = 10 + ysize = 10 + camera = "LsstCamSim" + det_name = "R22_S11" + # Manually create some off-detector photons in the base config to be drawn + # on the image. We expect a diagonal line in the resulting image. + photons = galsim.PhotonArray.fromArrays(x=np.arange(1., 11.), + y=np.arange(1., 11.), + flux=np.ones(10), + ) + # Assume that an existing image as well as off_detector_photons from the + # first pass have been read from file and are at this point available in the + # base config. + base = { + "image": { + "xsize": xsize, + "ysize": ysize, + "det_name": det_name, + "nobjects": 1, + "type": "LSST_FocalPlaneImage", + }, + "current_image": galsim.Image(xsize, ysize), + '_input_objs': { + 'off_detector_photons': [OffDetectorPhotons(camera, det_name, photons=photons)], + } + } + expected = np.diag(np.ones(10)) + image = galsim.config.BuildImage(base) + np.testing.assert_allclose(image.array, expected, rtol=1.e-10) if __name__ == "__main__": From 3dee97b56b6856cbc1b135777cfe7274271af2cd Mon Sep 17 00:00:00 2001 From: William Lucas Date: Tue, 24 Mar 2026 14:55:17 +0000 Subject: [PATCH 33/36] New test for LSST_FocalPlaneImage. --- tests/test_full_focal_plane.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_full_focal_plane.py b/tests/test_full_focal_plane.py index 8d0b8e69..ca6d68aa 100644 --- a/tests/test_full_focal_plane.py +++ b/tests/test_full_focal_plane.py @@ -6,7 +6,7 @@ import pytest from imsim.camera import get_camera -from imsim.full_focal_plane import gather_out_of_bounds_photons +from imsim.full_focal_plane import gather_out_of_bounds_photons, OffDetectorPhotons from imsim.utils import focal_to_pixel DATA_DIR = Path(__file__).parent / 'data' From d1c740bb6eb6d0e3c3b1c89831a663a4a04a9721 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Tue, 24 Mar 2026 14:58:48 +0000 Subject: [PATCH 34/36] Rework off-detector photons input to lazily read photons and also allow them to be manually set if desired, with what are hopefully appropriate warnings along the way. --- imsim/full_focal_plane.py | 68 ++++++++++++++++++++++++++++++--------- 1 file changed, 52 insertions(+), 16 deletions(-) diff --git a/imsim/full_focal_plane.py b/imsim/full_focal_plane.py index e07d0fd8..ddbb0b0b 100644 --- a/imsim/full_focal_plane.py +++ b/imsim/full_focal_plane.py @@ -88,7 +88,7 @@ class OffDetectorPhotons(object): then, and will now be read in this second pass to be accumulated on other sensors. """ - def __init__(self, camera, det_name, file_name=None, file_pattern=None, dir=None, logger=None): + def __init__(self, camera, det_name, file_name=None, file_pattern=None, photons=None, dir=None, logger=None): """ Initialize the off-detector photons input class. @@ -106,18 +106,27 @@ def __init__(self, camera, det_name, file_name=None, file_pattern=None, dir=None logger: logging.logger A logger object. (default: None) """ - self.logger = logger + logger = galsim.config.LoggerWrapper(logger) camera = get_camera(camera) + self.det_name = det_name self.det = camera[det_name] - # The Loader should have ensured we have one but not both of file_name or file_pattern. + # Lazy read the off-detector photons when they're accessed for the first time. + self.photons = None + # The Loader should have ensured we have one but not both of file_name + # or file_pattern from the config, but here we also allow the photons to + # be set directly during initialisation in case when instantiating in code. + # So, make sure again that two of these are None. + if sum([file_name is None, file_pattern is None, photons is None]) != 2: + raise galsim.GalSimIncompatibleValuesError("OffDetectorPhotons must be initialized with exactly one of file_name, file_pattern, or photons.",values={"file_name": file_name, "file_pattern": file_pattern, "photons": photons}) + if file_name is not None: # Only have a single file to read. file_names = [file_name] elif file_pattern is not None: # We want to read up to 188 files: potentially photons from every # detector except this one (as those photons were drawn normally in - # the first pass). We may want to allow det_num_start and - # det_num_end to be set in the config in a similar way to output + # the first pass). We may want change this to allow det_num_start + # and det_num_end to be set in the config in a similar way to output # where it's controlled with det_num.first, nfiles and only_dets. det_num_start = 0 det_num_end = 189 @@ -125,21 +134,45 @@ def __init__(self, camera, det_name, file_name=None, file_pattern=None, dir=None for i in range(det_num_start, det_num_end) if camera[i].getName() != det_name ] + else: + # An OffDetectorPhotons object can be be directly initialized in + # code with a photon array instead of reading from file. + file_names = [] + self.photons = photons if dir is not None: file_names = [os.path.join(dir, fname) for fname in file_names] # Allow for the possibility that not all the 188 files exist. self.file_names = [fname for fname in file_names if os.path.isfile(fname)] if len(self.file_names) == 0: - raise galsim.GalSimConfigValueError( - f"Detector {det_name} did not find any expected off detector photon files: {file_names}.") - logger.info(f"Detector {det_name} reading {len(self.file_names)} files for off-detector photons.") - # Would it be better to make photons a @property and read lazily when - # first accessed? - photon_storage = [] - for file_name in self.file_names: - photon_storage.append(galsim.PhotonArray.read(file_name)) - self.photons = galsim.PhotonArray.concatenate(photon_storage) - self.photons.x, self.photons.y = focal_to_pixel(self.photons.x, self.photons.y, self.det) + logger.warning(f"Detector {det_name} did not find any expected off detector photon files: {file_names}. No photons can be read from file.") + + @property + def photons(self): + if self._photons is None: + self._photons = self._read_photons() + return self._photons + + @photons.setter + def photons(self, photon_array): + # Allow the photons to be set manually from outside the class if necessary. + self._photons = photon_array + + def _read_photons(self, logger=None): + logger = galsim.config.LoggerWrapper(logger) + if len(self.file_names) > 0: + photon_storage = [] + logger.info(f"Detector {self.det_name} reading {len(self.file_names)} files for off-detector photons.") + for file_name in self.file_names: + photon_storage.append(galsim.PhotonArray.read(file_name)) + photons = galsim.PhotonArray.concatenate(photon_storage) + photons.x, photons.y = focal_to_pixel(photons.x, photons.y, self.det) + else: + # If we've reached this point, we've initialised the class without + # providing the names of any usable files and have now tried to + # access the photons without setting them from outside. + # Raise an error. + raise galsim.GalSimConfigError(f"No readable files provided for OffDetectorPhotons for {self.det_name}; photons cannot be read.") + return photons class OffDetectorPhotonsLoader(InputLoader): @@ -153,7 +186,10 @@ def getKwargs(self, config, base, logger): kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) if ('file_name' in kwargs) == ('file_pattern' in kwargs): - raise galsim.GalSimConfigValueError('off_detector_photons must contain one of either file_name or file_pattern') + raise galsim.GalSimConfigError('off_detector_photons must contain one of either file_name or file_pattern in config') + if 'file_pattern' in kwargs and ("DETNAME" not in kwargs['file_pattern'] or "DETNUM" not in kwargs['file_pattern']): + raise galsim.GalSimConfigError('file_pattern must contain both "DETNAME" and "DETNUM" placeholders to be replaced ' + 'with the detector name and number when reading off-detector photon files.') kwargs['logger'] = logger # Base assumption (for now) is that we have a per-detector set of inputs From e5b5134ae7dcdfa0bfabecb139d2d9079b5f1ea1 Mon Sep 17 00:00:00 2001 From: William Lucas Date: Thu, 26 Mar 2026 10:40:17 +0000 Subject: [PATCH 35/36] Minor tidying up in the second pass example config. --- examples/imsim-user-instcat-ffp-pass2.yaml | 26 ++++++++-------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/examples/imsim-user-instcat-ffp-pass2.yaml b/examples/imsim-user-instcat-ffp-pass2.yaml index 89d153c0..6b136b39 100644 --- a/examples/imsim-user-instcat-ffp-pass2.yaml +++ b/examples/imsim-user-instcat-ffp-pass2.yaml @@ -10,12 +10,6 @@ template: imsim-config-instcat # Make your changes below. ################################################################ -# Put your own commands that override the defaults below here. For example -# input.instance_catalog.file_name: ./imsim_cat_197356.txt -# input.instance_catalog.sort_mag: False -# input.tree_rings.only_dets: [R22_S11] -# image.nobjects: 5 - input.checkpoint: "" # input.instance_catalog: "" input.instance_catalog.file_name: $os.environ.get('IMSIM_HOME')+'/imSim/examples/example_instance_catalog_ffp.txt' @@ -25,6 +19,8 @@ input.instance_catalog.sort_mag: False input.tree_rings: "" input.atm_psf: "" +# In the second pass, we add off-detector photons to the image generated in the +# first pass from the the on-detector photons. Give the image here. input.initial_image: dir: output file_name: @@ -37,17 +33,13 @@ input.initial_image: - $det_name - "@output.det_num" -# Simple off detector photon input reading one file. -# input.off_detector_photons: -# file_name: "output/off_det_00398414-0-r-R22_S11-det094.fits" -# camera: $camera -# det_name: $det_name - -# If we're doing the preprocessing step then we've already reorganised the -# photons so they live in one file per destination detector, so we only need to -# read that one. But, we might feasibly also wish to directly specify a list of -# off-detector photon files, which presumably haven't been preprocessed, so I -# might want to add the capability to work from a file list too. +# Two methods can be used to read files containing off-detector photons. If a +# single file only is to be read, then give the file_name. Alternatively, we can +# give a file_pattern containing the fields DETNAME and DETNUM, which are +# wildcards for the detector name and number in an expected set of files +# matching the pattern. The loader will then attempt to read from all files it +# finds matching this patter, excluding the file with the same name/num as the +# current detector. input.off_detector_photons: dir: output # file_name: From f8baa872da00c7ecfd81396f0b970b1f18ba535f Mon Sep 17 00:00:00 2001 From: William Lucas Date: Thu, 26 Mar 2026 11:08:32 +0000 Subject: [PATCH 36/36] Add files in tests/data for new full focal plane tests. --- .../off_det_00398414-0-r-R22_S11-det094.fits | Bin 0 -> 8640 bytes .../off_det_00398414-0-r-R22_S12-det095.fits | Bin 0 -> 8640 bytes ...ff_det_multi_00398414-0-r-R22_S10-det093.fits | Bin 0 -> 8640 bytes ...ff_det_multi_00398414-0-r-R22_S12-det095.fits | Bin 0 -> 8640 bytes 4 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 tests/data/off_det_00398414-0-r-R22_S11-det094.fits create mode 100644 tests/data/off_det_00398414-0-r-R22_S12-det095.fits create mode 100644 tests/data/off_det_multi_00398414-0-r-R22_S10-det093.fits create mode 100644 tests/data/off_det_multi_00398414-0-r-R22_S12-det095.fits diff --git a/tests/data/off_det_00398414-0-r-R22_S11-det094.fits b/tests/data/off_det_00398414-0-r-R22_S11-det094.fits new file mode 100644 index 0000000000000000000000000000000000000000..28abc95347f1d13be8be587b6da9588894f2e55f GIT binary patch literal 8640 zcmeH~O=uHA6vzGkya?56ox{SGUTQuLUc%DYO<74dsXIc1Tqf?OE+o5Qvw_BYJzB5z zCPlo63SN5aQ1Bvx7isY#h!-Isg7)srY+@#`D=7&D;l0D|?7X-C`R{M%@WwJ$%0(Ta z9Eu5y7SMg)tNTH-jY1z43~ZrxNIaVac1%5+${VF445)^h@PQ1!+ z)?BO@7JH2KcZ2D9oo3AmkY67tyG_SyyS~>Ribx~9f^~B#K@UrXV-AJbn6*#_1b_e# z00KY&2mk>f@c$Fw?WARt%p6klhKU#RMO~G?Y5~<;j|5#5l3K$-&UzR%r&7u2@L2oj z^wIlg(T}I%`{y+8?;DP{7CzwpxIO5^5vxxMp_h)=OY#0bPLI_mMfq0guF3lNG0qEX zf!}GN770kx37sJRys-Q)z0?Rlp%g6ez><;o!3Ne=7W?wO_gs4mvXJ=n$VmjQtk{@zKDgi zEfJoSJ4=<%eDbOGthX_{v$`+k&Qax)FrmF*%G>+DrQCU{e5}qL?mf|7vg>AbM#^2F z%6nFSa9Mk^%{^^NxrypU;uX@00e*l5C8%|00;m9An=#KA5KYy A-R`fP^EBQ09gsN$s$|${r+DJhvgM*!7DX}Blk&+-5B*cPL zF(D8OV($#wsbWN|u$GmDrOe!&6LOKxC{-dPr1zFP-+TW2-Oo3?GmZUPMF%K?A}rX1 zhk@S=!gd640A<56A&PC^u|r2u&z4GtRWs@-dbgOK9fo!v96Pol?sr_Z%6OG**A0_B zD*fGLdVaUva6<^1BW17c`jHp-(O5*9>2*sdyD55>R5;~Sh)r3GWk>)CAOR$R1dsp{ zKmz|i0p3oUMwJvnD;dPvE>(0*+SMjBJl_uc5ZjHG3vMq?np3G{a{5pE=iJ5nXVH&m z)%|mh_xCN=KZ+mme%u*#;=t&)@jN zQ_ie%^yvB!C2vz=a50VnHFwvnWuq@Z;Hy4+mR6`KpKo3cb5*@gpUR zye2Q@E>n2YdH*%DMak0OZ13c{lv||m^!RY&{1YY1r)=^|$}LfNezHC|GbmZf6q$mQ zTc+^hGm9^jtcveSxhoW2a(^ZzYp+?nk#Z{(UU98&lw_Ed^~X|fmBQPN^V8!8lw9TG pp_E&r@Q#%azEg6I@8!0X8!^BU2_OL^fCP{L5F445)^h@PQ1!+ z)?BO@7JH2KcZ2D9oo3AmkY67tyG_SyyS~>Ribx~9f^~B#K@UrXV-AJbn6*#_1b_e# z00KY&2mk>f@c$Fw?WARt%p6klhKU#RMO~G?Y5~<;j|5#5l3K$-&UzR%r&7u2@L2oj z^wIlg(T}I%`{y+8?;DP{7CzwpxIO5^5vxxMp_h)=OY#0bPLI_mMfq0guF3lNG0qEX zf!}GN770kx37sJRys-Q)z0?Rlp%g6ez><;o!3Ne=7W?wO_gs4mvXJ=n$VmjQtk{@zKDgi zEfJoSJ4=<%eDbOGthX_{v$`+k&Qax)FrmF*%G>+DrQCU{e5}qL?mf|7vg>AbM#^2F z%6nFSa9Mk^%{^^NxrypU;uX@00e*l5C8%|00;m9An=#KA5KYy A))6eN}lwu`02P8~;TU+CkFHX)Q4ij#HnFph#a z>GVh%lZwM#(o03;iNrJX!>zKz8%{IF?*{mZBoPDRiBynu(DltOfZJ8xlYQNB{{S0VIF~kih>> zK#r5HQ?F5Cm7N;jDpzew`>J))3L}vWNGe)wpZJGqHk@iLlgocPKbOuwKdX7XXr7-- za=vf-(SG_!&g0%_5+~4~SIVzs_AAQy-sA`UdDXa8zgL4fev0$LeiHY)q$?89@l!uB zpBHxi=2x8Hr~J$=Y<@DI`hTs%x%kO^nID|uRL!1mL+%%AFc&|WFY}vUzD<=cX!jMu*}hli z%lzh&56{2DU4s`SfCP{L5;zxu1t>