1. INTRODUCTION

Computational imaging (CI) is a class of imaging systems that, starting from an imperfect physical measurement and prior knowledge about the class of objects or scenes being imaged, deliver estimates of a specific object or scene presented to the imaging system [1–7]. This is shown schematically in Fig. 1. The physical measurement, typically of light intensity on a digital camera, i.e., sampled on a pixel grid, we will refer to as raw measurement, or raw intensity image. The analytical relationship between the object and raw intensity image is the forward operator, whereas the analytical expression of prior knowledge is the regularizer, or simply prior. The premise on the existence of the regularizer is that, in general, any object that we measure in a natural scene obeys patterns that introduce correlations in the object’s spatial features [8]; these patterns, if discovered, should help to reduce uncertainty in the recovery of the object, i.e., to reduce ill posedness. As we will review in some detail in Section 2, in a CI system, the estimate of the object is generally obtained by solving an optimization problem. The functional to be optimized promotes competition between fitting the raw measurement according to the forward operator, and observing the prior. We refer to the optimization outcome as object estimate, scene estimate, solution to the inverse problem or, simply, image.

 

Fig. 1. General computational imaging system. The illumination source or source array is shaped by the condenser optics according to the operator $H_{\mathrm{i}}$ before reaching the object $f$. The radiation is subsequently shaped again by the imaging optics with collection operator $H_{\mathrm{c}}$, and the intensity $g$ is sampled by the digital camera. The signal $g$ is then processed by the CI algorithm, which takes into account prior knowledge $\mathrm{\Phi }$ about the class of objects being imaged (regularizer) in addition to the physical models $H_{\mathrm{i}}$, $H_{\mathrm{c}}$. The result of the computation is the image $\widehat{f}$. For detailed notation and description, please see Section 2.

Download Full Size | PDF

The motivation for employing computation instead of using the raw intensity image directly is that the latter may not be readily interpretable. The information may be there, but incomplete or hidden; so further processing is needed to complement and reveal it. The CI way of thinking also opens up new directions in optical design because it removes the traditional requirement that the raw image itself must observe some metric of spatial conformity and fidelity with respect to the object. Thus, it is no longer necessary to fully correct all the aberrations before forming the raw image; some of the correction may be assigned to the computational components of the system. Perhaps more importantly, with CI, it becomes possible to create optical systems producing raw images that are deliberately spatially non-conforming to their respective scenes. The intent would then be after computation to form representations that are not readily available in the traditional spatially conforming intensity images: tomography, i.e., reconstructing the interior of a volume from projections; and quantitative phase retrieval, i.e., extracting the phase of the optical field from the raw intensity, are both classical examples of this CI mode of operation.

Machine learning (ML) is a class of function interpolation algorithms that are informed by examples of the function being interpolated [9–13]. ML architectures are very general, e.g., they may be layered, recursive, etc. Later in this review, we will see that it is even possible to include explicit physical models as part of the computation. The specific architecture of interest here is based on the neural network (NN), a multilayered computational geometry. Each layer is composed of simple nonlinear processing units, also referred to as activation units (or elements); and each unit receives its inputs as weighted sums from the previous layer (except the very first layer, whose inputs are the quantities we wish the NN to process.) Until about two decades ago, students were advised to design NNs with up to three layers: the input layer, the hidden layer, and the output layer. Recent progress in ML has demonstrated the superiority of architectures with many more than three layers, referred to as deep NNs (DNNs) [14–17]. Figure 2 is a simplified schematic diagram of the multi-layered DNN architecture.

 

Fig. 2. Simplified schematic of a generic deep neural network. The input layer consists of the components of the input vector $u$. The dark circles denote the activation (nonlinear) elements, whereas $w_{kj}^{(l)}$ is the weight connecting the $j$-th activation element in layer $(l)$ with the $k$-th activation element in layer $(l+1)$. The last $(L)$-th layer produces an estimate of the output vector $v$. The quality of the approximation relating input to output depends on the training; please see Section 3.

Download Full Size | PDF

The process of adjusting the connection weights is known as training and is generally implemented as an optimization routine. In the supervised training mode, which is of most interest here, the weights are adjusted by minimizing the training loss function (TLF), a metric of the difference between the actual DNN output and its desired output. The former is determined by known DNN input–output pairs, the training examples. Thus, already the use of optimization emerges as a broad common ground between CI and ML. However, there is an important difference: in traditional CI, the optimization routine must be computed for each imaging operation, whereas in ML, the optimization is executed only during the training phase of the DNN. Imaging operations using ML, thus, are generally very fast because the trained DNN is a feed-forward computational architecture. Typical computation times may be from a few minutes to a few hours per frame for typical CI optimizations; several hours or days for DNN training; and milliseconds per frame for a trained DNN. Thus, as a rule, the ML approach yields images at rates approaching real time, even in common laptop computer or smartphone processors.

During the past few years, a number of researchers have shown convincingly that the ML formulation is not only computationally efficient, but it also yields high-quality solutions in several CI problems. In this approach, shown in Fig. 3, the raw intensity image is fed into a computational engine specifically incorporating ML components, i.e., multilayered structures as in Fig. 2 and trained from examples—taking the place of the generic computational engine in Fig. 1. CI problems so solved have included lensless imaging, imaging through scatter, bandwidth- or sampling-limited imaging (also referred to as “super-resolution”), and extremely noisy imaging, e.g., under the constraint of very low photon counts. With few early exceptions [18], the ML architectures used in these CI problems have relied on DNNs. The main purpose of this paper is to review these developments, providing a common framework for understanding and improving upon them as well as developing underlying design principles and intuition.

 

Fig. 3. CI architecture of Fig. 1 with an ML engine producing the image. As we see in Section 4, the ML engine generally includes a multilayered architecture such as the one shown in Fig. 2 and is informed on the physics of the illumination and collection optics, $H_{\mathrm{i}}$, $H_{\mathrm{c}}$, respectively, and the prior $\mathrm{\Phi }$. The three components $H_{\mathrm{i}}$, $H_{\mathrm{c}}$, $\mathrm{\Phi }$ or their combinations are incorporated in the ML engine either explicitly as approximants (Section 3.F) or implicitly through training with examples.

Download Full Size | PDF

It is important early on to make a clear distinction between image interpretation versus image formation (or image transformation, as it is referred to in the computer vision research community), and emphasize that the work presented here is exclusive to the latter. Methods for the former typically are agnostic to the physics of imaging, instead assuming that the input data to the ML algorithm are a vector or matrix of pixel values produced by a generic imaging system. These data are then used to arrive at a scalar or low-dimensional vector representation expressing classification or prediction outcomes e.g., identification of objects or faces in the scene, translating written text found in the scene, etc. The images used as input to the ML algorithm are assumed to meet the same spatial conformity and fidelity criteria that a human performing the image interpretation task would have demanded. Traditional image interpretation, in other words, aims to directly replace the human’s function with an algorithm, with the algorithm operating on input similar to what the human’s visual system would (or should) be receiving.

In contrast, the ML algorithms we are concerned with here are cognizant of imaging physics but agnostic to the contents of the object or scene, except to the degree that the contents provide regularizing priors. ML executes a regression operation: given an incomplete representation of the scene in the raw intensity measured by the camera, ML replaces the raw scene with one of higher spatial fidelity. In the process, information such as the phase of the field, which in the measured raw intensity may have been hidden or incomplete, is revealed. The user of the ML-produced images may well be a human demanding a guarantee that the images meet certain standards of fidelity to the true scene; or the images could be fed to another ML algorithm of the more common image interpretation type designed for a similar kind of spatial fidelity as the human would have required.

