Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Mechanical artifacts in optical projection tomography: classification and automatic calibration

Open Access Open Access

Abstract

Optical projection tomography (OPT) is a powerful tool for biomedical studies. It achieves 3D visualization of mesoscopic biological samples with high spatial resolution using conventional tomographic-reconstruction algorithms. However, various artifacts degrade the quality of the reconstructed images due to experimental imperfections in the OPT instruments. While many efforts have been made to characterize and correct for these artifacts, they focus on one specific type of artifacts, whereas a comprehensive catalog of all sorts of mechanical artifacts does not currently exist. In this work, we systematically document many mechanical artifacts. We rely on a 3D description of the imaging system that uses a set of angular and translational parameters. We provide a catalog of artifacts. It lists their cause, resulting effects, and existing correction methods. Then, we introduce an automatic calibration algorithm that is able to recover the unknown system parameters fed into the final 3D iterative reconstruction algorithm for a distortion-free volumetric image. Simulations with beads data and experimental results on a fluorescent textile fiber confirm that our algorithm successfully removes miscalibration artifacts in the reconstruction.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Conventional optical microscopy such as confocal microscopy is limited to the imaging of relatively thin samples. This limitation can be partially overcome with optical projection tomography (OPT) which was invented by Sharpe in 2002 [1]. Over the years, OPT has become a mature tool for the production of high-resolution 3D images of biological samples at mesoscopic scale [2,3] in a brightfield [4,5] or fluorescence [69] configuration. OPT is widely used in a variety of applications, such as the mapping of the distribution of proteins in embryos [1,10], the localization of metastases in lymph nodes [9], the display of vascular networks and amyloid depositions in the mouse brain model to study Alzheimer’s disease [2], and the imaging of the spatial arrangement of intestinal villi [3].

Since OPT is a tomographic-imaging technology, it falls into the same category as X-ray computed tomography (CT), single-photon electron computed tomography and electron tomography [1]. The contrast mechanism that these technologies rely upon is either the attenuation or the emission function of rays of light. For OPT, the rays follow optical straight lines that are geometrically only approximately straight. They project the 3D inner structure of the sample onto a 2D detector plane. The mathematical tool to describe such a straight-ray projection is the Radon transform [11]. The associated inverse problem is to reconstruct the 3D volume from the set of 2D projections acquired at various spatial positions of the sample. The discrete inversion formula of the Radon transform, the filtered backprojection (FBP) algorithm, is efficiently implemented and widely used in practice [3,9,12].

While OPT achieves high-resolution 3D imaging at a relatively low cost [13,14], the reconstructions often suffer from artifacts that impact the quality of the images due to several types of model mismatch. Among them, mechanical errors in the imaging system play a non-negligible role, especially in low-cost OPT systems [13]. The most common type of mechanical error is an offset of the rotation axis which results in a misaligned center of rotation during the experiment. This typically leads to distortions including double edges [7] or circles [15] in the reconstruction, depending on the type of sample. Many 2D methods have been proposed to address this issue, including maximum variance of a reconstructed slice under a set of guesses of the true rotation [7,9], sinogram unification of both fluorescent and bright-field OPT [15], center of mass or image registration [16], and total-variation regularization [17]. These 2D methods assume that the tilt of the rotation axis is negligible. Their extension to small tilt angles introduces a height dependence of the rotation axis [9,13]. For larger tilt angles, a method to account for tilt in 3D is still missing. Another type of mechanical error that has not been well studied is the angular errors of the rotation motor, which will result in seagull-shaped artifacts for point-like objects, as we show in Section 2. Artifacts may also arise due to optical effects. For example, the assumption that optical straight paths coincide with geometric ones need to be abandoned, thus requiring variations of the conventional Radon transform [6,1820]. Other types of model mismatch such as mismatches in the refractive index [2123], illumination fluctuations [7], spatial or temporal variations of the sensitivity of the detectors, their linearity, and background noise [7,21] can lead to various distortions as well, degrading the reconstruction quality.

In this paper, we first strive to provide a comprehensive catalog of many types of mechanical artifacts as a reference for OPT practitioners to assist them to calibrate their experiments. We rely on point-like objects (a very popular tool in the characterization of OPT setups) to demonstrate the cause and resulting appearance of each artifact. Then, we introduce a 3D auto-calibration algorithm to remove the artifacts due to mechanical errors listed in our catalog, for the convenience of OPT experimentalists.

This paper is organized as follows. In Section 2, we first describe the imaging geometry of OPT using a set of system parameters (angles and shifts). We point out typical places where the most critical mechanical errors occur. Then, we present a catalog of artifacts along with a detailed description of the three most common types of mechanical errors. In Section 3, we introduce a 3D forward model that characterizes all sorts of mechanical errors and a corresponding joint-reconstruction-calibration algorithm at coarse scale. It is able to recover the unknown set of characterizing angles including the rotation and tilt angles. The refined system parameters are then fed into the 3D reconstruction algorithm to achieve an image without miscalibration artifacts at a finer scale. Results on both simulated measurements of beads and experimental data of a fluorescent textile fiber are presented to validate our algorithms.

2. Characterization of mechanical artifacts

2.1 Common OPT geometry and reconstruction algorithms

