Image resolution and deconvolution in optical tomography

Jelle van der Horst; Jeroen Kalkman

doi:10.1364/OE.24.024460

1. Introduction

Three dimensional imaging of cells and small animals is of great importance in the biological and medical sciences as our knowledge on human diseases is based on the structure and function of cells and organisms. Optical techniques can provide cost effective, high-resolution, 3D images of the structure of cells and small animals. In addition, by using labelled molecules, fluorescence imaging can give information about the function of cells and tissue.

Currently, several optical techniques are used for whole cell and small animal imaging. Microscopy is a tool for imaging live cells. Light sheet microscopy uses a light illumination sheet and images the emission from perpendicular to the sheet. It is used for imaging sub-millimeter sized samples [1]. 3D imaging studies of millimeter sized samples have used optical projection tomography (OPT) [2]. OPT measures projections of transmission or emission and reconstructs the object from the projections. OPT is used for imaging small animals such as zebrafish and mouse embryos, whole organs taken from adult mice, and plant tissue. It can be used in-vivo, or ex-vivo in combination with optical clearing techniques that suppress light scattering. OPT has the advantage of being able to create 3D images in both transmission and fluorescent modes, hence OPT is useful to study gene expression, tissue morphology and locations of fluorescently labelled tissues. In transmission mode, OPT resembles X-ray CT, in fluorescent mode, OPT is more similar to single-photon emission computed tomography (SPECT) imaging. In both cases optical radiation is used instead of X-rays or gamma quanta.

In contrast to X-ray CT and SPECT imaging, the use of visible light in OPT allows for the use of lenses to relay and magnify the projection images and achieve micrometer spatial resolution. A drawback of the use of lenses is the limited depth sensitivity, related to the finite depth of field of the optical imaging system. Generally, the resolution of the images is inversely proportional to the numerical aperture (NA) of the collecting lens. However, high NA lenses have a small depth of field, meaning that objects are only in focus in a small region around the focal plane. Using high NA lenses in OPT causes some parts of the sample to have a high resolution in a projection, while other parts are out of focus and have a low resolution.

Several studies have identified this problem and proposed methods to reduce its effect. The depth of focus can be extended by focusing on a plane between the center of rotation and the edge of the sample and recording projections over the full range of 360 [3]. Alternatively, data from multiple focal planes can be combined by either simultaneous measurement [4] or scanning of the lens focus through the object [5]. While recording data from multiple focal planes solves the issue of the depth of field, the increased number of measurements requires longer acquisition times and increased complexity of the OPT system.

Several numerical approaches have been proposed to correct for resolution blurring in OPT. Based on the frequency distance relationship (FDR) [6], OPT tomograms can be corrected for the out of focus deterioration of resolution [7]. The effects of the axial intensity distribution and defocus on the point spread function (PSF) can also be compensated by using a weighted filtered backprojection [8]. More recently, the full modulation transfer function (MTF) of the imaging optics was included as an additional filter in the filtered backprojection process [9]. While obtaining a significant reduction in image background and artefacts using an MTF mask, the MTF correction with deconvolution did not fully correct for the tangential blurring observed in the system. McErlean et al. [10] investigated a possible spatial resolution improvement by image deconvolution with an experimentally determined PSF. However, they used a spatially-invariant PSF thereby obtaining less than optimum resolution.

The aim of this study is to quantify the effect of the imaging optics on the reconstructed images in optical tomography. In the first section we quantify this effect using a frequency domain analysis of the image resolution in a single projection and in the reconstructed tomographic image. In the second section we propose and verify our model for the spatially-variant resolution of the tomographic imaging system using OPT measurements of a sample consisting of fluorescent beads. Then, we use the derived point spread function to deconvolute the reconstructed image. A close to pixel-limited resolution after image deconvolution is demonstrated for the fluorescent bead sample. We also apply our method to a zebrafish sample. Finally, our results, and their implications for optical tomography, are discussed.

2. Theory

