1. Introduction

By utilizing high aperture angles, modern objective lenses provide lateral resolution at the diffraction limit. However, due to the asymmetry of an imaging system featuring just a single lens, the axial resolution is several times worse than the lateral resolution. By employing two opposing lenses and coherent illumination and/or detection from both sides, the family of 4Pi microscopes [1] doubles the aperture used for imaging, thus creating and/or detecting a total wavefront that comes closer to a complete spherical wavefront. Speculating about the use of the full solid angle in focusing, a paper of 1978 suggested that appropriate illumination from all directions would allow focusing the light virtually to a point, or at least much sharper than the diffraction limit, thus bringing about infinite spatial resolution in far-field microscopy [2]. This conjecture was based on the assumption that by holographically reconstructing the wavefronts originating from a nanosized light source (and reversing their propagation direction) a focal spot of the same nanosized dimension as the source could be created. This assumption clearly does not hold, because only the propagating part of the field but not the evanescent waves around the original source can be reconstructed. This is in fact the very essence of the diffraction limit: high spatial frequencies cannot be retrieved from light emitted by an object or encoded in the light being focused because they decrease exponentially with increasing distance from the light source regardless of the focusing geometry. The conjecture is thus invalid and in contrast, a 4Pi microscope as introduced [1] and known today relies on a total wavefront that is far from being completely spherical. Moreover, the goal is not to overcome the diffraction limit, but to provide a 3–7 fold increase in axial resolution that is still within the realm of diffraction. This important gain in axial resolution is achieved due to the alleviation of the asymmetry inherent to focusing from a single side only. By enlarging the total solid focusing angle as much as it is technically possible with aplanatic objective lenses, this arrangement pushes the far-field optical resolution in all three dimensions towards the diffraction limit.

The fact that the converging or collected total wavefront of the two-lens system is not spherical leads to side-lobes in both the illumination intensity pattern (I-PSF) of a 4Pi type A or C microscope and the detection PSF (D-PSF) of a 4Pi type B or C microscope and thus also in the effective PSF (E-PSF) of the system [3]. In the data, these lobes produce characteristic replica of the object which have to be removed from the final image by deconvolution.

Reliable and artifact-free image reconstruction is currently only possible if the lobe height is sufficiently low (< 50% of the main maximum for constructive interference) [4] and, because the lobe pattern depends on the relative phase of the opposing beam-paths, if the relative phase of the two beam paths is either constant [5–7] or if it can be estimated a-priori at every point of the object [8].

Importantly, the lobe height is readily controlled by employing e.g. two-photon excitation (TPE), confocal detection or using coherent detection in a type C arrangement [3]. In fact, employing an embedding medium of high refractive index, e.g. the water-miscible TDE (2,2-Thidodiethanol) [9], oil immersion and very high numerical apertures (NA) has resulted in virtually lobe-free TPE 4Pi microscopy [10]. Unfortunately only a sub-class of fixed samples allows this combination of measures and the phase-problem remains for all other applications. The a-priori estimation of an unknown phase relies on isolated, axially thin objects which have to be distributed throughout a sufficiently sparse sample [8]. In most application it is therefore currently mandatory to match the refractive indices of the specimen and the embedding medium and to either accurately monitor or compensate for the phase variation introduced by the axial scanning process [11]. While this combination of measures is often feasible, it is time-consuming, adds experimental complexity and restricts the choice of possible samples: For specimen which intrinsically feature a non-negligible and unknown variation of refractive index, the relative phase and thus the form of the E-PSF will always vary throughout the object rendering reliable image interpretation based on existing methods virtually impossible. We have therefore developed a parametric-blind deconvolution method that allows for simultaneous estimation of both, the unknown, spatially varying phase and the object, without prior knowledge of either. Based on a maximum a-posteriori approach, it allows for a space-variant relative phase function (RPF) ϕ(r) which is either assumed to be smooth or approximated by a linear combination of a few polynomial basis functions. Both approaches are tested on both synthetic and real data and compared. We verified the applicability of our method to one-photon confocal 4Pi type A microscopy, because this particular configuration is experimentally most simple, least demanding on sample brightness and stability and commercially available and because the data analysis is more challenging than in the case of TPE or coherent detection which result in lower side lobes and therefore less ill-posed problems.