The general OPT geometry is illustrated in Fig. 1. A 3D sample is placed in a rotating cylinder and imaged on a camera. This camera records 2D projections of either the absorption of the sample (transmission OPT) [18] or its emitted fluorescence (emission OPT) [6,24]. These tomographic projections, not considering more intricate models of light propagation, are enabled in reason of the very low numerical aperture of the OPT imaging system, in sharp contrast with conventional optical microscopy. The sample is rotated around the axis of the cylinder for tomographic reconstruction.

 figure: Fig. 1.

Fig. 1. Diagram of OPT. The light that passes through or is emitted by the sample gets imaged on a CCD camera using a lens system. In this ideal configuration, a slice of the sample perpendicular to the rotation axis corresponds to the information in the temporal sequence recorded by a row of pixels on the CCD camera.

Download Full Size | PDF

We place the sample in a 3D coordinate system ($O-xyz$) (Fig. 2). In an ideal setup [10], the center of the cylinder is at $O$ while the rotation axis would be perfectly parallel to the $z$-axis. In our refined model, the center of the sample cylinder is at $G$ and its orientation is described by the three angles $\theta =(\boldsymbol {\varphi },\psi _1,\psi _2)$ for a rigid-body rotation. The $P$-dimensional vector $\boldsymbol {\varphi } = (\varphi _1, \ldots, \varphi _P)$ represents the set of $P$ rotation angles around the rotation axis, one for each captured projection. The tilted rotation axis is described by the remaining two angles $\psi _1$ and $\psi _2$ which are scalars that represent the out-of-plane and in-plane tilt, respectively, common to all views. The detector plane is described by a 2D coordinate system ($O'-\xi \eta$) that is perpendicular to the optical axis.

 figure: Fig. 2.

Fig. 2. OPT imaging system with mechanical errors. (The lens system is omitted here for simplicity.) The cylinder represents the tilted sample that rotates around its orientation axis. The center of the cylinder is translated to $G$. The detector plane is described by a 2D coordinate system ($O'-\xi \eta$) perpendicular to the optical axis ($y$ axis). The projection of $O$ and $G$ in the detector plane are $O'$ and $G'$, respectively. The angle $\psi _1$ represents the out-of-plane tilt angle between the orientation axis of the sample and the $z$-axis direction (gray vertical arrow). The rotation angle $\psi _2$ represents the in-plane tilt angle of the detector plane around the optical axis.

Download Full Size | PDF

Many reconstruction algorithms are available for OPT, from FBP [19,25] to optimization-based methods [20,26]. They mainly come from the field of CT, a canonical tomography application using X-rays with a similar parallel-beam geometry [27]. In the simplistic model where the rotation axis is assumed to be aligned with the $z$-axis and when the in-plane tilt vanishes, one horizontal plane of the sample corresponds to one line of the camera. Many reconstruction algorithms thus operate on camera images line by line, reconstructing the object plane by plane.

A mismatch between the geometry of the setup and that assumed by the reconstruction algorithm obviously impact the quality of the final reconstruction. This induces artifacts that have been reported previously in the literature, along with dedicated procedures to correct for them.

2.2 Previously reported artifacts and correction techniques

It may occur that the center of rotation (COR) $G$ is off-center, which is a common issue for OPT experiments. The consequences depend on the direction of this translation mismatch. A constant shift in the $x$-axis direction (parallel to the detector plane) has been reported in both OPT [7,15] and X-ray CT [16]. It will result in a constant shift in the 2D projections generating ringing artifacts: each bright point in the sample will be imaged into a circle after OPT reconstruction. For example, Tang et al. reported circles in the imaging of a zebrafish embryo [15]. Due to the finite depth-of-focus in OPT, these circles may take different flavors. Walls et al. observed a so-called double-edge artifact while imaging the cardiac region of a mouse embryo [7] and Donath et al. observed tuning-fork artifacts of a point object [16]. A shift in the $y$ axis (the direction of light propagation) will result in out-of-focus blur artifacts due to limited depth-of-field that is unique to OPT [18].

Different methods are available to correct for an off-center COR in both OPT and the field of X-ray CT. Walls et al. proposed a heuristic method that calculates the variance of each reconstruction of a 2D slice using a range of assumed shift values in a sinogram: the shift value that produces the maximum variance serves as a close guess for the true shift [7]. This method is fast to implement and yields satisfactory reconstructions when the volume is not too sparse; it is commonly used as a preprocessing step by OPT practitioners [3]. The center-of-mass method is another popular method to determine the center position in OPT. It is based on the property that the center of mass of the object is projected onto the center of mass of the sinogram. Hence, a shift of the object will result in a corresponding shift of the center of mass in the sinogram. This constant shift is found by solving a linear system [16]. Other methods such as image registration [16] and methods that take advantage of the symmetrical structure of certain samples [14] are also used in OPT but with limitations because they require a priori information.

The tilt of the rotation axis in 3D can be described by a combination of the out-of-plane and in-plane tilt. Only in-plane tilts of a small angle have been studied. With this assumption, each transverse slice of the sample still approximately corresponds to a row of pixels on the camera; then, the compensation of the 3D tilt reduces to the 2D problem of finding the true center of rotation for each slice. This is done using the same techniques as for dealing with an off-center COR. For example, Torres et al. [9] used the maximum-variance method [7] to find the COR at two different heights. They then performed a linear regression to obtain the COR shift at any intermediate horizontal plane.

Previous techniques are dedicated to the correction of a single type of mechanical artifact. To the best of our knowledge, there has been no description of the effect of other tilt angles. Moreover, the tilt-correction strategy is limited to small tilt angles, and no robust correction method has been proposed yet.

2.3 Catalog of mechanical artifacts

The geometry of OPT is described by the three-dimensional coordinates $G (x_G, y_G, z_G)$ of the center of the cylinder center and the angles $\theta =(\boldsymbol {\varphi },\psi _1,\psi _2)$. In this parameterization, all parameters are fixed during an OPT acquisition except the rotation angles $\boldsymbol {\varphi }$. To characterize different mechanical imperfections, we introduce a catalog of the imaging artifacts that result from a mismatch between the “true” value and the value assumed at reconstruction for each of the parameters of out refined model. This catalog can be useful for experimentalists to identify artifacts they may encounter in the reconstructions, so as to correct for them either experimentally (re-calibration) or computationally. We summarize in Table 1 the results with references to previous works when relevant, accompanied with visualizations shown in Fig. 3 based on numerical simulations using point-like beads, a popular tool to characterize OPT systems.

 figure: Fig. 3.

Fig. 3. Artifacts in the reconstruction of simulated beads due to a controlled mismatch the geometric parameters of our refined model and the simplistic geometry assumed by traditional reconstruction methods. (a)-(b) Circle artifacts due to a constant shift of −4 pixels (a) and 8 pixels (b) of the COR in the $x$ direction. (c)-(d) Seagull artifacts due to a negative (c) and positive (d) random accumulative error with maximum amplitude 0.03 degrees in the rotation angles. (e)-(p) Various artifacts due to the tilt of the rotation axis. (e)-(h) Reconstruction at slice 192 (64 slices below the central slice). (i)-(l) Reconstruction at slice 256 (central slice in $z$-axis). (m)-(p) Reconstruction at slice 384 (128 slices above the central slice). The intensities of images are normalized and saturated for better visualization.

Download Full Size | PDF

Tables Icon

Table 1. Table of setup-related imaging artifacts and their dependence on the field of view (FOV) and on height (location along the $z$ axis).

The numerical simulation is done using Tomosipo, a convenient and versatile tomographic toolbox for 3D simulations and reconstructions of any geometric setup [30]. We use a ground-truth object of 10 scattered beads of size $6^3$ pixels in the FOV with their centers on the same plane perpendicular to the $z$ axis. This object is placed in a cube of size $512^3$ pixels. The number of rotation angles is 1200 over the range of 0 to 360 degrees to generate a set of 2D projections. The location of this plane was moved along the $z$ axis to observe the dependence of a certain type of mechanical artifacts on the $z$ location. We added different types of errors in the geometry of the forward model to simulate a realistic imperfect OPT system. Then, we used the standard FBP algorithm which corresponds to a generic imaging geometry with no correction of mechanical errors.

Shifts along the $x$-axis (orthogonal to the axis of rotation and light propagation) result in circle artifacts. The size of artifacts is uniform across the FOV; however, their severity, as is indicated by the radius of the circle, depend on the absolute value of the shift. The larger it is, the bigger the circle artifacts will be (Fig. 3(a) and 3(b)). This characterization can be applied to any object, by considering it as a sum of point-like elements. For continuous samples, their borders will appear as double edges as reported in previous OPT experiments [7]. Shifts along the $y$-axis (direction of light propagation) cause an out-of-focus blur when parts of the sample move out of the depth-of-field of the imaging system. Shifts along the $z$ axis (rotation axis) cause an overall shift of the reconstruction in the $z$ direction, but do not generate artifacts as long as the object is still within the field of view of the camera. Random vibrations may also blur the final reconstruction. They can be included in our framework as a projection-dependent translation error. This is different from the offset of the rotation axis that corresponds to a systematic error shared across all projections.

The tilt angles $\psi _1$ and $\psi _2$ generate artifacts of different shapes, especially when the mismatch is not small. To study the impact of the tilt of the rotation axis in different directions and levels of severity, as well as its dependence on the vertical location along $z$ axis, we moved the $z$ location of the beads to three different transverse slices: slice 256 (central slice), slice 192 (64 slices below the central slice), and slice 384 (128 slices above the central slice). We added errors in the two tilt directions described previously with two values: 5 and 10 degrees. The result is shown in Fig. 3(e) – 3(p).

The first two columns of images of the second to the fourth rows in Fig. 3 show that an out-of-plane tilt results in various artifacts that appear more severe close to the boundary of the image. The next two columns of images of the second to the fourth rows in Fig. 3 show that in-plane tilt in general leads to more severe artifacts of various shapes compared to the other direction. The dependence of the size of the artifacts on the location along the vertical direction is stronger than that of the tilt in the other direction. The central slice shows the least artifacts in the third row of Fig. 3, while slices far away from it show bigger artifacts in the second and last row of Fig. 3. As the tilt angle increases, the artifacts get stronger if we compare the first two columns or the last two columns of the tilt images in Fig. 3. We notice the circle artifacts in Fig. 3(g), 3(h) and 3(o) which indicate that, when the impact of the tilt is small, it can be approximated by a shifted COR problem. In the case of a relatively large unknown tilt, such an approximation does not hold anymore, as seen in Fig. 3(p). To the best of our knowledge, correction for the tilt in 3D has not been studied yet.

The last parameter describing our refined OPT geometry is the set of rotation angles $\boldsymbol \varphi$. On one hand, since each rotation angle is different for each camera image, a random mismatch on this angle does not induce artifacts but introduces noise that degrades the image quality. On the other hand, a miscalibration of the rotation motor may introduce a positive or negative drift in the rotation angle, which accumulates over time. We simulated this kind of model mismatch and observed seagull-shaped artifacts in the reconstruction (Fig. 3(c) and 3(d)). Depending on the sign of the random accumulative error, the seagull is either “flying” toward or outward the center of the image. The size of the seagull in each transverse slice depends on the distance of the bead relative to the boundary of the object. The closer to the boundary, the bigger the seagulls will be while the bead located at the exact center of the slice appears unaffected since it coincides with the rotation axis (Fig. 3(c) and 3(d)).

In Section 3, we introduce an automatic calibration algorithm to correct for the mechanical errors presented in Table 1. We start with a 3D forward model that fully characterizes our refined OPT geometry, including the 3D tilt that is out of reach of existing 2D methods. Then, we formulate the calibration of the system parameters as a multiscale joint reconstruction-calibration optimization problem. This multiscale scheme allows us to overcome the memory bottleneck of 3D models and helps us to accelerate the calibration. Once the system is calibrated, we are able to reconstruct an artifact-free volume.

3. Automatic calibration of mechanical artifacts

In this section, we propose a computational framework that automatically calibrates the parameters of our refined model of an OPT imaging system. This framework is able to detect the model mismatch between the simplistic forward model and the measurement data. It outputs a set of calibrated system parameters and allows us to improve the reconstruction of the 3D volume. We first show how to characterize the three types of mechanical errors mentioned in Section 2.3. Then, we present our multiscale joint reconstruction-calibration algorithm and show how to remove the artifacts in both simulated and experimental data.

3.1 Forward model and characterization of mechanical errors

The coordinate system to describe the OPT geometry is the same as described in Fig. 2. When there exist mechanical errors in the system, we describe the angular errors as the perturbation vector $\boldsymbol {\delta }=(\boldsymbol {\delta }_{\boldsymbol {\varphi }},\delta _{\psi _1},\delta _{\psi _2})$. The actual angles $\theta ={\theta }^*+\boldsymbol {\delta }$ can be expressed as a sum of the error vector $\boldsymbol {\delta }$ and the ideal angle vector ${\theta }^* = (\boldsymbol {\varphi }^*, \psi _1^*, \psi _2^*)$ that characterizes a simplistic OPT system: equidistant rotation angles $\boldsymbol {\varphi }^*$ between 0 and 360$^{\circ }$ and $\psi _1^* = \psi _2^* = 0$. In addition, translation error of the sample is described by a 2D vector $\boldsymbol {t}=(t_1,t_2)\in \mathbb {R}^2$ in the detector plane.

We omit optical effects and complex PSF models to keep the forward model computationally efficient, as is done in most OPT experiments. This omission is often acepted in the focal-sheet-scanning OPT setup [31]. To fully characterize the mechanical errors, we adopt the 3D X-ray transform $\mathcal {P}_{\theta }$ that provides a mathematical description of the straight-ray projections of a sample at any 3D pose [32] as

$$b^{\theta,\,\boldsymbol{t}}(\boldsymbol{v})=\mathcal{P}_{\theta}\{f\}(\boldsymbol{v}-\boldsymbol{t}) + n,$$
where the compactly supported function $f(\boldsymbol {u})\in L_2(\mathbb {R}^3)$, $\boldsymbol {u}=(x,y,z)$ represents the 3D sample to be reconstructed. The measurement in the detector plane is $b^{\theta,\,\boldsymbol {t}}(\boldsymbol {v})$ for a location $\boldsymbol {v}=(\xi,\eta )$, a given sample orientation $\theta$, and a shift vector $\boldsymbol {t}$. The additive random noise is $n$. To numerically implement the forward model, we discretize Eq. (1) under a sampling framework described in [33] and obtain the following linear system. The details of the derivation of Eq. (2) and the construction of the matrix $\textbf {H}(\boldsymbol {\theta }, \boldsymbol {t})$ are provided in Supplement 1.
$$\textbf{b}=\textbf{H}(\boldsymbol{\theta}, \boldsymbol{t})\textbf{c}+\textbf{n},$$
where $ {\textbf{b}, \;\textbf{n}}\in {\mathbb R}^{\mathbb{MP}}$ are the measurement and noise vectors of $P$ projections and where the system matrix $\textbf {H}(\boldsymbol {\theta }, \boldsymbol{t})\in \mathbb{R}^{MP\times N}$ is a function of system parameters $\boldsymbol {\theta }$ and $\boldsymbol {t}$. The 3D volume is represented by a finite-dimensional coefficient vector $\textbf{c}\in \mathbb{R}^{N}$ using the optimized Kaiser-Bessel window functions that are well suited for tomographic settings [34]. This choice of basis is convenient to compute analytic gradients for the tilt angles and shifts [32].

Tables Icon

Algorithm 1. Automatic Calibration

3.2 Multiscale calibration-reconstruction algorithm

The measurement data generated by a real OPT experiment usually has a very large size. For instance, current cameras can acquire OPT images up to size $2048^2$ pixels with pixel size less than 1$\mu$m per projection [3]. Moreover, the rotation motor can achieve a step angle as precise as 0.3 degrees, resulting in 1200 projections. This means the storage of the vector $\textbf {b}$ can take up to several tens of gigabytes. 3D reconstruction and calibration on such large datasets are either infeasible memory-wise on GPU or infeasible speed-wise on CPU. We thus propose the multiscale scheme illustrated in Fig. 4 to overcome the computational bottleneck. We first downsample the original measurement data to a computationally feasible scale and then run our automatic calibration algorithm at coarse scale. The inverse problem here is formulated as:

$$\textbf{c}^*,(\boldsymbol{\theta}^*,\boldsymbol{t}^*)\in\arg\underset{\textbf{c},\boldsymbol{\theta},\boldsymbol{t}}{\min}\left\{\frac{1}{2}\|\textbf{H}(\boldsymbol{\theta},\boldsymbol{t})\textbf{c}-\textbf{b}\|^2_2\right\}.$$

The reconstruction pipeline is detailed in Algorithm 1. It consists of the recovery of the coefficient vector $\textbf {c}$ by solving the optimization problem in Step 3 and calibrating the system parameters $\boldsymbol {\theta }$ and $\boldsymbol {t}$ by solving another optimization problem in Step 4 in alternating fashion. We start with an initial guess for the system parameter $\boldsymbol {\theta }^0$ and $\boldsymbol {t}^0=(0, 0)$ that corresponds to a perfect setup. Inspired by the observations of Section 2.3, which indicate that the radius of the circle artifacts is directly related to value of the shift, we shift the projections over the range of values estimated from the severity of the artifacts, apply the FBP algorithm to observe the evolution of the circle artifacts in the reconstruction, then choose the value such that the circle artifacts are sufficiently suppressed as the initial guess for the value of the shift. This allows us to attain a relatively close initial guess for $\boldsymbol {t}$ in a fast manner, which is crucial for the convergence of the calibration algorithm. Such a set of $\boldsymbol {\theta }$ and $\boldsymbol {t}$ serves as a good initial guess for the calibration algorithm, which is important to avoid local minima since 3 is a nonlinear optimization problem.

 figure: Fig. 4.

Fig. 4. Workflow of the multiscale calibration-reconstruction algorithm. We first perform a joint reconstruction-calibration at coarse scale. We then use the calibrated mechanical parameters for the final high-resolution reconstruction.

Download Full Size | PDF

The alternating process is repeated until the system parameters are well refined. The calibrated system parameters are then used to reconstruct an artifact-free 3D image by running an extra Step 3 in Algorithm 1 at a finer scale, as described in Fig. 4.

3.3 Results on simulated and experimental data

To validate Algorithm 1, we simulate the ground truth as a 3D cube of size $512^3$ pixels in which we have randomly inserted 150 fluorescent beads. A positive constant accumulative error of 0.05 degrees in the rotation angles and a constant shift of 8 pixels along the $x$ axis are added to the forward model to simulate a set of 300 projections of realistic OPT measurements. Then, we downsampled the measurements to a coarse scale of ($128\times 128\times 300$) pixels using the Matlab function imresize with bicubic interpolation and applied Algorithm 1. The initial guess was set as described in Section 3.2. The reconstruction step of the algorithm is implemented using the GlobalBioIm library, an open-source library that provides easy-to-use computational modules for solving inverse problems in computational bio-imaging [35]. The calibration step uses functions from the Cryo-Em-Refinement library [32].

In total, 10 global rounds of joint reconstruction-calibration were used, each of which composed of only 30 iterations of reconstruction and 6 iterations of calibration to avoid overfitting in the presence of model mismatch. After having obtained the calibrated system parameters, we reconstructed the 3D volume at a finer scale ($512^3$) using the FBP algorithm. In Fig. 5(a), we show that the reconstruction without any calibration or shift correction contains both the circle and seagull-shaped artifacts. The reconstruction with only naive shift correction (Fig. 5(b)) still suffers from seagull artifacts due to residual model mismatch. After applying our calibration algorithm, both the circle and the seagull artifacts are successfully removed all at once (Fig. 5(c)) and the reconstruction is very close to the ground-truth image in Fig. 5(d).

 figure: Fig. 5.

Fig. 5. (a) Reconstruction with circle and seagull artifacts at one slice (147) without calibration. (b) Reconstruction result at slice 147 after a naive shift correction. (c) Calibrated reconstruction at slice 147. Both the circle and seagull artifacts are successfully removed. (d) Ground truth at slice 147. The intensities of images are normalized and saturated for better visualization.

Download Full Size | PDF

Despite the best effort of the experimentalists to calibrate the hardware system, the rotation axis may still contain undesirable tilting which degrades the quality of the reconstructed OPT images. Similar to the procedures in the simulation, we follow the steps in the flowchart in Fig. 4 to further validate our algorithm on an experimental dataset of a fluorescent textile fiber that contains errors in the the tilt angles. The data are acquired using a focal-sheet-scanning OPT system [31]. It uses a lateral light-sheet illumination to reduce the out-of-focus blur in the images, enabling the assumption of straight-ray projections for our forward model. We use $720$ projections of size ($256 \times 256$) that we downsample to a coarse scale of ($64 \times 64 \times 720$) pixels. We used 6 global rounds of joint reconstruction-calibration, each composed of 30 iterations of reconstruction and 5 iterations of calibration.

The 3D visualization of the reconstruction result of the fiber with an out-of plane tilt angle of approximately 4 degrees is displayed in Fig. 6. Without any correction, the reconstruction shows multiple ghost artifacts due to a combination of misaligned COR and tilt of the rotation axis in Fig. 6(a). The severity of these ghost shadows are reduced after naive shift correction, as seen in Fig. 6(b) but the top and middle parts of the reconstructed fiber still suffer from ghost shadows. The 3D reconstruction using the calibrated system parameters output by Algorithm 1 effectively removes all the shadows and achieves an artifact-free image. The evolution of the cost function during the joint reconstruction-calibration is displayed in Fig. 7(a) and shows a 73% reduction of the cost compared to the uncalibrated configuration. After 6 global rounds, the calibration algorithm found an out-of-plane tilt angle of 3.7 degrees which is very close to the controlled 4 degrees. This further confirms that the calibration manages to reduce the model mismatch due to tilt. Additionally, we repeat this process for four different tilt angles and display the results of tilt calibration in Fig. 7(b). All four scenarios of different magnitudes and directions of tilt angles led to success as the calibrated angle is very close to the controlled true value.

 figure: Fig. 6.

Fig. 6. Textile fiber data. (a)–(b) Snapshots of the 3D visualization of the reconstruction results using the 2D FBP algorithm in a slice-by-slice fashion with no correction (a) and correcting only for the center of rotation (b). (c) 3D reconstruction result using the calibrated system parameters. The intensities of images are normalized and saturated for better visualization.

Download Full Size | PDF

 figure: Fig. 7.

Fig. 7. Convergence and accuracy. (a) Evolution of the cost function of the reconstruction calculated based on 3 during the calibration process with (solid) and without (dashed) calibration. (b) The calibrated in-plane tilt angle against the true in-plane tilt angle for 4 different values indicated by the triangles. The dotted line represents $y=x$. The closer the triangles are to this line, the closer the calibrated angles are to the true angles.

Download Full Size | PDF

4. Conclusion

We presented a comprehensive study of certain geometrical artifacts in a poorly calibrated OPT imaging system. We summarized our results in the form of a catalog by combining existing results on the shifted center of rotation and our own contribution on the 3D tilt of the rotation axis and the inaccurate rotation angles. In doing so, we are able to explain the various types of mechanical artifacts, its appearance, cause, and properties, as well as the associated correction methods. This catalog serves as a reference for OPT practitioners to gain insight into their experimental setup and help them better calibrate their hardware system. To fill the vacancy of a versatile computational method to account for all types of mechanical artifacts, we propose an automatic calibration algorithm. It is based on a refined 3D mathematical model of the OPT imaging geometry that is able to characterize mechanical errors in the system. Moreover, the algorithm adapts to the large-size measurement datasets of OPT by performing the calibration on a coarse scale to overcome the computational bottleneck, while the final reconstruction of the 3D volume is achieved on a finer scale. Our multiscale calibration scheme has first been validated on the synthesized bead data where we simulate the shifted center of rotation and the imprecise rotation angles then successfully remove the resulting artifacts. We have further applied our algorithm on an experimental dataset of a fluorescent textile fiber that suffers from a tilted rotation axis, managing to detect the model mismatch and recover the true tilt angle. In the visualization of the final reconstructed volume, the associated artifacts are taken care of and we obtain a clean 3D image. In practice, to further improve the quality of the image, regularization may be introduced in the final reconstruction. In principle, our automatic geometric calibration framework can be adapted to other tomographic modalities, for instance, transmission OPT, optical coherence refraction tomography [36] or optical diffraction tomography [37].

Funding

European Research Council (No. 692726 GlobalBioIm); Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (200020_179217, 206021_164022).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data and code underlying the textile-fiber results are available in [38].

Supplemental document

See Supplement 1 for supporting content.

References

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

2. B. W. Lindsey and J. Kaslin, “Optical projection tomography as a novel method to visualize and quantitate whole-brain patterns of cell proliferation in the adult zebrafish brain,” Zebrafish 14(6), 574–577 (2017). [CrossRef]  

3. C. Schmidt, A. L. Planchette, D. Nguyen, G. Giardina, Y. Neuenschwander, M. D. Franco, A. Mylonas, A. C. Descloux, E. Pomarico, A. Radenovic, and J. Extermann, “High resolution optical projection tomography platform for multispectral imaging of the mouse gut,” Biomed. Opt. Express 12(6), 3619–3629 (2021). [CrossRef]  

4. Y. Wang and R. K. Wang, “Optimization of image-forming optics for transmission optical projection tomography,” Appl. Opt. 46(27), 6815–6820 (2007). [CrossRef]  

5. A. Bassi, L. Fieramonti, C. D’Andrea, G. Valentini, and M. Mione, “In vivo label-free three-dimensional imaging of zebrafish vasculature with optical projection tomography,” J. Biomed. Opt. 16(10), 1–4 (2011). [CrossRef]  

6. 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]  

7. J. R. Walls, J. G. Sled, J. Sharpe, and R. M. Henkelman, “Correction of artefacts in optical projection tomography,” Phys. Med. Biol. 50(19), 4645–4665 (2005). [CrossRef]  

8. J. McGinty, H. B. Taylor, L. Chen, L. Bugeon, J. R. Lamb, M. J. Dallman, and P. M. W. French, “In vivo fluorescence lifetime optical projection tomography,” Biomed. Opt. Express 2(5), 1340–1350 (2011). [CrossRef]  

9. V. C. Torres, C. Li, W. Zhou, J. G. Brankov, and K. M. Tichauer, “Characterization of an angular domain fluorescence optical projection tomography system for mesoscopic lymph node imaging,” Appl. Opt. 60(1), 135–146 (2021). [CrossRef]  

10. 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]  

