A) Desiderata for basis and merit functions

There are a number of desirable properties or desiderata we seek for the full digital–optical representation implied by equation (1). (Of course it is unlikely that all can be fulfilled simultaneously; there will likely be tradeoffs.) These properties include:

• Mathematical completeness: The representation should be able to describe any optical system and any linear digital processing operation, and thus describe any linear transformation from external (optical) scene to final digital image.

• Ortho-normality: Elements in the basis function set should be linearly independent and therefore non-redundant. This property is necessary for a representation of a computational imaging system to be unique, and for transformations to be invertible – at least when the system is non-degenerate. (An overcomplete basis may be desirable, however, because it can enable the representation of typical computational systems to be sparse.) The basis elements should be normalized for ease of calculation and for direct interpretation and comparison.

• Interpretability: To the extent possible, the basis functions should be interpretable or at least easily visualized to facilitate understanding and scholarly presentation.

• Computationally efficient: The representation should admit efficient computation just as the computation of coefficients in an expansion, such as the DCT admit very fast implementations through the fast Fourier transform (FFT).

• Balanced complexity: The representation of “typical” or “representative” computational imaging systems should require a small number of terms (basis elements) in each domain, optical, and digital. In an (undesirable) unbalanced joint representation, one might find that a single basis element in one domain (e.g., optical) requires an infinite number of basis functions in the other domain (e.g., digital).

• Low condition number for digital signal processing: All computational imaging systems possess noise, such as Poisson photon shot noise, sensor read noise, analog- to-digital (A/D) conversion noise, and others. The condition number for a system matrix $\kappa \lpar {\bf R}\rpar = \vert \lambda_{max}/\lambda_{min}\vert \ge 1$ – the ratio of the matrix's largest to its smallest eigenvalue – should be small in order to compute accurately the final image.

There are a number of properties or desiderata we seek for end-to-end merit function. (Of course it is unlikely that all can be fulfilled simultaneously; there will likely be tradeoffs.) These properties include:

• Scalar: The merit function should be a scalar (as are RMS wavefront error, RMS point-size, digital image RMS error, etc.), so that system designers can rank order candidate systems, and so that architects of computational imaging design software can craft gradient-descent iterative or other optimization routines.

• Non-negative: The merit function should be non-negative so that when the merit function value is zero, the ideal or “perfect” digital–optical performance is obtained, just as in the limit a vanishing RMS point-spread-function spot size corresponds to a “perfect” optical image and a vanishing RMS pixel error corresponds to a “perfect” digital image.

• Scale and unit independent: The merit function should not depend upon overall physical scale of the optical component or number of pixels in the final digital representation, so as to facilitate comparisons of designs and imaging architectures independent of scale.

• Extendable to alternate image estimation tasks: The merit function should be naturally extendable/adaptable to image-based tasks other than faithfully reproducing a scene, for instance estimating visual motion, creating scene depth maps (appropriate for certain imaging architectures), and others.

We consider the current state of basis and merit function first in the optical domain (Section III.B), then in the two-dimensional digital signal-processing domain (Section III.C), then finally in the joint digital–optical domain (Section III.D).

B) Optical basis and merit functions

Current optical design (for axially symmetric systems) relies heavily on the circular Zernike polynomials, $Z_{n}^{m}\lpar \rho\comma \; \theta\rpar $ , defined as

(2) $$\matrix{ {Z_n^m (\rho ,\theta )} \hfill & { = R_n^m (\rho )\cos (m\theta )\quad ({\rm even}\;{\rm polynomials})} \hfill \cr {Z_n^{ - m} (\rho ,\theta )} \hfill & { = R_n^m (\rho )\sin (m\theta )\quad ({\rm odd}\;{\rm polynomials})} \hfill \cr } $$

where ρ is the normalized radial distance from the optical axis, θ the azimuthal angle, and n the order of the radial polynomial

(3) $$\eqalign{ R_n^m (\rho )\,{=}\,\sum_{k=0}^{n-m/2}\!{(-1)^k (n-k)! \over k! (({n+m / 2})- k)! ({(n-m / 2})- k)! } \rho^{n-2 k}\!{,}}$$