When the raw intensity data on the camera are not readily usable by a human or a traditional image interpretation algorithm for the purposes of recognition, an interesting dilemma occurs: does it make sense to first form an image of high spatial fidelity and then feed it to the image interpretation algorithm? Or should one design the image interpretation algorithm anew, operating directly on the incomplete raw intensity data? (For the human user, there is of course no choice: spatial fidelity is imperative.) The second approach has often been proven to be productive, e.g., for digital holographic particle localization and tracking [19,20] and diseased cell identification [21,22]; and recognition of handwritten characters through multimode fibers [23]. In this paper, we will not take a position on this dilemma or discuss it any further; yet, we wish to point out that, in our opinion, the dilemma is one of the least explored and most fascinating in the rapidly progressing field of merging ML into CI.

After these introductory remarks, in Section 2, we briefly summarize the principles of CI to establish notation and a common framework. With similar goals, we summarize ML architectures and design principles in Section 3, emphasizing those that have proven to be most useful for CI to-date and the principles for including physical knowledge of the forward operator into the ML computation. In Section 4, we first describe how ML principles apply specifically to CI, and then in Sections 4.A–4.D, respectively, we review the related literature on the following canonical CI problems: super-resolution, i.e., recovery of lost or suppressed spatial frequencies; lensless phase retrieval from intensity; imaging of dark scenes; and imaging through scatter. The use of ML in magnetic resonance imaging (MRI) and x-ray tomography [24–27] has been reviewed elsewhere in detail [28]; and overviews specific to the use of ML in computer vision-related reconstruction problems [29] and quantitative phase imaging [30] have also been published. Concluding remarks and further suggestions for future work are given in Section 5.

2. OVERVIEW OF COMPUTATIONAL IMAGING

A. General Formulation

Referring to Fig. 1, let $f$ denote the object or scene that the imaging system’s user wishes to retrieve. To avoid complications that are beyond the scope of this review, we will assume that even though objects are generally continuous, a discrete representation suffices [31–33]. Therefore, $f$ is a vector or matrix matching the spatial dimension where the object is sampled. Light–object interaction is denoted by the illumination operator $H_{\mathrm{i}}$, whereas the collection operator $H_{\mathrm{c}}$ models propagation through the rest of the optical system.

How $f$, $H_{\mathrm{i}}$, and $H_{\mathrm{c}}$ are constructed depends on the state of the illumination and how light is scattered by the object and optical elements in the system. Possible models for $H_{\mathrm{c}}$ are thin transparency, Born or Rytov expansions, beam propagation, transmission matrix, etc. If multiple scattering is present in the optical system, e.g., a strong diffuser with $f$ denoting the index of refraction distribution, then $H_{\mathrm{c}}$ becomes random and nonlinear. In Section 4.D, we discuss how using ML improves image quality in the latter most challenging case.

The illumination operator $H_{\mathrm{i}}$ is often designed to improve the condition of the inverse problem. Structured and coded illumination have been important for CI since the 1960s for improving resolution and sectioning contrast [34–51], and more recently for phase retrieval [22,52–54]. Another interesting instance is ghost imaging, whose ML implementation we discuss in Section 4.C. Similarly, coded aperture approaches are effective for designing the pupil function in $H_{\mathrm{c}}$ to improve performance [6,55–58].

If the thin object approximation is valid, $H_{\mathrm{i}}f$ is simply a multiplication whose outcome serves as input to the collection optics. For quasi-monochromatic spatially incoherent illumination, $f$ is the modulus-squared of the complex transmittance as a function of lateral position, and $H_{\mathrm{i}}f$ is the intensity immediately after the object plane. For monochromatic spatially coherent light, $f$ is the thin object’s complex transmittance, and $H_{\mathrm{i}}f$ is the optical field. Sections 4.A and 4.B describe the ML implementations of super-resolution and quantitative phase retrieval, respectively, where these assumptions hold.

The output of the collection optics is optical intensity $g$, sampled and digitized at the output (camera) plane. After aggregating the illumination and collection models into the forward operator $H=H_{\mathrm{c}}{H}_{\mathrm{i}}$, the noiseless measurement model is

Uncertainty in the measurements and/or the forward operator is the main challenge in inverse problems. Typically, an optical measurement is subject to signal-dependent Poisson statistics due to the random arrival of signal photons, and additive signal-independent statistics due to thermal electrons in the detector circuitry. Thus, the deterministic model (1) should be replaced by

B. Linear Inverse Problems, Regularization, and Sparsity

For linear forward operators $H$, the image is obtained by minimizing the Tikhonov [3,4] functional

Andrey Tikhonov [3,4] proposed the regularizer

Instead of (4), modern regularizers $\mathrm{\Phi }(·)$ promote sparsity in the object $f$ or its transformed version $Sf$. The linear transform $S$ that converts the object to a sparse representation is called the sparsifier. Sparse or compressive sensing and imaging schemes [61–65] are intuitively justified because we know that most signals, especially images and video, are compressible; therefore, the degrees of freedom expressed by most images of interest are far fewer than the raw number of pixels implies. Moreover, if the raw images are corrupted by noise, $S$ most likely will sparsify the signal but not the noise; as a result, the noise should be penalized and vulnerable to removal by the regularizer.

Compressed sensing is realized with either a fixed or learned set of basis functions. A common example of a fixed set is wavelet decomposition, which is highly redundant [66–70] and, therefore, typically results in a highly sparse representation of the signal. This is also the reason that wavelets are the basis of the JPEG compressed image format [71]. Another popular alternative is the class of nonlinear diffusion and Weikert operators [72–76]. The special case of

Early instances of the idea that sparse representations may instead be learned from examples were motivated by the power spectral density of natural images and analogies with the primate visual cortex [83–87]. This led to the concept of a dictionary operator $D$, ideally chosen such that the transformation

Sparsity-promoting regularizers typically involve non-differentiable metrics, e.g., $L^1$ [94], so their iterative optimization requires the use of proximal gradients [95]. Common algorithms include lasso [96], iterative shrinkage and thresholding algorithm (ISTA) [97] and its variants fast ISTA (FISTA) [98] and TwIST [99], message passing [100], coordinate descent (CoD) [101], alternating direction method of multipliers (ADMM) [102], and Adam [103]. Shortcomings in the iterative procedures specified by these algorithms, ISTA and CoD especially, inspired the idea of training the nonlinearities in a DNN to approximate the sparse regularizer with increased flexibility and computational efficiency [104–106]. We will pick up on this point in Section 3.F.

3. OVERVIEW OF NEURAL NETWORKS

A. Neural Network Fundamentals

The generic computational architecture of a NN is shown in Fig. 2. Its purpose is to learn and implement the map relating the input vector $u$ to the output vector $v$. The standard NN consists of layers, and each layer consists of several simple activation units that are connected to the units in the preceding and succeeding layers. The number of layers $L$ is the depth of the NN. The input vector is fed into the input layer $l=1$, and the output vector emerges from the output layer $l=L$; the counter $l$ is used to enumerate the layers. We will use the notation