2. Imaging model

Since imaging in the scanning 4Pi florescence microscope is incoherent, the relation between the expected image ḡ(r) and the object f(s) in the absence of noise can be described by the following equation

where h _eff is the 4Pi E-PSF, s the position in the sample and r is the position in the image back projected into the sample, i.e. the scanning position. For space-invariant PSFs where h _eff(r,s) = h _eff(r - s) Eq. (1) is a convolution and the forward operator in common restoration techniques can be efficiently implemented using fast Fourier transforms (FFT). Here, the space variant phase function precludes this by rendering the form of the PSF position dependent. Fortunately the special form of the 4Pi E-PSF allows its decomposition into a ϕ-dependent linear combination of a number of phase-independent and thus space-invariant convolution operators which can then be efficiently implemented using FFT algorithms. This is either achieved by using a convenient approximation of the PSF [8] or by deriving an adequate representation of the PSF using vectorial focusing theory [5]. Importantly, this approach also allows re-calculation of the forward operator for a given object but various phase functions without re-evaluation of the convolutions. In a previous publication we have used the latter technique to efficiently estimate both the object and a constant global phase directly from the image [5]. Here we extend the approach by replacing the constant global phase by the space-variant RPF. The form of the E-PSF of a 4Pi fluorescence microscope of type A using one photon excitation is given by

where E _1,2 denotes the counter-propagating focal fields, with E ₂(r) = ME ₁(Mr). The transformation matrix M accounts for the counter-propagation of the beams and h _det is the D-PSF. In the case of confocal detection, the D-PSF is given by the integral of the field created by a dipole in the image plane over the finite-sized pinhole [12] and E _1,2 are readily evaluated using standard vectorial focusing theory [13]. The forward operator H is thus defined by

Where the operator H(ϕ)f depends nonlinearly on the phase function ϕ but linearly on f for any given ϕ. Using Eq. (2) it can be written as

where we have also introduced its derivative with respect to the phase, H ^′. The convolution operators are independent of the RPF and given by

They are readily computed using the Fourier space multiplication technique requiring a few FFTs. A similar decomposition using five convolutions involving fourth powers of the fields is readily derived for two-photon excitation, where the I-PSF depends on the squared intensity. We note that the forward operator proposed by Baddeley and coworkers [8] requires the same number of convolutions but is based on an approximation of the E-PSF that becomes inaccurate for high numerical apertures.

Images are usually acquired on a regular 3-dimensional grid using photon counting devices. If we identify each voxel by its index n and denote the discretized object, expected image and RPF by the vectors f, ḡ and ϕ, we can write Eqs. (4) and (3) as

where the A _j are the convolution matrices associated with the convolution operators A_j. Here and in all subsequent equations multiplication and division of one vector by another is meant voxel by voxel. We introduced also a possible additional background contribution b. The measurement process is dominated by shot noise and count rates are usually in the range of zero to a few hundred photons per voxel. The measured value g(n) at a voxel n is therefore the realization of an independent Poisson random variable with its expectation value given by ḡ(n). During reconstruction we maximize the a-posteriori probability of observing the specimen f and RPF ϕ given the image g. It is given by

P(ϕ), P(g) and P(f) denote the a-priori probability distributions for the phase, for the image and for the specimen, respectively and the probability of observing g is given by

where poi(g∣λ) = λ^g exp(-λ)/g!. We assume that the random variables associated with the object and the RPF are Gibbs random fields and thus

where x = f,ϕ for object and RPF respectively, μ_x > 0 is the regularization parameter, Z is the partition function, and J_x is the penalization function, or simply prior. Maximization of the a-posteriori probability distribution with respect to f and ϕ is then equivalent to minimization of its negative logarithm, the so called cost function

where constant terms are omitted and the likelihood term J_o is given by