for n − m even, and 0 for n − m odd [Reference Born and Wolf6, p. 523–525]. These circular polynomials obey the ortho-normalization relation

(4) $$\eqalign{ \int\nolimits_{\theta = 0}^{2 \pi} \! \int\nolimits_{\rho = 0}^1 \!\! Z_n^m(\rho, \theta) Z_{n^\prime}^{m^\prime}\!(\rho, \theta)\ \rho\ d\rho d\theta = {\epsilon_m \pi \over 2(n+1)} \delta_{n, n^\prime}\delta_{m,m^\prime},}$$

where the Neumann factor ε _m is 2 if m = 0 and 1 otherwise and $\delta_{\cdot\comma \cdot}$ is the standard discrete Dirac delta or indicator function (Fig. 2). The difference between any wavefront on a disk (e.g., exit pupil) and a spherical reference wavefront can be represented as the linear sum of these basis functions (Fig. 3). The coefficients in such a representation reveal the class of aberration present and hence what optical corrections are needed in order to reduce RMS wavefront error or the RMS point-spread function spot size.

Fig. 2. The 21 lowest-order circular Zernike polynomials, $Z_{n}^{m}\lpar \rho \comma \; \theta\rpar $ (cf. equations (2)–(4), below): odd functions , even functions and circularly symmetric functions .

Fig. 3. Every complex wavefront (on a disk) can be represented as the sum of circular Zernike basis functions. Typically one represents the wavefront error – the difference between the actual wavefront and an ideal, spherical wavefront – as a weighted sum of circular Zernike polynomials.

A number of alternative bases sets for circular Zernike polynomials have been proposed, such as the Stephenson polynomials [Reference Stephenson10], shown in Fig. 4. Some of the benefits of this complete representation are that the mean-power error is naturally expressed as a term in a wavefront expansion. Moreover, if the wavefronts produced by two optical systems are very similar, their corresponding Stephenson coefficients will generally be more similar than will their corresponding Zernike coefficients. For this reason, Stephenson polynomials have found increasing use in searching lens design and optical patent databases, for analyzing optical manufacturing constraints, and so on. Finally, the coefficients in an expansion of the wavefront error in terms of the Stephenson basis are independent of pupil size and this property facilitates comparisons of optical systems at different scales. One drawback is that these basis functions do not relate to traditional aberrations as naturally as the Zernike polynomials do.

Fig. 4. The 15 lowest-order basis functions in the Stephenson basis set. Some of the lowest-order basis functions are equivalent to individual circular Zernike polynomials but all others (an infinite number) differ from Zernike polynomials.

For the last several centuries, the vast majority of optical systems – especially optical imaging systems – have employed circular apertures in large part because the lens and mirror elements are easier to manufacture than other apertures. As such, the circular Zernike polynomials have served the optics community well and their use is justifiably widespread. In the current era of computational imaging, however – where there is greater freedom in design and manufacture of optical components and where digital signal processing is an essential step in the data pipeline – the circular Zernike basis will be less appropriate and will likely present an impediment to progress. A basis set for square or rectangular optical systems, analogous to the circular Zernike polynomials, differs rather markedly from the circular Zernike polynomials (Fig. 5) [Reference Al-Hamdani and Hasan11]. It is unlikely that even this basis will be appropriate for many computational imaging systems. (Image-processing bases on a square or rectangular support are likely to be simple, however, because they may be separable into a product of one-dimensional bases, i.e., $G\lpar x\comma \; y\rpar = g\lpar x\rpar h\lpar y\rpar $ .)

Fig. 5. The first 20 basis functions in an extension of Zernike polynomials to imaging systems with a square aperture. These functions are derived from circular Zernike polynomials by orthonormality defined over a square aperture support.

There are a number of traditional optical merit functions, such as the RMS wavefront error, Strehl ratio, Rayleigh criterion, the size of the point-spread function, or a number of measures of the Modulation Transfer Function or Optical Transfer Function (Fig. 6). As mentioned above, these measures are inappropriate for computational imaging systems for a number of reasons. A small point-spread-function may have zeros in its modulation transfer function thereby losing spatial information that cannot be recovered through digital signal processing. More generally and importantly, such merit functions do not describe the full end-to-end performance of a digital–optical computational imager.