The number of units in each layer is the layer’s width. Unit width varies across the layers in various progressions. Classification tasks generally produce representations of much lower dimension than that of the input images; therefore, the width decreases progressively toward the output, following the contracting architecture in Fig. 4(a). Up-sampling tasks, as in the image super-resolution examples that we discuss in Section 4.A, require output dimension larger than the input, so expanding architectures such as Fig. 4(b) may be considered. The concatenation of the two is the encoder–decoder architecture in Fig. 4(c). The unit widths progressively decrease, forming a compressed (encoded) representation of the input near the waist of the structure, and then progressively increase again to produce the final reconstructed (decoded) image. In the encoder–decoder structure, skip connections are also used to transfer information directly between layers of the same width, bypassing the encoded channels. The U-net architecture [107] (see also Section 3.D and Fig. 7) was an early successful convolutional implementation of the encoder–decoder principle, and has been the most common implementation of ML approaches for CI to-date.

 

Fig. 4. Simplified schematics of layer width progression strategies for DNNs. (a) Contracting architecture. (b) Expanding architecture. (c) Encoder–decoder architecture.

Download Full Size | PDF

The layers are implemented as arrays of identical activation units. Each activation unit expresses a simple nonlinear function, typically the same for all units, which we will denote as $r(·)$. With the exception of the first layer, the inputs to the units are weighted sums of the outputs of the units in the previous layer. Let $s_k^{(l)}$ denote the output of the $k$-th unit in the $l$-th layer, for $l=2,\dots ,L$. This is obtained as

The nonlinearity $r(·)$ in modern NNs is most commonly implemented as the function

B. Training and Testing Neural Networks

The power of NNs to perform demanding computational tasks is drawn from the complex connectivity between very simple activation units. The training process determines the connectivity from examples, and can be supervised or unsupervised. The supervised mode has generally been used for CI tasks, though unsupervised training has also been proposed [111–113]. After training, performance is evaluated from test examples that were never presented during training.

The supervised training mode requires examples of inputs $u$ and the corresponding precisely known outputs $\tilde{v}$. In practice, one starts from a database of available examples and splits them to training examples, validation examples, and test examples. The training examples are used to specify the network weights; the validation examples are used to determine when to stop training; and the test examples are never to be used during the training process, only to evaluate it.

Suppose $N$ such pairs of examples $\{u_{(n)},{\tilde{v}}_{(n)}\},n=1,\dots ,N$ are available for training. The TLF is defined as the distance between the desired outputs ${\tilde{v}}_{(n)}$ and the actual outputs $v_{(n)}$ summed over all training examples as

Training specifies the network weights $W$ so as to minimize the TLF. Formally,

In Fig. 2 and Eq. (9), it is evident that the input–output relationship $\mathrm{NN}(·)$ involves the computation of nested (implicit) instances of the nonlinearity $r(·)$, and therefore, implementing the derivative in (15) will require multiple applications of the chain rule. This results in the computational procedure known as back propagation [115] (or back-prop.) To avoid excessive computational cost when dealing with large networks and training datasets, it is common to follow stochastic optimization approaches [116–120] using small, randomly chosen subsets of the dataset in each iteration. The optimizations are implemented in multi-dimensional numerical platforms such as TensorFlow [121] or Theano [122], and standard optimization libraries such as Adam [103] or ML-specific libraries, e.g., Caffe [123], CNTK [124], and Keras [125].

If $M$ pairs of test examples $\{u_{(m)},{\tilde{v}}_{(m)}\},m=1,\dots ,M$ are available, the test error is computed as

Even if the test metric is the same as the training metric, generally the two do not evolve in the same way during training. Recall that test examples are not supposed to be used in any way during training; however, the test error may be monitored and plotted as a function of training epoch $t$, and typically its evolution compared to the training error is as shown in Fig. 5. The reason test error begins to increase after a certain training duration is that overtraining results in overfitting: network function becomes so specific to the training examples that it damages generalization. It is tempting to use the test error evolution to determine the optimum training duration $t_{\mathrm{opt}}$; however, that is not permissible because it would contaminate the integrity of the test examples. This is the reason we set aside the third set of validation examples; their only purpose is to monitor error on them, and stop training just before this validation error starts to increase. Assuming that all three sets of training, test, and validation examples have been drawn so that they are statistically representative of the class of objects of interest, there is a reasonable guarantee that $t_{\mathrm{opt}}$ for validation and test error will be the same.

 

Fig. 5. Behavior of the training, test, and validation error as function of training epoch $t$.

Download Full Size | PDF

The generalization ability of a computational learning architecture depends on many factors, including: the network architecture, i.e., the depth and the width of each layer; the choice of nonlinearity (ReLU generally performs better than sigmoids, as we saw); the quality of examples, i.e., if they were drawn representatively enough from the distribution of interest; and on the accuracy of the numerical algorithm used to implement Eq. (14). The design of NNs involves a certain degree of intuition and art, and is beyond the scope of the present review. Detailed strategies are given in textbooks, e.g., [114,126–132].

C. Weight Regularization

Overtraining and overfitting relate to the complexity of the model being learned vis-à-vis the complexity of the NN. Here, we use the term complexity in the context of degrees of freedom in a computational model [133,134]. For learning models, in particular, model complexity is known as Vapnik–Chervonenkis (VC) dimension [135–138], and it should match the complexity of the computational task. Unfortunately, the VC dimension itself is seldom directly computable except for the simplest ML architectures.

In practice, degrees of freedom in deep learning models are controlled by placing restrictions on the weights during the training process. This strategy is called weight regularization, and it is motivated by arguments similar to those in Section 2.B in the context of inverse problems. For example, weight decay [139] modifies the training functional (14) as

More drastic is weight pruning: degrees of freedom are restricted by setting all weights below a certain threshold to zero, or by eliminating the activation units associated with weak weights [140,141]. Recent improvements are Dropout [142,143] and DropConnect [144], where instead multiple “thinned” versions of the network with shared weights are trained, each with subsets of activations or weights, respectively, randomly removed. The subnetworks are then recombined into a single network for testing. More complicated Bayesian strategies [145] sample models according to their posterior distributions given the training data. Alternatively, generalizability can be improved with a sparsity-promoting TLF, e.g., the minimum absolute error (MAE); see Section 3.E. One way to restrict weights appropriately for many CI circumstances is to exploit invariances in the object class, as in the convolutional model that we discuss in Section 3.D.

D. Convolutional Neural Networks

Certain tasks, such as speech and image processing, are naturally invariant to temporal and spatial shifts, respectively. This may be exploited to regularize the weights through convolutional architectures [146,147]. The convolutional NN (CNN) principle limits the spatial range on the next layer, i.e., the neighborhood where each unit is allowed to influence, and make the weights spatially repeating. Thus, (9) is replaced by

In the encoder–decoder DNN architecture in Fig. 4(c), convolutions may be used to reduce the dimension of signal representation on the encoder side, i.e., to the left of the DNN waist, and increase it on the decoder side. The convolutional encoding principle is shown in Fig. 6. The intermediate outputs ${\zeta }_j^{(l-1)}$ of the $(l-1)$-th layer are first fed into a pooling layer, which combines neighboring $\zeta$'s to reduce the size of the $l$-th layer, typically by a factor of two. To avoid overcomplicating the notation, we will describe only the 1D case here; notation in higher dimensions becomes tensorial without altering the concept. Common choices for the pooling function are

 

Fig. 6. Convolutional encoding principle between layers $(l-1)$ and $l$ of decreasing width. Shift-invariant weights of limited range feed into activation units, e.g., ReLU, followed by a pooling layer.

Download Full Size | PDF

