Camera characterization#
Estimating the gain#
This section shows how the gain of a camera can be empirically measured. Note that for a typical CMOS camera, the gain can be controlled via software, and therefore depends on camera settings and the selected “readout mode”. One of these settings might be called “GAIN”, but this is just a software parameter, typically between 0 and 100. In this context, the actual gain we are interested in is sometimes called “system gain”, but in the following we’ll just refer to it as gain.
Theory#
Recall that the gain \(g\) can be seen as the conversion factor between electrons (\(e\)) and ADUs (\(a\)), assuming that this gain is constant over the full dynamic range of the camera:
The main idea of the present characterization is that the number of electrons per pixel is drawn from a Poisson distribution. This distribution has the property that its variance equals its expectation value. In other words, if we would repeat an idealized observation of \(e\) in a given pixel a large number of times, illuminating this pixel with a constant source, the standard deviation \(\sigma_e\) of these \(e_i\) would be the square root of their average \(\langle e \rangle\):
As we will see, this statement does not hold when considering ADUs instead of electrons. Indeed the camera “arbitrarily” divides numbers of electrons by \(g\) to get ADUs. This operation also scales the standard deviation \(\sigma_a\) of these ADU values, so that
Putting everything together, we can write:
and equivalently
Therefore, comparing this variance of ADUs over repeated idealized observations of a given pixel to the expected number of ADUs would give us direct access to the gain of that pixel.
Let’s assume that all pixels have the same gain. In this case, instead of computing the variance and average over reapeated observations, it’s tempting to simply take the variance and average over all (or at least many) different pixels of a flat frame. This would however neglect that the pixels have different sensitivities (and illuminations), which increases their spatial variance compared to the pure shot noise. Ignoring this would lead to a biased gain characterization.
A simple trick to mitigate this issue is to measure the spatial variance on the difference between two flats taken in identical conditions. This largely cancels out the systematic variations from the flat, leaving only shot noise. And the variance of a difference (or sum) of uncorrelated random variables is simply the sum of the variances. We can therefore simply divide this variance of the difference by two before comparing it to the average ADU level.
In practice, we’ll perform these measurements on pairs of flats at many different light levels, covering the entire range of the camera. The gain can then be obtained via linear fit, as the inverse slope of a “signal-variance plot”. Of course, a good masterbias should be subtracted from the flat frames prior to the analysis. The dark current can often be ignored, if the exposures are short.
For a more detailed description of the procedure outlined above, see for example https://www.mirametrics.com/tech_note_ccdgain.php .
from pathlib import Path
import numpy as np
import astropy
from astropy import units as u
%matplotlib ipympl
import matplotlib
from matplotlib import pyplot as plt
import ccdproc
from astropy.modeling import models, fitting
from astropy.nddata.blocks import block_reduce
# Path to the raw data files (we follow good practice and won't modify these files):
data_dir = Path("/export/data1/fprak1/Data/2024-09-06")
# Path to a directory where intermediate files and results can be written:
work_dir = Path("/export/data1/fprak1/tests_mtewes/camera_charact/workdir")
work_dir.mkdir(exist_ok=True) # We create this directory in case it does not yet exist.
# Our FITS files have long headers.
# Just to make some summary tables easier to read, we give a list of the few important header keywords we care about:
ifc_header_keywords = ["imagetyp", "filter", "exptime", "object", "xbinning", "ybinning", "naxis1", "naxis2"]
# We'll ignore some astropy warnings that get raised as our FITS headers (from NINA) are not 100% standards compliant.
import warnings
warnings.simplefilter('ignore', category=astropy.wcs.FITSFixedWarning)
Creating a masterbias
bias_files = ccdproc.ImageFileCollection(data_dir / "BIAS", keywords=ifc_header_keywords)
bias_files.summary
# We check the level and noise in these files:
for ccd, filename in bias_files.ccds(ccd_kwargs={'unit': 'adu'}, return_fname=True):
print(f"{filename}: standard deviation {ccd.data.std():.2f} ADU, median level {np.median(ccd.data):.2f} ADU")
combiner = ccdproc.Combiner(bias_files.ccds(ccd_kwargs={"unit":"adu"}))
masterbias = combiner.average_combine()
masterbias.meta['combined'] = True
masterbias.data = masterbias.data.astype('float32') # Converts to float32 to save space
masterbias.write(work_dir / 'masterbias.fits', overwrite=True)
print(f"masterbias: standard deviation {masterbias.data.std():.2f} ADU, median level {np.median(masterbias.data):.2f} ADU")
flat_files = ccdproc.ImageFileCollection(data_dir / "FLAT", keywords=ifc_header_keywords)
flat_files.summary.show_in_notebook()
flat_filters = set(flat_files.summary["filter"])
flat_objects = set(flat_files.summary["object"])
print(f"Filters: {flat_filters}, objects: {flat_objects}")
We select some flats to use for the measurements
gain_filter = 'i'
gain_object = 'Panelflat Zenith'
gain_flat_files = flat_files.filter(object=gain_object, filter=gain_filter)
gain_exptimes = sorted(list(set(gain_flat_files.summary["exptime"])))
print(f"Exposure times of selected flat frames: {gain_exptimes}")
# And measure the variance and signal in different ways
# We prepare a table to hold the results
results = []
# Read the masterbias
masterbias = ccdproc.CCDData.read(work_dir / 'masterbias.fits')
# Loop over the flat frame exposure times
for exptime in gain_exptimes:
files_match_mask = np.isclose(gain_flat_files.summary["exptime"], exptime)
matched_flats = gain_flat_files.summary["file"][files_match_mask]
# We need (at least) two flats per exposure time
assert(len(matched_flats) > 1)
print(f"Processing exptime {exptime}")
# Read them
flat_1 = ccdproc.CCDData.read(matched_flats[0], unit='adu')
flat_2 = ccdproc.CCDData.read(matched_flats[1], unit='adu')
# Subtract the masterbias
flat_1 = ccdproc.subtract_bias(flat_1, masterbias)
flat_2 = ccdproc.subtract_bias(flat_2, masterbias)
# Compute the diff
flat_diff = flat_1.subtract(flat_2)
if False: # Used this to check for issues in the images
plt.figure(figsize=(8, 4))
plt.title(f"Exptime {exptime}")
cbar = plt.imshow(flat_diff, origin='lower', cmap='Greys_r', interpolation='nearest')
plt.colorbar(cbar, label="Pixel value [ADU]")
plt.xlabel("x [pixel]")
plt.ylabel("y [pixel]")
plt.tight_layout()
plt.show()
if False: # Used this to check that the difference is roughly Gaussian
plt.figure(figsize=(8, 4))
plt.title(f"Exptime {exptime}")
plt.hist(flat_diff.data[1000:1200,1000:1200].flatten(), bins=500)
plt.xlabel("Pixel value [ADU]")
plt.tight_layout()
plt.show()
# Compute stats:
res = {}
res["exptime"] = exptime
res["median"] = np.median(flat_1.data)
res["std"] = np.std(flat_1.data)
res["std_diff"] = np.std(flat_diff.data)
res["std_mad_diff"] = 1.4826 * np.median(np.fabs(flat_diff.data - np.median(flat_diff.data)))
res["var"] = res["std"]**2
res["var_diff"] = (res["std_diff"]**2)/2.0
res["var_mad_diff"] = (res["std_mad_diff"]**2)/2.0
results.append(res)
tab = astropy.table.Table(results)
# We robustly estimate the variance of this difference image.
# Robustly, as there could potentially be a few hot pixels or defects that might otherwise strongly influence the result.
#https://en.wikipedia.org/wiki/Median_absolute_deviation
# std_robust = 1.4826 * mad # Assuming Gaussian noise, but that should be ok.
#Median absolute deviation(x_i) = median(|x_i - median(x)|
# Fitting a line through the signal-variance plot
# initialize a linear fitter
fit = fitting.LinearLSQFitter()
# initialize a linear model
line_init = models.Linear1D()
# fit the data with the fitter
fitted_line = fit(line_init, tab["median"], tab["var_mad_diff"])
gain = 1.0 / fitted_line.slope
print(f"Gain: {gain:.3f} e/ADU")
# Manufacturer values:
# https://www.qhyccd.com/astronomical-camera-qhy600/
# Plotting the best fit line and the measurements in a signal-variance plot
plt.rcParams.update({
"text.usetex": True,
"font.family": "Helvetica"
})
plt.figure(figsize=(7, 4))
ax = plt.subplot()
ax.plot(tab["median"], tab["var_diff"], ls="None", marker="o", label="Variance computed using difference of two flats")
ax.plot(tab["median"], tab["var_mad_diff"], ls="None", marker="+", label="Same, but more robust esimation via MAD")
ax.plot(tab["median"], fitted_line(tab["median"]), label=rf"Fitted line: gain = 1/slope = {gain:.3f} e/ADU")
#ax.plot(tab["median"], tab["var"], ls="None", marker="o", label="Var on single flat")
ax.set_xlabel(r"Median flat level [ADU]")
ax.set_ylabel(r"Variance of shot noise [ADU$^2$]")
ax.set_title("Gain estimation (QHY600M, Readout Mode 0, settings GAIN=30 OFFSET=30)")
plt.grid()
plt.legend()
plt.show()
Visualizing the distortion model#
from pathlib import Path
import numpy as np
import astropy
from astropy import units as u
%matplotlib ipympl
import matplotlib
from matplotlib import pyplot as plt
import ccdproc
from astropy.modeling import models, fitting
from astropy.nddata.blocks import block_reduce
data_dir = Path("/export/data1/fprak1/Data/2024-08-27")
work_dir = Path("/export/data1/fprak1/tests_mtewes/camera_charact/workdir")
work_dir.mkdir(exist_ok=True) # We create this directory in case it does not yet exist.
# Our FITS files have long headers.
# Just to make some summary tables easier to read, we give a list of the few important header keywords we care about:
ifc_header_keywords = ["imagetyp", "filter", "exptime", "object", "xbinning", "ybinning", "naxis1", "naxis2"]
light_files = ccdproc.ImageFileCollection(data_dir / "LIGHT_WCS", keywords=ifc_header_keywords)
light_files.summary.show_in_notebook()
after_flip = ccdproc.CCDData.read(data_dir/"LIGHT_WCS"/"2024-08-28_01-33-12_r_-5.00_40.00s_0042.fits", unit="adu")
wcs = after_flip.wcs
h = after_flip.header
nx, ny = h["NAXIS1"], h["NAXIS2"]
print(f"nx: {nx}, ny: {ny}")
x, y = np.arange(0, h["NAXIS1"]), np.arange(0, h["NAXIS2"])
xx, yy = np.meshgrid(x, y)
#u, v = wcs.all_pix2world(xx, yy, 0)
skycoords = wcs.array_index_to_world(xx, yy)
u[0]
plt.subplot()
plt.scatter(u)
plt.grid(color='white', ls='solid')
plt.show()
Comparing different flat frames#
print(f"Filters: {flat_filters}, objects: {flat_objects}")
df_east = flat_files.filter(object='Domeflat East', exptime=0.12)
df_east_exptimes = sorted(list(set(df_east.summary["exptime"])))
print(df_east_exptimes)
df_west = flat_files.filter(object='Domeflat West')
df_west_exptimes = sorted(list(set(df_west.summary["exptime"])))
print(df_west_exptimes)
df_east_play = flat_files.filter(object='Domeflat East after filter wheel rotation')
df_east_play_exptimes = sorted(list(set(df_east_play.summary["exptime"])))
print(df_east_play_exptimes)
pf = flat_files.filter(object='Panelflat Zenith', exptime=1.3)
pf_exptimes = sorted(list(set(pf.summary["exptime"])))
print(pf_exptimes)
df_east_ccd = next(df_east.ccds(ccd_kwargs={'unit': 'adu'}))
print(np.median(df_east_ccd.data))
df_east_ccd.data = df_east_ccd.data / np.median(df_east_ccd.data)
df_west_ccd = next(df_west.ccds(ccd_kwargs={'unit': 'adu'}))
print(np.median(df_west_ccd.data))
df_west_ccd.data = df_west_ccd.data / np.median(df_west_ccd.data)
df_east_play_ccd = next(df_east_play.ccds(ccd_kwargs={'unit': 'adu'}))
print(np.median(df_east_play_ccd.data))
df_east_play_ccd.data = df_east_play_ccd.data / np.median(df_east_play_ccd.data)
pf_ccd = next(pf.ccds(ccd_kwargs={'unit': 'adu'}))
print(np.median(pf_ccd.data))
pf_ccd.data = pf_ccd.data / np.median(pf_ccd.data)
plt.figure(figsize=(7, 4))
cbar = plt.imshow(block_reduce(df_east_ccd.data/df_west_ccd.data, 5, func=np.mean), origin='lower', vmin=0.92, vmax=1.08, cmap='viridis', interpolation='nearest')
plt.colorbar(cbar)
plt.title("Dome flat East of pier / Dome flat West of pier (filter i)")
plt.tight_layout()
plt.gca().get_xaxis().set_visible(False)
plt.gca().get_yaxis().set_visible(False)
plt.show()
plt.figure(figsize=(7, 4))
cbar = plt.imshow(block_reduce(df_east_ccd.data/df_east_play_ccd.data, 100, func=np.mean), origin='lower', vmin=0.999, vmax=1.001, cmap='viridis', interpolation='nearest')
plt.colorbar(cbar)
plt.title("Dome flat before / after a few filter wheel rotations (filter i)")
plt.tight_layout()
plt.gca().get_xaxis().set_visible(False)
plt.gca().get_yaxis().set_visible(False)
plt.show()
plt.figure(figsize=(7, 4))
cbar = plt.imshow(block_reduce(df_east_ccd.data/pf_ccd.data, 5, func=np.mean), origin='lower', vmin=0.92, vmax=1.08, cmap='viridis', interpolation='nearest')
plt.colorbar(cbar)
plt.title("Dome flat East of pier / Panel flat at Zenith (filter i)")
plt.tight_layout()
plt.gca().get_xaxis().set_visible(False)
plt.gca().get_yaxis().set_visible(False)
plt.show()