The regularization parameters are usually adjusted empirically and the penalization functions J_f,J_ϕ, used to introduce a-priori knowledge about the object and RPF into the reconstruction have to be chosen carefully. This kind of maximum a-posteriori (MAP) approach is well established in microscopy and therefore, suitable object penalization functions for typical objects have been proposed in the literature (for a recent review see [14] and references therein). The most common ones are Tikhonov regularization J_f(f) = ∥f∥² ₂/2 which imposes an energy constraint on the restored object, the H ¹ case J_f(f) = ∥∣∇f∣∥² ₂/2 that favors smooth objects and total variation (TV) regularization J_f(f) = ∥∣∇f∣∥₁ for objects with sharp edges and areas of (nearly) constant intensity. Here, ∥ · ∥₁ and ∥ · ∥₂ denote the usual l ₁ and l ₂ norms. For the microtubule images analyzed in this work we found that classical Tikhonov penalization resulted in visually appealing reconstructed objects. As a side-benefit, its implementation has no significant impact on calculation time when compared to the unconstraint case J_f(f) = 0.

The key to successful simultaneous estimation of phase and object is the appropriate choice of the prior for the phase, J_ϕ. It has been shown that the RPF induced solely by refractive index mismatch between mounting medium and immersion medium depends linearly on the axial scanning position for moderate sample thickness [15]. While contributions from the object, e.g. small organelles and the nucleus break this simple dependence, the RPF usually remains smooth and never exhibits phase jumps [11,16]. In this work we therefore explore two approaches for the penalization function, J_ϕ, based on either allowing for an arbitrary but smooth RPF or for such relative phase functions that can be described by a linear combination of a few suitable polynomials.

3. Arbitrary, smooth phase function

Here, we allow for an arbitrary phase function and hence we call this approach FPBR (full parametric blind reconstruction). Regularization is achieved by favoring smooth RPFs using H ¹ penalization. To discretize the H ¹ penalization function we assume, for simplicity, that the N_x × N_y × N_z array ϕ is expanded by taking ϕ (N_x + 1, n_y, n_z) = ϕ (N_x, n_y, n_z), ϕ(n_x,N_y + 1, n_z) = ϕ(n_x,N_y,n_z) and ϕ(n_x,n_y,N_z + 1) = ϕ(n_x,n_y,N_z). Under this assumption we write the discrete version of the H ¹ penalization function in the following form

where η _n = {(n_x + 1,n_y,n_z), (n_x,n_y + 1,n_z), (n_x,n_y,n_z +1)} are the upper direct neighboring voxels of n, ϕ (n) ∈ ℝ and J_ϕ is calculated independently of the interpretation of ϕ as a phase during minimization. We note that RPF’s that differ by constant offsets of multiples of 2π are obviously equivalent and identify the same solution making this a proper prior. To minimize the resulting cost function J(f, ϕ∣g) simultaneously in f and ϕ we use a similar alternating minimization method (AMM) as in our previous work [5]. For the l-th cycle we have

where Eq. (16) adds non-negativity of the object as a constraint. Thus, during each cycle first the RPF and then the object is kept fixed while the other is updated by minimizing the cost function w.r.t. to it. Therefore, J(f _l, ϕ _l∣g) is guaranteed to decrease as the number of cycles l increases. It is important to point out that, while convex in f for fixed RPF, J(f, ϕ ∣g) is non-convex as a function of both variables and therefore, convergence of the AMM to a global minimum is not guaranteed. The quality of the local minimum found may therefore depend on the initial guess of the RPF. We will see that a suitable strategy can be devised to ensure visually appealing restored objects. For the solution of Eq. (16) we chose the split-gradient-method (SGM) [17] due to its robustness, the simplicity of its implementation and its capability to enforce the non-negativity constraint in a natural fashion. In the case of Poisson data and Tikhonov penalization the SGM iterations are given by (see appendix for a derivation):

For μ_f = 0 the SGM algorithm is equivalent to the well-known Richardson-Lucy [18,19] algorithm for the linear operator H(ϕ). We note that a similar iterative algorithm can be obtained for all the object penalization functions mentioned above. For the RPF update we adopt a quasi-Newton method (QNM) [20] defined by the iteration

The descent direction ρ is given by