A realistic implementation of an encoder–decoder type convolutional network from the first U-net realization [107] is shown in Fig. 7. The original purpose of U-net was image segmentation, but similar ML architectures have found wide use in CI problems, as we will see in Section 4. The layer width contracts through $3\times 3$ convolutions and pooling (max, in this case) in the earlier encoding stages of the network. Multiple channels of convolutions are added to ensure information is not throttled by the contraction. Starting from input image size $572\times 572$, by the deepest layer, the signal has contracted to $30\times 30$ but with 1024 parallel convolution channels. In the subsequent decoding stages, upsampling is achieved through $2\times 2$ convolutions and concatenating with same-width layers from the contracting path through the skip connections.

 

Fig. 7. Detailed implementation of the encoder–decoder architecture, Fig. 4(c), and the convolutional principle, Fig. 6, in the original U-net that was built for image segmentation (after [107], reprinted with permission from Springer).

Download Full Size | PDF

Subsequently, residual learning [149], or ResNet, was introduced to facilitate learning in very deep NNs, including U-nets. In the example illustrated in Fig. 8, let $s^{(l)}$ denote the vector of inputs (activation pattern) arriving at the $l$-th activation (ReLU) layer, and let $\mathcal{W}$ denote the (nonlinear) map between $s^{(l-1)}$ and $s^{(l+1)}$. Since $s^{(l-1)}$ is added through the unit-weight (dashed line) connection before the $(l+2)$-nd activation layer, $\mathcal{W}$ need learn only the difference $s^{(l+1)}-s^{(l-1)}$ rather than the complete map $s^{(l-1)}\to s^{(l+1)}$. With the residual method, it became possible for networks of depth in excess of $\sim 50$ to be trained and generalize well, whereas previously there had been a tendency in very deep networks to get “stuck” in local TLF minima with large training error.

 

Fig. 8. Simplified schematic of residual learning through bypass connections with unit weight (dashed lines) [149].

Download Full Size | PDF

E. Training Loss Functions

The most obvious TLF choices are the $L^2$ (minimum square error, MSE) and $L^1$ (MAE) metrics. They are defined, respectively, as

These metrics are pixel-wise, i.e., they do not take into account what happens in each pixel’s neighborhood. For images, especially of natural objects, local correlations are prevalent, and neglecting them may lead to unnatural results. Instead, correlative metrics such as the negative Pearson correlation coefficient (NPCC) and the structural similarity image measure (SSIM) take local structure explicitly into account.

Given two arbitrary signals $a$, $b$ of length $P$, let us define the mean

The (N)PCC may be viewed as a simplification of the more general SSIM index [157]:

Sparsity is also promoted by the cross-entropy metric, inspired by information theory–see, e.g., [158] Sections 2.4 and 8.1–8.7, and also known for its utility in solving combinatorial optimization problems [159–162]. Typically, cross-entropy in ML is computed by first binarizing the training ground truth $f$ and reconstruction data $\widehat{f}$ and then penalizing deviations from one of the probability $\chi$ of agreement between $f$ and $\widehat{f}$ according to [107,144,163]

A related metric is the generative adversarial gain (payoff) [165]. In this case, we assume that outputs $v$ and inputs $u$ are drawn from probability distributions $p_{\mathrm{data}}$ and $p_{\mathrm{model}}$, respectively. For CI, these should be interpreted as the prior distribution according to which objects $f$ are generated and the probability of the optical system generating raw intensity images $g$, respectively. Then, two DNNs are used: the generator $G$ and the discriminator $D$. The generator is given inputs $u$ and produces samples $v$, while the discriminator emits a probability $d(v)$ that $v$ belongs to the training set. The two networks are then set to compete against each other as the generator tries to minimize and the discriminator to maximize the discriminator’s payoff:

Unfortunately, it is often the case that while TLFs or their combinations decrease during training, the visual fidelity of reconstructed images becomes worse. This phenomenon is well known in inverse problems literature [5,166,167]. The frustrating disagreement between numerical metrics and image quality perception by human observers is not easy to resolve. Recently, computer vision researchers have started using Amazon’s Mechanical Turk platform for that purpose [168,169].

To make image evaluation content-aware, a perceptual metric [170,171] has been proposed. The principle is shown in Fig. 9. During training, each training object $f$ is accompanied by a content target $f_{\mathrm{c}}$ and/or a style target $f_{\mathrm{s}}$. The training architecture consists of two DNNs: ${\mathrm{DNN}}_1$ for the image reconstruction task, i.e., for producing an image $\widehat{f}$ given the raw measurement $g$; and ${\mathrm{DNN}}_2$, whose job is to compute the perceptual loss function $\mathcal{L}$. ${\mathrm{DNN}}_2$ is structured as 16-layer VGG [172]. The VGG-16 is pre-trained to recognize objects in images, e.g., ImageNet [173]; during the training of ${\mathrm{DNN}}_1$, ${\mathrm{DNN}}_2$ remains fixed.

 

Fig. 9. Training an image-forming ${\mathrm{DNN}}_1$ with a perceptual loss function generated by a content- and style-sensitive ${\mathrm{DNN}}_2$ (after [170]).

Download Full Size | PDF

For an image reconstruction task, such as super-resolution (Section 4.A) the content target to be used in Fig. 9 is the object itself, i.e., $f_{\mathrm{c}}=f$. Let $s^{(l)}(u)$ denote the vector of activation outputs of the $l$-th layer when the input to ${\mathrm{DNN}}_2$ is $u$. The vectors $s^{(l)}(f)$ are encoded representations of objects $f$; indeed, it has been shown that using a regularized optimization routine similar to (3) these encoded representations may be “inverted” to reproduce estimates of their respective inputs $f$ [174]. The key insight [170] is that these encoded representations also capture the perceptual difference between the ideal reconstruction $f$ and an imperfect reconstruction $\widehat{f}$ produced by ${\mathrm{DNN}}_1$. Accordingly, the perceptual feature loss is defined as

It is also common to use mixtures of error metrics for training, i.e., of the form

F. Physics Priors

Unlike abstract classification, e.g., face recognition and customer taste prediction, in CI, the input $g$ and intended output $\widehat{f}\approx f$ of the NN are related by the known physics of the imaging system, i.e., by the operator $H$. Physical knowledge should be useful; how then to best incorporate it into an ML engine for imaging?

One possibility is to not incorporate it at all, as depicted in Fig. 10(a). In this end-to-end engine, the raw intensity measurement $g$ is input directly to the DNN, and the image $\widehat{f}$ is the direct output; in the notation of Section 3.A, Eq. (8):

 

Fig. 10. Ways to implement the ML engine in Fig. 3, with or without physical priors. (a) End-to-end ML engine. (b) Recurrent physics-informed ML engine [176]. (c) Cascaded physics-informed ML engine [176]. (d) Single-pass physics-informed ML engine. Here, $H^{*}=H^T$ if the forward operator $H$ is linear; otherwise, even a crude approximation to the nonlinear inverse has often been sufficient.

Download Full Size | PDF

A definitive approach to incorporate the physical model (forward operator) into the ML engine originated from sparse representations, Section 2.B. Learned ISTA (LISTA), in particular, was an early proposal to incorporate a DNN into dictionary training for a Tikhonov–Wiener functional (3) [104–106,176,177]. Let us denote as $H^T$ the operator conjugate to $H$ and $P_{\mathrm{\Phi }}$ the proximal gradient operator corresponding to $\alpha \mathrm{\Phi }$. Following Mardani et al. [176], the proximal gradient descent approach is expressed as

