Field of view optimisation for airborne sensor performance¶
This notebook implements a variation of the performance model developed by T M Goss et al.. Field of View Selection for Optimal Airborne Imaging Sensor Performance, Proc SPIE, Vol 9071, doi:10.1117/12.2049911, 2014. We also use the material developed in Electro-Optical System Analysis and Design: A Radiometry Perspective, CJ Willers, SPIE, DOI: 10.1117/3.1001964, 2013. The first reference is identified as [G] and the second reference is identified as [W] in the text below. If no reference is given, [G] is implied.
The present form of the model can account for atmospheric transmittance and path radiance effects, but the present calculation does not include this. The results presented here therefore only apply to short-range observations.
The model is a much simplified version of reality, and could benefit from improvement in many areas, but even in its present form it is already quite useful. In particular the absence of atmospheric modelling is somewhat restrictive, but can easily be overcome. The model author [G] stated the following (reasonable) constraints:
The model is useful to compare similar imaging sensors. The optimization and comparison process assumes the imaging sensor is well designed and that the following design aspects are taken care of:
- The optical design is such that the exit pupil is placed close to the front lens/window and limited vignetting is permitted thereby ensuring that the real optical aperture diameter equal to the theoretical optical aperture.
- The optics is always optimally focused.
- The optical design is such that no significant stray light or glare issues exist.
- The electronic design is such that no significant additional noise sources exist i.e. other than those identified in this paper.
Similarly the performance of any thermal imaging sensor may be optimized within the constraints of the application using the model, for comparing thermal imagers. The optimization and comparison process assumes the thermal imaging sensor is well designed and that the following design aspects are taken care of:
- The same assumptions regarding the optical and electronic designs above apply except that vignetting is not permitted.
- The cold stop efficiency is 100% (achieved by imaging the entrance pupil onto the cold stop).
- A cold filter ensures no out-of-band radiation reaches the detector pixels.
- The detector temperature regulation is sufficient and no low frequency noise components compromise the thermal imaging sensor’s performance i.e. only white noise is present.
Prepare environment¶
This report is created from an IPython notebook which includes all the code to calculate the results.
# to prepare the Python environment
import numpy as np
import scipy as sp
import pandas as pd
import os.path
from scipy.optimize import curve_fit
from scipy import interpolate
from scipy import integrate
from scipy import signal
from scipy import ndimage
import matplotlib.pyplot as plt
import scipy.constants as const
import pickle
import scipy.constants as const
import collections
%matplotlib inline
# %reload_ext autoreload
# %autoreload 2
import pyradi.ryplot as ryplot
import pyradi.ryplanck as ryplanck
import pyradi.ryfiles as ryfiles
import pyradi.rymodtran as rymodtran
import pyradi.ryutils as ryutils
import pyradi.rypflux as rypflux
import pyradi.rystare as rystare
from IPython.display import HTML
from IPython.display import Image
from IPython.display import display
from IPython.display import FileLink, FileLinks
import matplotlib as mpl
mpl.rc("savefig", dpi=300)
mpl.rc('figure', figsize=(10,8))
# %config InlineBackend.figure_format = 'svg'
pim = ryplot.ProcessImage()
pd.set_option('display.max_columns', 80) # Number of columns in display
pd.set_option('display.width', 80) # display width - how many chars
pd.set_option('display.max_colwidth', 150)
Target Model¶
The target illumination model is described in 07-Optical-Sources.ipynb.
# to print the available scenarios and spectral bands in the dfPhotRates table
pf = rypflux.PFlux()
dfPhotRates = pf.lllPhotonrates()
print(dfPhotRates.columns)
print(dfPhotRates.index)
print(dfPhotRates)
Index([u'ColourTemp', u'FracPhotop', u'Irradiance-lm/m2',
u'Radiance-q/(s.m2.sr)-LWIR', u'Radiance-q/(s.m2.sr)-MWIR',
u'Radiance-q/(s.m2.sr)-NIR', u'Radiance-q/(s.m2.sr)-SWIR',
u'Radiance-q/(s.m2.sr)-VIS', u'k'],
dtype='object')
Index([u'Overcast night', u'Star light', u'Quarter moon', u'Full moon',
u'Deep twilight', u'Twilight', u'Very dark day', u'Overcast day',
u'Full sky light', u'Sun light'],
dtype='object')
ColourTemp FracPhotop Irradiance-lm/m2 \
Overcast night 5000.0 0.0 0.0001
Star light 5000.0 0.2 0.0011
Quarter moon 4150.0 0.4 0.0108
Full moon 4150.0 0.6 0.1080
Deep twilight 10000.0 0.8 1.0800
Twilight 10000.0 1.0 10.8000
Very dark day 7000.0 1.0 107.0000
Overcast day 6000.0 1.0 1075.0000
Full sky light 12000.0 1.0 10752.0000
Sun light 5700.0 1.0 107527.0000
Radiance-q/(s.m2.sr)-LWIR Radiance-q/(s.m2.sr)-MWIR \
Overcast night 1.098964e+10 4.126570e+10
Star light 1.352951e+11 5.080282e+11
Quarter moon 3.790946e+12 1.362625e+13
Full moon 4.296000e+13 1.544163e+14
Deep twilight 2.807479e+13 1.166828e+14
Twilight 3.819177e+14 1.587304e+15
Very dark day 8.245559e+15 3.284063e+16
Overcast day 1.306910e+17 5.080461e+17
Full sky light 2.859098e+17 1.207679e+18
Sun light 1.551574e+19 5.977276e+19
Radiance-q/(s.m2.sr)-NIR Radiance-q/(s.m2.sr)-SWIR \
Overcast night 1.320815e+11 2.811788e+11
Star light 1.626075e+12 3.461634e+12
Quarter moon 2.676633e+13 7.321109e+13
Full moon 3.033231e+14 8.296474e+14
Deep twilight 1.107208e+15 1.361028e+15
Twilight 1.506199e+16 1.851486e+16
Very dark day 1.998632e+17 3.068274e+17
Overcast day 2.380423e+18 4.172085e+18
Full sky light 1.347956e+19 1.529608e+19
Sun light 2.537434e+20 4.676317e+20
Radiance-q/(s.m2.sr)-VIS k
Overcast night 1.446768e+11 1.613462e-14
Star light 1.781137e+12 1.986357e-13
Quarter moon 2.198762e+13 6.956533e-12
Full moon 2.491695e+14 7.883327e-11
Deep twilight 2.488231e+15 1.887297e-11
Twilight 3.384886e+16 2.567401e-10
Very dark day 3.306346e+17 8.220586e-09
Overcast day 3.313916e+18 1.552340e-07
Full sky light 3.397997e+19 1.578701e-07
Sun light 3.315415e+20 1.955067e-05
Sensor Model¶
# to define the sensor parameters used in this analysis
foclen = 1. # normalised focal length in m
# read the excel spreadsheet containing the sensor data and extract the example data
dfS = pd.read_excel('data/FOV-optimisation-parameters.xlsx', 'Sheet1', index_col=None, na_values=['NA'])
dfS.set_index(keys='Sensor',inplace=True)
# get rid of empty lines in excel spreadsheet
dfS.dropna(inplace=True)
# get a list of the spectral band definitions from the data set
dfS.loc['Spectral band Mid'] = (dfS.loc['Spectral band Lo'] + dfS.loc['Spectral band Hi']) / 2.
specBands = dfS.loc['Band'].unique()
wlcs = {}
for specBand in specBands:
wlcs[specBand] = dfS.loc['Spectral band Mid'][dfS.loc['Band']==specBand][0]
# get list of MTF labels and matching masks
mtfLabels = dfS.loc['MTF label'].unique()
dOptics = {}
for mtfLabel in mtfLabels:
dOptics[mtfLabel] = dfS.loc['Pupil mask'][dfS.loc['MTF label']==mtfLabel][0]
# get max spatial frequency for sensor aperture and wavelength
dfS.loc['SpatFreqMax cy/m'] = dfS.loc['Optics diameter'] / dfS.loc['Spectral band Mid']
dfS = dfS.append(pd.DataFrame(index=['SpatFreq cy/m']))
for colname in dfS.columns.values:
dfS.loc['SpatFreq cy/m'][colname] = np.linspace(0, dfS.loc['SpatFreqMax cy/m'][colname], 100)
print(specBands)
print(wlcs)
print(mtfLabels)
print(dOptics)
print(dfS.drop('SpatFreq cy/m', axis=0))
[u'Visible' u'MWIR']
{u'Visible': 5.399999999999999e-07, u'MWIR': 4.350000000000001e-06}
[u'Central 0.25' u'No obscuration']
{u'No obscuration': u'data/mask00.png', u'Central 0.25': u'data/mask01.png'}
MWIR Sensor A MWIR Sensor B Vis Sensor A \
Sensor
Band MWIR MWIR Visible
Sensor Type MCT MCT CMOS
Spectral band Lo 3.7e-06 3.7e-06 4.8e-07
Spectral band Hi 5e-06 5e-06 6e-07
Pixels H 640 640 1280
Pixels V 512 512 720
Pixel size 1.5e-05 1.5e-05 5.5e-06
Fill factor 0.98 0.98 0.98
Quantum efficiency 0.75 0.75 0.25
Well capacity 4980000 4980000 13500
Optics diameter 0.2 0.16 0.2
Obscuration ratio 0.25 0 0.25
Pupil mask data/mask01.png data/mask00.png data/mask01.png
MTF label Central 0.25 No obscuration Central 0.25
Aberration size factor 0.23 0.23 0.23
Optical transmittance 0.74 0.8 0.5
Frame rate 100 100 50
RMS jitter rad 5e-06 5e-06 5e-06
Linear smear rad/s 0.0005 0.0005 0.0005
Scene temperature 298 298 2000
Optics temperature 298 298 298
F-number limit 4 4 4
Minimum modulation 0.1 0.1 0.1
Maximum HFOV deg 2.5 2.5 2.5
Spectral band Mid 4.35e-06 4.35e-06 5.4e-07
SpatFreqMax cy/m 45977 36781.6 370370
Vis Sensor B
Sensor
Band Visible
Sensor Type CMOS
Spectral band Lo 4.8e-07
Spectral band Hi 6e-07
Pixels H 1280
Pixels V 720
Pixel size 5.5e-06
Fill factor 0.98
Quantum efficiency 0.25
Well capacity 13500
Optics diameter 0.11
Obscuration ratio 0
Pupil mask data/mask00.png
MTF label No obscuration
Aberration size factor 0.23
Optical transmittance 0.7
Frame rate 50
RMS jitter rad 5e-06
Linear smear rad/s 0.0005
Scene temperature 2000
Optics temperature 298
F-number limit 4
Minimum modulation 0.1
Maximum HFOV deg 2.5
Spectral band Mid 5.4e-07
SpatFreqMax cy/m 203704
Signal model¶
The [G] model calculates the signal and noise components from primary sensor parameters, allowing for the sensor MTF. The model describes the noise equivalent temperature or noise equivalent reflectance in terms of spatial frequency and then solves for the noise equivalent terms in terms of spatial frequency.
The sensor's perceived signal to noise ratio when imaging a square wave bar target in the presence of sensor noise is given by
\begin{equation} \textrm{SNR} = \frac{\textrm{CTF}(\xi) \cdot [n_{\textrm{signal}} - n_{\textrm{background}}] }{n_{\textrm{noise}}} \end{equation}where $\textrm{SNR}$ is the signal to noise ratio, $\textrm{CTF}(\xi)$ is the sensor's contrast transfer function, $\xi$ is spatial frequency, $n_{\textrm{signal}}$ is the signal electron count in the detector, $n_{\textrm{background}}$ is the background electron count in the detector, and $n_{\textrm{noise}}$ is the sensorn noise at the sensor output. The signal as defined here is the signal in one pixel.
Using a variation of the camera equation (Equation 9.56 in [W]) the electron count in the detector from a source contrast radiance $\Delta L$ is given by the flux on the detector multiplied by the quantum efficiency and the integration time:
\begin{equation} \Delta n_{\textrm{det}\; \lambda} = \frac{\pi \eta_{\lambda} t_{\textrm{int}} P K_N(\alpha) A_d \tau_{o\lambda}\tau_{f\lambda} \tau_{a\lambda}\left[L_{\textrm{signal}\lambda}- L_{\textrm{background}\lambda}\right] \cos^4\alpha}{4F_\#^2(1+|m|)^2}, \end{equation}where $\eta_{\lambda}$ is the detector spectral quantum efficiency, $t_{\textrm{int}}$ is the detector integration time, $P$ an efficiency factor accounting for a central obscuration (if present), $K_N(\alpha)$ is the vignetting losses (if present), $A_d$ is the area of the detector, $\tau_{o\lambda}$ is the optics transmittance, $\tau_{f\lambda}$ is the spectral filter transmittance, $\tau_{a\lambda}$ is the atmospheric transmittance, $L_{\textrm{signal}\lambda} = \epsilon_\textrm{signal} L_\lambda(T_\textrm{signal})$ is the signal source radiance, $L_{\textrm{background}\lambda} = \epsilon_\textrm{background} L_\lambda(T_\textrm{background})$ is the background source radiance, $\cos^4\alpha$ is the cosine-to-the-fourth losses, $\alpha$ is the field angle, $F_\#$ is the sensor f-number, and $m$ is the magnification.
At inifinte conjugates $m$=0, with on-axis viewing $\alpha$=0, with no vignetting $K_N(\alpha)$=0, the camera equation for target contrast becomes
\begin{equation} \Delta n_{\textrm{det}} = \left(\frac{\pi t_{\textrm{int}} P A_d }{4F_\#^2}\right)\cdot \int_0^\infty \eta_{\lambda} \tau_{o\lambda}\tau_{f\lambda} \tau_{a\lambda}\left[L_{\textrm{signal}\lambda}- L_{\textrm{background}\lambda}\right] d\lambda \end{equation}Rejoining [G] we define the noise equivalent quanta as the the contrast quantum signal where the signal to noise is equal to one. In keeping with [W] notation noise equivalent quanta is written as NEQ. Then
\begin{equation} \textrm{NEQ} = \frac{n_{\textrm{noise}}}{\textrm{CTF}(\xi)} = \left(\frac{\pi t_{\textrm{int}} P A_d }{4F_\#^2}\right)\cdot \int_0^\infty \eta_{\lambda} \tau_{o\lambda}\tau_{f\lambda} \tau_{a\lambda}\left[L_{\textrm{signal}\lambda}- L_{\textrm{background}\lambda}\right] d\lambda \end{equation}From [G]: If we assume the square wave response of the sensor is band-limited and that for the higher spatial frequencies only the fundamental of the square wave has sufficient response, then the amplitude conversion from the square wave $\textrm{CTF}$ to the sinusoid modulation transfer function, $\textrm{MTF}$, is a factor $4/\pi$
\begin{equation} \textrm{CTF}(\xi) \approx \frac{4}{\pi}\textrm{MTF}(\xi) \end{equation}After manipulation we arrive at Equation 7 in [G]:
\begin{equation} \pi\int_0^\infty \eta_{\lambda} \tau_{o\lambda}\tau_{f\lambda} \tau_{a\lambda}\textrm{LEQ}_\lambda d\lambda = \frac{\pi F_\#^2 n_{\textrm{noise}}}{ t_{\textrm{int}} P A_d\textrm{MTF}(\xi)} \end{equation}with units of q/(s$\cdot$m$^2$), and where $\textrm{LEQ}_\lambda = L_{\textrm{signal}\lambda}- L_{\textrm{background}\lambda}$ is understood to be the contrast radiance between two pixels that provides a signal to noise ratio of one.
Thermal sources: Expanding on Section 9.5.4 in [W] that for small temperature differences, the noise equivalent temperature difference for a thermal source can be approximated by
\begin{equation} \frac{\textrm{LEQ}}{\textrm{NETD}} = \frac{dL}{dT} = \int_0^\infty \eta_{\lambda} \tau_{o\lambda}\tau_{f\lambda} \tau_{a\lambda} \left( \frac{L_\lambda(T_{\textrm{background}})}{dT} \right) d\lambda \end{equation}Hence for thermal sources a form similar to Equation 9 in [G] is derived
\begin{equation} \textrm{NETD} = \frac{ F_\#^2]\, n_{\textrm{noise}}}{ t_{\textrm{int}} P A_d\textrm{MTF}(\xi)\left[ \int_0^\infty \eta_{\lambda} \tau_{o\lambda}\tau_{f\lambda} \tau_{a\lambda} \left( \frac{L_\lambda(T_{\textrm{background}})}{dT} \right) d\lambda\right]} \end{equation}Reflective sources: Similarly to the thermal source derivation, for target reflection differences, the noise equivalent reflection difference for a reflective source can be approximated by
\begin{equation} \frac{\textrm{LEQ}}{ \textrm{NE}\rho} = \frac{dL}{d\rho} = \int_0^\infty \eta_{\lambda} \tau_{o\lambda}\tau_{f\lambda} \tau_{a\lambda} L_{\textrm{source}\lambda} d\lambda \end{equation}Hence for reflective sources we derive an equation similar to Equation 7 in [G]:
\begin{equation} \textrm{NE}\rho = \frac{ F_\#^2\, n_{\textrm{noise}}}{ t_{\textrm{int}} P A_d\textrm{MTF}(\xi) \int_0^\infty \eta_{\lambda} \tau_{o\lambda}\tau_{f\lambda} \tau_{a\lambda} L_{\textrm{source}\lambda} d\lambda} \end{equation}Staring array sensor model¶
See the notebook 09b-StaringArrayDetectors.ipynb for theory.
The functions to calculate the electrons in the charge well are contained in rystare.py.
# to evaluate typical electron count from different contributors
tauAtmo = 1
tauSun = 0.5
tauFilt = 1.
tauOpt = 0.8
quantEff = 0.7
rhoTarg = .3
cosTarg = 1.
inttime = 0.02
pfrac = 1.0
detarea = (10e-6)**2
fno = 2.
tmptr = 300.
tmptrOpt = 300.
emisscene = .8
scenarios = dfPhotRates.index.values
print(scenarios)
print('electrons in charge well and photon noise')
for wl, specBand in zip([np.linspace(0.3, 0.8, 100), np.linspace(3.5, 4.8, 100)],
[u'Radiance-q/(s.m2.sr)-VIS','Radiance-q/(s.m2.sr)-MWIR']):
print('\n{}{}'.format(34*'-',specBand[specBand.rfind('-')+1:]))
print('{:18s} {} {}'.format('Dataframe data ','elec', 'noise'))
if specBand in dfPhotRates.columns.values:
for scenario in scenarios:
ncnt = rystare.nEcntLLightDF(tauAtmo, tauFilt, tauOpt, quantEff, rhoTarg, cosTarg,
inttime, pfrac, detarea, fno, scenario, specBand, dfPhotRates)
print('{:18s} {:.2e} {:.2e}'.format(scenario,ncnt, np.sqrt(ncnt)))
print('{:18s} {} {}'.format('\nTheoretical data ','elec', 'noise'))
ncnt = rystare.nElecCntReflSun(wl, tauSun, tauAtmo, tauFilt, tauOpt, quantEff, rhoTarg, cosTarg,
inttime, pfrac, detarea, fno)
print('Direct sun {:.2e} {:.2e}'.format(ncnt, np.sqrt(ncnt)))
ncnt = rystare.nElecCntThermalScene(wl, tmptr, emisscene, tauAtmo, tauFilt, tauOpt, quantEff, inttime,
pfrac, detarea, fno)
print('Scene termal {:.2e} {:.2e}'.format(ncnt, np.sqrt(ncnt)))
ncnt = rystare.nEcntThermalOptics(wl, tmptrOpt, tauFilt, tauOpt, quantEff, inttime, pfrac,
detarea, fno)
print('Hot optics thermal {:.2e} {:.2e}'.format(ncnt, np.sqrt(ncnt)))
['Overcast night' 'Star light' 'Quarter moon' 'Full moon' 'Deep twilight' 'Twilight' 'Very dark day' 'Overcast day' 'Full sky light' 'Sun light'] electrons in charge well and photon noise ----------------------------------VIS Dataframe data elec noise Overcast night 9.54e-03 9.77e-02 Star light 1.18e-01 3.43e-01 Quarter moon 1.45e+00 1.20e+00 Full moon 1.64e+01 4.05e+00 Deep twilight 1.64e+02 1.28e+01 Twilight 2.23e+03 4.73e+01 Very dark day 2.18e+04 1.48e+02 Overcast day 2.19e+05 4.68e+02 Full sky light 2.24e+06 1.50e+03 Sun light 2.19e+07 4.68e+03 Theoretical data elec noise Direct sun 2.60e+07 5.10e+03 Scene termal 3.31e-14 1.82e-07 Hot optics thermal 1.03e-14 1.02e-07 ----------------------------------MWIR Dataframe data elec noise Overcast night 2.72e-03 5.22e-02 Star light 3.35e-02 1.83e-01 Quarter moon 8.99e-01 9.48e-01 Full moon 1.02e+01 3.19e+00 Deep twilight 7.70e+00 2.77e+00 Twilight 1.05e+02 1.02e+01 Very dark day 2.17e+03 4.65e+01 Overcast day 3.35e+04 1.83e+02 Full sky light 7.97e+04 2.82e+02 Sun light 3.94e+06 1.99e+03 Theoretical data elec noise Direct sun 2.50e+06 1.58e+03 Scene termal 5.10e+06 2.26e+03 Hot optics thermal 1.59e+06 1.26e+03
Dark current shot noise¶
The average dark current noise electron count is given by the model from here.
# to test dark current noise function
inttime = 0.01
detarea = (10e-6)**2
Egap = 1.11
for temptr in [300, 400]:
print('inttime={} s, det area={} m2, temperature={} K, bandgap={} eV'.format(inttime, detarea,temptr, Egap))
dcurnoise = rystare.darkcurrentnoise(inttime, detarea,temptr, Egap)
print('Dark current = {} e-, dark current noise = {} e-'.format(dcurnoise**2, dcurnoise))
print(' ')
inttime=0.01 s, det area=1e-10 m2, temperature=300 K, bandgap=1.11 eV Dark current = 3.14484481049e-08 e-, dark current noise = 0.000177337103013 e- inttime=0.01 s, det area=1e-10 m2, temperature=400 K, bandgap=1.11 eV Dark current = 1.03730276592e-05 e-, dark current noise = 0.00322071850045 e-
Fixed pattern dark current noise¶
Fixed pattern dark noise is not considered here because of the relatively short integration times (compared to astronomy sensors).
Read noise also known as kTC reset noise¶
The kTC electron count noise is given by model from here.
# to test kTC noise function
import itertools
gvs = [50e-6, 100e-6, 150e-6]
temps = [300, 400]
grid = itertools.product(gvs, temps)
df = pd.DataFrame()
for (gv,temp) in list(grid):
df = df.append(pd.DataFrame({'gv':[gv],'temp':[temp]}),ignore_index=True)
df['nc'] = rystare.kTCnoiseGv(df.temp, df.gv)
print('kTC noise in electrons:')
print(df.pivot_table(values='nc',index=['gv'],columns=['temp']))
kTC noise in electrons: temp 300 400 gv 0.00005 22.738510 26.256169 0.00010 16.078554 18.565915 0.00015 13.128085 15.159006
Source follower noise¶
Source follower noise is ignored in this study.
Quantisation noise¶
Quantisation follower noise is ignored in this study.
Total noise¶
The total noise is then given by
\begin{equation} n_{\textrm{noise}} = \sqrt{n_{\textrm{photnoise}}^2+ n_{\textrm{NU}}^2 + n^2_{\textrm{darkcurrentnoise}} + n^2_{\textrm{kTCnoise}}} \end{equation}For visual band sensors (corresponding to Equation 10 in [G]) \begin{equation} n_{\textrm{noise}} = \sqrt{n_{\textrm{photnoise}}^2+ (U n_{\textrm{photons}})^2 + n^2_{\textrm{darkcurrentnoise}} + n^2_{\textrm{kTCnoise}}} \end{equation}
For infrared sensors \begin{equation} n_{\textrm{noise}} = \sqrt{n_{\textrm{photnoise}}^2+ \left(\sqrt{n_{\textrm{photons}}}\right)^{-2} + n^2_{\textrm{darkcurrentnoise}} + n^2_{\textrm{kTCnoise}}} \end{equation} \begin{equation} n_{\textrm{noise}} = \sqrt{n_{\textrm{photons}}+ n_{\textrm{photons}}^{-1} + n^2_{\textrm{darkcurrentnoise}} + n^2_{\textrm{kTCnoise}}} \end{equation}
MTF components¶
Detector MTF¶
The detector MTF is modelled as a sinc function
where $\xi$ is spatial frequency in cy/rad, and $\alpha$ is the detector angular size in rad, $\beta$ is the detector size in meters, $f$ is the focal length in meters.
The dimensional analysis of the argument is not immediately evident: one detector subtends half a cycle and this must be accounted for in the dimensional analysis:
\begin{equation} \left[\frac{\textrm{rad}}{1}\right] \left[\frac{\textrm{cy}}{\textrm{rad}}\right] \left[\frac{\textrm{rad}}{\textrm{pix}}\right] \left[\frac{\textrm{pix}}{\textrm{cy}}\right] \rightarrow \left[\frac{\textrm{rad}}{1}\right] \left[\frac{\textrm{cy}}{\textrm{rad}}\right] \left[\frac{\textrm{m}}{\textrm{pix}}\right] \left[\frac{\textrm{pix}}{\textrm{cy}}\right] \left[\frac{\textrm{1}}{\textrm{m}}\right] \rightarrow \textrm{rad} \end{equation}Note: the numpy sinc function already includes a $\pi$ in the function code, don't add it in your code!!!
# to calculate detector MTF
def MTFdetector(freq, pixsize, foclen):
# freq in cy/rad
# pixsize in m
# foclen in m
return np.sinc(freq * pixsize / foclen)
# to plot detector MTF as function of detector pixel size
p = ryplot.Plotter (1 ,1 ,1 , figsize =(8 ,4) )
if 'MTF detector' not in dfS.index.values:
dfS = dfS.append(pd.DataFrame(index=['MTF detector']))
ffreqmax = np.max(dfS.loc['SpatFreqMax cy/m'])
for colname in dfS.columns.values:
ffreq = dfS.loc['SpatFreq cy/m'][colname]
pixsize = dfS.loc['Pixel size'][colname]
dfS.loc['MTF detector'][colname] = MTFdetector(ffreq,pixsize,foclen)
p.plot(1, ffreq/1000., np.abs(dfS.loc['MTF detector'][colname]),
r'Horizontal and vertical detector MTF','Frequency cy/mm','MTF',
label=[r'{} {} $\mu$m pixel size'.format(colname, pixsize*1e6)],
pltaxis =[0,ffreqmax/1000.,-0,1]);
Optics aberration MTF¶
From [G], the optics aberration can be modelled by the generic function
\begin{equation} MTF_\textrm{opt} = \exp \left[ -2\left( \frac{\pi \sigma_\textrm{d}\xi}{\xi_0} \right)^2\right] = \exp \left[ -2\left( \frac{\pi \sigma_\textrm{d}\xi\lambda}{D} \right)^2\right] \end{equation}where $\sigma_\textrm{d}$ is the spot size factor, $\xi$ is spatial frequency in cy/m, $\xi_0$ is cutof or critical spatial frequency in cy/m, $\lambda$ is the mean wavelength, and $D$ is the aperture diameter in m.
The units of the square root of the argument are
\begin{equation} \left[\frac{\textrm{rad}}{\textrm{1}}\right] \left[\frac{\textrm{m}}{\textrm{m}}\right] \left[\frac{\textrm{cy}}{\textrm{rad}}\right] \left[\frac{\textrm{rad}}{\textrm{cy}}\right] \rightarrow \textrm{rad} \rightarrow \left[\frac{\textrm{m}}{\textrm{m}}\right] \rightarrow \textrm{1} \end{equation}and
\begin{equation} \left[\frac{\textrm{rad}}{\textrm{1}}\right] \left[\frac{\textrm{m}}{\textrm{m}}\right] \left[\frac{\textrm{cy}}{\textrm{rad}}\right] \left[\frac{\textrm{m}_{\lambda}}{\textrm{1}} \right] \left[\frac{\textrm{1}}{\textrm{m}_\textrm{d}} \right] \rightarrow \left[\frac{\textrm{m}}{\textrm{m}}\right] \rightarrow \textrm{1} \end{equation}# to calculate aberration MTF
def MTFaberation(freq, sratio, freqo):
# freq in cy/rad
# spot size factor
# cutof freq in cy/mrad
return np.exp(-2. * (np.pi * sratio * freq / freqo) ** 2)
# to calculate the MTF of diffraction through a circular aperture
def MTFcirc(freq):
f = np.abs(freq / np.max(freq))
return 2. * ( np.arccos(f) - f * np.sqrt(1-f**2)) / np.pi
# to plot the aberration MTF as function of spot size factor
p = ryplot.Plotter (1 ,1 ,1 , figsize =(8,4) )
if 'MTF aberration' not in dfS.index.values:
dfS = dfS.append(pd.DataFrame(index=['MTF aberration']))
ffreqmax = np.max(dfS.loc['SpatFreqMax cy/m'])
for colname in dfS.columns.values:
ffreq = dfS.loc['SpatFreq cy/m'][colname]
sratio= dfS.loc['Aberration size factor'][colname]
dfS.loc['MTF aberration'][colname] = MTFaberation(ffreq, sratio, dfS.loc['SpatFreqMax cy/m'][colname])
p.plot(1, ffreq/1000.,dfS.loc['MTF aberration'][colname],
r'Aberration MTF approximation','Frequency cy/mm','MTF',
label =[r'{}: {} spot size factor'.format(colname,sratio)],
pltaxis =[0,ffreqmax/1000.,-0,1]) ;
Diffraction MTF¶
The optical transfer function is the Fourier transform of the point spread function, and the point spread function is the square absolute of the inverse Fourier transformed pupil function. The optical transfer function can also be calculated directly from the pupil function. From the convolution theorem it can be seen that the optical transfer function is the auto-correlation of the pupil function https://en.wikipedia.org/wiki/Optical_transfer_function. If an image of the pupil function is known the diffraction MTF of the pupil can be calculated numerically.
The fractional clear aperture and area of the sensor are calculated by simply counting the number of pixels in the clear aperture (accounting for the fact that the aperture is circular, not square as is the image). The MTF is calculated using the SciPy two-dimensional correlation function signal.correlate2d.
In this calculation the MTF for a clear aperture and centrally obscured aperture is calculated, both of which have closed form analytical equations. The advantage of calculating the MTF from the pupil function is that any shape aperture can be used; just create an image of the pupil. The calculation takes a long time to calculate for large images, so the image can be subsampled to the required degree of accuracy.
In the MTF calculations below the MTF of a round aperture is calculated in square image coordinates. When the frequency scale is calculated, the distance to the corner of the image determines the full scale frequency of the calculation. This full scale frequency is an artifact of the calculation, and is not the aperture cutof frequency that has optical meaning. The full scale frequency is given by
\begin{equation} \xi_\textrm{fs} = \frac{\sqrt{2} D}{\lambda f} \end{equation}where the $\sqrt{2}$ accounts for the distance to the corner (relative to the round aperture diameter), $D$ is the aperture diameter, $\lambda$ is the wavelength, and $f$ is the distance to the screen (e.g., focal length).
Then the optical cutof frequency or critical frequency is given by
\begin{equation} \xi_0 = \xi_\textrm{fs}/\sqrt{2} = \frac{ D}{\lambda f} \end{equation}The calculation below was performed for a screen/focal distance of 1 m. This essentially means that the units of $\xi_0$ is cycles/rad. The spatial frequency in cycles/m for any focal length $f$ can then be found as $\xi_\textrm{cycles/m} = \xi_\textrm{cycles/rad}/f$
# to calculate the autocorrelation and clear area and write to file
# this function may take a long time to execute
def pupilautocor(PIK, dOptics, sampling = 5):
fracclear = {}
areasens = {}
pupil = {}
corrn = {}
polar_corrn = {}
polar_radius = {}
meanmtf = {}
for i,mtflabel in enumerate(dOptics.keys()):
filename = dOptics[mtflabel]['Pupil']
print(mtflabel,filename)
#read the mask image
pup = ndimage.imread(filename).astype(float)
#make binary image
pup[...,-1] = np.where(pup[...,-1]>0, 1.0, 0.0)
#area-circle=pi*r^2, area-square=4*r^2, ratio=pi/4
fracclear[mtflabel] = (4. / np.pi) * np.sum(pup[...,-1])/(pup[...,-1].shape[0]*pup[...,-1].shape[1])
# slice, cut out one plane and exact size
pupil[mtflabel] = pup[::sampling,0:pup.shape[0]:sampling,3]
#http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.signal.correlate2d.html
corr = signal.correlate2d(pupil[mtflabel], pupil[mtflabel], boundary='fill', mode='full')
corrn[mtflabel] = corr / np.max(corr)
#do for autocorrelation
polar_c, _, _ = pim.reprojectImageIntoPolar(corrn[mtflabel].reshape(corrn[mtflabel].shape[0],corrn[mtflabel].shape[1],1), None, False,cval=0.)
polar_c = polar_c[:,:,0]
polar_corrn[mtflabel] = polar_c
#do for radius
polar_c, _, _ = pim.reprojectImageIntoPolar(pupil[mtflabel].reshape(pupil[mtflabel].shape[0],pupil[mtflabel].shape[1],1), None, False,cval=0.)
polar_c = np.flipud(polar_c[:,:,0])
polar_radius[mtflabel] = polar_c
meanmtf[mtflabel] = polar_corrn[mtflabel].mean(axis=1)
data = (pupil, corrn, pupfn, mtflabels, fracclear, areasens, polar_corrn, polar_radius, meanmtf)
with open(PIK, "wb") as f:
pickle.dump(data, f)
# to calculate the mtf
#NOTE!!!!!!!!
# The reprojectImageIntoPolar function maps radial to cartesian coords.
# The radial image is however presented in a cartesian grid, the corners have no meaning.
# The radial coordinates are mapped to the radius, not the corners.
# This means that in order to map corners, the frequency is scaled with sqrt(2),
# The corners are filled with the value specified in cval, which must be zero
mtfpclfilename = 'ApertureMTF.pcl'
freqfsm = {}
rho = {}
fcrit = {}
#calculate if not already existing
if not os.path.exists(mtfpclfilename):
pupilautocor(mtfpclfilename, dOptics, sampling = 5)
#load the correlation data from the pickle file
with open(mtfpclfilename, "rb") as f:
(pupil, corrn, pupfn, mtflabels, fracclear, areasens, polar_corrn, polar_radius, meanmtf) = pickle.load(f)
I = ryplot.Plotter(1,2,4,'', figsize=(8,4));
for i,mtflabel in enumerate(dOptics.keys()):
I.showImage(4*i+1, pupil[mtflabel], ptitle=mtflabel, cmap=ryplot.cubehelixcmap(), titlefsize=10,
cbarshow=False);
I.showImage(4*i+2, polar_radius[mtflabel], ptitle='Unfolded pupil', cmap=ryplot.cubehelixcmap(), titlefsize=10,
cbarshow=False);
I.showImage(4*i+3, corrn[mtflabel], ptitle='Pupil diffraction MTF', cmap=ryplot.cubehelixcmap(), titlefsize=10,
cbarshow=True, cbarorientation = 'vertical', cbarfontsize = 7);
I.showImage(4*i+4, polar_corrn[mtflabel], ptitle='Pupil diffraction MTF', cmap=ryplot.cubehelixcmap(), titlefsize=10,
cbarshow=True, cbarorientation = 'vertical', cbarfontsize = 7);
p = ryplot.Plotter(2,1,1,'', figsize=(10,4));
if 'MTF diffraction' not in dfS.index.values:
dfS = dfS.append(pd.DataFrame(index=['MTF diffraction']))
ffreqmax = np.max(dfS.loc['SpatFreqMax cy/m'])
for colname in dfS.columns.values:
freqfsm = np.sqrt(2.) * dfS.loc['SpatFreqMax cy/m'][colname]
mtflabel = dfS.loc['MTF label'][colname]
rho = np.linspace(0,freqfsm,polar_corrn[mtflabel].shape[0])
ffreq = dfS.loc['SpatFreq cy/m'][colname]
dfS.loc['MTF diffraction'][colname] = np.interp(ffreq, rho, meanmtf[mtflabel])
p.plot(1, ffreq/1000.,dfS.loc['MTF diffraction'][colname],
r'Diffraction MTF approximation','Frequency cy/mm','MTF',
label =[r'{}'.format(colname)],pltaxis =[0,ffreqmax/1000.,-0,1]) ;
Line of sight movement MTF¶
Line of sight movement has two effects: (1) linear motion smear of the image during integration time, and (2) random movement blur during integration time. The two degradation mechanisms have to different effect in the image.
Image degradation attributable to linear smear is modelled as a sinc function
where $\xi$ is spatial frequency in cy/rad, $\dot{\theta}$ is linear motion in rad/s, and $t_i$ is integration time in s
The units of the argument are \begin{equation} \left[\frac{\textrm{rad}}{\textrm{cy}}\right] \left[\frac{\textrm{cy}}{\textrm{rad}}\right] \left[\frac{\textrm{rad}}{\textrm{s}}\right] \left[\frac{\textrm{s}}{\textrm{1}}\right] \rightarrow \textrm{rad} \end{equation}
The [G] model does not account for linear smear movement. Experiments have shown that smear results in more degradation than does random noise blur, so perhaps it should be included in the calculation.
Random noise blur is calculated as the 1-$\sigma$ residual random movement of the projected pixel direction after linear smear is removed.
The blur MTF is given by
\begin{equation} MTF_{\textrm{blur}} = \exp\left[-2(\pi \sigma_b \xi)^2\right] \end{equation}where where $\xi$ is spatial frequency in cy/rad, and $\sigma_b$ is the 1-$\sigma$ jitter during integration time in rad.
# to calculate smear MTF
def MTFsmear(freq, smearRate, inttime):
# freq in cy/rad
# smearRate in rad/s
# inttime in s
return np.sinc(freq * smearRate * inttime)
# to calculate blur MTF
def MTFblur(freq, rmsblur):
# freq in cy/rad
# rmsblur in rad
return np.exp(-2 * (np.pi * freq * rmsblur) ** 2)
# to plot the smear and blus MTF as function of linear motion and blur.
p = ryplot.Plotter (1 ,1 ,1 , figsize =(8 ,4) )
if 'MTF smear motion' not in dfS.index.values:
dfS = dfS.append(pd.DataFrame(index=['MTF smear motion']))
if 'MTF blur motion' not in dfS.index.values:
dfS = dfS.append(pd.DataFrame(index=['MTF blur motion']))
ffreqmax = np.max(dfS.loc['SpatFreqMax cy/m'])
for colname in dfS.columns.values:
ffreq = dfS.loc['SpatFreq cy/m'][colname]
rmsblur= dfS.loc['RMS jitter rad'][colname]
smearRate= dfS.loc['Linear smear rad/s'][colname]
dfS.loc['MTF smear motion'][colname] = MTFsmear(ffreq, smearRate, inttime)
dfS.loc['MTF blur motion'][colname] = MTFblur(ffreq, rmsblur)
p.plot(1, ffreq/1000., dfS.loc['MTF blur motion'][colname],r'Line of sight motion MTF',
'Frequency cy/mm','MTF',label=[r'{} {} $\mu$m rms blur'.format(colname,rmsblur*1e6)]);
inttime = 0.01
p.plot(1, ffreq/1000., dfS.loc['MTF smear motion'][colname],r'Line of sight motion MTF',
'Frequency cy/mm','MTF',label=[r'{} {} mrad/s smear rate'.format(colname,smearRate*1e3)],
pltaxis =[0,ffreqmax/1000.,-0,1]);
# to plot the MTF components per sensor
p = ryplot.Plotter (1 ,4 ,1 , figsize =(8 ,12) )
if 'MTF effective' not in dfS.index.values:
dfS = dfS.append(pd.DataFrame(index=['MTF effective']))
for i,colname in enumerate(dfS.columns.values):
ffreq = dfS.loc['SpatFreq cy/m'][colname]
rmsblur = dfS.loc['RMS jitter rad'][colname]
smearRate = dfS.loc['Linear smear rad/s'][colname]
ffreqmax = dfS.loc['SpatFreqMax cy/m'][colname]
obscuration = dfS.loc['MTF label'][colname]
dfS.loc['MTF effective'][colname] = dfS.loc['MTF blur motion'][colname] \
* dfS.loc['MTF smear motion'][colname] \
* dfS.loc['MTF detector'][colname] \
* dfS.loc['MTF aberration'][colname] \
* dfS.loc['MTF diffraction'][colname]
p.resetPlotCol()
p.plot(1+i, ffreq/1000., dfS.loc['MTF blur motion'][colname],
r'{}'.format(colname),'Frequency cy/mm','MTF',
label=[r'{} $\mu$m rms blur'.format(rmsblur*1e6)]);
p.plot(1+i, ffreq/1000., dfS.loc['MTF smear motion'][colname],
r'{}'.format(colname),'Frequency cy/mm','MTF',
label=[r'{} mrad/s smear rate'.format(smearRate*1e3)]);
p.plot(1+i, ffreq/1000., np.abs(dfS.loc['MTF detector'][colname]),
r'{}'.format(colname),'Frequency cy/mm','MTF',
label=[r'{} $\mu$m pixel size'.format(pixsize*1e6)]);
p.plot(1+i, ffreq/1000.,dfS.loc['MTF aberration'][colname],
r'{}'.format(colname),'Frequency cy/mm','MTF',
label =[r'{} spot size factor'.format(sratio)]) ;
p.plot(1+i, ffreq/1000.,dfS.loc['MTF diffraction'][colname],
r'{}'.format(colname),'Frequency cy/mm','MTF',
label =[r'{} diffraction'.format(obscuration)], ) ;
p.plot(1+i, ffreq/1000.,dfS.loc['MTF effective'][colname],
r'{}'.format(colname),'Frequency cy/mm','MTF',
label =[r'Effective MTF'], plotCol='k',
pltaxis =[0,ffreqmax/1000.,-0,1]) ;
Methodology¶
The text for this section is take with minor changes from [G].
In the optimization examples the optical aperture $D$ is considered the primary constraint, and thereafter the following process is followed:
Using $D$, the pixel format and the pixel centre-to-centre dimension; the focal lengths $f$ are calculated for a series of FOVs.
For each FOV the f-number is calculated and if it is less than a practical acceptable value, then the f-number is set to a chosen practical value.
For the input source/background temperature the stare time $t_\textrm{int}$ is adjusted until the detector well fill of 50% is achieved or until it reaches a limit determined by the frame time of the imaging sensor. The absolute well fill is calculated as $W_f = 100 \, ( n_{\textrm{scene}} + n_{\textrm{optics}})\,/ \,W_c$ where $W_f$ is the actual well fill percentage, $n_{\textrm{scene}}$ is the scene electrons, $n_{\textrm{optics}}$ is the hot optics electrons, and $W_c$ is total well capacity of the detector integration capacitor.
The sensor and/or system designer is then required to input how much pixel-to-pixel minimal processable modulation (MPM) is required to process individual pixels with the integrity and robustness applicable to the application e.g. 5% modulation is required for video auto-tracking.
All other sensor variables are manually input.
The model is run and the result is plotted against the FOV.
Finally the spatial frequency at which the required pixel-to-pixel MPM is achieved is calculated for each FOV and the model results from above is plotted against the spatial frequency at which the MPM is achieved.
# to set up the run investigating field of view
for i,colname in enumerate(dfS.columns.values):
wn = np.linspace(1e4/dfS.loc['Spectral band Hi'][colname],
1e4/dfS.loc['Spectral band Lo'][colname],300).reshape(-1,)
wl = 1e4 / wn # convert from wavenumber to wavelength
maxFOV = dfS.loc['Maximum HFOV deg'][colname]
nFOVplot = 10
# step 1 calculate focal lengths
fovs = np.linspace(maxFOV/nFOVplot, maxFOV, nFOVplot)
f = dfS.loc['Pixels H'][colname] * dfS.loc['Pixel size'][colname] / \
(2. * np.tan(np.pi * fovs / (2. * 180.)))
fno = f / dfS.loc['Optics diameter'][colname]
# step 2 limit f-number
fno = np.where(fno < dfS.loc['F-number limit'][colname], 4., fno )
dfS.loc['Optics diameter'][colname] = f / fno
# step 3 set stare time
#set integration time to 1.0, normalised - to be scaled later.
#nEcntThermal(wl, tmptr, emis, tauAtmo, tauFilt, tauOpt, quantEff, inttime, pfrac, detarea, fno):
inttime = 1
nEscene = rystare.nElecCntThermalScene(wl*1e6, dfS.loc['Scene temperature'][colname],
1.0, 1.0, 1.0, dfS.loc['Optical transmittance'][colname],
dfS.loc['Quantum efficiency'][colname],
inttime,
(1-dfS.loc['Obscuration ratio'][colname]),
dfS.loc['Pixel size'][colname] ** 2, fno)
nEoptics = rystare.nElecCntThermalScene(wl*1e6, dfS.loc['Optics temperature'][colname],
(1.0-dfS.loc['Optical transmittance'][colname]), 1.0, 1.0, 1.0,
dfS.loc['Quantum efficiency'][colname],
inttime,
(1-dfS.loc['Obscuration ratio'][colname]),
dfS.loc['Pixel size'][colname] ** 2, fno)
print(nEscene)
print(nEoptics)
print(nEscene+nEoptics)
print(' ')
# print(dfS.loc['Optics diameter'][colname])
# print(' ')
[ -2.23684289e+07 -8.94745674e+07 -1.69184995e+08 -1.69184995e+08 -1.69184995e+08 -1.69184995e+08 -1.69184995e+08 -1.69184995e+08 -1.69184995e+08 -1.69184995e+08] [ -7859177.72295223 -31437010.15020002 -59443376.64395394 -59443376.64395394 -59443376.64395394 -59443376.64395394 -59443376.64395394 -59443376.64395394 -59443376.64395394 -59443376.64395394] [ -3.02276066e+07 -1.20911578e+08 -2.28628372e+08 -2.28628372e+08 -2.28628372e+08 -2.28628372e+08 -2.28628372e+08 -2.28628372e+08 -2.28628372e+08 -2.28628372e+08] [ -2.06353795e+07 -8.25423036e+07 -1.85723130e+08 -2.43870263e+08 -2.43870263e+08 -2.43870263e+08 -2.43870263e+08 -2.43870263e+08 -2.43870263e+08 -2.43870263e+08] [ -5158844.86429685 -20635575.89346462 -46430782.41222973 -60967565.7886707 -60967565.7886707 -60967565.7886707 -60967565.7886707 -60967565.7886707 -60967565.7886707 -60967565.7886707 ] [ -2.57942243e+07 -1.03177879e+08 -2.32153912e+08 -3.04837829e+08 -3.04837829e+08 -3.04837829e+08 -3.04837829e+08 -3.04837829e+08 -3.04837829e+08 -3.04837829e+08] [ -5.36896222e+07 -2.14760533e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08] [ -1.36774768e-23 -5.47104282e-23 -5.56333487e-23 -5.56333487e-23 -5.56333487e-23 -5.56333487e-23 -5.56333487e-23 -5.56333487e-23 -5.56333487e-23 -5.56333487e-23] [ -5.36896222e+07 -2.14760533e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08 -2.18383369e+08] [ -3.03167400e+07 -1.21268114e+08 -2.72857586e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08] [ -3.30994940e-24 -1.32399236e-23 -2.97903008e-23 -4.45066790e-23 -4.45066790e-23 -4.45066790e-23 -4.45066790e-23 -4.45066790e-23 -4.45066790e-23 -4.45066790e-23] [ -3.03167400e+07 -1.21268114e+08 -2.72857586e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08 -4.07648955e+08]
Python and module versions, and dates¶
# you only need to do this once
#!pip install --upgrade version_information
%load_ext version_information
%version_information numpy, scipy, matplotlib, pyradi
| Software | Version |
|---|---|
| Python | 2.7.12 64bit [MSC v.1500 64 bit (AMD64)] |
| IPython | 5.1.0 |
| OS | Windows 7 6.1.7601 SP1 |
| numpy | 1.11.1 |
| scipy | 0.18.1 |
| matplotlib | 1.5.3 |
| Mon Nov 28 17:39:48 2016 South Africa Standard Time | |