11. F. Natterer, The Mathematics of Computerized Tomography (Society for Industrial and Applied Mathematics, USA, 2001).

12. C. Vinegoni, L. Fexon, P. F. Feruglio, M. Pivovarov, J.-L. Figueiredo, M. Nahrendorf, A. Pozzo, A. Sbarbati, and R. Weissleder, “High throughput transmission optical projection tomography using low cost graphics processing unit,” Opt. Express 17(25), 22320–22332 (2009). [CrossRef]  

13. P. P. Vallejo Ramirez, J. Zammit, O. Vanderpoorten, F. Riche, F.-X. Blé, X.-H. Zhou, B. Spiridon, C. Valentine, S. E. Spasov, P. W. Oluwasanya, G. Goodfellow, M. J. Fantham, O. Siddiqui, F. Alimagham, M. Robbins, A. Stretton, D. Simatos, O. Hadeler, E. J. Rees, F. Ströhl, R. F. Laine, and C. F. Kaminski, “OptiJ: Open-source optical projection tomography of large organ samples,” Sci. Rep. 9(1), 15693 (2019). [CrossRef]  

14. H. Zhang, L. Waldmann, R. Manuel, H. Boije, T. Haitina, and A. Allalou, “zOPT: an open source optical projection tomography system and methods for rapid 3D zebrafish imaging,” Biomed. Opt. Express 11(8), 4290–4305 (2020). [CrossRef]  