In practice, instead of the recurrent architecture in Fig. 10(b), it is easier to implement the unfolded, or cascaded ML engine shown in Fig. 10(c). In this version, the measurement $g$ first produces an estimate ${\widehat{f}}^{[0]}$ according to the conjugate operator $H^T$ only. We will refer to ${\widehat{f}}^{[0]}$ as the zeroth-order approximant, or simply the approximant. This approximant result is correct within the conjugate operator’s range. To recover information outside the range, i.e., from the conjugate operator’s null space, the approximant is passed consecutively through $M$ cascaded DNNs. The output at the $m$-th stage is projected onto the null space by the operator $\mathcal{N}$, and the zeroth-order approximant is added [see Eq. (33)] to produce the $m$-th-order approximant. The final image is the $M$-th-order approximant ${\widehat{f}}^{[M]}$. Note also that in this architecture, there is no discriminator; the DNN cascade is trained on an appropriately chosen loss function (see Section 3.E.).

The advantage of the cascaded ML engine is that it is faithful to the original Tikhonov–Wiener formulation. The forward operator is explicitly included while the NNs, in a sense, are replacing the regularizer, and the prior knowledge is now learned from examples rather than explicitly defined. One disadvantage is that training $M$ DNNs increases the complexity of the learning task and the associated risks of undertraining (not enough training examples) or overtraining (too many degrees of freedom).

A compromise is the single-pass ML engine in Fig. 10(d). Here, an approximate inverse operator $H^{*}$ produces the single approximant $f^{[0]}$. The single DNN is trained to receive $f^{[0]}$ as input and produce the image $\widehat{f}$ as output directly, rather than its projection onto the null space. In practice, the single-pass approach has proven to be robust and reliable even for CI problems with high ill posedness or uncertainty, as we will see in Sections 4.A–4.C. Last, it is important to note that $H$ must be linear for Eqs. (32) and (33) to have a fixed point; only then are Figs. 10(b) and 10(c) valid. If this is the case, then the single-pass approximant should be $f^{[0]}=H^{*}g=H^{T}g$. For nonlinear inverse problems, such as phase retrieval and imaging through strong scatter, Sections 4.B and 4.D, respectively, we will see that end-to-end or single-pass engines work well even with crude approximations to $H^{*}$.

4. COMPUTATIONAL IMAGING REALIZATIONS WITH MACHINE LEARNING

The strategy for using ML for computational image formation is broadly described as follows:

Generating the physical database for CI tasks can be hardware intensive and time consuming. Moreover, the quality of training depends on the examples chosen; therefore, care must be taken that they be representative of the class of objects the ML architecture would be used on after training. For highly ill-posed and noisy tasks, it is generally better to restrict the class so that stronger priors are learned. In such cases, it is also a good strategy to adopt an ML engine including the forward operator explicitly, rather than end-to-end (Section 3.F and Fig. 10).

This section is organized according to the type of forward operator the ML algorithm is called upon to invert, and the severity of the ill-posed problem: super-resolution in photography, machine vision, and optical microscopy in Section 4.A; complex amplitude retrieval from intensity in Section 4.B; attenuation and phase retrieval under extremely low photon counts in Section 4.C; and imaging through diffusers in Section 4.D. We also wish to remind the reader of three recent and more focused reviews on tomography [28], computer-vision related reconstruction problems [29], and quantitative phase imaging [30].

A. Super-Resolution

The two-point resolution problem was first posed by Airy [179] and Lord Rayleigh [180]. In modern optical imaging systems, resolution is understood to be limited by mainly two factors: undersampling by the camera, whence super-resolution should be taken to mean upsampling; and blur by the optics or camera motion, in which case super-resolution means deblurring. Both situations or their combination lead to a singular or severely ill-posed inverse problem due to suppression or loss of entire spatial frequency bands; therefore, they have attracted significant research interest, including some of the earliest uses of ML in the CI context.

A comprehensive review of methods for super-resolution in the sense of upsampling, based on a single image, is in [181]. To our knowledge, the first-ever effort to use a DNN in the same context was by Dong et al. [182,183]. The key insight, as with LISTA (Section 3.F), was that dictionary-based sparse representations for upsampling [92,93] could equivalently be learned by DNNs. Both approaches similarly start by extracting compressed feature maps and then expanding these maps to a higher sampling rate. The difference is that sparse coding solvers are iterative; whereas, as we also pointed out in Section 1, with the ML approach, the iterative scheme takes place during training only; the trained ML engine operation is feed-forward and, thus, very fast. To combine super-resolution with motion compensation, a spatio-temporal CNN has been proposed, where, rather than simple images, the inputs are blocks consisting of multiple frames from video [184].

The ML approach to the super-resolution problem also served as motivation and testing ground for the perceptual TLF [170,171] (Section 3.E). The structure of the downsampling kernel was exploited in [177] using the cascaded ML engine architecture in Fig. 10(c) with $M=4$. Figure 11 is a representative result showing the evolution of the image estimates along the ML cascade, as well as their spatial spectra. It is interesting that, by the final stage, the ML engine has succeeded in both suppressing high-frequency artifacts due to undersampling and boosting low frequency components to make the reconstruction appear smooth.

 

Fig. 11. Evolution of the image estimate along the cascaded ML engine in Fig. 10(c) (reprinted with permission from [177], Fig. 10 ). Top row, left to right, and in the notation used in the present review: undersampled raw intensity image $g$, ${\widehat{f}}^{[0]}$, ${\widehat{f}}^{[1]}$, ${\widehat{f}}^{[2]}$, and $\widehat{f}\equiv {\widehat{f}}^{[3]}$. Bottom row: spatial Fourier transform magnitude of the corresponding images from the top row.

Download Full Size | PDF

Turning to inverse problems dominated by blur, early work [185] used a perceptron network with two hidden layers and a sigmoidal activation function to compensate for static blur caused by Gaussian and rectangular kernels, as well as motion blur [186]. Two years later, Sun Jiao et al. [187] showed that a CNN can learn to compensate even when the motion blur kernel across the image is non-uniform. This was accomplished by feeding the CNN with rotated patches containing simple object features, such that the network learned to predict the direction of motion.

In optical microscopy, blur is typically caused by aberrations and diffraction [188]. More than 100 years of research, tracing back to Airy and Rayleigh’s observations, have been oriented toward modifying the optical hardware—in our language, designing the illumination and collection operators—to compensate for the blur and obtain sharp images of objects down to sub-micrometer size. Thorough review of this literature is beyond the present scope; we just point out the culmination of optical super-resolution methods with the 2014 Nobel Prize in Chemistry [189–192]. Stochastic optical reconstruction microscopy (STORM) and fluorescence photoactivation localization microscopy (PALM) for single molecule imaging [193,194] and localization [195] are examples of co-designing the illumination operator $H_{\mathrm{i}}$ and the computational inverse to achieve performance vastly better than an unaided microscope could do.

Computationally, the blur kernel can be compensated for through iterative blind deconvolution [196,197] or learned from examples [198]. A DNN-based solution to the inverse problem was proposed for the first time, to our knowledge, by Rivenson et al. [199] in a wide-field microscope. The approach and results are summarized in Fig. 12. For training, the samples were imaged twice, once with a $40\times 0.95\mathrm{NA}$ objective lens and again with a $100\times 1.4\mathrm{NA}$ objective lens. The training goal was such that with the $40\times 0.95\mathrm{NA}$ raw images as input $g$, the DNN would produce estimates $\widehat{f}$ matching the $100\times 1.4\mathrm{NA}$ images, i.e., the latter were taken to approximate the true objects $f$. The number of pixels in the high-resolution images was ${(2.5)}^2\times$ the number of pixels in the low-resolution representation. Of course, the low-resolution images were also subject to stronger blur due to the lower-NA objective lens. Therefore, the inverse algorithm had to perform both upsampling and deblurring in this case. The ML engine was of the end-to-end type, as in Fig. 10(a), implemented as convolutional DNN with pyramidal progression for upsampling. The TLF was a mixture of the MSE metric (23) and a TV-like ${\partial }_{\mathrm{TV}}^2$ [Eq. (6)] penalty. Since then, ML has been shown to improve the resolution of fluorescence microscopy [200], as well as single-molecule STORM imaging [201] and 3D localization [202].

 