In an optical tomography system, the spatially-variant PSF is directly related to the PSF of the imaging system that makes every projection. Here we use a Fourier optics description of the PSF of the imaging system to derive the PSF of the tomography system. Our analysis is demonstrated for fluorescence tomography, but is equally valid for transmission optical tomography.

2.1. Image formation of a single projection

Consider an object with fluorescence distribution f(t, s, z). In fluorescent tomography the light emitted by the object is focussed onto a 2D detector to create projections of the object, as shown in Fig. 1. We assume that all fluorophores are excited at the same rate and emit isotropically. In addition, it is assumed that any emitted fluorescence can reach the detector unimpeded. Hence, effects of reabsorption, light attenuation and refraction of fluorescent light are not taken into account. The intensity in the image space is described by a convolution of the object function, f(t, s, z), with the incoherent PSF of the imaging optics |h(t, s, z)|², where h(t, s, z) describes the coherent PSF of the imaging system [11]. The system images the plane s = 0 in object space onto the plane s′ = 0 in image space. The intensity distribution in image space p(t′, s′, z′) is given by a 3D convolution with the PSF

p (t^{'}, s^{'}, z^{'}) = f (t, s, z) \otimes {| h (t, s, z) |}^{2} .

If a detector is placed in the focal plane s′ = 0, the intensity distribution on the detector is

p (t^{'}, z^{'}) = {\int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} f (t^{'} - t, s^{'} - s, z^{'} - z) {| h (t, s, z) |}^{2} d t d s d z |}_{s^{'} = 0} .

This equation shows the main difference between OPT and straight ray based tomography. If |h(t, s, z)|² = δ(t, z) Eq. 2 reduces to the line integral along s for the intensity of a single point in the projection. OPT deviates from straight ray based tomography as the measured projection points do not only sample a straight line from the source to the detector, but instead sample a complex volume of the object, described by the PSF |h(t, s, z)|².

Fig. 1 Schematic representation of the lens imaging system for the detection of a single projection

Download Full Size | PDF

2.2. Tomographic point spread function

In tomographic imaging a collection of projections, acquired at different angles θ are used to construct a tomographic image of the object. Figure 2(a) shows a top view of the object plane with the coordinate systems used in the derivation. In the following analysis we assume that the rotation axis coincides with the z-axis, i.e., the rotation axis is in the focal plane, and the detector is positioned at s′ = 0. In the reference frame of the detector (t′, s′, z′), a projection p at angle θ is given by Eq. (2). If the projection is taken at angle θ, this is equivalent to a rotation of the object f through angle −θ. The relation between the object coordinates in the rotating frame (x, y, z) and the stationary (detector) frame of reference (t, s, z) is given by

[\begin{matrix} t \\ s \\ z \end{matrix}] = [\begin{matrix} cos (- θ) & - sin (- θ) & 0 \\ sin (- θ) & cos (- θ) & 0 \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} x \\ y \\ z \end{matrix}] .

The resolution of the tomographic imaging system is calculated by determining the response to a point object function, i.e., f(x, y, z) = δ(x − x_o, y − y_o, z − z_o) = δ(t − t_o, s − s_o, z − z_o), as shown in Fig. 2(a). Without loss of generality we assume the point object to be placed in the plane z_o = 0 hence the PSF from Eq. (2) is

p (θ, t^{'}, z^{'}) = {| h (t^{'} - t_{o}, - s_{o}, z^{'}) |}^{2} .

Every individual horizontal slice of the object is reconstructed from horizontal sections of all the 2D projections. These horizontal sections correspond to one or several rows of detector pixels at the corresponding vertical position of the slice.

p (θ, t^{'}) = \int_{Δ z - δ z}^{Δ z + δ z} {| h (t^{'} - t_{o}, s_{o}, z) |}^{2} d z,

with Δz the vertical distance between the point source and the slice, and 2δz the slice thickness. The set of all 1D projections is called a sinogram.

