2.1.

Radiometry

The imaging chain starts with modeling the radiometry, i.e., the electromagnetic energy that is the essence of creating the image. Imaging is fundamentally the act of detecting and measuring electromagnetic radiation in an attempt to understand a remote object or collections of objects. The electromagnetic spectrum ranges from gamma rays with a wavelength $\lambda$ in the picometers to radio waves with wavelengths of hundreds of kilometers. The shorter wavelength electromagnetic spectrum is detected via photonic interactions with matter and is often visualized as packets of energy called photons, where the longer wavelengths are detected with the aid of antennas and are often visualized as waves. In this paper, we are concerned with the range of the electromagnetic spectrum that includes the visible spectrum ($\lambda =400$ to 800 nm) that can be redirected with refractive or reflective optics and detected with bandgap semiconductors. The visible spectrum is a subset of the optical spectrum, which involves frequencies from ultraviolet up to the long-wave infrared (about $\lambda =10\mathrm{nm}$ to 1 mm).

We are most interested in modeling the spectral radiance from the object of interest, $L_{\text{object}}(x,y,\lambda )$, but this is only one of the contributions to the spectral radiance of the entire scene, $L_{\text{scene}}(x,y,\lambda )$, that enters the camera aperture.³ Modeling the spectral radiance for the scene is very challenging given the complex variations and sources of the light energy that will arrive at the camera (Fig. 4). For example, the spectrum of the light changes dramatically in the imaging chain if the camera moves from sunlight illumination to indoor incandescent light. To further complicate the model for outdoor lighting, the radiometry is heavily dependent on the atmospheric conditions that the light travels though, which is constantly changing. The radiometry from indoor light is also complicated by the variety of light sources and the illumination condition from the light reflecting off walls and other objects. Probably, the best example of varying radiometric conditions for indoor lighting occurs when an image is captured with or without a camera flash (Fig. 5).

Fig. 4

Modeling the spectral radiance that enters the camera is more complicated than simply modeling the spectral radiance from the object of interest.

Fig. 5

The radiometry of the illumination source can dramatically change the appearance of the image, such as incandescent lighting (a) versus a camera flash (b).

Various software programs have been created to assist modeling the complex scene spectral radiance. Atmospheric models, such as moderate resolution atmospheric transmission (MODTRAN), produced by Spectral Sciences Inc., Burlington, Massachusetts and the U.S. Air Force, are usually employed to model the atmospheric propagation of the electromagnetic radiation. The scene can be modeled using a physics-based radiometric scene generation tool that builds up an object from smaller facets, with each facet comprising its own material properties. An example of a radiometric scene generating tool is the digital imaging and remote sensing image generation (DIRSIG) model developed at the Rochester Institute of Technology.

2.2.

Image Formation

The image formation element of the imaging chain starts with the spectral radiance from the scene, $L_{\text{scene}}(x,y,\lambda )$, entering the aperture of the camera. As the light waves propagate through the imaging elements to form the image on the sensor, we must consider the modifications to the spatial distribution of the light energy before we can calculate the resulting photons that strike each detector. If we assume that the optical system behaves as a linear shift-invariant system, then we can model the spatial modifications by understanding the effects on a single point of light, resulting in a point spread function (PSF).^4–6 The resulting spatial modification to the spectral radiance in the imaging plane, $L_{\text{image}}(x,y,\lambda )$, can then be modeled by simply convolving the PSF for the system, ${\mathrm{PSF}}_{\text{system}}$, with scene spectral radiance, i.e.,

It is also very helpful to understand how the imaging system is modifying the contrast of the spatial frequencies in the scene by calculating the Fourier transform of the system PSF to give us the system transfer function. In Fourier space, the convolution operation becomes a multiplication operation, so the transfer function for the entire camera system, $H_{\text{system}}$, is obtained by simply multiplying together the individual transfer functions from each element in the imaging chain that contribute to the system transfer function, i.e.,

The MTF shows us how the transfer function modifies the contrast of each spatial frequency and the PTF shows us how the transfer function modifies the phase of each spatial frequency.

A broader PSF translates into a transfer function that drops off more rapidly in the higher spatial frequencies, resulting in a blurrier image (Fig. 6). Thus, the PSF and its corresponding transfer function have a direct impact on the image quality, so it is very important to understand the elements in the imaging chain that contribute to them. For most camera systems, the biggest contributor to the system PSF is the optical PSF and its corresponding optical transfer function (OTF).

Fig. 6

An imaging system with a lower transfer function will produce an image with more blurring.

2.2.2.

Optical diffraction PSF and OTF

Let us go back to the entrance aperture and follow the spatially coherent section of a spherical field to the focal plane. Depending on the distance to the point source and, to a lesser extent, the size of the entrance aperture, the field at the entrance aperture may appear to be a plane wave, a parabolic wave, or a spherical wave. The optical system is designed to take the diverging wavefront at the entrance aperture and turn it into a converging spherical wave at the exit pupil. This converging spherical wave is called the Gaussian reference sphere. The center of the Gaussian reference sphere is on the focal plane and is called the Gaussian reference point. If the point sources in the object space are mapped to their corresponding Gaussian reference points, then we have what is called the Gaussian reference image. The Gaussian reference image is an imaginary perfect unaberrated image that does not have the blurring effect of the PSF and is the starting point for modeling the optical wavefronts in the imaging chain. Each point in the focal plane can be thought of as the center of a set of converging waves on a Gaussian reference sphere from an exit pupil, but each point may be the result of different views of the exit pupil. This effect can manifest itself in two ways, first as a warping of the Gaussian reference image and second as a spatially varying PSF. If the spatially varying PSF is small enough to be ignored, then the system is said to be a linear shift-invariant system. Another term that is often used is isoplanatic patch, which refers to the region of the image plan for which the spatial variability of the PSF is small enough to be considered constant.