Fig. 12. Deep learning microscopy (adapted from [199], Figs. 1 and 5, with permission). (a) Training of the end-to-end ML engine with low-resolution (downsampled and blurred) inputs to the DNN produced by a $40\times 0.95\mathrm{NA}$ objective lens and high-resolution outputs produced by a $100\times 1.4\mathrm{NA}$ objective lens. (b) Operation of the trained DNN engine, with the test outputs successfully upsampled and deblurred. (c) Resolution test demonstrating the DNN with inputs from the $40\times 0.95\mathrm{NA}$ objective lens achieving in testing performance comparable to the raw performance of the $100\times 1.4\mathrm{NA}$ objective lens the DNN was trained with.

Download Full Size | PDF

B. Quantitative Phase Retrieval and Lensless Imaging

The forward operator relating the complex amplitude of an object to the raw intensity image at the exit plane of an optical system is nonlinear. Classical iterative solutions are the Gerchberg–Saxton algorithm [203,204]; the input–output algorithm, originally proposed by Fienup [205] and subsequent variants [206–208]; and the gradient descent [209] or its variants, steepest descent and conjugate gradient [210]. This inverse problem has attracted considerable attention because of its importance in retrieving the shape or optical density of transparent samples with visible light [211,212] and x rays [213,214].

In the case of weak scattering, the problem may be linearized through a quasi-hydrodynamic approximation leading to the transport of intensity equation (TIE) formulation [215,216]. Alternatively, if a reference beam is provided in the optical system, the measurement may be interpreted as a digital hologram [217], and the object may be reconstructed by a computational back-propagation algorithm [218,219] (not to be confused with the back-propagation algorithm for NN training, Section 3.B.) Ptychography captures measurements effectively in the phase (Wigner) space, where the problem is linearized, by modulating the illumination with a quadratic phase and structuring it so that it is confined and scanned in either space [220–224] or angle [225–227]. Due to the difficulty of the phase retrieval inverse problem, compressive priors have often been used to regularize it in its various linear forms, including digital holography [228,229], TIE [82,230], and Wigner deconvolution ptychography [231,232].

When the linearization assumptions do not apply or regularization priors are not explicitly available, an ML engine may instead be applied directly on the nonlinear inverse problem. To our knowledge, this investigation was first attempted by Sinha et al. with binary pure phase objects [233], and subsequently with multi-level pure phase objects [234]. Representative results are shown in Fig. 13. The phase objects were displayed on a reflective SLM, and the light propagated in free space until intensity sampling by the camera. The ML engine of the end-to-end type [Fig. 10(a)] was of the convolutional DNN type with residuals. Training was carried out by drawing objects from standard databases, Faces-LFW, and ImageNet, converting each object’s grayscale intensity to a phase signal in the range $[0,\pi )$, and then displaying that signal on the SLM. Because of the relatively large range of phase modulation, linearizing assumptions would have been invalid in this arrangement.

 

Fig. 13. Phase extraction neural network (PhENN) (reprinted from [234], Fig. 5, with permission). Columns I and II are the ground truth pixel values driving the SLM and corresponding phase images $f$ produced by the SLM, according to calibration. Groups of columns III-V, VI-VIII, and IX-XI are different propagation distances $z=37.5\mathrm{cm}$, 67.5 cm, and 97.5 cm, respectively. Columns III, VI, IX are the raw images $g$ captured by the camera. Columns IV, VII, X are reconstructions $\widehat{f}$ when PhENN was trained with Faces-LFW [235], whereas columns V, VIII, XI are estimates $\widehat{f}$ when PhENN was trained with ImageNet [173]. The rows represent different databases the PhENN was tested on as a, Faces-LFW (disjoint from training); b, ImageNet (disjoint from training); c, Characters [236]; d, MNIST [237]; e, Faces-ATT [238]; f, CIFAR [239].

Download Full Size | PDF

Retrieval of the complex amplitude, i.e., of both the magnitude and phase, of biological samples using ML in the digital holography (DH) arrangement was reported by Rivenson et al. [240]; see Fig. 14. The samples used in the experiments were from breast tissue, Papanicolaou (Pap) smears, and blood smears. In this case, the ML engine used a single-pass physics-informed preprocessor, as in Fig. 10(d), with the approximant $H^{*}$ implemented as the (optical) backpropagation algorithm. The DNN was of the convolutional type. Training was carried out using up to eight holograms to produce accurate estimates of the samples’ phase profiles. After training, the ML engine was able, with a single hologram input, to match imaging quality, in terms of SSIM (Section 3.E) of traditional algorithms that would have required two to three times as many holograms, and was faster as well by a factor of three to four times.

 

Fig. 14. Retrieval of phase and amplitude with a physics-informed ML engine (reprinted with permission from Springer [240], Fig. 2 ). (a)–(h) Pap smear and (j)–(p) breast tissue section reconstructions. (a), (i) Zoom-ins to the (optical) backpropagation results from a single intensity image; (h), (p) corresponding bright-field microscopy images, shown for comparison. (b), (c), (j), (k) (Optical) backpropagation results showing artifacts due to twin image and self-interference effects; (d), (e), (l), (m) corresponding ML engine reconstructions from a single hologram, showing quality comparable to (f), (g), (n), (o) traditional reconstructions from eight holograms.

Download Full Size | PDF

The 3D image formation properties of phase retrieval have also been investigated as depth prediction with robust automatic focusing [241,242], image reconstruction with a single-shot in-line hologram [243], extended depth of field [244,245], and transparent 3D sample reconstruction from diffraction images obtained at multiple angles (phase tomography) [246]. In the results shown in Fig. 15 from some early work [244], the single-pass physics-assisted ML engine was trained with physical objects at various depths, and out-performed both single-height and multi-height (optical) back propagation in terms of depth invariance. Resolution improvements in digital holographic microscopy were obtained through aberration correction and background rejection [247] instead of the compressive approach in [79]. Adversarial payoff [165] was recently used in holographic microscopy with resolution below the nominal pixel size [248] and ptychography [249] using SSIM and MAE, respectively, as TLFs.

 

Fig. 15. Physics-informed DNN trained with samples at different propagation distances of up to $\pm 100\mathrm{\mu m}$ from the focal plane, exhibiting improved depth invariance compared to (optical) backpropagation and multi-height phase retrieval (MH-PR) (reprinted with permission [244], Fig. 3). The results are from a human breast tissue sample, captured at different propagation distances $\mathrm{d}z$. Distances marked in red exceed the training range of the DNN. $\text{Scale bar}=20\mathrm{\mu m}$.

Download Full Size | PDF

In the same context as strongly scattering phase objects, ML was used early on in an intriguing alternative to our general architecture in Fig. 3 where the object’s 3D refractive index distribution, instead of output by the ML engine, is stored as the weights in a multilayered deep network [250,251]. This is possible because of the formal analogy between the beam-propagation method [252] and the NN back-propagation weight training method, and is implemented by training the DNN to match the illumination fields from multiple angles as input and the corresponding raw intensity images as outputs. This technique recently inspired a similar solution to the simpler problem of a thin phase object with phase delays stored in a single layer and raw intensity measurements obtained according to the Fourier ptychography scheme [253].