Fig. 2 (a) Schematic representation of the coordinate systems in the object and camera frames of reference. The point object is represented by coordinates (x_o, y_o), (b) Frequency domain representation of the reconstructed image for a point object. The frequency content is composed of the frequency content of the individual projections, as indicated for a single projection at angle θ. (c) Reconstructed image resulting from a point object

Download Full Size | PDF

To obtain the frequency description of the projection we take the Fourier transform of Eq. (5) along the transverse coordinate t, to obtain

P (θ, f_{t}) = ℱ {p (θ, t^{'})} .

Assuming a parallel-beam geometry the projection slice theorem states that the 2D Fourier transform of the image of the object, denoted with OTF(f_x, f_y), is composed of the frequency content of the projections P(θ, f_t) [12]. Each projection P(θ, f_t) forms the radial cross-section of OTF(f_x, f_y) at angle θ, as shown in Fig. 2(b). The OTF is calculated from Eq. (6) by substitution of

s_{o} = s_{o} (θ) = r_{o} sin (θ_{o} - θ),

θ = {tan}^{- 1} (\frac{f_{y}}{f_{x}}),

r_{o}^{2} = x_{o}^{2} + y_{o}^{2},

and

f_{t} = \sqrt{f_{x}^{2} + f_{y}^{2}} .

Subsequent 2D inverse Fourier transform of OTF(f_x, f_y) yields the position dependent PSF, as illustrated in Fig. 2(c).

2.2.1. Gaussian beam PSF

The PSF of the tomographic imaging system can be calculated from the PSF of the imaging optics |h(t, s, z)|². In general, |h(t, s, z)|² can have a complicated shape. Here, we assume the imaging PSF is a focused Gaussian beam described by

{| h (t, s, z) |}^{2} = \frac{2}{π w^{2} (s)} exp (- \frac{2 (t^{2} + z^{2})}{w^{2} (s)}),

where w(s) is the beam waist, given by

w (s) = \sqrt{w_{0}^{2} + \frac{λ^{2} s^{2}}{π^{2} w_{0}^{2}}},

with λ the emission wavelength of the fluorophore and w₀ the Gaussian beam waist defined as the e⁻¹ value of the field amplitude in focus [13]. Substituting Eq. (11) into Eq. (5) and performing the integration in the z direction over the height of the detector row results in

p (θ, t^{'}) = \sqrt{2 π} \frac{A (θ)}{π w^{2} (s_{o} (θ))} exp (- \frac{2 {[t^{'} - t_{o}]}^{2}}{w^{2} (s_{o} (θ))}),

with

A (θ) = \frac{1}{2} [erf (\frac{\sqrt{2}}{w^{2} (s_{o} (θ))} (Δ z + δ z)) - erf (\frac{\sqrt{2}}{w^{2} (s_{o} (θ))} (Δ z - δ z))] .

Fourier transformation in the transversal t′ direction gives

P (θ, f_{t}) = A (θ) exp (- \frac{π^{2} w^{2} (s_{o} (θ)) f_{t}^{2}}{2}) exp (- 2 π i t_{o} f_{t}) .

When the slice is thick with respect to the spotsize of the point source on the detector, we can take the limit δz → ∞, resulting in A(θ) → 1. Using the relations (7)–(10), the OTF of the tomographic system is

OTF (f_{x}, f_{y}) = exp (- \frac{π^{2} (f_{x}^{2} + f_{y}^{2})}{2} [w_{0}^{2} + \frac{λ^{2} r_{o}^{2} {sin}^{2} (θ_{o} - θ)}{π^{2} w_{0}^{2}}]) exp (- 2 π i (x_{o} f_{x} + y_{o} f_{y}))

The last exponential term in Eq. (16) is a phase term resulting from the shift of the point object (x_o, y_o) in real space. Equation (16) can be simplified by introducing a shifted and rotated coordinate system (u, v), where u and v respectively represent the radial and tangential directions in the reconstructed image, as shown in Fig. 2(c) and given by