15. X. Tang, M. van’t Hoff, J. Hoogenboom, Y. Guo, F. Cai, G. Lamers, and F. Verbeek, “Fluorescence and bright-field 3D image fusion based on sinogram unification for optical projection tomography,” in 2016 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), (2016), pp. 403–410.

16. T. Donath, F. Beckmann, and A. Schreyer, “Automated determination of the center of rotation in tomography data,” J. Opt. Soc. Am. A 23(5), 1048–1057 (2006). [CrossRef]  

17. J. Michalek, “Total variation-based reduction of streak artifacts, ring artifacts and noise in 3D reconstruction from optical projection tomography,” Microsc. Microanal. 21(6), 1602–1615 (2015). [CrossRef]  

18. O. Koskela, T. Montonen, B. Belay, E. Figueiras, S. Pursiainen, and J. Hyttinen, “Gaussian light model in brightfield optical projection tomography,” Sci. Rep. 9(1), 13934 (2019). [CrossRef]  

19. V. Koljonen, O. Koskela, T. Montonen, A. Rezaei, B. Belay, E. Figueiras, J. Hyttinen, and S. Pursiainen, “A mathematical model and iterative inversion for fluorescent optical projection tomography,” Phys. Med. Biol. 64(4), 045017 (2019). [CrossRef]  