The problem of lateral resolution in phase retrieval algorithms is interesting, because performance depends on the optics as well as the inverse algorithm. In regularized algorithms, e.g., Gerchberg–Saxton, gradient descent, or compressed holography, the effectiveness of the computation depends on how faithful the prior is to the relevant class of phase objects. Similarly, in ML-based phase retrieval, one would expect qualities of the example database, such as the spatial power spectral density (PSD), to have a similar effect. In fact, due to the nonlinearities of the training and reconstruction processes, spatial frequencies exhibiting higher prevalence in the training database may become so dominant that they inhibit successful reconstruction of less prevalent spatial frequencies. Invariably, lower spatial frequencies are more popular, especially in datasets of natural images due to the well-known inverse-square PSD law [254]. The result is that the ML-reconstructed phase estimates are strongly low-pass filtered [156]. One way around that problem is to violate the training prior, e.g., by spectrally pre-modulating the examples to flatten the PSD. This results in resolution improvement, at the expense of some edge-enhancement artifacts, as shown in Fig. 16. Alternatively, the low- and high-frequency bands can be processed by two separate NNs [255] and recombined by a third NN, all three networks trained separately [256]. The idea can be thought of as a computational analogue to the HiLo super-resolution method [46,47,49,50], and in [256], it was applied to both super-resolution and phase retrieval.

 

Fig. 16. Resolution improvement by spectral premodulation of PhENN (from [156], Figs. 5, 8, and 9; reprinted with permission). (a) Two-phase point-object resolution target, implemented on the SLM; (b) PhENN image, not resolving phase point-objects separated by three pixels; (c) with spectral pre-modulation, the three pixel-separated phase-point objects are resolved; (d) sample test image from the ImageNet database; (e) PheNN reconstruction; (f) PhENN reconstruction with spectral pre-modulation, showing sharper features but also edge-enhancement artifacts.

Download Full Size | PDF

C. Imaging of Dark Scenes

The challenges associated with super-resolution and phase retrieval become much exacerbated when the photon budget is tight or other sources of noise are strong. This is because deconvolutions, in general, tend to amplify noise artifacts [5]. In standard photography, histogram equalization and gamma correction are automatically applied by modern high-end digital cameras and even in smartphones; however, “grainy” images and color distortion still occur. In more challenging situations, a variety of more sophisticated denoising algorithms utilizing compressed sensing and local feature representations have been investigated and benchmarked [257–262]. What these algorithms exploit, with varying success, is that natural images are characterized by the prior of strong correlation structure, which should persist even under noise fluctuations that much exceed the signal. Understood in this sense, ML presents itself as an attractive option to learn the correlation structures and then recover high-resolution content from the noisy raw images.

The first use of a CNN for monochrome Poisson denoising, to our knowledge, was by Remez et al. [263]. More recently, a convolutional network of the U-net type was trained to operate on all three color channels under illumination and exposure conditions that, to the naked eye, make the raw images appear entirely dark while histogram- and gamma-corrected reconstructions are severely color distorted [169]; see Fig. 17. The authors created a see-in-the-dark (SID) dataset of short-exposure images, coupled with their respective long-exposure images, for training and testing; and used Amazon’s Mechanical Turk platform for perceptual image evaluation by humans [168]. They also report that, unlike other related works, neither skip connections in U-net nor generative adversarial training led to any improvement in their reconstructions.

 

Fig. 17. ML applied to reconstruct a severely photon-limited scene (after [169], Fig. 5; reprinted with permission). (a) Signal captured by a Fujifilm X-T2 camera with ISO 800, aperture f/7.1, and exposure of 1/30 s, with illuminance of approximately 1 lux; (b) reconstruction using denoising, deblurring, and enhancement, which the authors refer to as “traditional image processing pipeline”; (d) reconstruction obtained using the radiational image processing pipeline on a Sony $\alpha 7\mathrm{SII}$ camera; (e) reconstruction obtained using BM3D [264], considered as benchmark for denoising; (c), (f) ML engine reconstructions of the respective dark images.

Download Full Size | PDF

ML has also been shown to be effective for phase retrieval from intensity data with extremely low photon counts [265]. The recurrent and cascaded engines in Figs. 10(b) and 10(c) are not, strictly speaking, applicable because this inverse problem is nonlinear. An end-to-end DNN, as in Fig. 10(a), yielded acceptable results only with moderate photon flux. For lower flux, the single-pass engine in Fig. 10(d) was used. The approximant $H^{*}$ was obtained by running the standard Gerchberg–Saxton algorithm up to a single iteration. With noisy data, as expected, the estimates ${\widehat{f}}^{[0]}$ produced by $H^{*}$ were rather poor, as seen in Figs. 18(c) and 18(g); yet with this as input, the DNN was capable of significantly improving its output; compare Figs. 18(d), 18(h) and 18(e), 18(i), respectively. It is also worth noting that running the Gerchberg–Saxton algorithm to more iterations produced discernible improvements neither to its own reconstructions nor as an alternative approximant to the DNN.

 

