Skip to content

Commit e950232

Browse files
Higgins00Higgins00
committed
Separate classes so it can be merged with main.
1 parent 1440417 commit e950232

2 files changed

Lines changed: 671 additions & 0 deletions

File tree

‎roman_imsim/sca.py‎

Lines changed: 300 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,304 @@
88
from galsim.image import Image
99
import gc
1010

11+
class RomanSCAImageBuilder(ScatteredImageBuilder):
12+
13+
def setup(self, config, base, image_num, obj_num, ignore, logger):
14+
"""Do the initialization and setup for building the image.
15+
16+
This figures out the size that the image will be, but doesn't actually build it yet.
17+
18+
Parameters:
19+
config: The configuration dict for the image field.
20+
base: The base configuration dict.
21+
image_num: The current image number.
22+
obj_num: The first object number in the image.
23+
ignore: A list of parameters that are allowed to be in config that we can
24+
ignore here. i.e. it won't be an error if these parameters are present.
25+
logger: If given, a logger object to log progress.
26+
27+
Returns:
28+
xsize, ysize
29+
"""
30+
# import os, psutil
31+
# process = psutil.Process()
32+
# print('sca setup 1',process.memory_info().rss)
33+
logger.debug(
34+
"image %d: Building RomanSCA: image, obj = %d,%d",
35+
image_num,
36+
image_num,
37+
obj_num,
38+
)
39+
40+
self.nobjects = self.getNObj(config, base, image_num, logger=logger)
41+
logger.debug("image %d: nobj = %d", image_num, self.nobjects)
42+
43+
# These are allowed for Scattered, but we don't use them here.
44+
extra_ignore = [
45+
"image_pos",
46+
"world_pos",
47+
"stamp_size",
48+
"stamp_xsize",
49+
"stamp_ysize",
50+
"nobjects",
51+
]
52+
req = {"SCA": int, "filter": str, "mjd": float, "exptime": float}
53+
opt = {
54+
"draw_method": str,
55+
"use_fft_bright": bool,
56+
"stray_light": bool,
57+
"thermal_background": bool,
58+
"reciprocity_failure": bool,
59+
"dark_current": bool,
60+
"nonlinearity": bool,
61+
"ipc": bool,
62+
"read_noise": bool,
63+
"sky_subtract": bool,
64+
"ignore_noise": bool,
65+
}
66+
params = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore + extra_ignore)[0]
67+
68+
self.sca = params["SCA"]
69+
base["SCA"] = self.sca
70+
self.filter = params["filter"]
71+
self.mjd = params["mjd"]
72+
self.exptime = params["exptime"]
73+
74+
self.ignore_noise = params.get("ignore_noise", False)
75+
# self.exptime = params.get('exptime', roman.exptime) # Default is roman standard exposure time.
76+
self.stray_light = params.get("stray_light", False)
77+
self.thermal_background = params.get("thermal_background", False)
78+
self.reciprocity_failure = params.get("reciprocity_failure", False)
79+
self.dark_current = params.get("dark_current", False)
80+
self.nonlinearity = params.get("nonlinearity", False)
81+
self.ipc = params.get("ipc", False)
82+
self.read_noise = params.get("read_noise", False)
83+
self.sky_subtract = params.get("sky_subtract", False)
84+
85+
# If draw_method isn't in image field, it may be in stamp. Check.
86+
self.draw_method = params.get("draw_method", base.get("stamp", {}).get("draw_method", "phot"))
87+
88+
# pointing = CelestialCoord(ra=params['ra'], dec=params['dec'])
89+
# wcs = roman.getWCS(world_pos = pointing,
90+
# PA = params['pa']*galsim.degrees,
91+
# date = params['date'],
92+
# SCAs = self.sca,
93+
# PA_is_FPA = True
94+
# )[self.sca]
95+
96+
# # GalSim expects a wcs in the image field.
97+
# config['wcs'] = wcs
98+
99+
# If user hasn't overridden the bandpass to use, get the standard one.
100+
if "bandpass" not in config:
101+
base["bandpass"] = galsim.config.BuildBandpass(base["image"], "bandpass", base, logger=logger)
102+
103+
return roman.n_pix, roman.n_pix
104+
105+
# def getBandpass(self, filter_name):
106+
# if not hasattr(self, 'all_roman_bp'):
107+
# self.all_roman_bp = roman.getBandpasses()
108+
# return self.all_roman_bp[filter_name]
109+
110+
def buildImage(self, config, base, image_num, obj_num, logger):
111+
"""Build an Image containing multiple objects placed at arbitrary locations.
112+
113+
Parameters:
114+
config: The configuration dict for the image field.
115+
base: The base configuration dict.
116+
image_num: The current image number.
117+
obj_num: The first object number in the image.
118+
logger: If given, a logger object to log progress.
119+
120+
Returns:
121+
the final image and the current noise variance in the image as a tuple
122+
"""
123+
full_xsize = base["image_xsize"]
124+
full_ysize = base["image_ysize"]
125+
wcs = base["wcs"]
126+
127+
full_image = Image(full_xsize, full_ysize, dtype=float)
128+
full_image.setOrigin(base["image_origin"])
129+
full_image.wcs = wcs
130+
full_image.setZero()
131+
132+
full_image.header = galsim.FitsHeader()
133+
full_image.header["EXPTIME"] = self.exptime
134+
full_image.header["MJD-OBS"] = self.mjd
135+
full_image.header["DATE-OBS"] = Time(self.mjd, format="mjd").datetime.isoformat()
136+
full_image.header["FILTER"] = self.filter
137+
full_image.header["ZPTMAG"] = 2.5 * np.log10(self.exptime * roman.collecting_area)
138+
139+
base["current_image"] = full_image
140+
141+
if "image_pos" in config and "world_pos" in config:
142+
raise galsim.GalSimConfigValueError(
143+
"Both image_pos and world_pos specified for Scattered image.",
144+
(config["image_pos"], config["world_pos"]),
145+
)
146+
147+
if "image_pos" not in config and "world_pos" not in config:
148+
xmin = base["image_origin"].x
149+
xmax = xmin + full_xsize - 1
150+
ymin = base["image_origin"].y
151+
ymax = ymin + full_ysize - 1
152+
config["image_pos"] = {
153+
"type": "XY",
154+
"x": {"type": "Random", "min": xmin, "max": xmax},
155+
"y": {"type": "Random", "min": ymin, "max": ymax},
156+
}
157+
158+
nbatch = self.nobjects // 1000 + 1
159+
for batch in range(nbatch):
160+
start_obj_num = self.nobjects * batch // nbatch
161+
end_obj_num = self.nobjects * (batch + 1) // nbatch
162+
nobj_batch = end_obj_num - start_obj_num
163+
if nbatch > 1:
164+
logger.warning(
165+
"Start batch %d/%d with %d objects [%d, %d)",
166+
batch + 1,
167+
nbatch,
168+
nobj_batch,
169+
start_obj_num,
170+
end_obj_num,
171+
)
172+
stamps, current_vars = galsim.config.BuildStamps(
173+
nobj_batch, base, logger=logger, obj_num=start_obj_num, do_noise=False
174+
)
175+
base["index_key"] = "image_num"
176+
177+
for k in range(nobj_batch):
178+
# This is our signal that the object was skipped.
179+
if stamps[k] is None:
180+
continue
181+
bounds = stamps[k].bounds & full_image.bounds
182+
if not bounds.isDefined(): # pragma: no cover
183+
# These noramlly show up as stamp==None, but technically it is possible
184+
# to get a stamp that is off the main image, so check for that here to
185+
# avoid an error. But this isn't covered in the imsim test suite.
186+
continue
187+
188+
logger.debug("image %d: full bounds = %s", image_num, str(full_image.bounds))
189+
logger.debug(
190+
"image %d: stamp %d bounds = %s",
191+
image_num,
192+
k + start_obj_num,
193+
str(stamps[k].bounds),
194+
)
195+
logger.debug("image %d: Overlap = %s", image_num, str(bounds))
196+
full_image[bounds] += stamps[k][bounds]
197+
stamps = None
198+
199+
# # Bring the image so far up to a flat noise variance
200+
# current_var = FlattenNoiseVariance(
201+
# base, full_image, stamps, current_vars, logger)
202+
203+
return full_image, None
204+
205+
def addNoise(self, image, config, base, image_num, obj_num, current_var, logger):
206+
"""Add the final noise to a Scattered image
207+
208+
Parameters:
209+
image: The image onto which to add the noise.
210+
config: The configuration dict for the image field.
211+
base: The base configuration dict.
212+
image_num: The current image number.
213+
obj_num: The first object number in the image.
214+
current_var: The current noise variance in each postage stamps.
215+
logger: If given, a logger object to log progress.
216+
"""
217+
# check ignore noise
218+
if self.ignore_noise:
219+
return
220+
221+
base["current_noise_image"] = base["current_image"]
222+
wcs = base["wcs"]
223+
bp = base["bandpass"]
224+
rng = galsim.config.GetRNG(config, base)
225+
logger.info("image %d: Start RomanSCA detector effects", base.get("image_num", 0))
226+
227+
# Things that will eventually be subtracted (if sky_subtract) will have their expectation
228+
# value added to sky_image. So technically, this includes things that aren't just sky.
229+
# E.g. includes dark_current and thermal backgrounds.
230+
sky_image = image.copy()
231+
sky_level = roman.getSkyLevel(bp, world_pos=wcs.toWorld(image.true_center))
232+
logger.debug("Adding sky_level = %s", sky_level)
233+
if self.stray_light:
234+
logger.debug("Stray light fraction = %s", roman.stray_light_fraction)
235+
sky_level *= 1.0 + roman.stray_light_fraction
236+
wcs.makeSkyImage(sky_image, sky_level)
237+
238+
# The other background is the expected thermal backgrounds in this band.
239+
# These are provided in e-/pix/s, so we have to multiply by the exposure time.
240+
if self.thermal_background:
241+
tb = roman.thermal_backgrounds[self.filter] * self.exptime
242+
logger.debug("Adding thermal background: %s", tb)
243+
sky_image += roman.thermal_backgrounds[self.filter] * self.exptime
244+
245+
# The image up to here is an expectation value.
246+
# Realize it as an integer number of photons.
247+
poisson_noise = galsim.noise.PoissonNoise(rng)
248+
if self.draw_method == "phot":
249+
logger.debug("Adding poisson noise to sky photons")
250+
sky_image1 = sky_image.copy()
251+
sky_image1.addNoise(poisson_noise)
252+
image.quantize() # In case any profiles used InterpolatedImage, in which case
253+
# the image won't necessarily be integers.
254+
image += sky_image1
255+
else:
256+
logger.debug("Adding poisson noise")
257+
image += sky_image
258+
image.addNoise(poisson_noise)
259+
260+
# Apply the detector effects here. Not all of these are "noise" per se, but they
261+
# happen interspersed with various noise effects, so apply them all in this step.
262+
263+
# Note: according to Gregory Mosby & Bernard J. Rauscher, the following effects all
264+
# happen "simultaneously" in the photo diodes: dark current, persistence,
265+
# reciprocity failure (aka CRNL), burn in, and nonlinearity (aka CNL).
266+
# Right now, we just do them in some order, but this could potentially be improved.
267+
# The order we chose is historical, matching previous recommendations, but Mosby and
268+
# Rauscher don't seem to think those recommendations are well-motivated.
269+
270+
# TODO: Add burn-in and persistence here.
271+
272+
if self.reciprocity_failure:
273+
logger.debug("Applying reciprocity failure")
274+
roman.addReciprocityFailure(image)
275+
276+
if self.dark_current:
277+
dc = roman.dark_current * self.exptime
278+
logger.debug("Adding dark current: %s", dc)
279+
sky_image += dc
280+
dark_noise = galsim.noise.DeviateNoise(galsim.random.PoissonDeviate(rng, dc))
281+
image.addNoise(dark_noise)
282+
283+
if self.nonlinearity:
284+
logger.debug("Applying classical nonlinearity")
285+
roman.applyNonlinearity(image)
286+
287+
# Mosby and Rauscher say there are two read noises. One happens before IPC, the other
288+
# one after.
289+
# TODO: Add read_noise1
290+
if self.ipc:
291+
logger.debug("Applying IPC")
292+
roman.applyIPC(image)
293+
294+
if self.read_noise:
295+
logger.debug("Adding read noise %s", roman.read_noise)
296+
image.addNoise(galsim.GaussianNoise(rng, sigma=roman.read_noise))
297+
298+
logger.debug("Applying gain %s", roman.gain)
299+
image /= roman.gain
300+
301+
# Make integer ADU now.
302+
image.quantize()
303+
304+
if self.sky_subtract:
305+
logger.debug("Subtracting sky image")
306+
sky_image /= roman.gain
307+
sky_image.quantize()
308+
image -= sky_image
11309

12310
class RomanSCAImageBuilderCMOS(ScatteredImageBuilder):
13311

@@ -358,3 +656,5 @@ def addNoise(self, image, config, base, image_num, obj_num, current_var, logger)
358656

359657
# Register this as a valid type
360658
RegisterImageType("roman_sca_cmos", RomanSCAImageBuilderCMOS())
659+
# Register this as a valid type
660+
RegisterImageType("roman_sca", RomanSCAImageBuilder())

0 commit comments

Comments
 (0)