20. A. K. Trull, J. van der Horst, W. J. Palenstijn, L. J. van Vliet, T. van Leeuwen, and J. Kalkman, “Point spread function based image reconstruction in optical projection tomography,” Phys. Med. Biol. 62(19), 7784–7797 (2017). [CrossRef]  

21. U. J. Birk, A. Darrell, N. Konstantinides, A. Sarasa-Renedo, and J. Ripoll, “Improved reconstructions and generalized filtered back projection for optical projection tomography,” Appl. Opt. 50(4), 392–398 (2011). [CrossRef]  

22. G. C. Antonopoulos, D. Pscheniza, R.-A. Lorbeer, M. Heidrich, K. Schwanke, R. Zweigerdt, T. Ripken, and H. Meyer, “Correction of image artifacts caused by refractive index gradients in scanning laser optical tomography,” in Three-Dimensional and Multidimensional Microscopy: Image Acquisition and Processing XXI, vol. 8949 (2014), p. 894907.

23. Y. Liu, J. Dong, C. Schmidt, A. Boquet-Pujadas, J. Extermann, and M. Unser, “Artifacts in optical projection tomography due to refractive-index mismatch: Model and correction,” Opt. Lett. 47(11), 2618–2621 (2022). [CrossRef]  

24. J. van der Horst and J. Kalkman, “Image resolution and deconvolution in optical tomography,” Opt. Express 24(21), 24460–24472 (2016). [CrossRef]  

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