where C _l,i is the Hessian matrix of the function J(f _l+1,ϕ∣g) w.r.t. ϕ. An analytical expression for the gradient is derived in the appendix while an approximation of C can be calculated using the Broyden, Fletcher, Goldfarb and Shanno (BFGS) formula [20]. However, due to the large size of ϕ on the order of a megavoxel, we implemented the limited-BFGS approach where C ^-1 is never explicity formed or stored. Rather, a direct estimate for the descent direction ρ is obtained by maintaining a history of past updates (< 10) of ϕ and the gradient ∇_ϕ [J(f _l+1,ϕ∣g)] [20]. The step size τ_l,i can be efficiently determined by the Armijo line-search method [20], because the several calculations of J(f _l+1, ϕ∣g) necessary for various RPFs do not require re-evaluation of the FFTs due to the decomposition of H. As outlined below, suitable stopping criteria for both the SGM and QNM parts of the cycle have to be defined.

4. Sparsity-based approach

Voxel-by-voxel estimation of the RPF has the disadvantage that the number of parameters to estimate (one parameter per voxel) is very large effectively slowing down the convergence. By using a sparsity-based approach, an emerging and powerful method in image processing [21], the number of free parameters can be significantly reduced. It is based on the observation, that many natural signals can be represented as a linear combination of a few (sparse) well chosen basis functions. We assume that given a suitable set of basis functions φ_p, any RPF is approximated well by a sparse linear combination of them. Accordingly we will refer to this approach as SPBR (sparse parametric blind reconstruction). We can write

with ∥ψ∥₀ ≪ N, N being the number of voxels. We call ψ the RPF’s coefficients vector and the operator D the dictionary and ∥ · ∥₀ is the l ₀ pseudo-norm, which counts the number of nonzero elements in a vector. The choice of dictionary plays a crucial role since it implies strong assumptions about the behavior of the RPF. With ϕ being a smooth function, several choices of smooth basis functions (which need not be orthogonal) seem reasonable, including (steer-able) wavelets and curvelets. In this work we selected products of polynomials in normalized cartesian coordinates r = (x,y,z)

with p ₁ + p ₂ + p ₃ ≤ N_ϕ and the normalization of the coordinates was chosen such that the longest edge has length 1. The optimum normalization and dictionary size, i.e. the choice of N_ϕ largely depends on prior knowledge about the specimen. Choosing a large dictionary may increase computational complexity and slow convergence while small dictionaries may not be able to model a realistic RPF. In the limiting case where the phase variation is solely induced by refractive index mismatch between mounting and immersion medium the RPF is known to be linear and we could set N_ϕ to 1. For the experimental data analyzed in this work N_ϕ ≤ 6 (i.e. a dictionary size of N_ψ ≤ 84) proved adequate.

The cost function is now given by J(f, ϕ∣g) = J(f,D ψ∣g) and has to be minimized w.r.t f and the component vector ψ. Sparsity could be enforce using the l ₀ pseudo-norm penalization function, i.e J_ϕ (ϕ) = ∥ ψ ∥₀ but its minimization results in a hard numeric problem. It was shown that for several problems, a sparse solution can also be obtained by replacing the l ₀ pseudonorm with the l ₁ norm, i.e. J_ϕ(ϕ) = ∥ψ∥₁ [21]. Motivated by this work, we implemented the minimization using the same AMM as above, simply replacing the large vector ϕ _l,i by the component vector ψ _{l, i}. An analytical expression of ∇_ψ[J(f _l+1,D ψ∣g)] is given in the appendix and the Hessian matrix C _l,i w.r.t. to the much smaller vector ψ can now be directly implemented using the BFGS method.

4.1. Stopping criteria and initialization

During each cycle, the SGM algorithm was stopped when the relative change during one iteration dropped below a certain threshold ε _SGM, i.e.

To stop the QNM algorithm it turns out to be much safer to terminate the iteration when the norm of the gradient is sufficiently small compared to the norm of the gradient at the initialization value [20]:

where x = ϕ, ψ for the FPBR and SPBR approach respectively. The choice of the two parameters ε _SGM and ε _QNM was based on the following consideration. At early AMM cycles the object estimate changes significantly from one SGM iteration to the next, and thus it is desirable to update the RPF after just a few updates to the object. However, because the object estimate is changing rapidly, the RPF estimate changes significantly every time it is updated, and thus making the tolerance ε _QNM too small at the early AMM cycles would waste computation time. Therefore, we choose the two parameters ε _SGM and ε _QNM as a function of the AMM cycle. We therefore start the AMM algorithm with somewhat loose tolerance values which are tightened in later cycles. From the test cases that we ran, we found that the combination of tolerance values ε _SGM = 10^{max(-6,-⌊l/25⌋-1)} and ε _QNM = 10^{max(-4,-⌊l/25⌋-2)} gave a good compromise between the computation time required and the quality of the results obtained. To stop the AMM algorithm we used a criterion based on the difference of the values of the cost function between two consecutive AMM cycles

A typical value for ε _AMM is 10^-8 for simulation and 10^-7 for real data. In all examples we initialized f ₀ with a constant vector and ϕ ₀ with a constant RPF.

5. Validation on simulated data

Both PBR approaches were tested on simulated data and an example is shown in Fig. 1. We used a 2-dimensional (xz) image in order to speed up evaluation and parameter optimization and created a test pattern that is reminiscent of different subcellular structures, i.e. microtubules, Golgi and mitochondria. Using a strongly varying cubic RPF (N_ϕ = 3) we obtained the simulated noise-free image using the forward operator from Eq. (9), added a constant background and subsequently applied a Poisson noise perturbation. For quantitative comparison of the quality of the restored objects, we use the Kullback-Leibler divergence (D_KL) of the estimate f _E from the original f

 

Fig. 1. Image restoration on simulated data (qualitative results): (a) phantom (136×136 pixels, pixel size 20×20nm in both lateral and axial direction); (b) cubic RPF; (c) simulated image obtained using the following parameters: objective NA 1.46 oil (RI 1.515), excitation wavelength λ_ex = 632 nm, emission wavelength λ_em = 680 nm, pinhole 0.5 AU, constant background of 2 photon counts and expected photon counts of 100 in the bright pixel; (d) restored object (the RPF is known a-priori); (e) restored object and (f) reconstructed RPF for FPBR; (g) restored object (h) and reconstructed RPF for SPBR.

Download Full Size | PDF

which is the best measure in the presence of a non-negativity constraint [23]. The quality of the reconstructed RPFs was assessed using a modified Euclidean distance (D _E) between the original RPF ϕ and the estimate ϕ _E

where the exponential takes into account for the phase wrapping and T is a mask setting all voxels to zero where the estimated object f _E is below a very low threshold of 10^-6 of its maximum. This ensures that only regions with non-negligible contribution to the image formation are used for the calculation of D _E.

The proposed algorithms were implemented in MATLAB. For comparison, we first tested the performance of the SGM algorithm separately from the AMM algorithm, i.e. we assumed that the RPF is known a-priori and we estimated only the object. Figures 1(a)–1(d) shows the results for μ_f = 10^-5 and demonstrates the well known ability of the SGM algorithm to generate artefact-free estimate of the object and to compensate for phase variation. Moreover, noise was suppressed and structures are well preserved. We then compared these results to the estimates obtained by using the proposed AMM algorithms to simultaneously reconstruct the RPF and the object. Keeping the same μ_f we used μ_ϕ = 50 for FPBR and μ_ϕ = 10^-3, N_ϕ = 4 for SPBR. Figures 1(e)–1(h), and 2 compare the results obtained using FPBR and SPBR: Both approaches generated suitable restored objects and reliably reconstructed RPFs. However, an accurate visual inspection of the restored objects spotlights a slightly superior performance of the SPBR respect to FPBR. Interference artifacts that were partially removed by FPBR were completely removed by SPBR and consequently this approach also reaches better values for the Kullback-Leibler divergence D _KL and for D _E. Because the cost function is convex in f for a given RPF the difference between the two results is rooted in their different ability to approach the global optimum for the RPF. Since we had assumed a polynomial RPF the superior performance of SPBR is not surprising given the fact that the number of free parameters is much lower and the model allows for a perfect representation of the RPF using the basis functions. Please also note that the algorithms cannot be expected to yield a good estimate of the RPF in regions far away from the object structures; a fact also reflected in the calculation of D _E.