If we calculate the Fresnel number to determine what fidelity of the diffraction theory is required to propagate the field at the exit pupil to the image plane, we would find that the Fresnel number is not much less than one, which is required in order to use the Fraunhofer diffraction theory. The diffraction calculation for a camera requires Fresnel diffraction but because the optical system has created a wavefront that is converging to a point in the image plane, the parabolic phase term in Fresnel diffraction is canceled with the parabolic phase of the converging wavefront. This, fortunately, leaves us with the mathematics of the Fraunhofer diffraction theory. It should be pointed out that our visualization is for converging spherical waterfronts, but our mathematics is for converging parabolic waves, which is a result of the paraxial approximation, i.e., considering the rays that are near the optical axis.

The converging field at the exit pupil has an amplitude that depends on the intensity of the point source in object space and on the complexities of the propagation from the source to the imaging system. Normally, the concept of radiometry is separated from the magnitude of the PSF. The field at the exit pupil is taken to be one such that the PSF from all point sources have the same magnitude. The radiometry of the scene is then applied by scaling the values of the Gaussian reference image. The converging field at the exit pupil also has phase aberrations relative to the perfect Gaussian reference sphere. These phase aberrations are called wavefront error (WFE).

As an example of the Fraunhofer diffraction PSF, let us look at a clear circular aperture of diameter $D$ with no WFE. The Fruanhofer diffraction pattern for light at wavelength $\lambda$ with incident electric field $E_o$ imaged with a lens of focal length $f$ will be

Eq. (6)

$${\mathrm{PSF}}_{\mathrm{diff}\text{-}\mathrm{circ}}(r,\lambda )={\left[\frac{E_o\pi D^2}{2\lambda f}\right]}^2{\left[\frac{2J_1\left(\frac{\pi Dr}{\lambda f}\right)}{\frac{\pi Dr}{\lambda f}}\right]}^2,$$

where $J_1$ is the first-order Bessel function and $r=\sqrt{x^2+y^2}$. This shape is often called an Airy pattern and consists of a pattern of concentric rings as shown in Fig. 7. The first ring of zeros occurs at

Eq. (7)