26. A. Darrell, H. Meyer, U. Birk, K. Marias, M. Brady, and J. Ripoll, “Maximum likelihood reconstruction for fluorescence optical projection tomography,” in 2008 8th IEEE International Conference on BioInformatics and BioEngineering, (2008), pp. 1–6.

27. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (Society for Industrial and Applied Mathematics, 2001).

28. 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]  

29. D. Ancora, D. D. Battista, G. Giasafaki, S. Psycharakis, E. Liapis, A. Zacharopoulos, and G. Zacharakis, “Optical projection tomography via phase retrieval algorithms for hidden three dimensional imaging,” in Quantitative Phase Imaging III, vol. 10074 (SPIE, 2017), p. 100741E.

30. A. A. Hendriksen, D. Schut, W. J. Palenstijn, N. Viganó, J. Kim, D. M. Pelt, T. van Leeuwen, and K. J. Batenburg, “Tomosipo: fast, flexible, and convenient 3d tomography for complex scanning geometries in python,” Opt. Express 29(24), 40494–40513 (2021). [CrossRef]  

31. F. Marelli and M. Liebling, “Optics versus computation: Influence of illumination and reconstruction model accuracy in focal-plane-scanning optical projection tomography,” in 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI), (2021), pp. 567–570.

32. M. Zehni, L. Donati, E. Soubies, Z. Zhao, and M. Unser, “Joint angular refinement and reconstruction for single-particle cryo-em,” IEEE Trans. on Image Process. 29, 6151–6163 (2020). [CrossRef]  