An important aspect is the difference of convergence rates for the two methods. As seen from Fig. 2 the reduction of free parameters in the SPBR approach leads to convergence after less cycles of the AMM as expected. For large data sets the calculation of the updated object will usually dominate due to the 𝒪(Nlog(N)) complexity of the FFTs and therefore this advantage directly translates into faster convergence in terms of calculation time. In our case the SPBR approach was usually twice as fast. However we note that since the two approaches implement different algorithms for the update of the RPF with different regularization parameters and stopping criteria, their speed is co-determined by the choice of parameters and the time needed for each RPF calculation. The latter is of complexity 𝒪(N) for APBR and 𝒪(N_ψN) for SBPR. For smaller images and/or (much) larger dictionaries this advantage will therefore become smaller and eventually vanish.

 

Fig. 2. Image restoration on simulated data (quantitative results): (a) Kullback-Leibler divergences between test phantom and reconstructed object as a function of the number of cycles for FPBR vs. SPBR; (b) Euclidean distances between reconstructed and assumed RPF as a function of the number of cycles for FPBR and SPBR.

Download Full Size | PDF

We also note that the local quality of the restoration may slightly vary with the RPF, an effect rooted in the dependence of the frequency content of the PSF on the local phase. This very small effect can be completely eliminated by obtaining several images with different relative phases [5].

6. Validation on experimental data

In order to test our method on real data we employed a standard beam scanning confocal microscope (Leica Microsystems, Mannheim, Germany) to which we attached a custom build 4Pi module to mount two opposing objective lenses (HCX PL APO 100×, 1.46 NA OIL CORR; Leica Microsystems). A driven-piezo mirror controls the difference in pathlength between both cavity arms. Three-dimensional imaging was performed on Vero cells. Cells were immunola-belled with a monoclonal mouse antiserum directed against the alpha-tubulin (Sigma-Aldrich St. Louis, MO, USA), the building-block of microtubules. The primary antibodies were detected with secondary antibodies (sheep anti-mouse; Jackson ImmunoResearch Laboratories, West Grove; PA; USA) custom labelled with the fluorophore KK114. Samples were mounted in the TDE (2,2’-Thidodiethanol, Sigma,-Aldrich, St. Louis, MO) water-miscible mounting medium, in order to adjust the refractive index of the sample to ~1.5 [9]. To provide a more challenging test situation, index matching was not performed completely, so that a mismatch between nucleus and medium remained. Additionally, a linear phase shift along z was introduced by changing the relative length of the two cavity arms while scanning axially. Parameters were choosen such that constructive interference on the top of the xz slice [Fig. 3(a), region of interest 1] and destructive interference on the bottom [Fig. 3(a), region of interest 2] were observed, rendering the RPF a combination of externally-induced and sample-induced phase shifts. To obtain a genuine 3D structure with a thickness of a few micrometers, imaging of microtubules wrapping the nucleus was performed [Fig. 3(a)]. Fluorescence after one-photon excitation of the fluorophore KK114 at 635nm was detected over the wavelength range from 660nm to 700nm by a photon-counting avalanche photodiode (APD) (PerkinElmer, Waltham, MA, USA). The confocal pinhole was set to 0.5 AU. We started by validating the method on a single 2D axial image. We used μ_f = 10^-4, μ_ϕ = 25 for FPBR and μ_ϕ = 10^-3, N_ϕ = 3 for SPBR. A constant background of 2 counts was estimated directly from the images.

 

Fig. 3. Image restoration on experimental data: (a) raw xz image of microtubules wrapping around the nucleus in a Vero cell (186 × 352 pixels); (b) deconvolved data assuming constant RPF; (c) restored object and (d) reconstructed RPF for FPBR; (e) restored object and (f) reconstructed RPF for SPBR. Scale bar 1μm.

Download Full Size | PDF

 

Fig. 4. Intensity profiles of the restored objects taken at the positions indicated by arrows in Fig. 3(a). (a) intensity profile at site 1; (b) intensity profile at site 2.

Download Full Size | PDF

