Image resolution and deconvolution in optical tomography.

We present a frequency domain analysis of the image resolution of optical tomography systems. The result of our analysis is a description of the spatially-variant resolution in optical tomographic image after reconstruction as a function of the properties of the imaging system geometry. We validate our model using optical projection tomography (OPT) measurements of fluorescent beads embedded in agarose gel. Our model correctly describes both the radial and tangential resolution of the measured images. In addition, we present a correction of the tomographic images for the spatially-varying resolution using a deconvolution algorithm. The resulting corrected tomographic reconstruction shows a homogeneous and isotropic pixel-limited resolution across the entire image. Our method is applied to OPT measurements of a zebrafish, showing improved resolution. Aside from allowing image correction and providing a resolution measure for OPT systems, our model provides a powerful tool for the design of optical tomographic systems.

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.

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.

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)| 2 , 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 If a detector is placed in the focal plane s = 0, the intensity distribution on the detector is This equation shows the main difference between OPT and straight ray based tomography. If |h(t, s, z)| 2 = δ(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)| 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 The resolution of the tomographic imaging system is calculated by determining the response to a point object function, i.e., f (x, 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 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.
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.
x y s To obtain the frequency description of the projection we take the Fourier transform of Eq. (5) along the transverse coordinate t, to obtain Assuming a parallel-beam geometry the projection slice theorem states that the 2D Fourier transform of the image of the object, denoted with OT F ( 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 OT F ( f x , f y ) at angle θ, as shown in Fig. 2 and Subsequent 2D inverse Fourier transform of OT F ( f x , f y ) yields the position dependent PSF, as illustrated in Fig. 2(c).

Gaussian beam PSF
The PSF of the tomographic imaging system can be calculated from the PSF of the imaging optics |h(t, s, z)| 2 . In general, |h(t, s, z)| 2 can have a complicated shape. Here, we assume the imaging PSF is a focused Gaussian beam described by where w(s) is the beam waist, given by with λ the emission wavelength of the fluorophore and w 0 the Gaussian beam waist defined as the e −1 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 with Fourier transformation in the transversal t direction gives 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 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 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 where θ u = θ o − θ now denotes the angle with the u axis. Further simplification using sin 2 Equation (19) is separable into two factors with dependence on u and v only, and is further simplified to with

Taking the inverse Fourier transform of
Eq. (20) we obtain the PSF of the tomographic imaging system which is the main theoretical result of this section. Figure 3  w 0 = 10 µm and λ = 515 nm. The PSF in the center of the image is isotropic and has a width w 0 in all directions. The effect of the lens system PSF away from the center results in a radial resolution of w 0 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.

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.

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 0 , which is subsequently used to predict the PSF in the reconstructed images. The measured beam waist in water at the focus is w 0 = (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.

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 refractiveindex 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.

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 −1 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.

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.

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.
The radial and tangential FWHM of all the selected beads in the reconstructed images are shown in Fig. 6(a)   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.
Radial distance from center of rotation (mm)  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.  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

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 2 = 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 0 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 0 the resolution becomes poor due to the increased beam waist w 0 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,o pt for each The optimal beam waist w 0,o pt represents the in-focus beam waist for which the tangential resolution at the edges of the object is minimal. 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.

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.