|Home | About | Journals | Submit | Contact Us | Français|
We describe an adaptive image deconvolution algorithm (AIDA) for myopic deconvolution of multi-frame and three-dimensional data acquired through astronomical and microscopic imaging. AIDA is a reimplementation and extension of the MISTRAL method developed by Mugnier and co-workers and shown to yield object reconstructions with excellent edge preservation and photometric precision [J. Opt. Soc. Am. A 21, 1841 (2004)]. Written in Numerical Python with calls to a robust constrained conjugate gradient method, AIDA has significantly improved run times over the original MISTRAL implementation. Included in AIDA is a scheme to automatically balance maximum-likelihood estimation and object regularization, which significantly decreases the amount of time and effort needed to generate satisfactory reconstructions. We validated AIDA using synthetic data spanning a broad range of signal-to-noise ratios and image types and demonstrated the algorithm to be effective for experimental data from adaptive optics–equipped telescope systems and wide-field microscopy.
Images acquired using any optical system are fundamentally limited in resolution by diffraction and corrupted by measurement noise. Aberrations intrinsic to the optical system and imaging medium result in further degradation and distortions of the observed images. In ground-based astronomical imaging, atmospheric turbulence is the primary source of aberrations. In microscopic and biological imaging, significant aberrations arise as a result of index-of-refraction inhomogeneities within the sample under study.
Aberration artifacts can be largely corrected using adaptive optics (AO) methods.1 Limited by the spatial and/or temporal response of AO hardware, however, such corrections remain imperfect. AO-corrected images are often contaminated by residual blurring that can significantly reduce the contrast of fine image details. Significant denoising and improved image contrast can be obtained using post-acquisition deconvolution techniques,2 implying that both hardware and software correction strategies are needed for optimal image recovery.
Deconvolution is an explicit attempt to model and computationally compensate for measurement nonidealities. Classic approaches presume that the imaging point-spread function (PSF) of the optical system is exactly known. In practice, however, the PSF is estimated either theoretically3,4 or by imaging a subresolution pointlike object (e.g., guide star or fluorescent bead).5,6 Such estimates may deviate significantly from the true PSF, yet no margin is given in classical methods for the PSF to adjust to a more appropriate estimate. Using a fixed, imperfect PSF thus inherently limits one’s ability to generate the most accurate and highest-resolution object reconstructions.
Myopic or blind deconvolution approaches allow an imprecise or unknown PSF estimate to adapt to a more correct form and thereby offer the possibility of improved object reconstructions over classical methods. The success of these myopic–blind methods, however, is dependent upon a priori constraints that compensate for the lack of information associated with having the PSF be variable.7–11
In this paper, we describe an adaptive image deconvolution algorithm (AIDA) for myopic deconvolution of two-dimensional (2D) and three-dimensional (3D) image data within a maximum a posteriori framework. AIDA is a de novo implementation and extension of the MISTRAL (Myopic Iterative STep-preserving Restoration ALgorithm) method, originally developed by Mugnier and co-workers12 to effectively deconvolve a broad range of astronomical targets with superior photometric restoration and sharp-edge feature preservation. We have significantly improved AIDA’s run-time performance over the original MISTRAL implementation and have developed a simple yet effective scheme to balance maximum-likelihood estimation with object regularization in the deconvolution process. Moreover, AIDA has capabilities to process multiple image frames simultaneously, thereby leveraging the information available through multiple observations.2,11 In Section 2, we present the deconvolution approach. In Section 3, we describe how AIDA was implemented and describe our automatic regularization scheme. In Section 4, we demonstrate AIDA’s effectiveness on both synthetic and experimental single-frame data. In Sections 5–7, we present the application of AIDA to multiple-image-frame data and 3D images. We conclude with a survey of possible algorithmic improvements and applications, offering AIDA as an open-source alternative to MISTRAL for further development.
Consider an image, i(r), of an object, o(r), observed through a telescope or microscope system and measured using a CCD detector array. This image may be viewed as a probabilistic mapping of the object’s brightness distribution to an intensity count distribution sampled over the discrete pixel/voxel position, r:o(r)i(r). Assuming that (i) image formation is linear and space invariant (isoplanatic approximation), (ii) the response of each CCD pixel element is equivalent and independent of all others, and (iii) signal-independent Gaussian and signal-dependent Poisson noise sources are present,13 the image formed can be described by the following equation:
where h(r) is the PSF, g(r) denotes the noise-free image, G(r) is a Gaussian random variable characterized by variance , and P(r) represents a stochastic Poisson process with variance . The operator, , denotes a convolution, and denotes a pixel-by-pixel operation. While the response of CCD pixel elements is rarely uniform in practice, we will assume that any nonuniformity can be accounted for through image flat fielding with negligible effect on the validity of Eq. (1). Moreover, we assume that if any constant image background is present, it can be subtracted from i(r) so that G(r) is zero centered.
When both Gaussian and Poisson noise sources are present and images are not photon-limited, a nonstationary but additive weighted-Gaussian noise model with variance
where capitalization denotes the Fourier transform of the variable, H(k) is the optical transfer function (OTF), and k is the conjugate spatial frequency. For brevity, the dependence on r and k will often be implicit hereafter.
The goal of deconvolution is to invert Eq. (1). Classical deconvolution approaches aim to find the best estimate, ô, of the true object given a single image frame, i, and an exactly known PSF convolution kernel, h. Such approaches are ill-posed (lacking a unique solution, or having a solution that is discontinuous with respect to the data) and ill-conditioned (numerically sensitive to small errors and thus unstable) for two reasons: (1) h is intrinsically band-limited by the resolution limit of the optical system, and (2) noise is present at frequencies beyond the band limit.8,15 This situation is further complicated in the case of myopic or blind deconvolution where the characteristics of the PSF kernel are poorly known, if at all. Because of ill-posedness, the quality of the deconvolution depends critically on the quantity and quality of a priori information that is incorporated into the inversion process.8,16 This a priori information can be divided into three classes related to , o, and h.
Owing to the presence of noise, deconvolution may be viewed as a problem of stochastic inversion. It is helpful to state the goal of deconvolution in Bayesian terms, namely, to maximize the a posteriori probability of observing the object, o, and PSF, h, given an image, i, and a set of model assumptions, a,
p(i | o, h, a) is the posterior probability density of observing an image, i, as expressed by the forward-imaging equation, Eq. (1). This term is the focus of maximum-likelihood methods, which aim to optimize the fidelity of the observed data to a set of parameters and subject to a particular noise model. p(o | a), p(h | a), and p(i | a) are the a priori probability distributions for the object, PSF, and image, respectively. These a priori distributions must be inferred based on the assumptions, a. In classical deconvolution methods for which the PSF is known, for example, p(h | a) is assumed to be a constant. In maximum-entropy deconvolution methods, p(o | a) is set implicitly by the definition of the entropy measure used.17 When the positivity of the variables o, h, and i can be assumed (e.g., under incoherent imaging conditions), the a priori probabilities for negative values can be set to zero.
where we have used the subscripts n to denote noise-model-related data fidelity terms, o to denote the terms arising from the a priori object distribution, h to denote the terms arising from the a priori assumptions for the PSF, and i to denote the terms arising from the a priori distribution of images. The mode or best estimate for both o and h can be found by maximizing Eq. (6) with respect to these variables or, equivalently, by minimizing the corresponding negative log-likelihood, J(o, h | i, a),
Since Ji(i | a) is formally independent of variables o and h given the set of assumptions a, we have dropped this term in Eq. (7). We have also dropped the constant term involving the ratio of partition functions, which embodies information on the relative normalization of the component probability distributions [cf. Eq. (6)]; we use the subscript Z to serve as a reminder of this.
Our goal is to minimize Eq. (7) subject to a specific set of model assumptions for Jn(i | o, h, a), Jo(o | a), and Jh(h | a). We follow the recommendations of Mugnier et al.12 in assigning functional forms to each of these component terms as detailed below.
Assuming the mixed-Gaussian noise model of Eq. (2), the fidelity of the reconstructed object ô and PSF ĥ with respect to the observed image i can be described by the following weighted maximum-likelihood term:
Deconvolution approaches that are based solely on this term often lead to noise amplification and severe ringing artifacts. The Landweber method and the Richardson–Lucy or expectation-maximization algorithm are examples of such approaches, which assume a stationary-Gaussian and Poisson noise model for w(r), respectively.8,17 To minimize noise amplification artifacts and find a unique and stable solution in practice, Eq. (8) must be regularized. In the aforementioned methods, regularization is accomplished empirically by limiting the number of deconvolution iterations.
Equation (8) may also be regularized through a quadratic penalty term based on an object’s spatial gradient.15,16 Quadratic regularization, however, often yields results that are oversmoothed and have compromised image contrast when applied uniformly to all object features. Using a roughness penalty that is instead subquadratic for regions of high contrast has been very successful in preserving edges and other sharp object features.15,20–22 The underlying assumption here is that large gradient discontinuities in the image arise from genuine object features and should be penalized comparatively less than small gradients due to noisy background features. We use the isotropic edge-preserving prior proposed by Mugnier et al.,12 which is based on the work of Brette and Idier23:
where ‖ô(r)‖ = [(xô(r))2 + (yô(r))2 + (zô(r))2]1/2 is the norm of the spatial gradient of the object, θr and λo are auxiliary parameters or hyperparameters of the object prior distribution, γ is a reduced gradient modulus, and Φ(γ) is called the clique potential. Φ(γ) is a function that characterizes the local object texture at a position r based on a subset or clique of neighboring pixels. This clique is defined in practice through the calculation of the gradient norm in Eq. (11). For large values of γ, Φ(γ) ≈ γ, whereas for small values of γ, Φ(γ) = γ − (γ − γ2/2 + ) ≈ γ2/2, resulting in so-called L1–L2 (linear–quadratic) behavior. Numerous L1–L2 regularization functionals have been suggested in the literature (e.g, see Teboul et al.22). The advantage of Eq. (10) over other forms is that it is convex and its derivative with respect to ô does not involve any transcendental or exponential functions, making cost function optimization easier and less expensive (see Subsection 3.C).
The scaling parameter λo plays an important role in balancing maximum-likelihood fidelity to the data, with the preservation of high-contrast features in the object estimate. The hyperparameter, θr, sets the width and shape of the Gibbs distribution in Eq. (5). It governs the point at which regularization transitions from being quadratic to being linear. In Mugnier et al.’s treatment,12 the same scalar pair of values (λo, θ) is applied to each pixel element of the object. We have found that using an inhomogeneous hyperparameter model as advocated by others,24–27 in which θr is pixel/voxel dependent (as indicated by the subscript) and adapted to the local object texture, results in better deconvolution results.
To myopically reconstruct the PSF, the following Fourier domain constraint is used:
where λH controls the degree of the OTF regularization constraint relative to the data fidelity term [Eq. (8)], Ĥ(k) is the true estimate of the OTF, (k) is a measured OTF, and the overbar denotes an average over l measured OTF samples. υ(k) is the OTF sampling variance or power spectral density defined as
υ(k) serves as a spring constant to harmonically constrain each OTF k component to a mean value, consistent with a set of measured OTFs. Equation (12) intrinsically handles band-limitedness of the OTF; frequencies beyond the optical system’s resolution are essentially ignored, since they are not represented in the measured samples. Conan and co-workers28,29 have shown that this harmonic OTF constraint performs noticeably better toward recovering the true OTF than a simple band-limited constraint typically used in blind deconvolution methods.7,30 An harmonic constraint for each spatial frequency, |k|, which is functionally equivalent to using a radially averaged υ(k), may be used, although we have found that using the less stringent constraint, Eq. (12), is sometimes more robust.
The focus thus far has been on a single image frame. One of our goals in developing AIDA was to combine the demonstrated strengths of MISTRAL with the multiple-frame synthesis capabilities available in a method such as IDAC, the Iterative Deconvolution Algorithm in C.2,30,31 Christou et al.31 have argued that the use of multiple observations can serve as an additional deconvolution constraint: the ratio of unknown variables to measured quantities being reduced from 2:1 for a single image frame to (M + 1):M for M image frame observations. The simultaneous analysis of multiple observations implicitly accounts for correlations that may exist among variables as well as between variables and the data.32 Consequently, multiple-frame deconvolution should result in systematically lower error bounds with more reliable results than when individual image frames are deconvolved separately or when multiple frames are merged into an averaged “shift-and-added” image (i.e., an image generated by averaging the image frames after appropriate pixel shifts are made to maximize image correlation) and then deconvolved.2,11,33–36
The extension to multi-frame deconvolution is straightforward. For multiple-image observations, Eq. (1) may be expressed generally in vector form:
where ⊗̈ specifies a convolution performed over appropriate oj:hj pairs and we have assumed the noise model of Eq. (2). In general, for Mi measured images, there may be Mo unique objects and Mh unique PSFs: Mo ≤ Mi ≥ Mh. In addition to mono-frame data sets where Mi = Mo = Mh = 1, we consider two multi-frame data-set types in this work: (i) multi-PSF data sets where Mi = Mh and Mo = 1 and (ii) multi-object data sets where Mi = Mo and Mh = 1. Multi-PSF deconvolution may be used to process AO images for which there is a common target object but a variable PSF per image observation. Multi-object deconvolution may be used to process time-lapsed microscopy images for which a single common PSF does not change significantly between frames.
The cost function to be minimized for multi-PSF deconvolution is given by
and for multi-object deconvolution by
where α and β are used to index multiple objects and PSFs, respectively.
We implemented AIDA using Numerical Python–Numarray,37 with calls to a specialized C++ conjugate gradient (CG) optimizer (see Subsection 3.B), which were handled by code generated using the Simplified Wrapper and Interface Generator38,39 (SWIG). Fast Fourier transforms were computed using the fftw (version 2.1.5) subroutine library40 (see also http://www.fftw.org) in lieu of the standard Numarray fftpack library, resulting in about a factor of 2 improvement in the overall speed of the algorithm. A schematic of the algorithm is shown in Fig. 1.
AIDA begins with a preprocessing stage to estimate data fidelity weights, w (see below, Subsection 3.C), and to calculate the mean OTF, , and OTF variance, υ. It is assumed that all the images supplied have been properly flat fielded and optionally background subtracted. In cases where the image does not have negative pixels following background subtraction (as is the case for an image without true dark areas), the user must supply either a value for σG or a dark image from which it can be estimated.
The present version of AIDA expects images of reference PSFs (e.g., of a guide star or subdiffraction-sized bead), which are normalized to 1 and used to compute and υ. If only one PSF image is supplied, υ is calculated based on the noise statistics of the image as for w. AIDA is equipped with an optional clean-up module to remove hot–dark pixels from these PSF images and remove noise according to some user-defined threshold. An option to use a radially averaged OTF variance is provided to enable a more stringent harmonic constraint of spatial frequencies (see Subsection 2.C.3).
The default mode for AIDA uses automatic hyperparameter settings as described below in Subsection 3.D. The option to directly specify hyperparameter values or a scale factor by which to multiply the automatic estimates is available for fine-tuning purposes. For mono-frame deconvolutions, AIDA is also capable of performing unsupervised deconvolutions over a grid of λo and θr hyperparameter values centered about automatic estimates or user-defined centers.
Although it is possible to simultaneously estimate both sets of objects, ô, and PSFs, ĥ, by stacking them into a single variable to be optimized [see Eq. (7)], doing so could result in slower convergence, since significant differences in magnitude between ô and ĥ can result in a skewed optimization landscape and ill-conditioning.41 Although variable renormalization could solve this issue, we have chosen instead to alternate between the minimization of ô and ĥ in the current version of AIDA, as advocated by Mugnier et al.12
For nonquadratic cost functions, solution convergence can often be improved by periodically restarting the CG minimization after a defined number of steps so as to interlace steepest-descent steps with CG steps. We have found this partial conjugate gradient (PCG) approach41 to be more effective than a simple CG approach in minimizing the quasi-quadratic cost functions Eqs. (15) and (16), consistent with the findings of Mugnier et al.12
Starting with each PSF in ĥ set to the mean of the sampled PSFs (−1), each object in ô is optimized via a PCG approach. CG optimization is capped by a set number of iterations, ζ (typically 25), constituting a CG block and repeated for πo PCG iterations. The resulting estimate for ô is then fixed, and each PSF in ĥ is optimized via πh PCG iterations. The multi-frame estimates ô and ĥ are alternatively optimized, with each pair of estimations constituting one AIDA optimization round. The number of PCG iterations per optimization round for ô and ĥ is typically increased progressively, with the possibility of separate PCG iteration plans for ô and ĥ. By default, the number of PCG iterations executed per optimization round is given by PCG[j] = 2(j − 1) + 1, where j is the optimization round, from 1 to η, the maximum default number of optimization rounds (typically 8). Progressively increasing the number of PCG iterations in this manner ensures that the optimization of the current variable (e.g., ô) does not get fixed too quickly relative to the other variable (e.g., ĥ), which may yet be suboptimal. Multi-frame optimization of ô and ĥ is continued until the fraction of individual ôj and ĥj frame estimates that have converged is greater than some tolerance, ξ (typically >0.9), or until a specified maximum number of optimization rounds is reached. The convergence of each ôj or ĥj frame optimization is achieved when the root-mean-square deviation between two consecutive PCG iteration estimates falls below a specified tolerance for at least three times within one optimization round. We have found these default settings sufficient for processing most data sets; stricter convergence–stopping criteria typically do not yield significantly improved results.
AIDA’s quasi-quadratic cost function was minimized using a constrained CG algorithm developed by Goodman and co-workers42 and is freely available as part of the EDEN Holographic Method package.43,44 This algorithm incorporates three significant advances over the conventional CG method.45 First, to ensure that solutions are positive (or within a user-specified bound), a projected gradient or active sets approach is used.41 Johnston et al.46 have shown that such an approach is superior to maintaining solution positivity via reparametrization, since reparametrization often leads to the creation of spurious minima that can complicate the optimization process. Second, to prevent zig-zagging behavior that can arise when using an active sets approach or minimizing nonquadratic functions, an adaptive bending line search is used to set the most effective conjugate direction step size (typically called α). Third, to better preserve conjugacy between successive directions, the CG deflection parameter (typically called β) is computed using the Hestenes–Stiefel formula instead of the standard Fletcher–Reeves or Polak–Ribiere formula.41
To facilitate modification and future developments of AIDA, the calculation of the cost function was written in an extensible manner in which cost function terms may be turned on or off. For computational efficiency, only terms that are dependent upon the variable being estimated are computed (e.g., for ôj, data fidelity and object regularization terms, but not the OTF constraint, are computed).
The first term accounts for Gaussian detection–electronic readout noise, , which can be estimated using the average over all negative pixels in the image. For images of extended objects that do not have any negative-pixel areas (common in microscopy), a separate dark image is required from which can be computed directly. The second term in Eq. (17) accounts for Poisson photonic noise, ; this term is derived from the fact that the variance equals the mean and the mode for a Poisson distribution. Although this term should technically be determined using a noise-free image estimate, , we did not observe a significant improvement in deconvolution quality to merit using this more accurate though algorithmically complicated approach.
The estimates for the variances in Eq. (17) implicitly assume that i has been properly background subtracted so as to lead to a properly centered and sampled Gaussian distribution for readout noise. Only noise arising from the image formation is accounted for here. “Scientific noise” (e.g., cellular autofluorescence in microscopy imaging), which may be irrelevant to image features of scientific interest, are not accounted for here explicitly but treated as an optically genuine component of the object under observation.
The clique potential [see Eq. (10)] used for edge-preserving object regularization requires that effective spatial gradients of the object estimate be computed. This can be done efficiently by convolving the object estimate with a gradient mask:
where Gr, is a 3 × 3 matrix operator corresponding to the gradient of interest in the direction r and χ is a scaling normalization factor. Many different gradient masks that have been developed for image segmentation may be used.17,47 We prefer masks based on the work of Frei and Chen,48 since it is equally effective on horizontal, vertical, and diagonal edges, and we have found these operators to be more effective in recovering subtle object features than traditional nearest-neighbor finite-difference approximations (see, e.g., Press et al., Section 5.745). In two dimensions this is given by
and in three dimensions it is given by
where and κ is a z-resolution compensation factor. In 3D microscopic imaging, the OTF support in the axial direction is significantly smaller than in the radial direction. This leads to a greater loss of information and thus increased blurring in the z direction relative to the x or y direction; κ is used compensate for a more diffuse and uncertain gradient observed in the z direction of the image stack so that axial and lateral gradient information are on equal footing. Given the lateral and axial resolutions of a microscope, rxy and rz, κ can be estimated as κ ~ rxy/rz. If we define optical resolution as the distance between the central maximum and the first minimum of the lateral or axial component of a PSF Airy disk, the lateral and axial resolutions of a microscope are given by rxy = 0.6λem/NA and rz = 2λemn/NA2, where λem is the wavelength of light, n is the index of refraction of the sample, and NA is the numerical aperture of the microscope objective lens.49 Thus,
and, using values typical in microscopic imaging (n ≈ 1.33, NA ≈ 1.4), κ ≈ 3.
Minimizing the AIDA cost function [Eq. (15) or (16)] with the CG method requires analytical derivatives with respect to both object and PSF estimates. These can be determined through functional differentiation50 and are given by
where denotes a correlation. In practice, the terms in curly brackets are computed in the Fourier domain, in accordance with the convolution- and correlation-Fourier theorems.45,51 We assume that the arrays (or region-of-interest subarrays) used in Fourier calculations are sufficiently padded so that boundary aliasing problems can be ignored. In computing the derivative of the OTF constraint with respect to h [rightmost term in Eq. (23)], we have used the property of the discrete Fourier transform, [x*] = Nd−1[x], where x* is the conjugate of x.
The spatial Laplacian of the object in Eq. (22) may be computed by convolving the spatial object gradient with a gradient mask [cf. Eq. (18)] as proposed by Mugnier et al.52 Alternatively, the object may be convolved directly with the following Laplacian operator mask, which we find to be faster and yield finer results:
where in two dimensions
and in three dimensions
where κ once again compensates for the relative loss in resolution in the z versus xy directions (typically ~3).
Methods to estimate the hyperparameters that tune object regularization terms such as Eq. (9) have been a subject of considerable attention.24–27,53–59 A number of approaches have been advocated including L-curve analysis and generalized cross validation.54,55 These heuristic methods are computationally expensive, essentially requiring that multiple deconvolutions be performed over a grid of λo values for each image to be processed. Other more advanced and theoretically rigorous approaches attempt to optimize hyperparameters jointly with object reconstruction.54,58,59 These methods aim to maximize the marginal likelihood of observing the measured image given an incomplete data set over the space of hyperparameters: (r, o) = arg maxθr,λop(i | θr, λo); this is functionally equivalent to maximizing the ratio of partition functions, Z/ZoZn [cf. Eq. (6)], with respect to the hyperparameter variables.27,59 In practice, these methods require nontrivial Monte Carlo expectation-maximization sampling steps prior to object reconstruction, which increases the computational expense of a deconvolution considerably.24,57 In contrast to all of these methods, our AIDA approach directly calculates hyperparameter estimates using a semiempirically-based scheme, forgoing any stochastic sampling steps or comprehensive grid searching.
Our initial efforts to derive an automatic scheme were founded upon a large collection of deconvolution results generated over a grid of θr and λo values spanning several orders of magnitude. We used a variety of different 2D object types and natural scenes to build a reference set of images covering a broad range of signal-to-noise ratios. A subset of these reference objects is shown in Fig. 2. These reference images were used to assess deconvolution quality as a function of hyperparameter pairs. From a grid search over hyperparameters, a plane of acceptable (θr, λo) solutions (determined by visual inspection) was found to exist, in agreement with observations by Jalobeanu et al.26 This finding implies that one hyperparameter may be defined while the other hyperparameter is optimally adjusted to balance data fidelity with object regularization. Within the AIDA cost function framework for a single image frame, we found a balance can be achieved by setting θr according to
and computing λo directly via the approach detailed below. The form of θr was motivated by general trends observed in the aforementioned set of grid search results as well the desire for a simple scalar form for λo (see below).
From Eqs. (8) and (9), the following partition function-like integrals may be defined over the distribution of possible data-model variations, δ i − o h, and the distribution of possible gradient norm values for each pixel element:
A convenient relation linking θr and λo can be obtained by equating these integrals:
where the approximation holds for λo 10. The element-by-element equivalence of these integrals essentially assumes that the behavior of each pixel/voxel element can be decoupled and that the Gibbs distribution (and thus partition function Z) of Eq. (5) can be represented as a product of separable functions (i.e., a mean-field approximation).59 Equating these integrals effectively defines the balance of maximum-likelihood estimation with edge-preserving regularization: it is achieved by properly normalizing the probability distributions for data fidelity and object gradient norms with respect to one another. In more rigorous marginal likelihood–based hyperparameter estimation approaches,24,54,57–59 partition functions over primitive model variable(s) (e.g., i or o) are used, which lead to nonanalytical equalities that require expectation-maximation sampling in order to be solved. Our scheme estimates the sum over all states using conglomerate variables instead [Eqs. (28) and (29)], leading to the approximate though analytical relation of Eq. (30). Solving for λo in expression (30):
This definition, along with the vector definition of θr, Eq. (27), leads to a simple, pixel-independent scalar expression for λo:
From Eq. (17) and given the quantized nature of real, noisy data, σG is guaranteed to be such that θr and λo are well defined by Eqs. (27) and (32). Using w(r) as defined in Eq. (17) and object gradients and Laplacians calculated according to expressions (18)–(26), this estimation scheme is quite robust for data with PSFs of compact spatial extent (effective FWHM 8 pixels). For imaging data with spatially extended or oversampled PSFs, the pixel-by-pixel integral equivalence approximation used in Eq. (30) breaks down and can lead to somewhat overregularized results. In such cases, scaling the single scalar hyperparameter estimate, λo, down by typically no more than a factor of 10–100 is sufficient to generate optimal reconstructions. It is important to note that careful estimates of σG and w(r) in accordance with Eq. (17) are important for the success of this estimation scheme.
For the OTF constraint, a quadratic term in real space [Eq. (8)] must be balanced with a quadratic term in Fourier space [Eq. (12)]. Consistent with the fast Fourier transform40 normalization scheme used in our algorithm, we have found that this balance can be approximately achieved by setting
where Nd, the number of pixel/voxel elements, is assumed. The heuristic motivation for this comes from the power conservation relation of Parseval’s theorem for discrete Fourier, transforms, in which .
In Fig. 3, we present classical deconvolution results for one of our synthesized data sets to demonstrate the effectiveness of the automatic estimation scheme. The brain object (256 × 256 pixels) shown in Fig. 3(A) is from a magnetic-resonance imaging (MRI) scan available from the Computer Vision Group at the University of Granada.60 This object was convolved with a Gaussian PSF of FWHM of 4 pixels and normalized to a maximum intensity of 1000. This noise-free image, g(r), was subjected to a Poisson noise transformation. Varying amounts of Gaussian noise were subsequently added (mimicking CCD detector readout noise) according to a predetermined image signal-to-noise ratio (SNR), which we define as
where var[g(r)] is the variance of the noise-free image.8
Significant denoising can be observed after deconvolution [Fig. 3(B)] with a contrast enhancement of about 50%. Average contrast improvement was computed by multiple (N ≥ 6) comparisons of average intensities over an area of 3 × 3 pixels within a region of interest (IROI) versus over an adjacent background region (Ibackground) (separated by at least 4 pixels, the FWHM of the PSF):
Using the definition
we see signal-to-noise improvements of 6.2, 4.2, and 2.4 dB for the deconvolution results of SNR = 0, 10, and 20 dB images, respectively.
Figure 4 shows the deconvolution results for the SNR = 20 dB image of Fig. 3 over a grid of λo or θr values that are 20 times larger or smaller than those automatically estimated. Using the estimated hyperparameters (Fig. 4, center) gave the best visual results and balance between data fidelity and regularization. Using the estimated o and a value of θr= r/20 also gave acceptable results (though contrast was slightly compromised). In general, the deconvolution results were generally less sensitive to changes in θr than λo over the range of values examined. Although not shown, we note that AIDA’s hyperparameter estimation scheme works equally well for a range of maximum intensity scalings (i.e., images for which the maximum intensity of the noise-free image is 100 or 10,000). Deconvolution results were typically generated within 30–90 s per (256 × 256) image pixels on a 2.8 GHz Intel Xeon Linux machine.
In Figs. 5 and and6,6, we demonstrate the capabilities of the myopic deconvolution approach with a synthetic phantom composed of pointlike, line, and smooth extended elements. The object in Fig. 5(A) was convolved with a true PSF [Fig. 6(B), left] taken from a set of aberrated PSFs generated using pupil functions with random Zernike polynomial phase components of up to order 15 (Gaussian-distributed amplitudes with mean = 0 and standard deviation = 0.1). The resulting noise-free image was normalized, and Poisson and Gaussian noise was added as described above for a combined image SNR of 17 dB [Fig. 5(B)].
Classical deconvolution of this image using a fixed average PSF () results in significant denoising and contrast enhancement (ΔSNR = 1.7 dB), although artifacts can be seen in the reconstructed object [Fig. 5(C), bottom]. Allowing the PSF to relax through myopic deconvolution [Fig. 5(D)] helps to remove these artifacts and further improves image contrast (ΔSNR = 2.9 dB). Object recovery is not perfect, however, as highlighted in the bottom panel of Fig. 5(D): (1) dotlike features are larger than in the true object, and two out of the three dots shown are not fully resolved; (2) some residual haze surrounds the two intersecting line elements, and the square-on-square feature is slightly compressed in the lateral direction. The diameter of the dots can be reduced, and the remaining haze around the line elements can be removed by scaling the estimated o hyperparameter down by a factor of 2 (Fig. 5(E); ΔSNR = 4.2 dB). With slightly lower regularization, however, the square-on-square feature becomes less smooth, highlighting the intrinsic balance between noise suppression and edge preservation. For comparison, classical deconvolution results using the true PSF and the scaled λo hyperparameter value are shown in Fig. 5(F) (ΔSNR = 3.8 dB). The two lower dot features (separated peak to peak by ~3 pixels) remain unresolvable, although this is consistent with the resolution limitations of the simulated PSFs (FWHM of 3–4 pixels). Stricter a priori constraints that assume pointlike objects may lead to improved separation of these features.29,61 Owing to imperfect noise suppression, the edges of the square-to-square feature are more jagged in the classical result versus myopic deconvolution result, in which the PSF is allowed to relax. The relaxation of the PSF also results in better noise suppression and fewer noise speckles in the myopic deconvolution result versus the classical result; this leads to an improved ΔSNR for the myopic result over the classical result. Artifactual lateral compression of the square features is not seen in the classical result as it is in the myopic result, however.
Photometric comparisons with the true phantom object are shown in Table 1 for each of the highlighted features in Fig. 5. With the exception of dotlike features, myopic deconvolution using automatic hyperparameter estimates can recover intensity values to within ~10%; this is only slightly improved by λo scaling. However, using the true PSF or scaling down λo can dramatically improve the photometric recovery over the dotlike features by 15%–30%.
Displayed in Fig. 6(B) are the true PSF (htrue), the average PSF used as the initial guess in myopic deconvolution (), the myopically recovered PSF using AIDA (ĥ), and the myopically recovered PSF using a band-limited OTF constraint (ĥband-limited). Myopic deconvolution using a simple band-limited constraint results in an expected delta-function-like solution for the recovered PSF. Myopic deconvolution using a harmonic frequency constraint based on sampled PSFs prevents a delta-function-like solution and leads to the recovery of the Airy ring around the core of the true PSF that is only faintly visible in the average PSF. Artifactual line elements of the object are also present in the recovered PSF, however. This leads to a more laterally extended PSF and gives rise to the slight compression observed for the myopically reconstructed object [Figs. 5(D) and 5(E)]. Given the highly variable bounds of the sampled PSFs [Fig. 6(A)], complete separation of object and PSF features in myopic deconvolution is unlikely without further constraints.
Below, we demonstrate the effectiveness of AIDA in myopically deconvolving real imaging data for two astronomical targets, Io and Titan.
Io is the innermost Galilean satellite of Jupiter with a diameter similar to Earth’s moon (~3600 km) and is known to be volcanically active. To understand the origin of Io’s volcanism, its time evolution, and relationship to tidal heating, its volcanic activity needs to be monitored over a large time baseline. With the demise of the Galileo spacecraft that was in orbit around the Jovian system until 2003, the monitoring of Io volcanism now lies in the hands of ground-based observers.
When Io is closest to earth, its angular size is ~1.2 arcsec, which is very close to the natural angular resolution (seeing) provided by ground-based telescopes. Because of its brightness (apparent visual magnitude, mυ ~ 5), Io is ideally suited for observation by adaptive optics (AO) systems. Volcanism on Io has been monitored regularly in the near infrared (NIR) between 1 and 5 µm by one of us (F. Marchis) using the Keck 10m telescope AO system.62–64 The angular resolution provided by AO varies with the wavelength range of observations from 55 milli-arcsec (mas) in the Kc band (centered at 2.2 µm) to 100 mas in the Ms band (4.7 µm), corresponding, respectively, to ~170 and ~305 km on the surface of the satellite. Such spatial resolution is comparable with that of the Galileo observations of Io in the same wavelength range.65
Marchis and co-workers62,64 used MISTRAL to process the first high-resolution AO images of Io volcanic activity. We compared the performance of AIDA (with automatic hyperparameter estimation) with that of MISTRAL with a set of Io images acquired in 2003. The deconvolution results for three different broadband filter observations are shown in Fig. 7. Each basic-processed filtered image was a shift-and-added synthesis of five observations (<5min each; background subtracted and flat fielded). The improvement in image contrast after deconvolution is obvious. In the Kc band, the surface reflectance or albedo markings including dark paterae and bright frost areas are visible on the surface of Io. The general features of Io are in excellent agreement with those of Galileo/Voyager maps shown in Fig. 8. AIDA and MISTRAL deconvolution results are extremely similar, with a correlation coefficient of 99.4% when calculated over the area of the satellite.
For a single, 512 × 512 image, our AIDA implementation was 15–20 times faster than the original MISTRAL implementation (e.g., ~25 min versus ~7 h on a 1.8 GHz iMac G5 computer running Mac OS X 10.3). In practice, multiple MISTRAL deconvolutions must typically be performed to hone in on hyperparameter values that yield the best results. This is often a time-consuming and laborious process: between 10 and 20 MISTRAL deconvolution runs are usually necessary to locate an optimal (θr, λo) pair. Thus, the practical gain in processing time of AIDA compared with MISTRAL is >100-fold.
The image of Io in the Ms band is radically different than for the Kc band, being dominated by the localized thermal emission of the volcanoes. In the Lp band (intermediate wavelength, ~3.8µm), large-scale albedo features on the surface are visible as are the thermal emissions of the active centers. After deconvolution, several additional hot spots were revealed on the hemisphere of Io. Most of them can be found in the basic-processed image upon more careful scrutiny. The Lp band result generated with AIDA using automatic hyperparameters is noticeably different (more diffuse bright spot and some slight ringing) compared with that of MISTRAL, although these differences can be reduced by manually adjusting the hyperparameters (data not shown).
The accurate recovery of image intensities from which the temperature and emission areas of these hot spots can be determined (e.g., assuming a blackbody emission law) is also of interest. Hot-spot flux was measured using aperture photometry on the deconvolved image, assuming that most of the flux is gathered in an area slightly larger than the angular resolution on the image.66 This is a good approximation for hot spots with a peak contrast lower than 20%, since the intensity of the first Airy ring is negligible compared with the variation of brightness on the surface. For the extremely bright hot spot (outburst) on the Ms band image, a prominent Airy ring remains after deconvolution. This residual artifact may be explained by the fact that the Keck PSF is hexagonal in shape67 and that its orientation changes with the position of the telescope; optimizing the rotation of the sampled PSFs (and thus the mean PSF to which the PSF estimate is constrained) would likely minimize this artifact. Since this problem would not significantly affect the scientific analysis of the image, we have not pursued this matter further. The hot spot can be seen on the basic-processed image with a very good SNR, and therefore its integrated intensity can be easily measured after comparison with the PSF. Overall, the deconvolution of Io images with AIDA provides excellent reconstructions, which can be used to analyze surface changes on Io and to detect the faintest active centers and quantify their intensities.
Titan, Saturn’s largest moon, was largely a mystery until very recently. Observations collected by the Voyager spacecraft in 198168 showed that Titan is obscured by a dense and opaque atmosphere consisting mainly of nitrogen. The surface of this 0.9″ angular-sized satellite, however, can be probed in the NIR through methane windows using such high-resolution techniques as speckle imaging69 and AO.70 Recent AO observations of its atmosphere revealed the presence of clouds and a complex structure with seasonal variability. The NASA-ESA Cassini–Huygens probe in orbit within the Saturnian system and an intensive campaign of observations using AO systems available on the Keck 10m telescope (Mauna Kea, Hawaii) and the ESO-8 m Very Large Telescope (Cerro Paranal, Chile) are in place to help understand this complex satellite.
In Fig. 9(A), we show a ground-based observation of Titan taken on January 15, 2005, one day after the Huygens probe landed on its surface. Titan was observed with the Keck AO using the NIRC-2 camera with a pixel scale of 9.94 mas through a narrowband He filter (2.06 ± 0.03 µm). At this wavelength, the atmosphere is nearly transparent, and most of the structures visible on the image are larger than 330 km (corresponding to 55 mas). A remarkable gain in image contrast is obtained after AIDA deconvolution, as shown in Fig. 9(B). This imaged hemisphere contains the landing site of the Huygens probe and was regularly observed by the Cassini spacecraft [Figs. 9(C) and 9(D)]. The similarity between the Imaging Science Subsystem images (with a slight rotation of Titan) is striking. The smallest albedo structures detected after deconvolution have clear equivalents in the higher-resolution image71 (see arrow markings). This comparison validates the efficiency of our algorithm and demonstrates the absence of significant artifacts on the deconvolved image. A full scientific analysis of this and numerous other Titan observations and deconvolution results is presented elsewhere.71
When multiple AO images of a common object are acquired, they are often simply combined into a single shift-and-added image, which is then deconvolved. This practice has been demonstrated by others to be suboptimal; a more effective data reduction strategy would be to deconvolve the set of images in a global fashion, linking common variables while maintaining the distinctiveness of each observation. Extending the MISTRAL approach to simultaneously deconvolve multiple image frames is another feature of AIDA. Below, we present deconvolution results for two different multi-frame data sets. The first consists of AO images of Uranus’s atmosphere and is used to demonstrate AIDA’s multi-PSF deconvolution capabilities, in which there is a common object but a variable PSF. The second data set consists of time-lapsed fluorescence microscopy images of yeast microtubule dynamics and is used to demonstrate AIDA’s multi-object mode, in which there is a common PSF but different objects between frames.
Since the Voyager spacecraft encounter of the planet Uranus in 1986, interest in this planet has been revitalized with the discovery that its atmosphere is considerably active.72 High-angular-resolution imaging, however, is necessary to detect cloud motions,73 faint rings, and small satellite systems.74,75 The extended disk (diameter ~3.6 arcsec) of the planet (integrated apparent visual magnitude, mυ ~ 6) is bright enough to be used as a reference for wavefront sensor analysis on most AO systems. However, since the position of the centroid on the wavefront is not well determined in the case of a quad-cell aperture wavefront sensor for such an extended object, the atmospheric correction is degraded in the final image, and artifacts may appear in several frames.75 We tested AIDA on observations of Uranus taken on October 3, 2003, with the Keck AO system and its NIRC-2 camera, using a broadband filter centered at 1.6 µm (H band). Five 30 s frames recorded in less than 8 min were processed using standard near-infrared data reduction techniques (flat-field, sky subtraction, and bad pixel removal). To estimate the PSF for myopic deconvolution, we imaged Puck, a bright satellite of Uranus located 2.4″ away from the center of the planet and whose motion was negligible during the exposure time. Given the large imaged size of Uranus and size of the image frames (1024 × 1024 pixels), using MISTRAL for deconvolution would not have been practical due to the long processing time needed (~23 h/deconvolution on a Sun Ultra 10 computer), especially since we would have needed to run multiple deconvolutions to determine a good choice of regularization parameters. Deconvolution using AIDA with automatic hyperparameter estimation was significantly faster (45 min for mono-frame deconvolution and 1.5 h for multi-PSF deconvolution on a 2.8 GHz Intel Xeon Linux machine) with the possibility of analyzing all AO data frames simultaneously.
Deconvolution results in significant image sharpening (Fig. 10), with a gain in contrast of ~2–3 on the cloud features. A layered structure of the northern haze and some faint clouds at ~40° latitude are revealed, and the structure of the large clouds on the southern hemisphere is clearer after deconvolution. A ghost outer ring artifact present in previous observations using the same Keck AO system75 is visible in several of the individual AO-corrected image frames [Fig. 10(C)]. This artifact remains in the mono-frame deconvolution of the shift-and-added combined image but is half as intense in the multi-frame deconvolution result [cf. Figs. 10(D) and 10(E)]. The glare of Uranus (e.g., see area near the innermost ringlet) is also further reduced in the multi-frame deconvolution result than in the mono-frame deconvolution result. Overall, we find that simultaneous deconvolution of multiple-frame data is better able to restore low SNR features and minimize artifacts than the deconvolution of a single shift-and-added representation of the multiple-frame data.
Microtubules are hollow cylindrical polymers that radiate from near the nucleus of a cell and serve as tracks upon which cellular components are transported. Roughly 25 nm in diameter, these microtubules are formed from the stochastic polymerization and depolymerization of α- and β-tubulin proteins. The regulation of microtubule dynamics has been a topic of investigation for many years in cell biology, aided greatly by the direct observation of microtubules using time-lapsed video fluorescence microscopy.76
We used AIDA in multi-object deconvolution mode to process time-series images of microtubule dynamics in the fission yeast, Schizosaccharomyces pombe. Using the OMX wide-field fluorescence microscope system developed recently in our lab at the University of California, San Francisco (UCSF), a yeast cell whose microtubules were fluorescently labeled using the green fluorescence protein fused to α-tubulin was imaged every second. Each image was formed by physically sweeping the microscope focus (by linearly moving the sample stage) through the entire z depth of the cell (~4 µm in 50 ms) every second. Using estimates of the PSF based on a set of three images of a 0.1 µm fluorescent bead acquired under similar conditions, these time-series data were myopically deconvolved assuming a common (time-invariant) PSF for the whole data set and assuming each image was simply a snapshot of a distinct object.
In Fig. 11(A), we show the results of standard myopic deconvolution and multi-object deconvolution with automatic hyperparameter estimates for a single representative time slice. In Fig. 11(B), the corresponding kymograph plots—1D maximum intensity projections of each image as a function of time—are shown for these data. These kymograph plots provide a better perspective on the time-dependent features of microtubule growth and shrinkage. The mono-frame deconvolution results are significantly denoised with improved microtubule contrast. The multi-object deconvolution results have even better contrast enhancement, exhibiting thinner microtubule fibers and a more textured background within the cell cytoplasm. It is unclear how much of this texturing may be artifactual. However, the fact that each image slice was deconvolved independently with respect to the time axis and that a number of cell background features are temporally persistent in the kymograph suggest that some of these grainy features are genuine.
One main advance of AIDA is the extension of the MISTRAL method to deconvolved 3D data commonly encountered in biological imaging. Unlike the 2D PSFs encountered in low-numerical-aperture astronomical imaging, the PSFs in optical microscopy are more diffuse, with significant axial (z-dimensional) blurring on the order of three times the lateral blur. Deconvolution is expected to dramatically sharpen image data subject to such out-of-focus blur. Recently, Chenegros et al.77 demonstrated the effectiveness of MISTRAL’s edge-preserving regularization term in deconvolving synthetic 3D retinal images. Here, we show myopic deconvolution results for two 3D data sets, one synthesized from magnetic-resonance imaging (MRI) data of a frog and another of real, wide-field fluorescence microscopy data of chromosomes within cells undergoing cell division.
We constructed synthetic 3D frog images (128 × 256 × 256 pixels) by convolving a MRI volume data set from The Whole Frog Project (Lawrence Berkeley National Laboratory)78 with a PSF derived from microscopic imaging of a subresolution (0.1 µm) fluorescent bead; Poisson and Gaussian noise was added to the convolved image as described earlier. The PSF used had a FWHM in the lateral direction of ~3 pixels and an effective resolution loss in the z direction (κ) of ~3 [see Eq. (21)]. Using an ensemble of similarly acquired experimental PSFs, these frog images were myopically deconvolved using automatic hyperparameter estimates (~6 h on a 2.8 GHz Intel Xeon Linux machine).
Additive 2D volume projections for the raw and deconvolved 3D image stacks for image SNRs of 0 and 20 dB are shown in Figs. 12(A) (en face) and 12(B) (side view). The denoising and object reconstructions for these data are striking. The quality of the deconvolution results conveyed by these 2D projections is comparable to that seen from a comparison of individual 2D slices. Representative slices through the 3D volume stack of the original object, 20 dB SNR image, and deconvolution result are shown in Fig. 13; also shown are intensity line profiles (denoted by an asterisk) through the eye region of the 2D frog slices. Deconvolution with AIDA leads to substantial photometric restoration of the original frog data, with a signal-to-noise improvement (ΔSNR) of 5.7 and 5.1 dB for image SNRs of 0 and 20 dB, respectively.
Nearly 50 years since the atomic structure of DNA was elucidated, the higher-order structural organization of DNA within chromosomes of cells remains poorly understood. With recent advances in high-resolution microscopic imaging and fluorescent labeling technology, however, discerning the mesoscopic arrangements of DNA within living cells is becoming more of a reality. A primary interest of ours is to better understand the detailed structural changes of chromosomes as a cell divides in a process called mitosis. During mitosis, a cell’s chromosomes are unraveled, condensed, and separated; defects in chromosome structure during any of these mechanical steps could have devastating consequences on the fidelity of genetic transmission to daughter cells.79
Drosophila melanogaster (fruit fly) embryos offer a unique opportunity to study chromosome structural changes during mitosis. Cells in early embryos (within 2–3 h) undergo multiple rounds of cell division in an ordered and highly reproducible manner. Using the OMX microscope system mentioned earlier (Subsection 5.B), a 3D image stack (32 × 512 × 512 pixels) was acquired of a cell-cycle-10 D. melanogaster embryo fixed in 10% formaldehyde and mounted in glycerol. Cells in this embryo were stained with the DNA-specific dye, DAPI, and captured undergoing anaphase, the stage of mitosis in which chromosomes separate. This image stack was deconvolved myopically using a PSF derived from an image of a 170 nm fluorescent bead under similar imaging settings. Image pixel spacing was 80 nm in xy and 150 nm in z, for a total image stack thickness of 4.8 µm. ζ was set to 3.2 based on the extent of a measured OTF in the lateral versus axial directions.
Shown in Fig. 14 are 2D maximum intensity projections of representative portions of the original 3D image stack and the result after myopic deconvolution. Although the original data shown are of especially good quality so that most chromosome arms can be distinguished in Fig. 14(A), chromosome boundaries are significantly more demarcated in the deconvolution result. The benefits of deconvolution are even more pronounced in Fig. 14(B) in which there is greater blurring in the axial versus lateral directions: finer structures and corrugated banding patterns of the chromosome arms become noticeable; the arrows highlight a few representative areas showing improved contrast in fine image features. Some residual hour-glass PSF blur remains after deconvolution, however, and appears to become more prominent with increasing z depth (see, e.g., lower left of deconvolution result, Fig. 14(B)). This blur may be attributed to greater index-of-refraction aberrations between the microscope objective lens and the sample as one focuses deeper into the embryo. The true PSF in this case is thus likely to be depth dependent, although space-invariant PSFs are assumed in the current AIDA deconvolution framework.
To achieve the nonblurry, visually balanced deconvolution result of Fig. 14, we found it necessary to scale the automatic hyperparameter estimate, o, down by a factor of 10. Inaccurate hyperparameter estimation is likely due to at least one of four possible causes. First, since only a single PSF estimate was available for these data [in which case the OTF constraint is based simply on the photonic-noise variance (see Subsection 3.A)], the calculated OTF statistics may not be sufficient to guide the myopic deconvolution toward a more correct OTF. A lower o likely compensates for imprecise OTF statistics. Second, as alluded to above, depth-dependent variations of the true PSF are not accounted for in our imaging model and may lead to compromised object reconstructions. Third, there may be noise sources (e.g., out-of-focus, scattered background light) that are not accounted for by the assumed noise model; the effectiveness of the hyperparameter estimation scheme is predicated upon good estimates for the Gaussian and Poisson noise statistics (as discussed in Subsection 3.D). Fourth, out-of-focus contributions to the image stack from areas of the embryo outside the image stack are not accounted for in the current imaging framework. The effects of these factors on deconvolution outcome and strategies to compensate for them are currently being explored by our group.
We have reimplemented and extended the MISTRAL approach12 to myopically deconvolve, as far as we know for the first time, multiple-image-frame data and 3D image stacks. Unlike MISTRAL, which is implemented using the commercial Interactive Data Language (Research Systems, Inc., Boulder, Colorado) and has proprietary source code, our adaptive image deconvolution algorithm, AIDA, was implemented using freely available Numerical Python and is intended for open-source development. AIDA runs at least 15 times faster than the original MISTRAL implementation. Importantly, AIDA incorporates a simple yet robust scheme to estimate regularization hyperparameters, which greatly simplifies the tedious and delicate though necessary task of balancing maximum-likelihood estimation with object regularization and noise suppression. Object reconstructions can be generated using AIDA that are comparable with those of MISTRAL, with high photometric precision and good edge preservation and without the need to sample (typically 10–20) different hyperparameter settings in order optimize the degree of regularization. This results in a practical efficiency gain of AIDA over MISTRAL of greater than 100-fold.
Multiple image observations are commonly acquired in adaptive optics imaging, although they are often combined into a single averaged image before deconvolution. Deconvolving these images simultaneously, however, is a more effective data reduction strategy.11,31–33 The multi-frame deconvolution results for the Uranus AO observations show that leveraging invariable aspects of the data while retaining the unique variations between distinct observations leads to object reconstructions with crisper details than the corresponding mono-frame deconvolution result.
AIDA’s multi-frame deconvolution capabilities are currently limited to data with a single object and multiple variable PSFs (Mo = 1; Mh > 1) or a single PSF and multiple variable objects (Mh = 1; Mo > 1). It would be straightforward to extend the algorithm to handle data sets in which multiple objects are imaged using different though known transformations of a fundamental PSF describing the optical system. This is relevant, for example, to multi-wavelength imaging in astronomy34 and microscopy in which the PSF characteristics as a function of wavelength are well established and can be predicted. Such an approach could also be applied to process tomographic imaging data in which the dependence of the transfer function is known and parametrizable as a function of tilt angle. Using such a multi-object–multi-linked-PSF approach, our group is currently exploring the application of AIDA to deconvolve electron microscopy (EM) images, with the goal of improving 3D object reconstructions from EM tomographic data.
AIDA is equally effective in deconvolving 3D image data and 2D data, and deconvolution times scale linearly with the size of the image data. In the current AIDA framework, each image pixel element is treated as a variable to be optimized, leading to substantial computational demands as image size increases. Work in our group is in progress to recast the optimization of the PSF in terms of the more computationally compact pupil function that characterizes the optical wavefront at the exit pupil of an imaging system arising from a point source.80,81 In addition to greater computational efficiency for larger image data sets, myopic deconvolution using the pupil function could provide explicit insight into the inherent or dynamic aberration modes of an optical system (e.g., by Zernike mode decomposition). The ease with which the pupil function can be modified to account for aberrations also makes it particularly amenable to use in cases where the PSF is space variant80,81 (e.g., with depth-dependent index-of-refraction variations in microscopy or anisoplanatic imaging in astronomy). Moreover, use of the pupil function could help bridge the synthesis of wavefront-sensing data from AO and imaging data in the deconvolution process.52
At least four issues merit further development and exploration. First, the reasons for the success of our automatic hyperparameter estimation scheme. While this semiempirical scheme is effective in deconvolving a broad range of image data, the theoretical foundations for its robustness deserve future study. The assumption of quasi-independent pixel/voxel prior distributions and the assumption that the balance of maximum-likelihood estimation and object regularization is best achieved by normalizing these prior distributions with respect to one another should be explored in relation to the partition functions of Eq. (6) and other more rigorous marginal likelihood approaches. Second, the development of a multi-object deconvolution mode more specifically tailored for time-series data. In deconvolving the microtubule dynamics data in subsection 5.B, the temporal independence of each object in the time series was assumed. While this was helpful in highlighting common, persistent features between time frames, incorporating a cost function term or procedure within the deconvolution algorithm to maximize the temporal correlation between adjacent time slices may help reinforce object features that are self-similar and suppress temporally uncorrelated noise artifacts. Third, as image data sets become larger and/or deviations from the assumed noise model become more pronounced, the time-to-optimization convergence may become seriously compromised. Convergence might be improved by toggling between a weighted least-squares (L2-norm) form for the data fidelity term [Eq. (8)] and a robust L1-norm form that is computationally simpler and less sensitive to noise model mismatch and data outliers.82–84 Deconvolution efficiency might also be improved by a reparametrization of the object, for example, using wavelets,17 and by incorporating aspects of multi-resolution/hierarchical scaling into the deconvolution algorithm.17,85–87 Finally, it would be interesting to see how the myopic capabilities and edge-preserving noise suppression advantages of AIDA deconvolution could improve the processing of data from such superresolution imaging modalities as multi-frame mosaicing84,88 and structured illumination microscopy.89,90
We thank members of the Sedat Laboratory Z. Kam, and M. Gustafsson for discussions and comments on the manuscript. We thank Eugene Ingerman [University of California, Davis, and Lawrence Livermore National Laboratory (LLNL)] and Stephen Lane (LLNL) for making the conjugate gradient code available to us. We thank the authors of MISTRAL for allowing F. Marchis privileged use of their program and a reviewer for bringing to our attention the work of Chenegros et al.77 The astronomical data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and NASA. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. Work by E. F. Y. Hom, S. Haase, and J. W. Sedat was supported by NIH grant GM25101-26. F. Marchis was supported in part by NASA grant PAST04-0000-0025 and by the National Science Foundation Science and Technology Center for Adaptive Optics, managed by the University of California at Santa Cruz under cooperative agreement AST-9876783. T. K. Lee was supported by Howard Hughes Medical Institute summer fellowships. D. A. Agard was supported by NIH grant GM31627 and is a Howard Hughes Medical Institute investigator.
OCIS codes: 100.1830, 100.3020, 100.3190, 010.1080, 180.0180, 180.6900.
AIDA source code, news, and updates will be available at http://msg.ucsf.edu/AIDA.
Erik F. Y. Hom, Graduate Group in Biophysics and Department of Biochemistry and Biophysics, University of California, San Francisco, Genentech Hall, 600 16th Street, San Francisco, California 94143-2240, USA.
Franck Marchis, Department of Astronomy, University of California, Berkeley, 601 Campbell Hall, Berkeley, California 94720, USA.
Timothy K. Lee, Department of Biochemistry and Biophysics, University of California, San Francisco, Genentech Hall, 600 16th Street, San Francisco, California 94143-2240, USA.
Sebastian Haase, Department of Biochemistry and Biophysics, University of California, San Francisco, Genentech Hall, 600 16th Street, San Francisco, California 94143-2240, USA.
David A. Agard, Howard Hughes Medical Institute and Department of Biochemistry and Biophysics, University of California, San Francisco, Genentech Hall, 600 16th Street, San Francisco, California 94143-2240, USA.
John W. Sedat, Department of Biochemistry and Biophysics, University of California, San Francisco, Genentech Hall, 600 16th Street, San Francisco, California 94143-2240, USA.