[\begin{matrix} u \\ v \end{matrix}] = [\begin{matrix} cos (- θ_{o}) & - sin (- θ_{o}) \\ sin (- θ_{o}) & cos (- θ_{o}) \end{matrix}] [\begin{matrix} x - x_{o} \\ y - y_{o} \end{matrix}],

with the frequency domain counterpart (f_u, f_v). This centers the point source on the origin in the real domain and removes the Fourier-shift term in the OTF, giving

OTF (f_{u}, f_{v}) = exp (- \frac{1}{2} π^{2} (f_{u}^{2} + f_{v}^{2}) [w_{0}^{2} + \frac{λ^{2}}{π^{2} w_{0}^{2}} (r_{o}^{2} {sin}^{2} (θ_{u}))]),

where θ_u = θ_o − θ now denotes the angle with the u axis. Further simplification using

{sin}^{2} (θ_{u}) = \frac{v^{2}}{u^{2} + v^{2}} = \frac{f_{v}^{2}}{f_{u}^{2} + f_{v}^{2}}

reduces the expression to

OTF (f_{u}, f_{v}) = exp (- \frac{1}{2} π^{2} [(f_{u}^{2} + f_{v}^{2}) w_{0}^{2} + \frac{λ^{2}}{π^{2} w_{0}^{2}} (r_{o}^{2} f_{v}^{2})]),

Equation (19) is separable into two factors with dependence on u and v only, and is further simplified to

OTF (f_{u}, f_{v}) = \frac{1}{2 π w_{0}^{2}} exp (- π^{2} (f_{u}^{2} a_{u} + f_{v}^{2} a_{v})),

with

a_{u} = w_{0}^{2} / 2

and

a_{v} = (w_{0}^{2} + \frac{λ^{2} r_{o}^{2}}{π^{2} w_{0}^{2}}) / 2 = w^{2} (r_{o}) / 2

. Taking the inverse Fourier transform of Eq. (20) we obtain the PSF of the tomographic imaging system

PSF (u, v) = \sqrt{\frac{1}{π^{2} a_{u} a_{v}}} exp (- [\frac{u^{2}}{a_{u}} + \frac{v^{2}}{a_{v}}]),

which is the main theoretical result of this section. Figure 3(a) shows the modelled optical tomography PSF at several positions in the reconstructed image. The example is calculated with w₀ = 10 μm and λ = 515 nm. The PSF in the center of the image is isotropic and has a width w₀ in all directions. The effect of the lens system PSF away from the center results in a radial resolution of w₀ independent of radial position, however the tangential resolution deteriorates as the radial position of the object increases. Figure 3(b) shows a quantification of this effect as the tangential resolution increases with the radial distance from the center.

Fig. 3 (a) Effect of the PSF on the tomographic image for point sources located at various positions (x_o, y_o), indicated at the top and right. (b) Theoretical FWHM resolution in the radial and tangential direction as function of the radial distance from the center of rotation.

Download Full Size | PDF

3. Methods

3.1. Experimental setup

A schematic of the fluorescence OPT setup is shown in Fig. 4. A collimated beam of light with a wavelength of 488 nm is created using an Argon laser (150m Select, Laser Physics). The output power of the laser is adjusted using a neutral density filter wheel (NDC-100C04M, Thorlabs). The beam is expanded with a 4f lens system consisting of two lenses with a focal length of 10 mm (LB1157, Thorlabs), and 1000 mm (LB1859, Thorlabs). The sample is placed in a cuvette (Hellma, HELL704001-30-10) with inner dimensions 30×30×30 mm (length × width × height). After passing through a color filter, the fluorescent light is detected using a camera (ORCA-ER, Hamamatsu) with an imaging lens assembly (Optem Fusion, Qioptic). The imaging lens assembly consist of a 1.67× objective lens (35-00-04-000, Qioptic) and a 1.0× tube lens (35-08-06-000, Qioptic). The numerical aperture of the lens assembly can be adjusted using an adjustable aperture stop (35-07-25-000, Qioptic) placed between the objective and the tube lens. The sample is mounted on a stage assembly to perform x, y and z translation (8MT167M-25LS, Standa) as well as a θ rotation (8MR151, Standa). An x − y translation mount (SCP05, Thorlabs), and a tip-tilt mount (KC05-T/M, Thorlabs) enable fine tuning of the position of the rotation axis and the sample. The agarose sample is attached to the lower end of the tip-tilt mount using a custom metal cylinder and glue.