Fig. 18. ML applied to quantitative phase retrieval on a severely photon-limited signal (after [265], Fig. 2; reprinted with permission). (a) Ground truth $f$ for a test example; (b), (f) raw intensity signal at photon flux of 1,050 photons/pixel; (c), (g) corresponding outputs ${\widehat{f}}^{[0]}$of the approximant $H^{*}$ of Fig. 10(d), implemented as a single-iteration Gerchberg–Saxton algorithm; (d), (h) corresponding reconstructions by the end-to-end ML engine of Fig. 10(a); (e), (i) corresponding reconstructions by physics-informed single-pass ML engine according to Fig. 10(d). In (a) and reconstructions (c)–(e, (g)–(i), grayscale tone represents phase delay $[\mathrm{0,2}\pi )$.

Download Full Size | PDF

A radically different approach to imaging under extremely low light flux is computational ghost imaging. Interferometric ghost imaging originated in quantum optics [266], but it was soon discovered theoretically [267,268] and experimentally [269,270] that it has a classical analogue with pseudo-thermal light. The raw intensity image is obtained not by a camera but by a spatially integrating (bucket) detector on one arm of the interferometer, while an SLM displays random phase patterns, and a pinhole detector is scanned laterally on the other arm. The image is obtained by correlating the two measurements and removing a background term. Following a dispute about the explanation of pseudo-thermal ghost imaging [271,272], it was shown that in the classical case, one arm of the interferometer can be removed and replaced by computation [273]. Because the bucket detector collects all available photons, ghost imaging methods are attractive for applications where the available illumination is extremely weak, much less than a few photons per pixel [274] and in harsh environments [275]. However, to reconstruct images with a good signal-to-noise-ratio (SNR), usually a long acquisition time is needed [276]. Therefore, the sampling ratio

Lyu et al. [279] used deep learning with the single-pass physics-informed engine [Fig. 10(d)] and approximant $H^{*}$computed according to the original computational ghost imaging [273]. Due to the low sampling rate and the noisy nature of the raw measurements, the approximant reconstructions ${\widehat{f}}^{[0]}$ were corrupted and unrecognizable. However, when these ${\widehat{f}}^{[0]}$were used as input to the DNN, high-quality final estimates $\widehat{f}$ were obtained even with sampling rates $\beta$ as low as 5%, as shown in Fig. 19.

 

Fig. 19. Experimental computational ghost image reconstructions with sampling ratio $\beta =0.1$ (left-hand side group) and $\beta =0.05$ (right-hand side group) (after [279]; reprinted with permission from Springer Nature). Top row: ground truth; second-row: basic ghost reconstruction according to [273]; third row: deep learning ghost reconstructions using the images in the second row as approximants ${\widehat{f}}^{[0]}$ to the physics-informed ML engine of Fig. 10(d); last row: compressive ghost reconstructions according to [278].

Download Full Size | PDF

D. Imaging in the Presence of Strong Scattering

Imaging through diffuse media [280,281] is a classical challenging inverse problem with significant practical applications ranging from non-invasive medical imaging through tissue to autonomous navigation of vehicles in foggy conditions. The noisy statistical inverse model formulation (2) must now be reinterpreted with the forward operator $H$ itself becoming random. When $f$ is the index of refraction of the strongly scattering medium itself, then $H$ is also nonlinear. Not surprisingly, this topic has attracted considerable attention in the literature, with most attempts generally belonging to one of two categories. The first is to characterize the diffuse medium $H$, assuming it is accessible and static, through (incomplete) measurement of the operator $H$, which in this context is referred to as transmission matrix [282–284]. The alternative is to characterize statistical similarities between moments of $H$. The second-order moment, or speckle correlations, are known as the memory effect. The idea originated in the context of electron propagation in disordered conductors [285] and of course is also valid for the analogous problem of optical disordered media [286–290].

The first-ever, to our knowledge, use of a non-dictionary (Section 2.B) ML algorithm for imaging through diffusers was by Horisaki et al. [18]. The authors created an experimental configuration explicitly non-suited to the transmission matrix formulation by sandwiching two SLMs, both displaying the same image, between three acrylic diffuser plates. The optical system was illuminated coherently, and it was lensless, i.e., without imaging optics on the beam path. Examples for training were drawn from the Caltech Computer Vision Database [291]. The ML engine was not implemented as a DNN, but rather as an extension of support vector machines [292] known as support vector regression (SVR) [293].

In SVR, the image is expressed in linear form as

 

Fig. 20. Learning-based imaging through diffusers using support vector regression (SVR) (after [18], Fig. 3; reprinted with permission). (a) Training examples from a database of faces and (b) their corresponding speckle patterns; (c) test face examples and their corresponding (d) speckle patterns, (e) test reconstructions using SVR, and (f) test reconstructions using pattern matching.

Download Full Size | PDF

 

Fig. 21. Generalization study of the SVR method (after [18], Fig. 5; reprinted with permission). (a) Test non-face examples and their corresponding (b) speckle patterns, (c) test reconstructions using SVR, and (d) test reconstructions using pattern matching.

Download Full Size | PDF

Deep learning solutions to the problem were first presented in [299] and [155], using end-to-end fully connected and residual-convolutional (CNN) architectures, respectively. Results are shown in Figs. 22 and 23. The fully connected solution [299] is motivated by the physical fact that when light propagates through a strongly scattering medium, every object pixel influences every raw image pixel in shift non-invariant fashion. However, the large number of connections creates risks of undertraining and overfitting, and limits the space-bandwidth product (SBP) of the reconstructions due to limited computational resources. On the other hand, the CNN trained with NPCC loss function [155,300], despite being designed for situations when limited range of influence and shift invariance constraints are valid, Section 3.D, does a surprisingly good job at learning shift variance—through the ReLU nonlinearities and pooling operations, presumably—and achieves larger SBP. Both methods work well with spatially sparse datasets, e.g., handwritten numerical digits, and Latin and Chinese characters. Compared to Horisaki et al. [18], the deep architectures perform comparably well with spatially dense datasets of restricted content, e.g., faces, and also hallucinate when tested outside their learned priors.

 

Fig. 22. Learning-based imaging through optically thick diffusers using a fully connected DNN with MAE TLF (after [299], Fig. 3; reprinted with permission). First row: speckle patterns input to the DNN; second row: corresponding ground-truth images selected among a database of English letters; third row: corresponding reconstructions by the DNN.

Download Full Size | PDF

 

Fig. 23. Learning-based imaging through diffusers using IDiffNet: a residual-convolutional DNN with NPCC TLF (after [155], Fig. 6; reprinted with permission). Columns I, II; ground truth pixel values and fields modulated by the SLM, after calibration; III-VI: results with 600-grit diffuser. III are the raw images, IV-VI test reconstructions from IDiffNet trained with Faces-LFW [235], ImageNet [173], and MNIST [237] datasets, respectively. Rows correspond to the dataset the test image is drawn from, as (a) Faces-LFW, (b) ImageNet, (c) Characters [236], (d) MNIST, (e) Faces-ATT [238], (f) CIFAR [239].

Download Full Size | PDF

Approaches [18,155,299,300] are all meant to be used with the specific diffuser they were trained with. It was recently shown that for spatially sparse datasets, it is possible to train a single DNN to image through arbitrary diffusers belonging to a class of similar statistics, by adopting a cross-entropy TLF (27), Section 3.E, and with the training procedure including samples of the same training object imaged with randomly chosen diffusers from among the class [163]. It has also been shown that imaging through multi-core fibers can be relatively tolerant to variability in the transmission matrix, in terms of both image reconstruction quality and recognition of simple patterns, e.g., characters from the MNIST database [23].

Non-line-of-sight (NLOS) imaging, recognition, and tracking belong to a related class of problems, because capturing details about objects in such cases must rely on scattering, typically of light pulses [301–309] or spatially incoherent light [310–313]. Convolutional DNNs have been found to be useful for improving gesture classification [314], and person identification and three-dimensional localization [315]; in the latter case even with a single-photon, single-pixel detector only.

5. CONCLUDING REMARKS

The diverse collection of ML flavors adopted and problems tackled by the CI community in a relatively brief time period, mostly since $\sim 2010$ [104], indicate that the basic idea of doing at least partially the job of Tikhonov–Wiener optimization by DNN holds much promise. A significant increase in the rate of related publications is evident—we had trouble keeping up while crafting the present review—and is likely to accelerate, at least in the near future. As we saw in Section 4, in many cases, ML algorithms have been discovered to offer new insights or substantial performance improvements on previous CI approaches, mostly compressive sensing based, whereas in other cases, particular challenges associated with acute CI problems have prompted innovations in ML architectures themselves. This productive interplay is likely to benefit both disciplines in the long run, especially because of the strong connection they share through optimization theory and practice.

Beyond these bidirectional improvements, we wish to single out and return to the open question of image interpretation versus image formation that we touched upon in Section 1. It is possible that image-forming instruments aided by ML will continue doing just that, producing high-quality images for consumption by humans and image interpretation algorithms that require high-quality images as their input. It is also possible that, to some degree, image formation and interpretation will fuse, resulting in algorithms that interpret the raw data directly, without bothering to form high-quality images; or other algorithms might turn out to be useful by falling somewhere in between these two extremes. We saw already that DH and imaging through multimode fiber bundles have emerged as fertile ground for this kind of cross-over [19–23]. In our view, further developments in that direction will carry significant fundamental value, as existing algorithms and insights, originating as they have from either side of formation or interpretation, will likely have to be at least partially reconsidered. Practical implications are also to be expected if pattern recognition and image interpretation, rather than relying mostly on human-recognizable optical intensity patterns as it does presently, become more directly linked to the complete electromagnetic representation of the optical field and its interaction with matter, i.e., the broadest definition of CI.