$$r\approx 1.22\left(\frac{\lambda f}{D}\right)=1.22(\lambda f_{/\#}),$$

where

Eq. (8)

$$f_{/\#}=\frac{f}{D}.$$

Fig. 7

The left column shows shape of the exit aperture and the right column shows the shape of the point spread function (PSF) on the focal plane, often called an Airy pattern. The top row is a cross-section of the functions and the bottom row shows an intensity plot of the functions. The first zero of the PSF occurs when $1.22(\lambda f/D)$.

The Fourier transform of the diffraction PSF gives the same result as autocorrelating of the aperture pattern to give us the diffraction OTF.

Eq. (9)

$$H_{\mathrm{diff}\text{-}\mathrm{circ}}(\rho ,\lambda )=\{\begin{array}{ll}\frac{2}{\pi }[{\cos }^{-1}(\frac{\rho }{{\rho }_c})-\frac{\rho }{{\rho }_c}\sqrt{1-{(\frac{\rho }{{\rho }_c})}^2}]& \text{for}\frac{\rho }{{\rho }_{\mathrm{c}}}\le 1\\ 0& \text{for}\frac{\rho }{{\rho }_{\mathrm{c}}}>1\end{array},$$

where

Eq. (10)

$${\rho }_{\mathrm{c}}=\frac{D}{\lambda f}=\frac{1}{\lambda f_{/\#}},$$

and $\rho =\sqrt{{\xi }_x^2+{\xi }_y^2}$. Note that the diffraction OTF, as shown in Fig. 8, has a cutoff frequency ${\rho }_{\mathrm{c}}$, where the contrast of the higher spatial frequencies are zero; thus, the diffraction from the aperture imposes a limit to the resolution of the optical system that improves with increasing aperture size. The optical cutoff frequency, therefore, imposes a fundamental resolution limit on the entire camera system because the system transfer function will be zero for all spatial frequencies higher than ${\rho }_{\mathrm{c}}$, even when all of the individual transfer functions are multiplied together.

Fig. 8

The optical transfer function for an optical system that is linear in intensity and has an exit pupil as shown in Fig. 7. A system that is linear in intensity is often called an incoherent system. The optical transfer function (OTF) has a cut-off at ${\rho }_{\mathrm{c}}=D/(\lambda f)$.

2.2.3.

WFE and OTF

Although many simple imaging chain models only consider the aperture diffraction effects on the OTF, there are many other factors that will alter the functional form of the OTF. As discussed earlier, aberrations in the imaging elements will cause departures from the ideal spherical wavefronts in the optical system. Accurately modeling the WFE caused by these aberrations can be very complicated, but is essential in order to properly model the image quality of the camera. We will return to this point as the focus of Sec. 3.

The addition of the aberrations will result in a MTF that is equal to or lower than the diffraction MTF, i.e.,

Eq. (11)

$${\mathrm{M}\mathrm{T}\mathrm{F}}_{\text{aberrated}}({\xi }_x,{\xi }_y)\le {\mathrm{M}\mathrm{T}\mathrm{F}}_{\mathrm{diff}}({\xi }_x,{\xi }_y).$$

The WFE caused by aberrations will, therefore, reduce the OTF and degrade the image quality with additional blurring (Fig. 9). It is common to characterize the WFE of an optical system by the root-mean-square (RMS) WFE, ${\mathrm{WFE}}_{\mathrm{RMS}}$, i.e., the statistical variance of the optical path difference over the entire wavefront. When the aberrations are not known, such as early in the design phase of an optical system, a random phase error can be used to model a nominal WFE produced by the aberrations and the optics manufacturing process. The optics degradation is modeled as an optical quality factor (OQF) transfer function that is multiplied with the diffraction transfer function to get the OTF. Hufnagel modeled the OQF transfer function for aberrations and high-frequency surface roughness as⁹

Eq. (12)

$$H_{\mathrm{OQF}\text{-}\text{Hufnagel}}(\rho )=e^{-4{\pi }^2{\mathrm{WFE}}_{\mathrm{RMS}}^2(1-e^{-4\frac{{\rho }^2}{l^2}})},$$

where $l$ is the spatial frequency of the correlation length but can be used as a tuning parameter to fit measured optical data if they are available. It should be pointed out that some of the original work on what is now called the Hufnagel model was done by O’Neill¹⁰ and developed further by Barakat.¹¹ One key idea that is often overlooked is that the term $e^{-4({\rho }^2/l^2)}$ represents a correlation function that may not be a Gaussian.

Fig. 9

The top left plot shows a contour plot of an example wavefront error on a circular aperture optic. The wavefront error is made up of a random set of aberrations with a high spatial frequency (HSF) wavefront error with a correlation length about a tenth the radius of the optic. The plot in the top right shows the resulting modulation transfer function (MTF) when the wavefront error is scaled to 0.05, 0.1, and 0.2 root-mean-square (RMS) wavefront error (WFE) at 670 nm. The dashed line is for zero wavefront error.

Although there is a correlation between ${\mathrm{WFE}}_{\mathrm{RMS}}$ and the reduction of the OTF, different aberrations at the same ${\mathrm{WFE}}_{\mathrm{RMS}}$ can have different effects on the resulting image quality, so it is very important to model the WFE from the actual aberrations if they are known. As an example, let us look at aberrations that cause a low spatial frequency WFE compared to aberrations that cause a high spatial frequency WFE. Figure 10 shows the distinct differences of the low-frequency and high-frequency aberration optics PTF and PSF even though they both have ${\mathrm{WFE}}_{\mathrm{RMS}}=0.2\lambda$. It is also important to note that different weightings between the types of optical aberrations, e.g., spherical aberration and astigmatism, can occur, which will result in the same RMS WFE but result in different optics MTFs and PSFs. This variability as well as the accuracy of the Hufnagel approximation is tighter for a high spatial frequency WFE than it is for a low spatial frequency WFE, as illustrated in Fig. 11. Figure 12 shows the different effects that the aberrations have on image quality even though they have the same ${\mathrm{WFE}}_{\mathrm{RMS}}$. It should be noted that these image examples did not have noise added or any image restoration processes applied to them.

Fig. 10

The left column shows WFE, phase transfer function (PTF), and the PSF for low spatial frequency (LSF) WFE. The right column shows the same plots for high spatial WFE. Both have the same RMS WFE of $0.2\lambda$ at 670 nm. A key observation is that low spatial WFE disrupts the core of the PSF, but the high spatial WFE preserves the PSF core but creates a very broad base. The numbering of the PSFs relates to the PSFs illustrated in Fig. 11

Fig. 11

These MTFs and PSFs are all different even though they all have ${\mathrm{WFE}}_{\mathrm{RMS}}=0.2\lambda$. Note that the variation in possible transfer functions with the same RMS WFE is larger for LSF WFE than for HSF WFE. Note also that the Hufnagel approximation is a reasonable fit to the average of the range of MTFs for both the LSF and the HSF case when the spatial frequency of the correlation length $l$ is modified for each case.

Fig. 12

The top plots show the MTFs for three instates of low and high spatial WFE at two different RMS WFE values. In both cases, the MTF’s for the high spatial WFE drop quickly. The low spatial WFE has more variability between instances with the same RMS WFE. The top row of image simulations is for the high spatial WFE and the bottom row of image simulations is for the low spatial WFE. Notice the middle row of images look hazy, while the bottom row looks blurrier as one would expect by looking at the PSF’s shown in Fig. 10.

2.3.

Motion

The next consideration in the imaging chain is the relative motion between the camera and the scene that will cause additional blurring to the image. General motion transfer functions can be derived from

Fig. 13

Smear and jitter are the most common motion blurs encountered with cameras.

Smear is directional blurring caused by a constant linear motion, such as snapping a long exposure image from the window of a fast-moving car. If we let ${\overrightarrow{x}}_o(t)=vt$, where $v$ is the relative velocity, and evaluate Eq. (13), then we find the transfer function to be

The phase term $e^{-\pi i{d}_{\text{smear}}\xi }$ shifts the centroid of the image by half the smear distance. This phase term is normally ignored because it only represents the shift in the image pixels and does not affect image quality; hence, usually only the smear MTF as a sinc function is used in the modeling. If the smear direction is not along an axis, then a more general smear transfer function can be created by replacing $d_{\text{smear}}\xi$ with ${\overrightarrow{d}}_{\text{smear}}·\overrightarrow{\xi }$. The smear PSF along the $x$-direction is modeled as

Note that convolving the smear PSF with the scene performs a uniform integration of the scene radiance over the local distance $d_{\text{smear}}\xi$.

Jitter is the result of very high-frequency vibrations [$f_{\text{jitter}}\gg (1/t_{\exp })$] of the camera system. Jitter can be derived by visualizing the integration time being broken up into many very short subintegration times $\delta t$ such that the camera motion during the time $\delta t$ can be assumed stationary. These individual points and subintegration times can be placed into Eq. (13), which will return the transfer function that is the average of all of the spatial shifts. However, if the probability distribution is known or assumed for the spatial location of the points, then one can find the expected value of the average of the transfer function. Because of the form of Eq. (13), the expected value of the average of the transfer function is the Fourier transform of the probability distribution. This allows us to visualize the probability distribution as the PSF for jitter. Often jitter is modeled as a Gaussian probability distribution given by

2.4.

Image Sensing

Modeling the sensor requires us to understand the impact that the sensor has on the spatial properties of the image as well as the electromagnetic energy reaching the detector. Let us first consider how a digital sensor alters the spatial properties of the image. The predominant effect in most sensors is the integration of the incident light across the aperture shape of each detector. For a sensor with rectangular detectors of size $d_x$ by $d_y$, this integration will cause a blurring modeled as

Other blurring effects need to be incorporated into the sensor model, such as carrier diffusion, but these will depend on the specific sensor design and architecture and will not be considered here for our simple imaging chain model.

The other predominant effect that needs to be modeled for a digital sensor is the sampling. The simplest sensor model can be given by the blurring of the optical image by the detector PSF, ${\mathrm{PSF}}_{\det }(x,y)$, followed by a sampling at each detector location (Fig. 14). For an $M\times N$ detector array with each detector a distance $p_x$ and $p_y$ from the next, we can model the detector blur and sampling on the incident radiance $L(x,y)$ as

Fig. 14

A simple model for an imaging sensor is a blur operation followed by a sampling operation.

Mathematically, $g_{\text{detector}}(x,y)$ is a distribution function because of the Dirac delta distribution in its definition. The operation of sampling is where our visualization transitions from continuous to discrete. An alternative representation of $g_{\text{detector}}(x,y)$ can be purely discrete, i.e.,

The distances $p_x$ and $p_y$ are called the sampling pitch of the detector and can have a profound effect on the image quality. All spatial frequencies higher than the Nyquist frequency, defined by

Fig. 15

High frequency line patterns appear as lower frequency patterns as the sampling distance between detectors increases.

Let us now consider the light energy that reaches each detector. For a detector with a spectral bandpass of ${\lambda }_{\min }$ to ${\lambda }_{\max }$, the radiant flux (watts) on the detector is given by

The digital imaging sensor model converts the spatially distributed light energy of the image into the digital counts that record the image. The most common digital sensors are composed of individual detectors that generate electrons from the incident light through the photoelectric effect; thus, the average number of photoelectrons generated at the detector during the exposure time $t_{\exp }$ is given by

Unfortunately, the digital count value of the image does not have a one-to-one relationship with the scene radiance, $L_{\text{scene}}(x,y,\lambda )$, due to the blurring caused by the PSF, the sampling, and the noise that is introduced. The noise introduces an undesired randomness to the brightness of the digital count value for each pixel with the primary noise contributors in most digital cameras coming from the photons’ noise, dark noise, and quantization noise. The photons arriving at the sensor do not arrive in a steady stream but in random fluctuations that follow a Poisson distribution, giving rise to the photon noise that increases with the square root of the signal intensity. The dark noise appears as fluctuations in the digital counts even when no signal is present. There are many causes for the dark noise, such as error in the analog-to-digital converter that creates the digital count value, so this value in the imaging chain model is captured from measured data for the sensor being modeled. Finally, we have quantization noise that arises from quantizing the signal electrons into integer digital count values, creating a randomness that follows a uniform distribution. Putting these noise sources together gives us the root-sum-square noise for the sensor electrons as

An interesting detail regarding the photon noise is that the arriving photons arrive following a Poisson distribution but then drive a binomial probability distribution that has a probability of success given by the quantum efficiency. It just so happens that when a Poisson process drives a binomial process, the result is also a Poisson process, so it is, therefore, still valid to use a Poisson distribution for the photoelectrons.

2.5.

Additional Considerations When Putting the Imaging Chain Elements Together

Figure 16 illustrates the fundamental steps of putting together all of the elements of a simple imaging chain model to simulate the image capture portion of the imaging chain. Image processing, such as calibration, is generally applied to the output of the image capture portion of the imaging chain to generate the final image produced by the camera. The image simulations can then be processed to enhance the image quality, displayed, and evaluated by image users to complete the imaging chain for assessing the image quality of the camera design.

Fig. 16

The fundamental steps of the imaging chain model.

It is easy for many developers to oversimplify the imaging chain models, but these simplifications can lead to misleading results when the evaluations are conducted using the image simulations. For example, many simple image chain models will only consider the PSF of the Airy pattern, shown in Fig. 7, in their system. Although it is the optics diffraction that sets the fundamental resolution limit of the system PSF, using the Airy pattern as the system PSF overlooks the blurring effects caused by the other elements of the imaging chain that can significantly impact the resulting image quality. Also, the Airy pattern is the circular aperture diffraction PSF for a single wavelength of light and ignores the integration of the PSF of the spectral bandpass of the imaging system. If the camera integrates the light over the spectral range ${\lambda }_{\min }$ to ${\lambda }_{\max }$, then the radiant flux, ${\varphi }_{\text{image}}(x,y)$, reaching the image plane is a function of $L_{\text{scene}}(x,y,\lambda )$ by

As an illustration, let us consider a gray world where the spectral radiance of the scene, $L_{\text{scene}}(x,y,\lambda )$, is flat across the spectral range ${\lambda }_{\min }$ to ${\lambda }_{\max }$ and the spectral response of the system is flat across the spectrum as well. In this special case, the integration over the spectral bandpass will simply average the spatial scaling of the PSF over the range ${\lambda }_{\min }$ to ${\lambda }_{\max }$. Even if we only look at the diffraction PSF and corresponding OTF for a circular aperture over the visible spectral bandpass ${\lambda }_{\min }=0.4\mu \mathrm{m}$ and ${\lambda }_{\max }=0.8\mu \mathrm{m}$, as shown in Fig. 17, we see differences from the single wavelength model. Note that the side rings of the Airy pattern smooth out after the PSF is integrated over the spectral bandpass.

Fig. 17

Both plots show the different wavelengths as dashed lines and the spectral averaged as the solid line. The left plot shows the PSF and the right shows the OTF. The inset on the left is a zoomed scale of the PSF that appears in the gray region on the main plot.

When building an imaging chain model for a digital imaging system, particular attention also needs to be paid to the relationship between the optics and the digital sensor. We saw earlier that each imposes a fundamental limit on the detail that can be imaged with a digital camera and the image quality is dependent on the relationship between these limiting factors. The diffraction from the optical aperture imposes a limit on the highest spatial frequency that can be captured [Eq. (10)] and likewise the detector sampling imposes a limit from the Nyquist frequency [Eq. (24)]. The ratio of the sampling frequency and the optical cutoff frequency is a design parameter $Q$ defined as¹²

Historically, $Q$ has been defined as the ratio of the sampling frequency to the optical cutoff frequency rather than the Nyquist frequency to the optical cutoff frequency, so the detector and optics resolution limits are equal when ${\xi }_N={\rho }_{\mathrm{c}}$ and $Q=2$. The resolutions of digital cameras that are designed with $Q>2$ are fundamentally limited by the optics diffraction, whereas digital cameras designed with $Q<2$ are fundamentally limited by the detector sampling.

It may seem intuitive that digital cameras should all be designed with $Q=2$, but other factors, such as signal-to-noise ratio and motion blurring, will influence the final image quality and $Q=2$ may not be the best design when all image quality factors are considered (Fig. 18). Most digital cameras will produce brighter, sharper images when $Q<2$, with $Q$ typically in the range $0.5<Q<1.5$. It is important to note that $Q$ only compares the resolution limits between the optics and the detector sampling and does not take into consideration the system transfer function. The impact of the OTF is a significant factor in the system transfer function when $Q\gtrsim 0.1$ (Fig. 19), so for most imaging systems, we must understand the effect of the optical WFE on the system transfer function, which will influence the image quality.

Fig. 18

The value of $Q$ will influence the final system image quality and must be optimized for each digital camera design.

Fig. 19

Comparing the diffraction-limited optics MTF with the detector aperture MTF shows that the detector aperture MTF dominates for $Q\lesssim 0.1$, but the OTF is a significant factor when $Q\gtrsim 0.1$.

3.1.

Problem Motivation

In the previous section, we only considered a very simple form of the OTF that included the diffraction from a clear circular aperture with a generalized WFE that can be modeled as an OQF, although the importance of properly modeling the WFE was illustrated. If we simply extend the aperture model to include a circular central obscuration of diameter $D_{\mathrm{obs}}$, as we would see in a Cassegrain telescope, then the aperture diffraction PSF will now have the form

The diffraction OTF for a circular aperture with a central obscuration is given by

Note that adding the central obscuration to the clear circular aperture has the effect of moving the concentric rings in the PSF closer toward the center while more energy is moved from the central peak to the outer rings (Fig. 20). Also note that the contrast of the mid-spatial frequencies in the diffraction OTF is reduced as the central obscuration increases in size, but the cutoff frequency remains the same.

Fig. 20

(a) The PSF for a system with different obscuration sizes and (b) the OTF. The inset in (a) is a zoomed-in section of what is shown is a zoomed-in section of what is shown in the gray region. The online version is color and accompanied with an animation (Video 1, MOV, 6.4 MB) [URL: http://dx.doi.org/10.1117/1.OE.53.8.083103.1].

Although these models for the diffraction OTF are useful for simple imaging chain models, they are still inadequate for more complex optical designs, especially when we need to consider all of the potential sources for WFEs.

The primary driver for the shape of the PSF is the shape of the exit pupil. Historically and in many elementary texts on optics, the shape of the exit pupil is an unobstructed circle. However, in reality, the exit pupil has obstructions and in more complex systems, like segmented arrays and multitelescope systems, the overall shape of the exit pupil is very complex. A secondary effect on the shape of the PSF is that of the WFE. WFE is the spatial variation about the optimal optical path length over the exit pupil. Again, historically and in many elementary texts on optics, WFE would be modeled with a linear combination of Zernike polynomials, which are an infinite set of orthogonal functions over the unit disk. Obviously, this orthogonality is lost on exit pupils with complex shapes. One option is to find a new set of orthogonal functions for the new shapes.¹³ This works well for a small yet common set of shapes. Another option is to accept the loss of orthogonality and develop a method that allows the usage of any collection of basis functions, including the influence functions of adaptive optics on an arbitrary exit pupil geometry.

The method presented below is a generalized method for tracking WFE in a consistent manner that is applicable to arbitrary pupil geometries. WFE will be visualized as a vector in an infinite-dimensional vector space; however, calculations will be done in a finite-dimensional vector space. This vector space visualization will conveniently allow for the separation of WFE into multiple subspaces that can be aligned with the specifications for a future system. Calculations for RMS WFE, line-of-sight (LOS), and other macroscopic parameters are reduced to fast matrix calculations.

3.2.

OTF Generation

The impulse response function $h(\overrightarrow{x})$, which is a precursor to the PSF, can be calculated by evaluating the Fraunhofer diffraction integral.

Before we go much farther, it is important that we become explicitly clear about our choice of sign conventions. In general, sign conventions and handedness of axes are completely arbitrary, but we must be consistent and when obtaining data from an external source, one must take great care to understand the hidden assumptions with sign conventions.

We start off with the definition of the optical axis that runs along the center line of the optical system in the direction of propagation. This will be the $z$ axis of our coordinate system. In other words, the $z$ axis points from the exit pupil to the focal plane. We will use a right-handed coordinate system, which will dictate the order of the $x$ and $y$ axes in the focal plane, as shown in Fig. 21. The labels $u$ and $v$ have been used for the axes in the exit pupil. The polar angle $\theta$ will be positive when measured from the $x$ axis in the direction the fingers of the right hand point when the thumb of the right hand is pointing in the positive $z$ direction. This is the normal convention for the polar angle $\theta$.

Fig. 21

A right-handed coordinate system is used with the $z$ axis running in the direction of the propagation of light. The diagram shows the location of the PSF in the focal plane as the top of the wavefront tips toward the focal plane. Positive WFE $\varphi >0$ is colored bluish and is defined as moving in the direction of the $z$ axis.

The complex scalar field at the Gaussian reference sphere at the exit pupil is called the pupil function, $P(\overrightarrow{u})$, and is used in Eq. (40) to calculate the impulse response function. The pupil function can be broken into two real functions: the aperture function, $0\le A(\overrightarrow{u})\le 1$, and the WFE function, $\varphi (\overrightarrow{u})\in \mathit{R}$. The aperture function, which is often a binary function but, in principle, can take on values between 0 and 1, represents the relative spatial distribution of the transmission of the exit pupil. The WFE function represents deviations in the wavefront from the ideal Gaussian wavefront.

In this documentation, we have picked $\varphi (\overrightarrow{u})>0$ to be in the direction along positive $z$. This choice forces us to place a minus sign in the exponent of Eq. (40). This choice supports the visualization that the normal vector of the tip-tilt WFE shown in Fig. 21 points toward the center of the PSF.

The Fraunhofer diffraction integral of Eq. (40) has the same form as a Fourier transform integral, provided that we define the forwarded transform as

In practice, one will be creating sampled images of the exit pupil and then using a discrete Fourier transform (DFT) algorithm to calculate the PSF and OTF with the DFT implemented using a fast Fourier transform (FFT).

To avoid aliasing, one must sample the pupil function with sufficient zero padding. For this discussion, we will assume that a square grid is used and we will only talk about the linear scaling along one side. Let $w$ be the physical width of the sampling grid and $R_{\text{in}}$ be the enclosing radius of the pupil function. For a simple round optic, the enclosing radius is the radius of the optic, but for a complex multisegmented system, it is the radius of the smallest circle that can contain the full system. The enclosing radius must satisfy the following constraint:

If the square grid is sampled with an $n\times n$ array, the sample spacing is $w/n$. The left most axis on Fig. 22 shows the sampling of one side of this grid. The pupil function, $P(x_i)=P_i$, has been sampled at

Fig. 22

The left axis is the sampling of the pupil function. The second axis shows the scaling resulting from converting the Fraunhofer diffraction integral into a Fourier transform. The next axis shows the width and sampling of the PSF on the focal plane. The last axis shows the sampling of the OTF in the frequency domain.

An arbitrary choice has been made regarding the sampling that is done from $-(w/2)$ to $(w/2)-(w/n)$. We could have easily sampled from $-(w/2)+(w/n)$ to $w/2$, which would give us consistent results. The current choice is driven by using a computer programming language that references arrays with a zero-based index system. In other words, $\mathtt{A}[\mathtt{0}]$ is the first element of the array and $\mathtt{A}[\mathtt{n}-\mathtt{1}]$ is the last. For an array with an even number of points and length $n$, the center value of the function is at array index $\mathtt{A}[\mathtt{n}/\mathtt{2}]$.

Another requirement for sampling the pupil function is that the smallest structure of the pupil function is adequately sampled. This includes the physical structure in the aperture function and the variations in the WFE function. If the size of the resulting grid gets too large for the computing resources of the day, one can use the techniques of Eikenberry et al. or Ransom et al.^14,15

Let $\mathrm{FFT}(P_{\mathtt{i}})=h_{\mathtt{j}}$ represent the action of a DFT on the $P_{\mathtt{i}}$ data to create the $h_{\mathtt{j}}$ data. One should keep in mind that every data value $h_{\mathtt{j}}$ is a function of all the input data values $P_{\mathtt{i}}$. The PSF data can be calculated from the sampled pupil function, $P_{\mathtt{i}}$, as follows:

When modeling ${\mathrm{OTF}}_{\mathtt{k}}$, it is often convenient to consider the MTF and PTF, which are related as

The magnitude of the OTF, is the MTF which is scaled such that the zero frequency is equal to one, ${\mathrm{MTF}}_{n/2}=1$. The PSF and OTF are a Fourier Transform pair,

For a simple round optic of radius $R$ with WFE, one can expand the phase of the pupil function out into a set of orthogonal functions. A common choice for these orthogonal functions are Zernike polynomials.

Two key benefits of using a set of orthogonal functions to describe WFE are that it is easy to describe the key aberrations with just a few terms and the orthogonality relationship of the basis functions makes it easy to calculate macroscopic parameters with just the coefficients of the basis functions. The orthogonality relationship is as follows:

The coefficients that describe a given WFE, $\varphi (\overrightarrow{u})$, can be calculated by taking advantage of the orthogonality relationship.

Once the coefficients $(a_0,a_1,a_2,,a_n)$ have been calculated, we can easily calculate the following integral:

Therefore, if we know the coefficients $(a_0,a_1,a_2,,a_n)$, calculating the value for $\langle {\varphi }^2\rangle$ is as simple as summing the square of the coefficients. Another very important way of looking at this is if we think of the WFE as a vector made up of the coefficients

The rest of this work describes how to do these types of integrals when the optic is arbitrarily complex and the orthogonality relationship is no longer valid. This is done by careful book keeping of the WFE vector spaces and the calculation of a metric we will call the eta matrix, $\mathit{\eta }$. This matrix will allow us to keep the visualization that WFE is a vector in a vector space and will be used to calculate the new dot product in this vector space.

4.1.

Aperture Nesting Rules

In general, each aperture has an enclosing disk or other enclosing shape that contains the unobscured domain of the basis functions. For this discussion, we will assume that the enclosing shape is a disk with radius of $R_{\overrightarrow{q}}$, where $\overrightarrow{q}$ is the vector-index designator for the given aperture. The enclosing disk is the domain over which the WFE basis functions are defined. If one is using Zernike polynomials, then the basis functions are also orthogonal over the domain of the disk. However, if the domain of the basis functions is not the shape of a disk, then one can think of the radius $R_{\overrightarrow{q}}$ as a scale parameter of the enclosing shape. For example, if the basis functions are square image files, then the enclosing shape is a square and the scaling parameter $R_{\overrightarrow{q}}$ can be the distance from the center to one of the corners.

The domain over which the $\overrightarrow{q}$’th aperture is defined will be referred to as ${\mathcal{D}}_{\overrightarrow{q}}$. This region is shown in Fig. 23 as the gray regions within the supporting circles. The gray regions get darker as there are more nested apertures. The $\overrightarrow{q}$ vector indexing method helps visualize the tree structure of the nested apertures. Figure 24 shows the nesting hierarchy for the general hierarchy of nested apertures shown in Fig. 23. The $\overrightarrow{q}$ vector is as long as the deepest trail on the graph from the root to the leaf. Each level corresponds to a dimension in the $\overrightarrow{q}$ vector. The $\overrightarrow{q}$ vector is the trail from the root to the node in question with zero padding. The union of the domains of the leafs of the aperture tree corresponds to the final pupil mask as shown in the inset of Fig. 23.

Fig. 24

The tree structure shown here is for the aperture system shown in Fig. 23. The union of the domain of the gray boxes; the leaf nodes make up the pupil mask, which is the inset in Fig. 23.

Apertures can be combined into a hierarchal arrangement to aid in modeling a complex system. Any aperture can have children apertures, which can be used to build up a complex hierarchy. The following rules apply to hierarchical aperture arrangements:

4.2.

Calculation of the $\mathit{\eta }$ Matrix

The goal of the $\mathit{\eta }$ matrix is to replace time-consuming numerical integrations of the product of two functions defined over the domain of the pupil with a quick matrix multiplication.

This matrix multiplication can be thought of as the dot product between the two vectors $\overrightarrow{\varphi }$ and $\overrightarrow{\psi }$ in a nonorthogonal space described by the metric $\mathit{\eta }$, where the elements of $\overrightarrow{\varphi }$ and $\overrightarrow{\psi }$ are the expansion coefficients of the functions $\varphi (\overrightarrow{r})$ and $\psi (\overrightarrow{r})$. If we first look at a simple single aperture system where the WFE is expanded on the following basis functions, $B_i(\overrightarrow{x})$, then Eq. (60) can be rewritten as

It is now possible to switch the order of the integration and summations to arrive at

This seemingly trivial algebraic operation has profound computational effects. The integral in Eq. (62) can be calculated for all combinations of $\mathtt{i}$ and $\mathtt{j}$ with just the knowledge of the geometry of the aperture, $\mathcal{D}$. This will become the $\mathit{\eta }$ matrix. The instance of the WFE is contained in the expansion coefficients $a_{\mathtt{i}}$ and $b_{\mathtt{j}}$, which become the WFE vectors $\overrightarrow{\varphi }$ and $\overrightarrow{\psi }$. Basically, we have separated the geometry of the aperture from the instance of the wavefront error.

In the most general case, the $\mathit{\eta }$ matrix is the overlap matrix between every combination of two basis functions for every combination of two apertures over the domain these two apertures make in the exit pupil. The $\mathit{\eta }$ matrix is calculated by evaluating the following integral for each matrix element:

$B_{{\overrightarrow{q}}_a:\mathtt{i}}(\overrightarrow{r}/R_{{\overrightarrow{q}}_a})$ in Eq. (63) is the $i$’th basis function on the ${\overrightarrow{q}}_a$’th aperture that is used to describe the WFE of the optical system.

As a result of there being no restrictions on the basis functions used and the fact that the same type of basis functions can be used on parent apertures along with children apertures, the $\mathit{\eta }$ matrix will not have full rank and there will be redundancy in how a given WFE can be represented on the basis functions. However, along the diagonal, there are blocks that have full rank.

4.3.

Derivation of Support Matrices

We need to derive some supporting matrices to assist in the calculation of the system RMS WFE, global pointing plane (GPP), and LOS. To do this, one needs to calculate vector representations for the simple functions 1, $x$, and $y$.

As a result of using multiple complete sets of basis functions over the same region of space, there is some freedom in how one can expand any given function on the basis functions. When modeling a system, one has the freedom to pick appropriate basis functions. If the basis functions are polynomial, then the lowest-order basis functions are most likely of the form $c_0$, $c_{1}x$, and $c_{2}y$, or some linear combination thereof, where $c_i$ are constants. In this situation, the vector expressions for the functions 1, $x$, and $y$ will evaluate to a finite number of nonzero terms. If the basis functions are not polynomial, then the expansion of the functions 1, $x$, and $y$ may not terminate. Due to this condition, it is suggested that at least three basis functions are polynomial and of the form $c_0$, $c_{1}x$, and $c_{2}y$. The remainder of the basis functions can be anything else that best matches the needs of the modeling the WFE of the system.

There is no loss in generality by requiring three basis functions for each leaf node to be the first three Zernike polynomials. If one chooses to use a different set for the first three basis functions, one is required to derive the final expressions shown in this section.

An arbitrary function $f(\overrightarrow{r})$ can be written in terms of a linear combination of basis functions $B_{{\overrightarrow{q}}_a:\mathtt{i}}(\overrightarrow{r})$ over the domain of the ${\overrightarrow{q}}_a$’th aperture as

When solving for the expansion coefficients $b_{{\overrightarrow{q}}_a:\mathtt{i}}$ by multiplying by the $\mathtt{j}$’th basis function, one can integrate over the full domain, not just the domain as masked by the ${\overrightarrow{q}}_a$ aperture. This, of course, assumes that the function $f(\overrightarrow{r})$ is known outside of the pupil domain and everywhere within the domain of the basis functions.

If Zernike polynomials are used for the basis functions, the expansion coefficients can be calculated with

The point $(x_{{\overrightarrow{q}}_a},y_{{\overrightarrow{q}}_a})$ is the center of the ${\overrightarrow{q}}_a$’th aperture and $R_{{\overrightarrow{q}}_a}$ is the radius of the disk that encloses the aperture. Because the $\mathit{\eta }$ matrix of Eq. (63) contains the information for all the apertures and basis functions, the ${\mathbf{1}}_{{\overrightarrow{q}}_a}$, ${\mathit{X}}_{{\overrightarrow{q}}_a}$, and ${\mathit{Y}}_{{\overrightarrow{q}}_a}$ vectors are combined to form full system vectors, for example,

4.4.

Projecting Out the Mean and the GPP

Given an arbitrary WFE function $\varphi (\overrightarrow{r})$ defined over the union of all of the domains of the leafs of the pupil mask tree, the mean of the WFE can be found by integrating $\varphi (\overrightarrow{r})$ over the valid domain and then dividing by the area of the domain, giving

Our goal is to find a matrix operation that accomplishes this task. By using Eq. (60) and the vector $\mathbf{1}$, we can write down the integration as a matrix operation.

This leads us to two matrices, one that calculates the mean, ${\mathit{K}}_{\mu }$, and one that projects the mean out of the original WFE vector, ${\mathit{P}}_{\mu }$.

In addition to projecting out the mean, another common operation is calculating and projecting out the best-fit plane to the WFE. This best-fit plane is called the GPP. A normal vector to the GPP that intersects the optical axis points to the location of the maximum of the PSF in the image plane. This location will be called the LOS. GPP is the global piston, tip, tilt (PTT) of a system.

We first start with a least-squares fit cost function for the GPP.

We proceed in the normal fashion by setting the partial derivatives to zero and then solving for $\alpha$, $\beta$, and $\gamma$. During the mathematical manipulation, one will come across terms of the following form, which can be rewritten in a matrix form.

It is now possible to write down the matrices for calculating GPP and projecting out the GPP.

Finally, the ${\mathit{P}}_{\mathrm{\Gamma }}$ matrix, which is used to project out the GPP, can be calculated with

The coefficients for the GPP can be calculated by applying the ${\mathit{K}}_{\mathrm{\Gamma }}$ matrix on the WFE vector, $\overrightarrow{\varphi }$.

The WFE about the GPP is obtained by subtracting the GPP from the full WFE function, or in matrix terms, the WFE vector can have the GPP projected out, where we let ${\overrightarrow{\varphi }}_{\mathrm{\Gamma }}$ be the WFE vector with the GPP projected out, i.e.,

4.5.

Generalized Projection Operators

In general, one can project one subspace into another provided that there is some overlap between the subspaces. An excellent example is the case of adaptive optics, where a set of influence functions are used to back out a measured WFE. Let ${\widehat{\overrightarrow{\varphi }}}_1$ of length $n$ be the measured WFE, which is expanded out on the basis functions $B_1(\overrightarrow{x})$. This subspace will be called ${\mathcal{S}}_1$. Let ${\widehat{\overrightarrow{\varphi }}}_2$ of length $m$ be the induced WFE from the influence functions $B_2(\overrightarrow{x})$ in subspace ${\mathcal{S}}_2$. The hat above the vectors is used to denote that these vectors are not the full length that corresponds to the dimensions of the $\mathit{\eta }$ matrix. The influence functions are just another set of basis functions. The full WFE vectors are