Fig. 4 Schematic representation of the experimental OPT setup. LS: Light source, NDF: neutral density filter, PH: pinhole, C: cuvette, TS: translation stage assembly, F: filter, L1 L2 and L3 Lenses, CA, camera, XY: manual x and y translation mount, TT: tip-tilt mount, θ: rotation stage, X: x translation stage, Y: y translation stage, Z: z translation stage.

Download Full Size | PDF

A total of 360 projections are recorded at 1 intervals. Measurements are performed with a camera integration time of 5.0 seconds and a laser output power of 28.13 mW. The measurement time of a complete sinogram is approximately 35 minutes.

3.2. Calibration and performance

The PSF of the imaging lens system (L3) is characterized by its beam waist as function of the distance from the focal plane w(s). The PSF is measured by imaging a straight knife edge inside the water-filled cuvette at multiple distances from the focal plane using the translation stage assembly. The measured PSF is fitted using a Gaussian function for all the measured stage translations to determine w(s) and w₀, which is subsequently used to predict the PSF in the reconstructed images. The measured beam waist in water at the focus is w₀ = (6.2 ± 0.1) μm.

The axis of rotation of the sample is aligned with the middle pixel column of the camera and positioned in the focal plane using the stage assembly and the tip-tilt adjuster. The center of the sample is positioned at the center of rotation using the x − y adjuster.

3.3. Sample preparation

Green fluorescent silica particles (KI-PSI-G4.0, Kisker) with a diameter of 4.0 μm are suspended in 1%, low melting point agarose gel (H26417.14, VWR). The fluorescent particles have peak excitation and emission wavelengths of 485 nm and 510 nm respectively. During measurements the agarose phantom is submersed in a cuvette filled with demineralized water to reduce refractive-index differences at the borders of the sample.

A ten day old transgenic zebrafish larvae expressing membrane bound green fluorescent protein (GFP) from a beta-actin promoter is mounted in low melting point agarose and imaged in the OPT system. Before mounting, the zebrafish larvae is euthanized in ice water in the Erasmus Medical Center, Rotterdam according to animal welfare regulations. Animal experiments are approved by the Animal Experimentation Committee of the Erasmus MC, Rotterdam. Imaging of the zebrafish is performed using an additional 7 : 1 manual zoom section (35-31-10-000, Qioptic) in the lens assembly L3. Projections are measured at 1 intervals over the full 360 range with an integration time of 1 s. Ten vertical slices are averaged before applying the deconvolution algorithm for 4 iterations.

3.4. Data analysis

After aquisition, sinograms are created by the summation of 100 pixel rows. This ensures that the slice is thick with respect to the spotsize on the detector such that Eq. (16) is valid. Subsequently, a constant background is removed from the sinograms. After correction for photobleaching with a characteristic e⁻¹ timescale of 798 seconds, the center pixel row in the sinograms is aligned with the center of rotation. Reconstruction of tomographic images is performed from all recorded projections using the MATLAB function iradon.

After reconstruction, the fluorescent beads are manually identified in the tomographic images. All visible beads are included in the analysis, except beads on the agarose-water interface as there might be clustering and refraction effects at the boundaries. A 0.31 × 0.31 mm region of interest (ROI) is selected for each bead and fitted with the elliptical Gaussian of Eq. (21) using the MATLAB function fit.

3.5. Image deconvolution