Fig. 6. (Left) An Airy disk with red circle indicating a spot size of half maximum. Traditional optical merit functions favor small spot sizes. (Right) A one- dimensional Modulation Transfer Function (MTF) for a diffraction-limited imaging system, represented as a function of the normalized spatial frequency ν. Aberrated systems have MTFs that are lower than such a physical limit.

C) Signal-processing basis and merit functions

There are a number of basis function sets that have been used in image processing and image compression. The simplest mathematically and computationally – and most widely used – is based on the DCT, a separable product of cosine functions (equation (5)). One can index pixels by their discrete coordinates in orthogonal directions, n ₁ and n ₂, and then the basis functions can be represented as:

(5) $$\matrix{ {G_{k_1,k_2}(n_1,n_2)} \hfill & { = \cos \left[ {\displaystyle{\pi \over {N_2}}\left( {n_2 + \displaystyle{1 \over 2}} \right)k_2} \right]} \hfill \cr {} \hfill & {\quad \times \cos \left[ {\displaystyle{\pi \over {N_1}}\left( {n_1 + \displaystyle{1 \over 2}} \right)k_1} \right],} \hfill \cr } $$

where k ₁ and k ₂ serve as components of a two-dimensional discrete wave vector, and N ₁ and N ₂ denote the integral size of the the image block being transformed. Figure 7 shows the 64 basis functions for a DCT appropriate for representing an 8-by-8-pixel image block. The harmonic nature and separability of the DCT permits this transform to be easily computed using the FFT, either on uni-processors or on digital signal processing hardware.

Fig. 7. The 64 8-by-8-pixel basis functions of the DCT for the case N ₁ = N ₂ = 8 and {k ₁, k ₂} ∈ [0 − 7] × [0 − 7] in equation (5).

A widely recognized limitation of the DCT (and its JPEG instantiation) is the appearance of “blockies,” or discontinuities at the boundaries of the component image blocks [Reference Shen and Kuo12]. A limitation of particular relevance to computational imaging is that the DCT incorporates no priors or other information about the properties of aberrated and non-standard optics. For instance the fact that an optical subsystem may commonly be out of optical focus (and thus produce optical images with certain statistical characteristics) is not reflected in the image-processing basis set. Moreover, the DCT is not radially symmetric, which is likely appropriate when processing images produced by radially symmetric optical subsystems. These limitations make the DCT a very poor candidate for a component in representations for computational imaging systems.

A number of image representations based on wavelets have been used to overcome the image-processing and compression limitations inherent in the DCT, many of which have been based on Daubechies filters [Reference Prasad and Iyengar13]. and exploited in compression methods such as the JPEG 2000 standard. This standard has a number of both lossless and lossy modes, is reversible and admits a number of efficient computational implementations based on the hierarchical discrete wavelet transform (DWT). Another class of image-processing basis functions comprises Gabor wavelets of the form

(6) $$\eqalign{ G({\bf x}; {\bf k}, \phi_0, {\bf \Sigma})\,{=}\,e^{- 2 \pi i({\bf k}^t {\bf x} + \phi_0)} {1 \over 2 \pi \vert {\bf \Sigma} \vert^{-1/2}} e^{- {1 / 2}({\bf x} - {\bf x}_0)^t {\bf \Sigma}^{-1} ({\bf x} - {\bf x}_0)}\!{,} }$$

where ${\bf x} = \lpar x\comma \; y\rpar $ is a two-dimensional position in the sensed image, ${\bf k} = \lpar k_{x}\comma \; k_{y}\rpar $ is a two-dimensional vector describing the sinusoidal carrier and ϕ₀ its phase, ${\bf \Sigma}$ is a positive-definite matrix, t denotes transpose, and ${\bf x}_{0} = \lpar x_{0}\comma \; y_{0}\rpar $ is the geometrical center of the basis function. A full basis set consists of different orientations and phase, typically in quadrature phase for mathematical compactness (Fig. 8).

Fig. 8. A subset of the Gabor wavelet basis functions for Σ proportional to the identity matrix (circularly symmetric), and different orientations described by the vector k and phases described by ϕ₀ = 0 or π/2 in equation (6). Here red and blue denote positive and negative response values. Gabor wavelets have been used in radially symmetric image representations, and thus may serve as guidance for a new class of orientation-selective image-processing bases for joint digital-optical system representations.

