Open Access
Add to BibSonomy
Abstract
To achieve approximately parallel projection geometry, traditional optical projection tomography (OPT) requires the use of low numerical aperture (NA) objectives, which have a long depth-of-field at the expense of poor lateral resolution. Particularly promising methods to improve spatial resolution include ad-hoc post-processing filters that limit the effect of the system’s MTF and focal-plane-scanning OPT (FPS-OPT), an alternative acquisition procedure that allows the use of higher NA objectives by limiting the effect of their shallow depth of field yet still assumes parallel projection rays during reconstruction. Here, we provide a detailed derivation that establishes the existence of a direct inversion formula for FPS-OPT. Based on this formula, we propose a point spread function-aware algorithm that is similar in form and complexity to traditional filtered backprojection (FBP). With simulations, we demonstrate that our point-spread-function aware FBP for FPS-OPT leads to more accurate images than both traditional OPT with deconvolution and FPS-OPT with naive FBP reconstruction. We further illustrate the technique on experimental zebrafish data, which shows that our approach reduces out-of-focus blur compared to a direct FBP reconstruction with FPS-OPT.
© 2017 Optical Society of America
1. Introduction
Optical projection tomography (OPT) is a 3D microscopy technique for imaging small transparent animals (up to a few millimeters) using visible light [1,2]. Similar to x-ray computed tomography (CT), OPT involves the acquisition of multiple 2D projections through a 3D sample, with each projection taken from a different angle of rotation. However, unlike CT, which uses x-rays that travel in approximately straight lines, OPT uses visible light and microscope optics that accept light rays from a range of angles over a cone. In order to apply traditional tomographic reconstruction techniques to OPT, researchers use low numerical aperture (NA) objectives with OPT to reduce the acceptance angle of the system and achieve approximately straight-line projections. However, low NA objectives have worse lateral resolution than high NA objectives, thus limiting the lateral resolution of OPT systems.
The problem with high NA objectives is that they have a shallow depth of field, making them unsuitable for traditional OPT imaging of thick samples. With high NA objectives, part of the sample lying on the focal plane will be in focus, but out-of-focus light from above and below will also contribute to the resulting image. To adapt OPT for use with high NA objectives, previous studies have shown that one can scan the focal plane through the entire sample to create pseudoprojections with an extended depth of field [3,4]. Another previous study has also shown that one can achieve high NA OPT imaging with a selective plane illumination microscope using high-pass filtering and weighted averaging of multiple focal slices to create an extended depth of field projection [5]. In both cases, these pseudoprojections are then used to reconstruct a 3D image using filtered backprojection (FBP) [6].
While the image formation process through convolution with an extended PSF (rather than line integrals) has been acknowledged in OPT and methods have been proposed to mitigate its effect [7–10], they focus on fixed focal plane OPT, which precludes the use of high NA objectives. Other methods for fixed focal plane OPT also include a well-defined model of image formation, but require computationally intensive iterative algorithms for reconstruction [11]. Here, we instead consider the acquisition procedure in [3,4], which we will refer to as focal-plane-scanning OPT (FPS-OPT). For focal-plane-scanning OPT, we derive an analytic inversion formula that fully incorporates the system’s 3D optical point-spread-function. We then show that our inversion can be implemented computationally inexpensively with a series of 2D deconvolutions followed by the traditional filtered backprojection algorithm.
This paper extends on our original work in [12]. While our previous work described our method only in 2D, here we describe our problem formation and reconstruction algorithm for 3D volumetric data. Furthermore, in this paper, we compare our method to other existing methods for OPT reconstruction, and we provide software code and an experimental dataset for reproduction. This paper is organized as follows. In Section 2, we present the image formation process in optical projection tomography. In Section 3, we present our proposed image acquisition and reconstruction approach. In Section 4.1, we characterize our method with a 2D phantom and point-spread-functions of various NAs. In Section 4.2, we demonstrate our method on 3D data acquired with FPS-OPT of a Tg(fli1a:EGFP) zebrafish embryo.
2. Problem formulation
Assuming minimal optical attenuation and scattering, the ideal image formation process in optical projection tomography can be described by the X-ray transform (or Radon transform in 2D) of the underlying object. For a 3D object f(x, y, z), assuming the axis of rotation coincides with the z-axis, the observed projection at a particular angle θ is the set of line integrals (Fig. 1(a)),
where δ(·, ·) is the 2D Dirac delta function defined as δ(a, b) = δ(a) · δ(b), and p(u, v, θ) is the 2D projection of f(x, y, z) along the direction θ. However, this model is only accurate when (i) the entire sample is contained in the system’s depth of field, and (ii) the system’s optical point spread function (PSF) is infinitely narrow in the lateral (x–y) direction. In OPT, this model is inaccurate, particularly the second condition. The image formation model in OPT is more accurately described by a single focal plane sampled from the result of the convolution between the object and a rotated version of the system’s PSF (Fig. 1(b)), where Tθ{·} is a transformation operator that rotates a function by an angle θ about the z-axis, h is the point-spread-function of the system, and u and v are distances from the origin (which we assume to lie on the focal plane). In such a case, the standard filtered backprojection (FBP) algorithm can no longer be used to accurately recover f(x, y, z).
Fig. 1 (a) In ideal tomography, acquisitions are assumed to be straight line projections through the entire object. Here, a projection in the direction parallel to a vector θ̂ = shown for an OPT system with an axis of rotation parallel to the z-axis. (b) In traditional OPT, acquisitions are better described as the result of the 3D convolution between the object and the PSF, sampled on a single focal plane. (c) In FPS-OPT, each acquisition is obtained by summing up images obtained while scanning through all focal planes.
3. Proposed method
3.1. Increasing the depth of field
To address the first problem when the sample is not completely contained in the system’s depth of field, we follow the procedure described in [4] to create a pseudoprojection by scanning the focal plane through the entire sample (FPS-OPT). This can be done either by linearly sweeping the focal plane through the sample during a single camera integration period, or by acquiring a focal stack and retrospectively integrating along the projection direction by digitally summing the slices. By scanning the focal plane during imaging, we manually perform an X-ray transform as in Fig. 1(c),
However, the imaging model still differs from the ideal model in Fig. 1(a), because for each angle θ, the X-ray transform is taken of a differently blurred underlying image due to the rotated point-spread-function Tθ{h}.3.2. Deblurring with PSF-Aware filtered backprojection
Applying a standard filtered backprojection to Equation (3) will not allow us to recover f(x, y, z). Rather, it will reconstruct a blurred version of f(x, y, z) due to the effect of the optical point-spread-function h. However, we recall that according to the Fourier slice theorem, the M-dimensional Fourier transform of a projection (X-ray transform) of an N-dimensional function onto M dimensions is equivalent to an M-dimensional slice of that function’s N-dimensional Fourier transform [13],
where ℱM denotes an M-dimensional Fourier transform, denotes a projection (X-ray transform) at an angle θ from N dimensions to M dimensions, denotes an M-dimensional slice normal to the direction of θ from N dimensions, and ℱN denotes an N-dimensional Fourier transform. For consistency with Fig. 1, the remainder of this section considers the case where M = 2 and N = 3 for 3D OPT with a single component of rotation (θ) about the z-axis. In this scenario, the projection operator is a transformation from a three-dimensional x-y-z coordinate system to a two-dimensional u–v coordinate system.As the X-ray transform in Equation (3) is a projection operator , by the Fourier-slice theorem, its Fourier transform is equivalent to a slice from the Fourier transform of the PSF-blurred image (as illustrated in Fig. 2),
Since the convolution between two functions translates to the product of their Fourier transforms in the frequency domain, and since the point-wise product can be performed after the slicing operator, we can write: Applying the Fourier-slice theorem once more to the PSF-term in the product allows us to separate the Fourier-slice of the desired, underlying function via filtering of the observed projection with the inverse of the PSF’s X-ray transform:At this point, we recall that the filtered back projection (FBP) algorithm [6] can be used to invert the X-ray transform in Equation (3) by expressing the reconstructed image as:
where Qθk (u, v) is the projection at an angle θk after Fourier domain filtering with a ramp filter W(ωu, ωv) = |ωu|, Substituting (7) in the above, each filtered backprojection can be expressed in terms of the FPS-OPT-measured blurred projection p̃(u, v, θ). Specifically, we use where W̃(ωu, ωv) is the product of the ramp filter W(ωu, ωv) and a (regularized) version of the expression in Equation (7), for where λ is the regularization weighting parameter and r is a high-pass regularization filter. This regularization term in the denominator stabilizes the inverse filter and prevents it from growing large when the optical transfer function is close to zero. We use a 2D Laplacian as the regularization filter, which can be implemented in practice as a discrete 2D filter utilizing finite differences. While similar in form to a Wiener filter, this filter does not depend on knowledge of the signal or noise power spectrum. With this result, we can now apply the modified FBP in (8), to estimate the underlying image f(x, y, z) from the blurred projections p̃(u, v, θ).Fig. 2 The Fourier transform of an object’s projection is equivalent to a slice from the underlying object’s Fourier transform. The inclusion of the system’s point spread function introduces an extra convolution in space, or equivalently, an extra multiplication in the Fourier domain.
4. Experiments
While the modified filtered backprojection in Equation (10) can be implemented with a single filter W̃(ωu, ωv), in practice we separate this filter into the traditional ramp filter and the PSF inverse filter, as shown in Fig. 3. This allows us to take advantage of existing FBP packages that use traditional ramp filtering. In our experiments, we use the spline-based filtered backprojection in [14] to implement the inverse X-ray transform.
Fig. 3 In practice, our PSF-aware filtered backprojection algorithm is implemented as two-step process consisting of an initial PSF inverse filtering followed by traditional filtered backprojection.
4.1. Simulations with Shepp-Logan phantom
To evaluate our method, we simulated an optical projection tomography system in Matlab with a 2D Shepp-Logan phantom as a test image. In this simulation, we acquired 180 uniformly spaced projections over a total of 180° degrees (for an angular spacing of 1°). To model the image formation process, we used 2D point-spread-functions obtained by extracting the central y–z plane from 3D PSFs generated with the Born & Wolf model [15,16]. We compared image reconstruction quality in terms of peak signal-to-noise ratio, PSNR = 10 log10(max f̃)2/MSE where f̃ is the reconstructed image, and MSE is the mean squared error between f̃ and the original ground truth image.
We compared two approaches for data acquisition, traditional fixed-focal-plane OPT and focal-plane-scanning OPT (FPS-OPT). As our baseline, we reconstructed the acquired projections with the standard filtered backprojection (FBP) algorithm in both scenarios. Next, we reconstructed the acquired OPT and FPS-OPT projections with algorithms specifically designed for each acquisition procedure: the MTF deconvolution method [8] for traditional OPT, and our proposed PSF-aware FBP method for FPS-OPT. For the MTF deconvolution reconstruction, we implemented the Wiener filter deconvolution as described in [8], with the exception of selecting Su = 0.5, which was empirically chosen to produce a reconstruction with a high PSNR. For our PSF-aware filtered backprojection, we performed each simulation with multiple values of λ and selected the λ that produced the reconstruction with the best PSNR.
Fig. 4 Using point spread functions generated with numerical apertures of (a,f) 0.3 and (k,p) 0.5, we compared the following reconstruction algorithms: (b,g,l,q) traditional FBP under fixed-plane OPT, (c,h,m,r) MTF deconvolution [8] under fixed-plane OPT, (d,i,n,s) traditional FBP under FPS-OPT, and (e,j,o,t) our proposed PSF-aware FBP under FPS-OPT. (b–e) Reconstructions under the 0.3 NA PSF in (a,f) with a small phantom and no noise. (g–j) Reconstructions under the 0.3 NA PSF in (a,f) with a phantom larger than the depth-of-field and no noise. (l–o) Reconstructions under the 0.5 NA PSF in (k,p) with a phantom much larger than the depth-of-field and no noise. (q–t) Reconstructions under the 0.5 NA PSF in (k,p) with a phantom much larger than the depth-of-field and with projections corrupted by shot noise. Specifically, for each observed projection, the measured value y at any given detector followed a Poisson distribution y ∼ Poisson(x), where x was the deterministic (noise-free) projection value at that detector, scaled so all deterministic projection values were in the interval [0, 104].
We repeated the acquisition and reconstruction procedures in simulated OPT systems with various numerical apertures and object sizes (one per row of Fig. 4 and Table 1). We used the same sampling resolution (magnification factor) for each numerical aperture, which we chose to be 100 nm per pixel in both directions. We evaluated our PSF-aware filtered backprojection under this scenario without noise (Fig. 4(a–o)), and in each case, our method (which assumes focal-plane-scanning OPT) outperformed standard FBP under both fixed-plane and focal-plane-scanning OPT, as well as the fixed-plane MTF deconvolution method in [8] (Table 1). When ignoring the circular background halo artifact, the MTF deconvolution approach outperformed simple FBP when the phantom was smaller than the PSF’s axial width. However, as the depth-of-field decreases with increasing NA, this approach is no longer able to produce an accurate reconstruction. In contrast, our approach, which combines deconvolution with focal-plane-scanning geometry, is able to produce a good reconstruction in all scenarios. We also repeated this comparison for projections acquired with Poisson-distributed shot noise (Fig. 4(p–t)). In the presence of noise, our approach is still able to produce a good reconstruction. However, as the noise variance increases, we must choose a larger value of the regularization parameter λ to stabilize the inverse filter and avoid noise amplification, leading to a weaker overall deblurring.
Table 1. PSNR comparisons between traditional FBP, the MTF deconvolution reconstruction in [8], and our PSF-aware FBP with different numerical apertures. We computed both the overall PSNR of the entire image, as well as the PSNR for the Shepp-Logan phantom itself (excluding the background pixels). We performed all simulations using the full size phantom as in Fig. 4 (f–o), except where indicated by †, where we used the smaller phantom as in Fig. 4 (a–e). In addition, we performed all simulations assuming noise-free conditions, except where indicated by ★, where we assumed the projections were taken with Poisson-distributed shot noise as in Fig. 4 (q–t). To generate data corrupted by shot noise, we assumed that the measured value y at any given detector followed a Poisson distribution y ∼ Poisson(x), where x was the deterministic (noise-free) projection value at that detector, scaled so all deterministic projection values were in the interval [0, 104].
4.2. Application to zebrafish imaging
To demonstrate our approach in practice, we imaged the head of a 62 hpf (hours post-fertilization) Tg(fli1a:eGFP) zebrafish that expresses green fluorescence in its blood vessels following a protocol approved by the UCSB Institutional Animal Care and Use Committee. We acquired 1600 projections over 360° (with 8 focal slices per projection, over a total depth of 1000 μm) using a custom-built rotational stage consisting of a stepper motor connected to a fluorinated ethylene propylene (FEP) tube that has a refractive index close to that of water. We mounted the larval zebrafish inside the FEP tube in a 1% low melting-point agar solution. We placed this rotational stage on a Leica DMI6000B inverted widefield microscope and imaged in fluorescence with a 10 × /0.3 dry objective and a Hamamatsu C9100-13 EM-CCD camera. Figure 5 shows a sketch of the image acquisition setup and rotation orientation.
Fig. 5 Diagram of the OPT rotational acquisition procedure for zebrafish. For each projection angle θ, the focal plane is scanned through the entire zebrafish to create a full projection, even with a shallow depth-of-field objective.
At each angle, we acquired a scanned projection of the zebrafish fluorescence emission using FPS-OPT, and we computed the projection’s center of mass to determine the position of its rotational axis [17]. More advanced methods may also be used to correct for any sample movement or rotational drifts [18]. However, we found this unnecessary for our experimental data. From these projections, we reconstructed a 3D volume using both traditional FBP and our PSF-aware FBP (Fig. 6). Our PSF-aware FBP reconstruction contains significantly less out-of-focus blur compared to the traditional FBP reconstruction.
Fig. 6 We used focal plane scanning OPT (FPS-OPT) to image the head of a Tg(fli1a:eGFP) zebrafish in 3D fluorescence with a 10 ×/0.3 NA air objective. Under such conditions, single-plane OPT would be unable to produce an acceptable reconstruction due to the large sample thickness and shallow depth-of-field. With FPS-OPT, we compare 3D reconstructions from (a) standard filtered backprojection (FBP) and (b) our proposed PSF-aware FBP. Our proposed PSF-aware FBP algorithm reconstructs an image with less out-of-focus blur. Scalebar is 100 μm.
5. Discussion and conclusion
For the same magnification factor, a higher NA system allows for better resolution than a lower NA system due to the fact that the lateral width of the system’s point-spread-function, as defined by its full width at half maximum (FWHM), is inversely proportional to the system’s NA. However, since the axial width of a point-spread-function is inversely proportional to the square of the system’s NA, a high NA system also has a much shallower depth-of-field than a low NA system. As a result, fixed-plane OPT with a high NA is incompatible with thick samples when the depth-of-field is shallower than the sample thickness (Fig. 4). When the phantom size is larger than the PSF’s axial width, such as with an NA of 0.5 in our simulation, fixed-plane OPT is unable to produce an accurate reconstruction. Also, as the PSF’s lateral width increases with lower NA, the inverse filter also becomes more numerically unstable, and the reconstructed image contains more visible ringing artifacts. While we use direct regularized inverse filtering to implement deconvolution, more advanced iterative methods may produce better results with fewer artifacts, though at the expense of greater computational complexity and longer processing time [19].
In conclusion, we have derived a direct inversion algorithm for focal plane scanning optical projection tomography (FPS-OPT) that can be implemented using a modified filtered backprojection algorithm that incorporates knowledge of the system’s point-spread-function. Through simulations with a 2D phantom, we showed that our modified filtered backprojection offers a noticeable improvement in PSNR compared to both traditional FBP as well as the MTF deconvolution method for single focal-plane OPT [8], especially when the object size is large relative to the PSF’s depth-of-field. To demonstrate our method in practice, we used fluorescence FPS-OPT to image the head of a larval zebrafish, and we showed that our PSF-aware FBP algorithm is able to reconstruct an 3D volume with significantly reduced out-of-focus blur compared to a traditional FBP reconstruction. To allow reproducibility, we provide both our FPS-OPT dataset, consisting of 1600 scanned projections of a Tg(fli1a:EGFP) zebrafish head over 360°, as well as a Matlab implementation of our reconstruction algorithm (available at http://sybil.ece.ucsb.edu/pages/fpsopt/).
Funding
Swiss National Science Foundation: 200021_159227 and 206021_164022.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References and links
1. J. Sharpe, U. Ahlgren, P. Perry, B. Hill, A. Ross, J. Hecksher-Sørensen, R. Baldock, and D. Davidson, “Optical projection tomography as a tool for 3D microscopy and gene expression studies,” Science 296(5567), 541–545 (2002). [CrossRef] [PubMed]
2. R. J. Bryson-Richardson, S. Berger, T. F. Schilling, T. E. Hall, N. J. Cole, A. J. Gibson, J. Sharpe, and P. D. Currie, “FishNet: an online database of zebrafish anatomy,” BMC Biol. 5(1) 34 (2007). [CrossRef] [PubMed]
3. M. Fauver, E. J. Seibel, J. R. Rahn, M. G. Meyer, F. W. Patten, T. Neumann, and A. C. Nelson, “Three-dimensional imaging of single isolated cell nuclei using optical projection tomography,” Opt. Express 13(11), 4210–4223 (2005). [CrossRef] [PubMed]
4. Q. Miao, J. Hayenga, M. G. Meyer, T. Neumann, A. C. Nelson, and E. J. Seibel, “Resolution improvement in optical projection tomography by the focal scanning method,” Opt. Letters 35(20), 3363–3365 (2010). [CrossRef]
5. A. Bassi, B. Schmid, and J. Huisken, “Optical tomography complements light sheet microscopy for in toto imaging of zebrafish development,” Development 142(5), 1016–1020 (2015). [CrossRef] [PubMed]
6. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE1988).
7. J. R. Walls, J. G. Sled, J. Sharpe, and R. M. Henkelman, “Resolution improvement in emission optical projection tomography,” Phys. Med. Biol. 52(10), 2775–2790 (2007). [CrossRef] [PubMed]
8. L. Chen, J. McGinty, H. B. Taylor, L. Bugeon, J. R. Lamb, M. J. Dallman, and P. M. W. French, “Incorporation of an experimentally determined MTF for spatial frequency filtering and deconvolution during optical projection tomography reconstruction,” Opt. Express 20(7), 7323–7337 (2012). [CrossRef] [PubMed]
9. J. Guo, Y. Yang, D. Dong, L. Shi, H. Hui, M. Xu, J. Tian, and X. Liu, “A projection selection method to improve image quality in optical projection tomography,” in International Conference of the IEEE Engineering in Medicine and Biology Society (IEEE, 2014), pp. 206–209.
10. J. van der Horst and J. Kalkman, “Image resolution and deconvolution in optical tomography,” Opt. Express 24(21), 24460–24472 (2016). [CrossRef] [PubMed]
11. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comput. Imag. 2(1), 59–70 (2016). [CrossRef]
12. K. G. Chan and M. Liebling, “A point-spread-function-aware filtered backprojection algorithm for focal-plane-scanning optical projection tomography,” International Symposium on Biomedical Imaging (IEEE, 2016), pp. 253–256.
13. Z. P. Liang and D. C. Munson Jr., “Partial Radon transforms,” IEEE Trans. Image Process. 6(10), 1467–1469 (1997). [CrossRef] [PubMed]
14. S. Horbelt, M. Liebling, and M. Unser, “Discretization of the Radon transform and of its inverse by spline convolution,” IEEE Trans. Med. Imag. 21(4) 363–376 (2002). [CrossRef]
15. M. Born and E. Wolf, Principles of Optics (Cambridge Univ. Press, 2003), Chap. 8.
16. H. Kirshner, F. Aguet, D. Sage, and M. Unser, “3-D PSF fitting for fluorescence microscopy: implementation and localization application,” J. Microsc. 249(1), 13–25 (2013). [CrossRef]
17. S. G. Azevedo, D. J. Schneberk, J. Fitch, and H. E. Martz, “Calculation of the rotational centers in computerized tomography sinograms,” IEEE Trans. Nucl. Sci. 37(4), 1525–1540 (1990). [CrossRef]
18. U. J. Birk, M. Rieckher, N. Konstantinides, A. Darrell, A. Sarasa-Renedo, H. Meyer, N. Tavernarakis, and J. Ripoll, “Correction for specimen movement and rotation errors for in-vivo Optical Projection Tomography,” Biomed. Opt. Express 1(1), 87–96 (2010). [CrossRef]
19. D. Sage, L. Donati, F. Soulez, D. Fortun, G. Schmit, A. Seitz, R. Guiet, C. Vonesch, and M. Unser, “DeconvolutionLab2: An open-source software for deconvolution microscopy,” Methods 115, 28–41 (2017). [CrossRef] [PubMed]