The position-dependent PSF of the OPT system is used to deblur the reconstructed images. After standard reconstruction with an inverse Radon transformation the image is transformed into cylindrical coordinates (r, θ), which ensures a space-invariant PSF in the r direction. The PSF in the θ direction only changes with r. After subsequent 1D deconvolutions in the θ and r directions the image is transformed back into Cartesian coordinates to obtain a deblurred reconstruction image. The deconvolution is performed using the Lucy-Richardson method for 100 iterations implemented using the MATLAB function deconvlucy, The coordinate transformations and image deconvolution combined take several minutes to process on a standard desktop computer.

4. Results

A single slice of an OPT reconstruction of the sample with fluorescent beads is shown in Fig. 5(a). Several individual beads are visible in the image, showing the characteristic blurring in the tangential direction, as predicted by Eq. (21). Also the contours of the agarose phantom are visible due to some fluorescent beads sticking to this boundary. An example of a single-bead ROI, Fig. 5(b), and the corresponding fitted Gaussian Fig. 5(c) are also shown.

Fig. 5 (a) Reconstruced slice of the phantom showing single fluorescent beads. A ROI (b) around the bead is selected and the fluorescence emission is fitted (c).

Download Full Size | PDF

The radial and tangential FWHM of all the selected beads in the reconstructed images are shown in Fig. 6(a) for both the radial and tangential directions. The measured resolutions are compared to the theoretical model of Eq. (21) based only on the measured evolution of the Gaussian beam waist with displacement. The measured FWHM resolution agrees well with the theoretically predicted values. The resolution in the radial direction remains constant throughout the image, while the tangential resolution becomes worse for larger distances from the center of rotation.

Fig. 6 (a) Radial (blue) and tangential (red) FWHM resolution as function of the radial distance from the center of rotation. Theoretical curves (solid lines) are compared to measured FWHM resolution in radial (open blue squares) and tangential (open red circles) directions. In most cases the markersize exeeds the errormargin. (b) Radial (open blue squares) and tangential (open red circles) FWHM resolution after deconvolution. Theoretical curves plotted for comparison.

Download Full Size | PDF

Reconstructed images can be enhanced by deconvolution with the spatially-varying PSF of the OPT system based on our theory. The overall effect of the image deconvolution is shown in Fig. 6(b), where the radial and tangential FWHM of all the analyzed beads after deconvolution are plotted as function of their radial distance from the center of rotation. The theoretically predicted FWHM resolutions without performing deconvolution are plotted for comparison. The graph shows that the deconvolution algorithm reduces both the radial and tangential resolution for most of the beads to below the radial resolution limit.

Figure 7 shows the deconvolution results for two selected beads in the OPT image. An image of the beads before deconvolution (Figs. 7(b) and 7(f)) is compared to an image after deconvolution (Figs. 7(c) and 7(g)). Fluorescence cross sections in radial and tangential direction (Figs. 7(d) and (h)) are significantly wider in the original image compared to the corresponding cross sections in the deconvoluted image (Figs. 7(e) and 7(i)).

Fig. 7 (a) Standard reconstructed slice of the fluorescence phantom. The image of a typical bead (b) and (f) show significant sharpening after the deconvolution (c) and (g). Radial (blue) and tangential (red) cross-sections of the initial and deconvoluted images of the bead are shown in (d) and (h), and (e) and (i) respectively.

Download Full Size | PDF

The deconvolution algorithm is applied to an OPT scan of a GFP labelled zebrafish larvae to illustrate the performance on a biological sample. Figure 8 shows the resulting reconstruction of a single slice before and after deconvolution. Figure 8(a) shows a single projection of the zebrafish. Comparison of Figs. 8(b) and 8(c) shows a clear improvement in image quality is obtained by application of the spatially-variant deconvolution algorithm.

Fig. 8 (a) Projection of a zebrafish with a reconstructed slice before (b) and after (c) deconvolution. The reconstruction corresponds to the vertical position indicated by the white line in (a). The inset compares the intensity profiles allong the indicated lines.

Download Full Size | PDF