33. M. Unser, “Sampling—50 Years after Shannon,” Proc. IEEE 88(4), 569–587 (2000). [CrossRef]  

34. M. Nilchian, J. P. Ward, C. Vonesch, and M. Unser, “Optimized Kaiser–Bessel window functions for computed tomography,” IEEE Trans. on Image Process. 24(11), 3826–3833 (2015). [CrossRef]  

35. E. Soubies, F. Soulez, M. T. McCann, T.-a. Pham, L. Donati, T. Debarre, D. Sage, and M. Unser, “Pocket guide to solve inverse problems with GlobalBioIm,” Inverse Probl. 35(10), 104006 (2019). [CrossRef]  

36. K. Zhou, R. Qian, S. Degan, S. Farsiu, and J. Izatt, “Optical coherence refraction tomography,” Nat. Photonics 13(11), 794–802 (2019). [CrossRef]  

37. Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express 17(1), 266–277 (2009). [CrossRef]  

38. https://doi.org/10.5281/zenodo.7315108.

Supplementary Material (1)

NameDescription
Supplement 1       supplementary material on the derivation of the forward matrix

Data availability

Data and code underlying the textile-fiber results are available in [38].

38. https://doi.org/10.5281/zenodo.7315108.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1.
Fig. 1. Diagram of OPT. The light that passes through or is emitted by the sample gets imaged on a CCD camera using a lens system. In this ideal configuration, a slice of the sample perpendicular to the rotation axis corresponds to the information in the temporal sequence recorded by a row of pixels on the CCD camera.
Fig. 2.
Fig. 2. OPT imaging system with mechanical errors. (The lens system is omitted here for simplicity.) The cylinder represents the tilted sample that rotates around its orientation axis. The center of the cylinder is translated to $G$. The detector plane is described by a 2D coordinate system ($O'-\xi \eta$) perpendicular to the optical axis ($y$ axis). The projection of $O$ and $G$ in the detector plane are $O'$ and $G'$, respectively. The angle $\psi _1$ represents the out-of-plane tilt angle between the orientation axis of the sample and the $z$-axis direction (gray vertical arrow). The rotation angle $\psi _2$ represents the in-plane tilt angle of the detector plane around the optical axis.
Fig. 3.
Fig. 3. Artifacts in the reconstruction of simulated beads due to a controlled mismatch the geometric parameters of our refined model and the simplistic geometry assumed by traditional reconstruction methods. (a)-(b) Circle artifacts due to a constant shift of −4 pixels (a) and 8 pixels (b) of the COR in the $x$ direction. (c)-(d) Seagull artifacts due to a negative (c) and positive (d) random accumulative error with maximum amplitude 0.03 degrees in the rotation angles. (e)-(p) Various artifacts due to the tilt of the rotation axis. (e)-(h) Reconstruction at slice 192 (64 slices below the central slice). (i)-(l) Reconstruction at slice 256 (central slice in $z$-axis). (m)-(p) Reconstruction at slice 384 (128 slices above the central slice). The intensities of images are normalized and saturated for better visualization.
Fig. 4.
Fig. 4. Workflow of the multiscale calibration-reconstruction algorithm. We first perform a joint reconstruction-calibration at coarse scale. We then use the calibrated mechanical parameters for the final high-resolution reconstruction.
Fig. 5.
Fig. 5. (a) Reconstruction with circle and seagull artifacts at one slice (147) without calibration. (b) Reconstruction result at slice 147 after a naive shift correction. (c) Calibrated reconstruction at slice 147. Both the circle and seagull artifacts are successfully removed. (d) Ground truth at slice 147. The intensities of images are normalized and saturated for better visualization.
Fig. 6.
Fig. 6. Textile fiber data. (a)–(b) Snapshots of the 3D visualization of the reconstruction results using the 2D FBP algorithm in a slice-by-slice fashion with no correction (a) and correcting only for the center of rotation (b). (c) 3D reconstruction result using the calibrated system parameters. The intensities of images are normalized and saturated for better visualization.
Fig. 7.
Fig. 7. Convergence and accuracy. (a) Evolution of the cost function of the reconstruction calculated based on 3 during the calibration process with (solid) and without (dashed) calibration. (b) The calibrated in-plane tilt angle against the true in-plane tilt angle for 4 different values indicated by the triangles. The dotted line represents $y=x$. The closer the triangles are to this line, the closer the calibrated angles are to the true angles.

Tables (2)

Tables Icon

Table 1. Table of setup-related imaging artifacts and their dependence on the field of view (FOV) and on height (location along the z axis).

Tables Icon

Algorithm 1. Automatic Calibration

Equations (3)

Equations on this page are rendered with MathJax. Learn more.

b θ , t ( v ) = P θ { f } ( v t ) + n ,
b = H ( θ , t ) c + n ,
c , ( θ , t ) arg min c , θ , t { 1 2 H ( θ , t ) c b 2 2 } .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.