Loading CpicImgSim/camera.py 0 → 100644 +633 −0 Original line number Diff line number Diff line import os import yaml import numpy as np import scipy.ndimage as nd from astropy.io import fits from .config import cpism_refdata, solar_spectrum, MAG_SYSTEM from .utils import region_replace, random_seed_select from .io import log from .optics import filter_throughput def sky_frame_maker(band, skybg, platescale, shape): """ generate a sky background frame. Parameters ---------- band : str The band of the sky background. skybg : str The sky background file name. platescale : float The platescale of the camera in arcsec/pixel. shape : tuple The shape of the output frame. (y, x) Returns ------- sky_bkg_frame : numpy.ndarray The sky background frame. """ filter = filter_throughput(band) sk_spec = solar_spectrum.renorm(skybg, MAG_SYSTEM, filter) sky_bkg_frame = np.zeros(shape) sky_bkg_frame += (sk_spec * filter).integrate() * platescale**2 return sky_bkg_frame class CRobj(object): """ Cosmic ray object. Attributes ---------- flux : float The flux of the cosmic ray. angle : float The angle of the cosmic ray. sigma : float The width of the cosmic ray. length : int The length of the cosmic ray. """ def __init__(self, flux, angle, sigma, length) -> None: self.flux = flux self.angle = angle self.sigma = sigma self.length = length class CosmicRayFrameMaker(object): """ Cosmic ray frame maker. Parameters ---------- depth : float The depth of the camera pixel in um. pitch : float The pitch of the camera pixel in um. cr_rate : float The cosmic ray rate per second per cm2. """ def __init__(self) -> None: self.tmp_size = [7, 101] self.freq_power = -0.9 self.trail_std = 0.1 self.depth = 10 # um self.pitch = 13 # um self.cr_rate = 1 # particle per s per cm2 from Miles et al. 2021 def make_CR(self, length, sigma, seed=-1): """ make a image of cosmic ray with given length and sigma. Parameters ---------- length : int The length of the cosmic ray in pixel. sigma : float The width of the cosmic ray in pixel. Returns ------- output : numpy.ndarray The image of cosmic ray. """ h = self.tmp_size[0] w = self.tmp_size[1] freq = ((w-1)/2-np.abs(np.arange(w)-(w-1)/2)+1)**(self.freq_power) x = np.arange(w) - (w-1)/2 hl = (length-1)/2 x_wing = np.exp(-(np.abs(x)-hl)**2/sigma/sigma/2) x_wing[np.abs(x) < hl] = 1 cr = np.zeros([h, w]) center = (h-1)/2 for i in range(h): phase = np.random.rand(w)*2*np.pi trail0 = abs(np.fft.fft(freq*np.sin(phase) + 1j*x*np.cos(phase))) # TODO maybe somthing wrong trail_norm = (trail0 - trail0.mean())/trail0.std() cr[i, :] = np.exp(-(i - center)**2/sigma/sigma/2) \ * (trail_norm * self.trail_std + 1) * x_wing output = np.zeros([w, w]) d = (w-h)//2 output[d:d+h, :] = cr return output def _length_rand(self, N, seed=-1): """ randomly generate N cosmic ray length. """ len_out = [] seed = random_seed_select(seed=seed) log.debug(f"cr length seed: {seed}") for i in range(N): x2y2 = 2 while x2y2 > 1: lx = 1 - 2 * np.random.rand() ly = 1 - 2 * np.random.rand() x2y2 = lx * lx + ly * ly z = 1 - 2 * x2y2 r = 2 * np.sqrt(x2y2 * (1 - x2y2)) length = abs(r / z * self.depth) pitch = self.pitch len_out.append(int(length/pitch)) return np.array(len_out) def _number_rand(self, expt, pixsize, random=False, seed=-1): """ randomly generate the number of cosmic rays. """ area = (self.pitch / 1e4)**2 * pixsize[0] * pixsize[1] ncr = area * expt * self.cr_rate if random: seed = random_seed_select(seed=seed) log.debug(f"cr count: {seed}") ncr = np.random.poisson(ncr) else: ncr = int(ncr) self.ncr = ncr return ncr def _sigma_rand(self, N, seed=-1): """ randomly generate N cosmic ray sigma. """ sig_sig = 0.5 # asuming the sigma of size of cosmic ray is 0.5 seed = random_seed_select(seed=seed) log.debug(f"cr width seed: {seed}") sig = abs(np.random.randn(N))*sig_sig + 1/np.sqrt(12) * 1.2 # assume sigma is 1.2 times of pictch sig return sig def _flux_rand(self, N, seed=-1): """ randomly generate N cosmic ray flux. """ seed = random_seed_select(seed=seed) log.debug(f"cr flux seed: {seed}") u = np.random.rand(N) S0 = 800 lam = 0.57 S = (-np.log(1-u)/lam + S0**0.25)**4 return S def random_CR_parameter(self, expt, pixsize): """ randomly generate cosmic ray parameters, including number, length, flux, sigma and angle. Parameters ---------- expt : float The exposure time in second. pixsize : list The size of the image in pixel. Returns ------- CRs : list A list of cosmic ray objects. """ N = self._number_rand(expt, pixsize) log.debug(f"cr count: {N}") length = self._length_rand(N) if N > 0: log.debug(f"cr length, max: {length.max()}, min: {length.min()}") flux = self._flux_rand(N) log.debug(f"cr flux, max: {flux.max()}, min: {flux.min()}") sig = self._sigma_rand(N) log.debug(f"cr width, max: {sig.max()}, min: {sig.min()}") seed = random_seed_select(seed=-1) log.debug(f"cr angle seed: {seed}") angle = np.random.rand(N) * 180 CRs = [] for i in range(N): CRs.append(CRobj(flux[i], angle[i], sig[i], length[i])) return CRs def make_cr_frame(self, shape, expt, seed=-1): """ make a cosmic ray frame. Parameters ---------- shape : list The size of the image in pixel. expt : float The exposure time in second. seed : int, optional The random seed. The default is -1. If seed is -1, the seed will be randomly selected. Returns ------- image : numpy.ndarray The cosmic ray frame. """ image = np.zeros(shape) sz = shape cr_array = self.random_CR_parameter(expt, shape) cr_center = (self.tmp_size[1] - 1)/2 seed = random_seed_select(seed=seed) log.debug(f"cr position seed: {seed}") for i in range(len(cr_array)): cr = cr_array[i] x = np.random.rand() * sz[1] y = np.random.rand() * sz[0] cr_img = self.make_CR(cr.length, cr.sigma) cr_img *= cr.flux cr_img = abs(nd.rotate(cr_img, cr.angle, reshape=False)) if i == 0: pdin = False else: pdin = True if i == len(cr_array) - 1: pdout = False else: pdout = True image = region_replace( image, cr_img, [y-cr_center, x-cr_center], padded_in=pdin, padded_out=pdout, subpix=True ) image = np.maximum(image, 0) log.debug(f"cr image max: {image.max()}, min: {image.min()}") return image class EMCCD(object): """ EMCCD camera class Parameters ---------- config_file : str config file name Attributes ---------- switch : dict switch for each camera effects, including: - 'flat': bool, - 'dark': bool, - 'stripe': bool, - 'cic': bool, - 'cte': bool, - 'badcolumn': bool, - 'nonlinear': bool, - 'cosmicray': bool, - 'blooming': bool, """ def __init__(self, config_file="emccd_config.yaml"): self.plszx = 1024 self.plszy = 1024 self.pscan1 = 8 self.pscan2 = 0 self.oscan1 = 16 self.oscan2 = 18 self.udark = 6 self.bdark = 2 self.ldark = 16 self.rdark = 16 self.fullwell = 80_000 self.em_fullwell = 780_000 # if config file exists, load it, otherwise use default values config_file = cpism_refdata + '/camera/' + config_file log.debug(f"emccd config file: {config_file}") if os.path.exists(config_file): self.load_config(config_file) else: # pragma: no cover # set default values for EMCCD # note: these values are default values, you can change them by load_config() # ↓↓↓↓↓↓↓start default values setting↓↓↓↓↓↓ self.readout_noise = 40 self.ph_per_adu = 8 self.bias_level = 30 self.max_adu = 16_383 self.switch = { 'flat': True, 'dark': True, 'stripe': True, 'cic': False, 'cte': False, 'badcolumn': True, 'nonlinear': False, 'cosmicray': True, 'blooming': False, } self.dark_file = cpism_refdata + '/camera/emccd_dark_current.fits' self.flat_file = cpism_refdata + '/camera/emccd_flat_field.fits' self.cic_file = cpism_refdata + '/camera/emcid_cic.fits' self.bad_col_file = cpism_refdata + '/camera/emccd_bad_columns.fits' # ↑↑↑↑↑↑↑end default values setting↑↑↑↑↑↑ # note: these values are default values, you can change them by load_config() self.flat_shape = [self.plszy, self.plszx] darksz_x = self.plszx + self.rdark + self.ldark darksz_y = self.plszy + self.udark + self.bdark self.dark_shape = [darksz_y, darksz_x] biassz_x = darksz_x + self.pscan1 + self.oscan1 biassz_y = darksz_y + self.pscan2 + self.oscan2 self.image_shape = [biassz_y, biassz_x] self.flat = fits.getdata(self.flat_file) self.cic = fits.getdata(self.cic_file) self.dark = fits.getdata(self.dark_file) self.bad_col = fits.getdata(self.bad_col_file) def load_config(self, config_file): """ load config file. Only for internal use. """ with open(config_file, 'r') as f: config = yaml.load(f, Loader=yaml.FullLoader) self.switch = config['switch'] self.readout_noise = config['readout_noise'] self.ph_per_adu = config['ph_per_adu'] self.bias_level = config['bias_level'] self.max_adu = config['max_adu'] self.dark_file = cpism_refdata + "/camera/" + config['dark_file'] self.flat_file = cpism_refdata + "/camera/" + config['flat_file'] self.cic_file = cpism_refdata + "/camera/" + config['cic_file'] self.bad_col_file = cpism_refdata + \ "/camera/" + config['bad_col_file'] def vertical_blooming(self, image): """ vertical blooming effect """ fullwell = self.fullwell line = np.arange(image.shape[0]) yp, xp = np.where(image > fullwell) n_saturated = len(xp) log.debug(f"{len(xp)} pixels are saturated!") if n_saturated > 5000: log.warning(f"More than 5000({len(xp)}) pixels are saturated!") img0 = image.copy() for x, y in zip(xp, yp): image[:, x] += np.exp(-(line-y)**2/20**2) * img0[y, x] * 0.2 return np.minimum(image, fullwell) def nonlinear_effect(self, image): """ nonlinear effect """ fullwell = self.fullwell nonlinear_coefficient = 0.1 log.debug( f"nonlinear effect added with coefficient {nonlinear_coefficient}") image += (image / fullwell)**2 * nonlinear_coefficient * fullwell return image def emregester_blooming(self, image, max_iteration=5): """ emregester blooming effect """ line = image.flatten().copy() curve_x = np.arange(1300)+2 curve_y = np.exp(11*curve_x**(-0.19)-11) curve_y[0] = 0 curve_y /= curve_y.sum() over_limit_coe = 0.999 saturated = image > self.em_fullwell n_saturated = saturated.sum() if n_saturated > 0: log.debug(f"{n_saturated} pixels are saturated during EM process.") if n_saturated > 2000: log.warning( f"More than 2000 ({n_saturated}) pixels are saturated during EM process!") for index in range(max_iteration): over_limit = np.maximum( line - self.em_fullwell * over_limit_coe, 0) line = np.minimum(line, self.em_fullwell * over_limit_coe) blooming = np.convolve(over_limit, curve_y, mode='full')[ :len(line)] line = line + blooming n_over = (line > self.em_fullwell).sum() if n_over <= 0: break log.debug( f'{index}/{max_iteration} loop: saturated pixel number: {n_over}') return line.reshape(image.shape) def cte(self, image): """ cte effect """ image = self.emregester_blooming(image, max_iteration=5) return image def readout(self, image_focal, emgain, expt, image_cosmic_ray=None): """ emccd readout. Generate a image with emccd readout effect. Parameters ---------- image_focal : numpy.ndarray image at focal plane. Unit: electron/second emgain : float emgain of emccd expt : float exposure time. Unit: second image_cosmic_ray : numpy.ndarray, optional cosmic ray image. Unit: electron/second, by default None Returns ------- numpy.ndarray image with emccd readout effect. Unit: ADU Notes ----- 1. effects include: dark, flat, cte, blooming, nonlinear, etc. Can be turned on/off by switch. 2. size of input image_focal must be 1024x1024 3. size of output image is 1080x1056 (including overscan and dark reference region) 4. Q.E is not included in this function. It should be included in image_focal. See optics.py for details. """ log.debug( fr"EMCCD readout: {emgain=:}, {expt=:}, image_comic_ray:{'None' if image_cosmic_ray is None else 'Not None'}") log.debug( f"camera effects switch={self.switch}" ) image = image_focal * expt if self.switch['flat']: image = image * self.flat if self.switch['nonlinear']: image = self.nonlinear_effect(image) darksz_x = self.plszx + self.rdark + self.ldark darksz_y = self.plszy + self.udark + self.bdark img_dark = np.zeros((darksz_y, darksz_x)) img_dark[ self.bdark:self.plszy+self.bdark, self.ldark:self.ldark+self.plszx ] = image image = img_dark if self.switch['dark']: image += self.dark * expt if self.switch['cic']: image += self.cic if image_cosmic_ray is not None and self.switch['cosmicray']: image += image_cosmic_ray if self.switch['blooming']: image = self.vertical_blooming(image) if self.switch['badcolumn']: for i in range(self.bad_col.shape[1]): deadpix_x = self.bad_col[0, i] deadpix_y = self.bad_col[1, i] image[deadpix_y:, deadpix_x] = 0 biassz_x = darksz_x + self.pscan1 + self.oscan1 biassz_y = darksz_y + self.pscan2 + self.oscan2 img_bias = np.zeros((biassz_y, biassz_x), dtype=int) seed = random_seed_select() log.debug(f"photon noise seed: {seed}") img_bias[ self.pscan2:self.pscan2+darksz_y, self.pscan1:self.pscan1+darksz_x ] = np.random.poisson(image) image = img_bias if self.switch['cte']: image = self.cte(image * emgain) / emgain seed = random_seed_select() log.debug(f"gamma noise seed: {seed}") if emgain != 1: image = np.random.gamma(image, emgain) image = np.minimum(image, self.em_fullwell) seed = random_seed_select() log.debug(f"readout noise seed: {seed}") image += np.random.randn(biassz_y, biassz_x) * self.readout_noise image = image / self.ph_per_adu + self.bias_level if self.switch['stripe']: image += self.add_stripe_effect(image) image = np.minimum(image, self.max_adu) image = np.maximum(image, 0) return image.astype(np.uint16) def add_stripe_effect(self, image): """ add stripe effect """ shape = image.shape v_width = 1 v_amplitude = 30 v_limit = 0.01 v_base = 10 h_width = 20 h_amplitude = 10 h_limit = 0.9 h_base = 20 index = np.linspace(0, np.pi, shape[0] * shape[1]) def stripe(width, limit, amplitude, base, axis=0): seed = random_seed_select() log.debug(f"stripe noise seed: {seed}") dim_axis = shape[axis] dim_other = shape[0] * shape[1] // shape[axis] value = np.sin(index / width * dim_axis + np.pi * dim_axis / width * np.random.randint(1024)) value = np.maximum(value, -limit) value = np.minimum(value, limit) value = (value / limit + limit) / 2 * amplitude + base value = value.reshape(dim_axis, dim_other) if axis == 1: value = value.T return value output = stripe(v_width, v_limit, v_amplitude, v_base, axis=1) output += stripe(h_width, h_limit, h_amplitude, h_base, axis=0) return output # # plt.plot(horizontal_index, horizontal_value) # # # plt.xlim([0, 6.28]) # # plt.show() # fits.writeto('horizontal_value.fits', output, overwrite=True) # if __name__ == '__main__': # import matplotlib.pyplot as plt # emccd = EMCCD() # image_focal = np.zeros((emccd.plszy, emccd.plszx)) + 1000 # image_focal[100:105, 100:105] = 10_000_000 # after_cte = emccd.emregester_blooming(image_focal, max_iteration=100) # print(after_cte.sum(), image_focal.sum()) # fits.writeto('after_cte.fits', after_cte, overwrite=True) # # darksz_x = emccd.plszx + emccd.rdark + emccd.ldark # # darksz_y = emccd.plszy + emccd.udark + emccd.bdark # # iamge_cosmic_ray = np.zeros((darksz_y, darksz_x)) # # emgain = 10 # # expt = 10 # # image = emccd.readout(image_focal, emgain, expt, iamge_cosmic_ray) # # fits.writeto('test.fits', image, overwrite=True) # image = np.zeros((1000, 1000)) # make_cosmic_ray_frame = CosmicRayFrameMaker() # crimage = make_cosmic_ray_frame(image.shape, 3000) # fits.writeto('crimage.fits', crimage, overwrite=True) # # emccd.add_stripe_effect(image) Loading
CpicImgSim/camera.py 0 → 100644 +633 −0 Original line number Diff line number Diff line import os import yaml import numpy as np import scipy.ndimage as nd from astropy.io import fits from .config import cpism_refdata, solar_spectrum, MAG_SYSTEM from .utils import region_replace, random_seed_select from .io import log from .optics import filter_throughput def sky_frame_maker(band, skybg, platescale, shape): """ generate a sky background frame. Parameters ---------- band : str The band of the sky background. skybg : str The sky background file name. platescale : float The platescale of the camera in arcsec/pixel. shape : tuple The shape of the output frame. (y, x) Returns ------- sky_bkg_frame : numpy.ndarray The sky background frame. """ filter = filter_throughput(band) sk_spec = solar_spectrum.renorm(skybg, MAG_SYSTEM, filter) sky_bkg_frame = np.zeros(shape) sky_bkg_frame += (sk_spec * filter).integrate() * platescale**2 return sky_bkg_frame class CRobj(object): """ Cosmic ray object. Attributes ---------- flux : float The flux of the cosmic ray. angle : float The angle of the cosmic ray. sigma : float The width of the cosmic ray. length : int The length of the cosmic ray. """ def __init__(self, flux, angle, sigma, length) -> None: self.flux = flux self.angle = angle self.sigma = sigma self.length = length class CosmicRayFrameMaker(object): """ Cosmic ray frame maker. Parameters ---------- depth : float The depth of the camera pixel in um. pitch : float The pitch of the camera pixel in um. cr_rate : float The cosmic ray rate per second per cm2. """ def __init__(self) -> None: self.tmp_size = [7, 101] self.freq_power = -0.9 self.trail_std = 0.1 self.depth = 10 # um self.pitch = 13 # um self.cr_rate = 1 # particle per s per cm2 from Miles et al. 2021 def make_CR(self, length, sigma, seed=-1): """ make a image of cosmic ray with given length and sigma. Parameters ---------- length : int The length of the cosmic ray in pixel. sigma : float The width of the cosmic ray in pixel. Returns ------- output : numpy.ndarray The image of cosmic ray. """ h = self.tmp_size[0] w = self.tmp_size[1] freq = ((w-1)/2-np.abs(np.arange(w)-(w-1)/2)+1)**(self.freq_power) x = np.arange(w) - (w-1)/2 hl = (length-1)/2 x_wing = np.exp(-(np.abs(x)-hl)**2/sigma/sigma/2) x_wing[np.abs(x) < hl] = 1 cr = np.zeros([h, w]) center = (h-1)/2 for i in range(h): phase = np.random.rand(w)*2*np.pi trail0 = abs(np.fft.fft(freq*np.sin(phase) + 1j*x*np.cos(phase))) # TODO maybe somthing wrong trail_norm = (trail0 - trail0.mean())/trail0.std() cr[i, :] = np.exp(-(i - center)**2/sigma/sigma/2) \ * (trail_norm * self.trail_std + 1) * x_wing output = np.zeros([w, w]) d = (w-h)//2 output[d:d+h, :] = cr return output def _length_rand(self, N, seed=-1): """ randomly generate N cosmic ray length. """ len_out = [] seed = random_seed_select(seed=seed) log.debug(f"cr length seed: {seed}") for i in range(N): x2y2 = 2 while x2y2 > 1: lx = 1 - 2 * np.random.rand() ly = 1 - 2 * np.random.rand() x2y2 = lx * lx + ly * ly z = 1 - 2 * x2y2 r = 2 * np.sqrt(x2y2 * (1 - x2y2)) length = abs(r / z * self.depth) pitch = self.pitch len_out.append(int(length/pitch)) return np.array(len_out) def _number_rand(self, expt, pixsize, random=False, seed=-1): """ randomly generate the number of cosmic rays. """ area = (self.pitch / 1e4)**2 * pixsize[0] * pixsize[1] ncr = area * expt * self.cr_rate if random: seed = random_seed_select(seed=seed) log.debug(f"cr count: {seed}") ncr = np.random.poisson(ncr) else: ncr = int(ncr) self.ncr = ncr return ncr def _sigma_rand(self, N, seed=-1): """ randomly generate N cosmic ray sigma. """ sig_sig = 0.5 # asuming the sigma of size of cosmic ray is 0.5 seed = random_seed_select(seed=seed) log.debug(f"cr width seed: {seed}") sig = abs(np.random.randn(N))*sig_sig + 1/np.sqrt(12) * 1.2 # assume sigma is 1.2 times of pictch sig return sig def _flux_rand(self, N, seed=-1): """ randomly generate N cosmic ray flux. """ seed = random_seed_select(seed=seed) log.debug(f"cr flux seed: {seed}") u = np.random.rand(N) S0 = 800 lam = 0.57 S = (-np.log(1-u)/lam + S0**0.25)**4 return S def random_CR_parameter(self, expt, pixsize): """ randomly generate cosmic ray parameters, including number, length, flux, sigma and angle. Parameters ---------- expt : float The exposure time in second. pixsize : list The size of the image in pixel. Returns ------- CRs : list A list of cosmic ray objects. """ N = self._number_rand(expt, pixsize) log.debug(f"cr count: {N}") length = self._length_rand(N) if N > 0: log.debug(f"cr length, max: {length.max()}, min: {length.min()}") flux = self._flux_rand(N) log.debug(f"cr flux, max: {flux.max()}, min: {flux.min()}") sig = self._sigma_rand(N) log.debug(f"cr width, max: {sig.max()}, min: {sig.min()}") seed = random_seed_select(seed=-1) log.debug(f"cr angle seed: {seed}") angle = np.random.rand(N) * 180 CRs = [] for i in range(N): CRs.append(CRobj(flux[i], angle[i], sig[i], length[i])) return CRs def make_cr_frame(self, shape, expt, seed=-1): """ make a cosmic ray frame. Parameters ---------- shape : list The size of the image in pixel. expt : float The exposure time in second. seed : int, optional The random seed. The default is -1. If seed is -1, the seed will be randomly selected. Returns ------- image : numpy.ndarray The cosmic ray frame. """ image = np.zeros(shape) sz = shape cr_array = self.random_CR_parameter(expt, shape) cr_center = (self.tmp_size[1] - 1)/2 seed = random_seed_select(seed=seed) log.debug(f"cr position seed: {seed}") for i in range(len(cr_array)): cr = cr_array[i] x = np.random.rand() * sz[1] y = np.random.rand() * sz[0] cr_img = self.make_CR(cr.length, cr.sigma) cr_img *= cr.flux cr_img = abs(nd.rotate(cr_img, cr.angle, reshape=False)) if i == 0: pdin = False else: pdin = True if i == len(cr_array) - 1: pdout = False else: pdout = True image = region_replace( image, cr_img, [y-cr_center, x-cr_center], padded_in=pdin, padded_out=pdout, subpix=True ) image = np.maximum(image, 0) log.debug(f"cr image max: {image.max()}, min: {image.min()}") return image class EMCCD(object): """ EMCCD camera class Parameters ---------- config_file : str config file name Attributes ---------- switch : dict switch for each camera effects, including: - 'flat': bool, - 'dark': bool, - 'stripe': bool, - 'cic': bool, - 'cte': bool, - 'badcolumn': bool, - 'nonlinear': bool, - 'cosmicray': bool, - 'blooming': bool, """ def __init__(self, config_file="emccd_config.yaml"): self.plszx = 1024 self.plszy = 1024 self.pscan1 = 8 self.pscan2 = 0 self.oscan1 = 16 self.oscan2 = 18 self.udark = 6 self.bdark = 2 self.ldark = 16 self.rdark = 16 self.fullwell = 80_000 self.em_fullwell = 780_000 # if config file exists, load it, otherwise use default values config_file = cpism_refdata + '/camera/' + config_file log.debug(f"emccd config file: {config_file}") if os.path.exists(config_file): self.load_config(config_file) else: # pragma: no cover # set default values for EMCCD # note: these values are default values, you can change them by load_config() # ↓↓↓↓↓↓↓start default values setting↓↓↓↓↓↓ self.readout_noise = 40 self.ph_per_adu = 8 self.bias_level = 30 self.max_adu = 16_383 self.switch = { 'flat': True, 'dark': True, 'stripe': True, 'cic': False, 'cte': False, 'badcolumn': True, 'nonlinear': False, 'cosmicray': True, 'blooming': False, } self.dark_file = cpism_refdata + '/camera/emccd_dark_current.fits' self.flat_file = cpism_refdata + '/camera/emccd_flat_field.fits' self.cic_file = cpism_refdata + '/camera/emcid_cic.fits' self.bad_col_file = cpism_refdata + '/camera/emccd_bad_columns.fits' # ↑↑↑↑↑↑↑end default values setting↑↑↑↑↑↑ # note: these values are default values, you can change them by load_config() self.flat_shape = [self.plszy, self.plszx] darksz_x = self.plszx + self.rdark + self.ldark darksz_y = self.plszy + self.udark + self.bdark self.dark_shape = [darksz_y, darksz_x] biassz_x = darksz_x + self.pscan1 + self.oscan1 biassz_y = darksz_y + self.pscan2 + self.oscan2 self.image_shape = [biassz_y, biassz_x] self.flat = fits.getdata(self.flat_file) self.cic = fits.getdata(self.cic_file) self.dark = fits.getdata(self.dark_file) self.bad_col = fits.getdata(self.bad_col_file) def load_config(self, config_file): """ load config file. Only for internal use. """ with open(config_file, 'r') as f: config = yaml.load(f, Loader=yaml.FullLoader) self.switch = config['switch'] self.readout_noise = config['readout_noise'] self.ph_per_adu = config['ph_per_adu'] self.bias_level = config['bias_level'] self.max_adu = config['max_adu'] self.dark_file = cpism_refdata + "/camera/" + config['dark_file'] self.flat_file = cpism_refdata + "/camera/" + config['flat_file'] self.cic_file = cpism_refdata + "/camera/" + config['cic_file'] self.bad_col_file = cpism_refdata + \ "/camera/" + config['bad_col_file'] def vertical_blooming(self, image): """ vertical blooming effect """ fullwell = self.fullwell line = np.arange(image.shape[0]) yp, xp = np.where(image > fullwell) n_saturated = len(xp) log.debug(f"{len(xp)} pixels are saturated!") if n_saturated > 5000: log.warning(f"More than 5000({len(xp)}) pixels are saturated!") img0 = image.copy() for x, y in zip(xp, yp): image[:, x] += np.exp(-(line-y)**2/20**2) * img0[y, x] * 0.2 return np.minimum(image, fullwell) def nonlinear_effect(self, image): """ nonlinear effect """ fullwell = self.fullwell nonlinear_coefficient = 0.1 log.debug( f"nonlinear effect added with coefficient {nonlinear_coefficient}") image += (image / fullwell)**2 * nonlinear_coefficient * fullwell return image def emregester_blooming(self, image, max_iteration=5): """ emregester blooming effect """ line = image.flatten().copy() curve_x = np.arange(1300)+2 curve_y = np.exp(11*curve_x**(-0.19)-11) curve_y[0] = 0 curve_y /= curve_y.sum() over_limit_coe = 0.999 saturated = image > self.em_fullwell n_saturated = saturated.sum() if n_saturated > 0: log.debug(f"{n_saturated} pixels are saturated during EM process.") if n_saturated > 2000: log.warning( f"More than 2000 ({n_saturated}) pixels are saturated during EM process!") for index in range(max_iteration): over_limit = np.maximum( line - self.em_fullwell * over_limit_coe, 0) line = np.minimum(line, self.em_fullwell * over_limit_coe) blooming = np.convolve(over_limit, curve_y, mode='full')[ :len(line)] line = line + blooming n_over = (line > self.em_fullwell).sum() if n_over <= 0: break log.debug( f'{index}/{max_iteration} loop: saturated pixel number: {n_over}') return line.reshape(image.shape) def cte(self, image): """ cte effect """ image = self.emregester_blooming(image, max_iteration=5) return image def readout(self, image_focal, emgain, expt, image_cosmic_ray=None): """ emccd readout. Generate a image with emccd readout effect. Parameters ---------- image_focal : numpy.ndarray image at focal plane. Unit: electron/second emgain : float emgain of emccd expt : float exposure time. Unit: second image_cosmic_ray : numpy.ndarray, optional cosmic ray image. Unit: electron/second, by default None Returns ------- numpy.ndarray image with emccd readout effect. Unit: ADU Notes ----- 1. effects include: dark, flat, cte, blooming, nonlinear, etc. Can be turned on/off by switch. 2. size of input image_focal must be 1024x1024 3. size of output image is 1080x1056 (including overscan and dark reference region) 4. Q.E is not included in this function. It should be included in image_focal. See optics.py for details. """ log.debug( fr"EMCCD readout: {emgain=:}, {expt=:}, image_comic_ray:{'None' if image_cosmic_ray is None else 'Not None'}") log.debug( f"camera effects switch={self.switch}" ) image = image_focal * expt if self.switch['flat']: image = image * self.flat if self.switch['nonlinear']: image = self.nonlinear_effect(image) darksz_x = self.plszx + self.rdark + self.ldark darksz_y = self.plszy + self.udark + self.bdark img_dark = np.zeros((darksz_y, darksz_x)) img_dark[ self.bdark:self.plszy+self.bdark, self.ldark:self.ldark+self.plszx ] = image image = img_dark if self.switch['dark']: image += self.dark * expt if self.switch['cic']: image += self.cic if image_cosmic_ray is not None and self.switch['cosmicray']: image += image_cosmic_ray if self.switch['blooming']: image = self.vertical_blooming(image) if self.switch['badcolumn']: for i in range(self.bad_col.shape[1]): deadpix_x = self.bad_col[0, i] deadpix_y = self.bad_col[1, i] image[deadpix_y:, deadpix_x] = 0 biassz_x = darksz_x + self.pscan1 + self.oscan1 biassz_y = darksz_y + self.pscan2 + self.oscan2 img_bias = np.zeros((biassz_y, biassz_x), dtype=int) seed = random_seed_select() log.debug(f"photon noise seed: {seed}") img_bias[ self.pscan2:self.pscan2+darksz_y, self.pscan1:self.pscan1+darksz_x ] = np.random.poisson(image) image = img_bias if self.switch['cte']: image = self.cte(image * emgain) / emgain seed = random_seed_select() log.debug(f"gamma noise seed: {seed}") if emgain != 1: image = np.random.gamma(image, emgain) image = np.minimum(image, self.em_fullwell) seed = random_seed_select() log.debug(f"readout noise seed: {seed}") image += np.random.randn(biassz_y, biassz_x) * self.readout_noise image = image / self.ph_per_adu + self.bias_level if self.switch['stripe']: image += self.add_stripe_effect(image) image = np.minimum(image, self.max_adu) image = np.maximum(image, 0) return image.astype(np.uint16) def add_stripe_effect(self, image): """ add stripe effect """ shape = image.shape v_width = 1 v_amplitude = 30 v_limit = 0.01 v_base = 10 h_width = 20 h_amplitude = 10 h_limit = 0.9 h_base = 20 index = np.linspace(0, np.pi, shape[0] * shape[1]) def stripe(width, limit, amplitude, base, axis=0): seed = random_seed_select() log.debug(f"stripe noise seed: {seed}") dim_axis = shape[axis] dim_other = shape[0] * shape[1] // shape[axis] value = np.sin(index / width * dim_axis + np.pi * dim_axis / width * np.random.randint(1024)) value = np.maximum(value, -limit) value = np.minimum(value, limit) value = (value / limit + limit) / 2 * amplitude + base value = value.reshape(dim_axis, dim_other) if axis == 1: value = value.T return value output = stripe(v_width, v_limit, v_amplitude, v_base, axis=1) output += stripe(h_width, h_limit, h_amplitude, h_base, axis=0) return output # # plt.plot(horizontal_index, horizontal_value) # # # plt.xlim([0, 6.28]) # # plt.show() # fits.writeto('horizontal_value.fits', output, overwrite=True) # if __name__ == '__main__': # import matplotlib.pyplot as plt # emccd = EMCCD() # image_focal = np.zeros((emccd.plszy, emccd.plszx)) + 1000 # image_focal[100:105, 100:105] = 10_000_000 # after_cte = emccd.emregester_blooming(image_focal, max_iteration=100) # print(after_cte.sum(), image_focal.sum()) # fits.writeto('after_cte.fits', after_cte, overwrite=True) # # darksz_x = emccd.plszx + emccd.rdark + emccd.ldark # # darksz_y = emccd.plszy + emccd.udark + emccd.bdark # # iamge_cosmic_ray = np.zeros((darksz_y, darksz_x)) # # emgain = 10 # # expt = 10 # # image = emccd.readout(image_focal, emgain, expt, iamge_cosmic_ray) # # fits.writeto('test.fits', image, overwrite=True) # image = np.zeros((1000, 1000)) # make_cosmic_ray_frame = CosmicRayFrameMaker() # crimage = make_cosmic_ray_frame(image.shape, 3000) # fits.writeto('crimage.fits', crimage, overwrite=True) # # emccd.add_stripe_effect(image)