5. Discussion

In this paper, we have presented a theoretical model that describes the spatially-variant resolution of optical tomography systems. The theoretical model relates the size of the PSF in the reconstructed tomographic image to the PSF of the imaging optics. In the derivation of the model we have assumed that the reconstructed slice is thick with respect to the spotsize of a point emitter on the detector. This approximation reduces the model to 2D and neglects any 3D imaging effects. If the slice is thin then the factor A(θ) in Eq. (15) should be taken into account. An analytical expression of this effect proved to be difficult to calculate analytically. Consequently, taking 3D effects into account leads to complex PSF shapes in the reconstruction domain that are challenging to deconvolute.

The theoretical model described in this paper assumes that the PSF of the imaging optics can be described by a Gaussian beam. Although the PSF of a typical lens system is not necessarily a Gaussian beam, for the lens system used in our experiments this was a good approximation. A fit of the measured PSF with a Gaussian beam-shape had a value of R-squared of R² = 0.9962.

The measured resolution in both radial and tangential directions agrees well with the theoretical model. We attribute the spread in the measured dimensions of the spots to small imperfections in the alignment of the sample and the axis of rotation. The summation of 100 rows in the z-direction does not guaranty the absence of 3D effects as some beads will be close to edges of a section. This influences the resolution as well as reducing the validity of the fit-function, causing variations from the predicted resolution.

Knowledge on the spatially-variant resolution in an optical tomography setup allows for specific guidelines for the design of such systems if certain resolution targets are to be achieved. Our model describes how the trade-off between resolution and depth of focus of the imaging optics can be optimally balanced for a certain sample size. A possible design criterion might be to minimize the worst resolution in the entire reconstructed image. Since the radial resolution is constant, this amounts to limiting the tangential resolution at the largest distance from the center of rotation. Figure 9 shows the FWHM tangential resolution as a function of the in-focus Gaussian beam waist w₀ for several distances from the center of rotation. At a particular radial distance r a small beam waist results in a large Gaussian beam divergence, giving rise to a poor tangential resolution. As the beam waist increases the beam divergence decreases and the resolution improves. At large w₀ the resolution becomes poor due to the increased beam waist w₀ in-focus. From this result it can be seen that one can sacrifice resolution in the center to improve resolution at the edges of the image. There is clearly an optimal beam waist w_0,opt for each distance from the center, given by

w_{0, opt} = \sqrt{\frac{λ r}{π}} .

The optimal beam waist w_0,opt represents the in-focus beam waist for which the tangential resolution at the edges of the object is minimal.

Fig. 9 Tangential FWHM resolution as function of the beam waist, plotted for several distances from the rotation center.

Download Full Size | PDF

The theoretically derived PSF in the reconstructed images allows for image de-blurring. The Lucy-Richardson deconvolution algorithm applied in this work provides a possible implementation of such a method. It shows that significant improvement can be obtained in image resolution, although the resolution improvement comes at the cost of some increased noise and enhancement of artefacts in the image. We observed that the number of iterations in the deconvolution has a large effect on this trade-off. The use of 100 iterations in this paper was motivated by the increase in resolution mainly. Using fewer iterations degrades the final resolution, but also shows a smaller increase in noise and artefacts. The deconvolution results shown here are a proof of concept, further optimization of our implementation is the topic of current research.

While our model is derived for a Gaussian-beam shaped PSF of the imaging optics with the rotation axes in the focal plane, it can be extended to other beam shapes and focus arrangements with relative ease. Some authors position the focal plane between the center of rotation and the edge of the sample [3]. In addition, other imaging beams can be implemented using the proposed theoretical framework. The depth-of-field of the imaging optics can be extended by implementing a Bessel-beam shaped focus. The resolution in the tomographic reconstructed image obtained by these methods can be modelled and compared in a similar manner as described in this paper.

In a broader scope, the presented model can be useful in other fields that apply tomographic imaging techniques such as terahertz tomography, tomographic phase microscopy, electron microscopy, and PET/SPECT imaging.