Figure 3(b) shows the results obtained using the parametric blind deconvolution (PBD) assuming a single, constant phase as devised in our earlier work [5]. Since the algorithm cannot account for the varying phase, it is not able to completely remove the interference artifacts. On the other hand, both PBR methods generate interference artifact-free object estimates. As expected, the reconstructed RPFs show a change from destructive to constructive interference along the axial coordinate. To better appreciate the ability of the method to remove sidelobes Fig. 4 shows intensity profiles from the the sites indicated by arrows in Fig. 3. The artifacts remaining after constant phase PBR are clearly removed by both variable phase PBR approaches. While visual inspection shows some difference between FPBR and SPBR, most prominently in their reconstructed RPFs, the differences in the restored objects are not significative. As before, the main difference between the two methods is the rate of convergence. While SPRB requires ~ 190 cycles to reach the stopping criterion, FPBR requires ~ 450 cycles and is 3 times slower. We also performed SPBR with larger dictionaries (N_ϕ ≤ 8) but did not observe significant improvements.

Finally, we applied our methods to 3D data-stacks. Since SBPR proved to converge faster and yielded comparable or better results in our previous analysis we restricted ourselves to this method. We also used a simple trick to speed up reconstruction in the more complex 3D case. We first applied the SBPR algorithm to a dense xz slice of the 3D data-stack using ε _AMM = 10^-5 and used the resulting coefficient vector ψ to initialize the 3D reconstruction Fig. 5 shows the raw data and the restored object along with the reconstructed RPF. Importantly, also on experimental data, the AMM algorithm converges to appealing results even though the cost function is not convex in ϕ or ψ.

7. Discussion and conclusion

While they feature the highest axial resolution in diffraction limited imaging, microscopes that use the nearly spherical wavefront afforded by two opposing lenses usually suffer from side-lobe artifacts in their raw image data. Reliable image restoration is therefore mandatory for their routine application. In 4Pi-microscopy imaging conditions can often be tuned such that the PSF is space-invariant or its spatial dependence is known, allowing for the application of presently available techniques. However, such adjustments may be a major obstacle or even impossible in many real-life applications. To date the data acquired in such experiments could either not be used at all or had to be reconstructed by trying to manually estimate the RPF and repeatedly performing deconvolutions. Usually the most appealing recovery will then be chosen for further analysis introducing a subjective component into the interpretation of the results. We have addressed this problem by developing a parametric blind restoration method based on a maximum a-posteriori approach which uses as much of the available a-priori information about the noise, the object and the RPF as possible. While our method is not guaranteed to converge to a global-optimum due to the non-convex nature of the resulting cost function our tests have shown that the algorithm gracefully converges to visual appealing restored objects. By assuming a phase function represented by a sparse linear combination of polynomials rather than an arbitrary one we were able to speed up convergence while yielding similar or better results, both on experimental and simulated data. The method proposed here might be improved or accelerated by the use of new minimization algorithms, not using the AMM approach or by initializing the RPF with a more educated guess before starting the iteration. The problem of possible convergence to local minima could be addressed by systematically re-starting the algorithm with different RPFs or by implementing a scheme of simulated annealing. Also, more suitable penalization functions or dictionaries for the relative phase function could emerge. One could also imagine extending the 4Pi microscope to simultaneously map fluorescence and changes in the refractive index, e.g. by detecting differential interference contrast. Using an AMM approach very similar to the one presented here, this additional data could help to significantly improve the reconstruction of the RPF e.g. by imposing local restraints on its derivative. However, even in its present implementation the algorithm performs well on experimental 4Pi data and should be an important step towards more routine application of this method also to live cells by making data acquisition less demanding and faster. Finally a similar approach might prove useful when deconvolving raw data from other microscopy modalities when their E-PSF is position dependent, e.g. in the presence of sample induced abberations. To this end a decomposition of the forward operator could be based on approximation of the E-PSF by a Fourier series to enable fast re-calculation with various forms of the PSF.

 

Fig. 5. 3D Parametric blind restoration: (a) raw data (228 × 334 × 40 voxels); (b) restored object and (c) reconstructed RPF for SPBR. The dotted lines indicate the position of the xz- and yz-slices shown for each of the stacks. Scale bar 1μm.

Download Full Size | PDF