Figure 9 shows a hierarchical wavelet representation and synthesis of an image based on Symlets – Daubechies wavelets modified to have spatial symmetry [Reference Swee and Elangovan14]. The display at the left is scale-space, where the upper-left corner displays the small number of coefficients of large basis functions and squares to the lower-right display the large number of coefficients of small basis functions. Squares that lie off the diagonal represent coefficients for spatially oriented basis functions: vertical bases at the upper right and horizontal bases at the lower left.

Fig. 9. A wavelet synthesis of an image based on Symlets – a slight modification of Daubechies wavelets to ensure spatial symmetry. The left figure shows the wavelet coefficients in a scale-space display, where the few coefficients for spatially large basis functions are displayed in the small square at the upper left, and the more numerous coefficients for successively smaller basis functions are displayed in the other squares. Vertical bases to the right, horizontal to the left, joint along the diagonal of the figure.

Another general class of computational optical imagers exploits the technique of compressive sensing [Reference He and Carin15, Reference Takhar, Miller and Pollak16]. Such imagers are typically based on optical subsystem that produces a high-quality image on a sensor array (even if that “image” is a mere point on a single photodetector) filtered by a spatial light modulator (SLM). Each component image corresponds to a projection onto a subspace of the high-dimensional space of all possible images. If the patterns presented in the SLM have certain desirable properties of statistical independence, then a high-resolution image can be computed from these projections, for instance by Tikhonov regularization or related methods. In practice, the SLM patterns have been based on random partial Fourier matrices, Hadamard matrices, Gaussian, Bernoulli, multi-wavelet matrices, chirp matrices, all of which have provably good measurement properties. More complex matrices which also have good statistical independence properties have been generated by bipartite expander graphs [Reference Berinde, Gilbert, Indyk, Karloff and Strauss17]. Representation in such a basis can be viewed as spatial filtering. What if the optical image produced by the lens system is degraded or aberrated? For example, what set of spatial filters, when paired with a lens system with severe coma aberration, can lead to a high-quality digital image? This class of problems at the foundation of computational imaging have not yet been addressed.

D) Joint electro-optical basis and merit functions

As mentioned above, currently optical design methodology and software (such as Zemax, Code V, and OSLO) are dominated by function representations based on Zernike and Seidel polynomials and merit functions based on RMS wavefront error, RMS spot size, Strehl ratio, and measures of the modulation transfer function. Currently image-processing methodology and software (such as Matlab, Open CV, etc.) is dominated by function representations based on cosines and wavelets and merit functions are based on mean-squared error (RMSE) and its variants, such as peak RMS. The discipline of computational imaging, and the software tools that serve it, lack an appropriate end-to-end merit function.

Robinson and Stork [Reference Robinson, Stork, Howard and Koshel4] proposed an end-to-end merit function base on the predicted MSE and used it for iterative joint design of the optics and signal processing in simple computational imagers:

(7) $$MSE_{P} = \displaystyle{1 \over n}{\rm Tr }\left[ {({\bf RH}(\Theta ) - {\bf I}){\bf C}_s{({\bf RH}(\Theta ) - {\bf I})}^t + {\bf R}{\bf C}_n{\bf R}^t} \right],$$

where n is the number of pixels in the final digital image, R is the digital image-processing filter, ${\bf H}\lpar \Theta \rpar $ is the optical system matrix determined by the values of the optical variables Θ (focal lengths, lens separations, etc.), I is the n × n identity matrix, C _s is the point-to-point source correlation matrix, C _n is the noise correlation matrix, and Tr denotes the trace operation. This merit function was sufficient for gradient descent design in optical design parameters such as lens curvatures and spacings and digital-processing parameters such as the values of filter taps. The value of MSE _P given in equation (7) is a limit, which in practice is well approximated the actual MSE. It is not known which end-to-end merit function might have properties superior to that in equation (7), for instance in satisfying the desiderata and constraints listed in Section III.A.

In summary, despite some preliminary explorations and acceptable candidates for end-to-end merit functions, there remains no true joint basis function set or widely accepted end-to-end merit function derived especially for the needs of computational imaging.