6. Conclusion

We have presented and validated an analysis of the image resolution in optical tomography systems. The model describes the spatially-variant PSF in the reconstructed images as a function of the properties of the imaging optics. The presented model provides users with a description of the system resolution, provides guidelines for system design, and can be used for image restoration using deconvolution algorithms. The model can be easily adapted for various tomography applications.

Funding

This research is supported by the Dutch Technology Foundation STW, which is part of the Netherlands Organisation for Scientific Research (NWO), and which is partly funded by the Ministry of Economic Affairs.

Acknowledgments

The authors would like to thank L. J. van Vliet, L. N. van Opstal, and A. K. Trull for useful discussions, and T. J. van Ham for his assistance with the zebrafish experiments.

References and links

1. P. J. Keller, A. D. Schmidt, J. Wittbrodt, and E. H. Stelzer, “Reconstruction of zebrafish early embryonic development by scanned light sheet microscopy,” Science 322, 1065–1069 (2008). [CrossRef] [PubMed]

2. 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, 541–545 (2002). [CrossRef] [PubMed]

3. J. Sharpe, “Optical projection tomography,” Annu. Rev. Biomed. Eng. 6, 209–228 (2004). [CrossRef] [PubMed]

4. L. Chen, N. Andrews, S. Kumar, P. Frankel, J. McGinty, and P. M. French, “Simultaneous angular multiplexing optical projection tomography at shifted focal planes,” Opt. Lett. 38, 851–853 (2013). [CrossRef] [PubMed]

5. 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. Lett. 35, 3363–3365 (2010). [CrossRef] [PubMed]

6. W. Xia, R. M. Lewitt, and P. R. Edholm, “Fourier correction for spatially variant collimator blurring in SPECT,” IEEE Trans. Med. Imaging 14, 100–115 (1995). [CrossRef] [PubMed]

7. J. R. Walls, J. G. Sled, J. Sharpe, and R. M. Henkelman, “Resolution improvement in emission optical projection tomography,” Phys. Med. Biol. 52, 2775 (2007). [CrossRef] [PubMed]

8. A. Darrell, H. Meyer, K. Marias, M. Brady, and J. Ripoll, “Weighted filtered backprojection for quantitative fluorescence optical projection tomography,” Phys. Med. Biol. 53, 3863 (2008). [CrossRef] [PubMed]

9. L. Chen, J. McGinty, H. B. Taylor, L. Bugeon, J. R. Lamb, M. J. Dallman, and P. M. French, “Incorporation of an experimentally determined MTF for spatial frequency filtering and deconvolution during optical projection tomography reconstruction,” Opt. Express 20, 7323–7337 (2012). [CrossRef] [PubMed]

10. C. M. McErlean, E. Bräuer-Krisch, J. Adamovics, and S. J. Doran, “Assessment of optical CT as a future QA tool for synchrotron X-ray microbeam therapy,” Phys. Med. Biol. 61, 320 (2015). [CrossRef] [PubMed]

11. M. Gu, Advanced optical imaging theory, vol. 75 (Springer Science and Business Media, 2000). [CrossRef]

12. J. Hsieh, Computed Tomography: Principles, Design, Artifacts, and Recent Advances (SPIE, 2009), Chap. 3.

13. H. Kogelnik and L. Tingye, “Laser beams and resonators,” Appl. Opt. 5, 1550–1567 (1966). [CrossRef] [PubMed]

Image resolution and deconvolution in optical tomography

Abstract

1. Introduction

2. Theory

2.1. Image formation of a single projection

2.2. Tomographic point spread function

2.2.1. Gaussian beam PSF

3. Methods

3.1. Experimental setup

3.2. Calibration and performance

3.3. Sample preparation

3.4. Data analysis

3.5. Image deconvolution

4. Results

5. Discussion

6. Conclusion

Funding

Acknowledgments

References and links

Cited By

Figures (9)

Equations (22)

Optics Express