Accelerated Stimulated Raman Projection Tomography by Sparse Reconstruction from Sparse-View data

A. Forward modeling of SRPT

The SRPT imaging system has been described in our previous work [37]. Briefly, two Bessel beams generated by a pair of axicons and an objective lens were utilized as the pump and Stokes beams. Due to the considerable focal lengths of the Bessel beams [44], the detectable SRP signal may be an integration of SRS generation along with the Bessel beam Rayleigh length (or of the interaction length between the beams and the sample) [37].

in which S_SRP is the detected SRP signal by photodiode; d is the diameter of the Bessel beam’s central lobe; r₀ is the center of the Bessel beam’s central lobe; J₀(·) is the 0^th order Bessel function; B, β_p and β_S are constants which can be found in existing literature [37]; L₂ − L₁ indicates the thickness of the sample (within the Rayleigh length); C(r₀, z) is the distribution of the imaginary element of the third-order nonlinear susceptibility χ⁽³⁾ at the center of the Bessel beam’s central lobe along the beam propagation direction of z; P_p(r₀, z) and P_S(r₀, z) are the powers of the pump and Stokes beams respectively. In SRPT, because the sizes of the pump and Stokes beams are very small, the powers can be regarded as the average power on the focusing plane. Moreover, the central lobe intensity of the Bessel beams passing through the sample can be constant if the sample is always within the depth of focus of the beams. In this case, P_p(r₀, z) and P_S(r₀, z) can be recorded as P_p(r₀) and P_S(r₀).

In SRPT, both single point excitation and single point detection were employed, while the connection between the excitation and detection points were approximated as being perpendicular to the sample plane. Therefore, the SRPT can be considered as parallel beam imaging. By scanning the Bessel beams two-dimensionally on the lateral plane, we can generate a projection image at a certain rotation angle. This generated projection image can be calculated according to the following equation by changing the position of r:

where S(r) is the SRP signal at position of $\mathbf{r};J(\mathbf{r})=P_p(\mathbf{r})P_S(\mathbf{r}){\int }_0^{\frac{d}{2}}{\mathbf{r}}^{\prime }{J}_0^2\left({\beta }_p({\mathbf{r}}^{\prime }-\mathbf{r})\right)J_0^2\left({\beta }_S({\mathbf{r}}^{\prime }-\mathbf{r})\right)d{\mathbf{r}}^{\prime };C(\mathbf{r},\mathbf{z})$ is the imaginary element of the third nonlinear susceptibility at position of (r, z); (r, z) are coordinates for the lateral and propagation directions of the Bessel beams.

By rotating the sample with a step of a small angle, a number of angular-dependent projection images can be found and utilized to reconstruct the SRPT image:

where $S(\mathbf{r};\theta ),J(\mathbf{r};\theta ),L_1(\theta )\text{and}{L}_2(\theta )$ are the corresponding values at different projection angles.

Equation (3) demonstrates the forward model of SRPT. Using this equation, we can determine angular-dependent projection images as long as we are aware of all the relevant information about the laser beams and the sample. The reconstruction of SRPT recovers the distribution of the imaginary element of the third-order nonlinear susceptibility (C(r, z)) from the measured angular-dependent projection images (S(r; θ)). As can be seen in the following simulations, 180 angular-dependent projection images that were utilized for SRPT reconstruction were generated by projecting the sample to a specific plane while rotating it consecutively at one-degree steps along an axis, using Eq. (3).

B. TV-regularized SART algorithm

Eq. (3) denotes a line integral of C(r, z) from L₁ to L₂ along the beam propagation direction. Such line integral represents the Radon transform of the unknown or reconstructed quantity. Thus, the forward model of SRPT can be described by the following matrix equation:

where S is the measured SRP image; R is the Radon transform; C is the (unknown) distribution of the imaginary element of the third-order nonlinear susceptibility.

If we acquire a complete data set, the recovery of C from the noisy measurement S can be achieved by applying the FBP algorithm. If the measurements are incomplete, the SRPT images can be reconstructed by solving the following optimization problem:

where $\tilde{\text{C}}$ is the reconstructed distribution of the imaginary element of χ⁽³⁾; λ is the regularization parameter and ${\Phi }_{TV}(\text{C})$ is the TV regularization term [38]:

where the subscripts s, t are the pixel indices in matrix C.

The minimization problem in Eq. (6) can be solved by using various optimization algorithms, such as the steepest descent method and the conjugate gradient method. Here, we used the steepest descent method to determine the steepest descent direction for the TV based minimization problem in Eq. (6) [38]:

where э is a very small positive number. The steepest descent direction can be expressed by the vector d and its unit vector is denoted as $\widehat{\text{d}}$.

The solution to the minimization problem in Eq. (5) can be obtained using the SART, which is iterated as [45]:

where I_φ represents the set of index of rays passing through the jth pixel of the reconstructed image when projecting from angle $\phi ;r_{ij}$ denotes the weight of the jth pixel contributing to the ith ray; N is the total number of pixels in the reconstructed image; R_i is the ith row of the projector; k is the number of corresponding iterations; α is the relaxation factor.

The step-by-step procedure of the TV-regularized SART algorithm can be summarized as [42] [46] [47]:

1. Initialization: C = 0, C₀ =C;

2. SART iteration reconstruction using a projection model;

3. Positivity constraint is imposed on the reconstructed image: C_j = 0 when C_j < 0

4. Compute the residual: d_img = |C − C₀|;

5. Compute TV minimization term by repeating the following operation N_grad times:

   • Computing d and $\widehat{\text{d}}$ by using Eq. (7);
   • Update: C ← C − $\lambda d_{img}\widehat{\text{d}};$

6. C₀ = C;

7. Repeat Step 2 to Step 6 until there is no considerable change between two adjacent iterations.

There are several parameters which are both crucial and essential, including regularization parameter λ, relaxation factor α, iteration number N of SART, and iteration number N_grad of TV. Existing research shows that we can empirically determine these parameters [38][39][46][47][48][49][50]. To conduct the comparisons under the same condition, parameters were set as: λ = 0.08, α = l/(l + 0.5(n − 1)) (n is the current iterative number of SART), N=50, and N_grad = 10 in the simulations and experiments which follow.

C. Pixel vertex driven projection model

In the TV-regularized SART algorithm, we needed to repeatedly calculate projectors during iterative process, as stated in Step 2. The accuracy and computational burden of such projectors are major challenges for iterative reconstruction. The pixel vertex driven model (PVDM) was deduced from the separable footprint (SF) projector by assuming that the parallel beam projection was performed in SRPT. The SF projector has been reported to be more accurate and efficient than that of the distance-driven model (DDM), as can be seen in Supplementary Figure S1, whose capability has been demonstrated in the cone-beam computed tomography [43]. Here, we made some modifications on the SF projector so that it can be applied to parallel beam SRPT. Such PVDM was used to determine the projector of R in Eq. (5) or (8).

Figure 1 shows the diagram of the PVDM projection model for the parallel beam illumination, where the red rectangle frames the projection of one pixel on the detector. In the forward projector, projecting four vertices of one pixel onto the detector along the beam illumination direction provides the ladder-shaped footprint. In some particular cases, the shape of the footprint can be a rectangle at the projection angles of 0°, 90°, 180°, and 270°, or it can be a triangle at the angles of 45°, 135°, 225°, and 315°. Such diverse projection shapes make the PVDM best approximate the actual physical process of parallel beam projection imaging. The footprint intersects with the detector, and the intersection area denotes the weight of the pixel to the detector. The whole area of one footprint is set to be unity. The weight and the coordinates of the pixel are then recorded for the subsequent back projection. There are ten kinds of intersections for the determination of weight values, as shown in Fig. 2. The calculation of weight values for the ten intersections are detailed in the Appendix.

D. Assessment of image quality

Two evaluation indices, structural similarity (SSIM) and root mean square error (RMSE), were used to evaluate our sparsely reconstructed results. SSIM compares the local patterns of normalized pixel intensities in terms of luminance and contrast [51], [52]. The closer the SSIM is to unity, the better the sparsely reconstructed result is. The SSIM can be calculated as:

where μ_R and μ_A are the average values of the reconstructed and reference images, respectively; σ_R and σ_A are the standard deviations of the reconstructed and reference images, respectively; σ_RA is the covariance between them; c₁ and c₂ are the constants related to the dynamic ranges of the images.

The RMSE is a type of numerical index that can be used to measure the accuracy of an image restoration. The smaller the RMSE, the more accurate the sparsely reconstructed image. RMSE can be calculated as follows [53]:

in which C_xy and $C_{x,y}^{*}$ are the reconstructed and reference values of the imaginary element of χ⁽³⁾ on each pixel and N is the pixel number. Here, the original image was used as the reference image for simulation comparisons, and the image reconstructed by using 180 projections was selected as the reference image for the single cell experiment.

A. Simulation verification

The validity of the proposed TV-regularized SART algorithm was evaluated with two groups of numerical simulations. In the evaluation, the data sets employed for sparse reconstruction were constructed by equally sampled complete projections in a half circle, consisting of 180, 90, 60, 45, 36, 30, 15, 9, and 6 projections respectively. For comparisons, the DDM based TV-regularized SART algorithm and the FBP algorithm were applied as the references. A Ramp filter was utilized in the FBP algorithm [54]. For a fair comparison, the same data set was utilized for all three algorithms. The data set was generated by using discrete form of Eq. (3), whose process was very similar to ray-driven model. The generation process can be described as follows. First, by fixing the position of r and the projection angle of θ, the SRP signal can be calculated by integrating Eq. (3) along a light ray. Second, by changing the position of r two-dimensionally, a SRP image can be generated by calculating the SRP signal at each position. Last, a data set of angular-dependent SRP images can be obtained by changing the projection angle of θ and generating the SRP image at each projection angle. To ensure the accuracy of forward projection data, we discretized the detector into 1024 pixels.

First, a modified Shepp-Logan model was used to perform the evaluation simulations [17]. Figure 3(a) shows the original image of the modified Shepp-Logan model, and the reconstructed images by three algorithms using 180, 60, 30, 15, and 9 projections are presented in Fig. 4. Figure 4(a1-a5) depicts the reconstructed images of the PVDM based TV-regularized SART algorithm, Fig. 4(b1-b5) shows those of the DDM based TV-regularized SART algorithm, and Fig. 4(c1-c5) displays those of the FBP algorithm. We observed that almost the same performance was obtained by the PVDM and the DDM based TV-regularized SART algorithms, which were much better than that obtained by the FBP algorithm. Furthermore, we also found that good quality of the reconstructed images could be maintained even when the projection number was decreased to 30. When the number of projections was reduced to 15, there was some blur observed in the reconstructed images and small details could not be distinguished. Even worse, no acceptable images could be reconstructed when the number of projections was reduced to 9. To quantitatively evaluate the quality of the reconstructed images, we calculated the SSIM and RMSE values. Figure 5(a) and 5(b) show the changes of the SSIM and RMSE values as a function of the projection number respectively. Almost the same performance conclusions were drawn from the SSIM and RMSE curves. The SSIM and RMSE values changed slightly when the projection number was no smaller than 30 for iterative algorithms. However, it decreased (SSIM) or increased (RMSE) dramatically when the number of projections was smaller than 15. Importantly, the proposed PVDM based TV-regularized SART algorithm provided better results than the DDM based one in the sparse views case, for example the projection number was not larger than 60 in this simulation.

Second, an Arabidopsis silique model was utilized as a physical model in order to generate forward projection images and measure the performance of the proposed algorithm further. Figure 3(b) shows one cross-section of the original image of Arabidopsis silique model. Figure 6 presents the images reconstructed from the sparsely sampled data sets consisting of, respectively, 180, 60, 30, 15, and 9 projections, where the left column shows the results of the PVDM based TV-regularized SART algorithm, the middle column displays the DDM based results, and the right column shows those obtained from the FBP algorithm. Once again, better performance was obtained by applying the iterative algorithms than the FBP. The distribution and shape of the targeted objects could be well recovered by using as little as 30 projections. Even when the projection number was reduced to 15, the reconstructed image still showed little blur. A little distorted or displaying artifact interference was aroused when the number of projections was reduced to 9. The calculated SSIM and RMSE values shown in Fig. 5(c) and 5(d) also gave the same conclusion. It should be noted that the PVDM based TV-regularized SART algorithm provided better results than the DDM based one when the projection number was smaller than 30, showing the good ability in sparse reconstruction from sparse-view data. These results collectively indicate that our TV-regularized SART algorithm can provide satisfactory or acceptable reconstructed images (SSIM>0.9) even when the number of projections is reduced to 15, leading to a 12-fold improvement in imaging speed compared to the current SRPT scheme.

Although the TV-regularized iterative algorithms performed better than the FBP algorithm especially in the case of sparse- view data, they have over-smoothing effects on tiny details. We can see that the images reconstructed by the FBP algorithm look quite sharp in the few-view cases. Even with 15 or nine views, when the images reconstructed by the TV-regularized SART algorithm have been blurred, the edges of the structures in the FBP images are still visible, despite the conspicuous streaks. In the modified Shepp-Logan phantom, there are six dots near to the right side of the skull that are still visible in 15- and nine-view reconstructions, as shown in Figs. 4(c4) and 4(c5). In order to observe this phenomenon more intuitively, we extracted a line across two of six dots in the modified Shepp-Logan phantom, as labeled in Fig. 7(a). Figures 7(b)-7(f) plot the extracted lines of the three reconstruction algorithms by using 180, 60, 30, 15, and 9 projections respectively. The same conclusion can be drawn; which is that the FBP algorithm had a better detail resolution ability in few-view cases. This may be caused by the use of TV regularization, which may over-smooth edges and details [42] [55] [56]. This smoothing effect was influenced by the parameters of the TV regularization. We investigated this potential influence by changing the parameter of N_grad. We discovered that the larger the N_grad, the more smoothing the effect would be. If the N_grad was too small, the quality of the reconstructed image would be worse (see Supplementary Figure S2).

B. Applicability-potential demonstration

To further explore the potential of the proposed fast SRPT for live cell imaging, we three-dimensionally imaged lipid droplets (LDs) in adipose cells. Here, differentiated 3T3-L1 cells were used for testing [37]. As a mouse derived cell line, 3T3-L1 is widely used for adipogenesis studies [57], and can differentiate into adipose cells, which accumulate LDs in the cytosol under certain conditions. After the cells were cultured for some time, they were first mixed with 1.5% (w/v) agarose gel in PBS solution. The gel mixture was allowed to set in a cylindrical capillary tube with an inner diameter of ~50 microns. The tube was mounted onto a homemade sample holder for imaging during sample rotation [37]. The rotation was performed at 1^o per step for acquisition of 180 projection images. The image pixel dwell time was set to 10 μs and each image contained 300 × 300 pixels. Other imaging parameters were setup as described in our previous publication [37].

The collected projection images were first pre-processed for de-noising and compensation, and then used for reconstructing the distribution of lipids in the cell. By comparing the simulation results with the original image, we had found that the image reconstructed by the iterative algorithms had much better quality than that reconstructed by the FBP algorithm, with the SSIM values being higher even when using 180 projections. Therefore, the image reconstructed by the PVDM based SART algorithm from 180 projections was used as the reference to evaluate the performance of the reconstruction from the sparsely sampled data. Figure 8(a) shows the image reconstructed by the FBP algorithm with 180 projection images, and Fig. 8(b)-(h) show the reconstructed cell and distribution of lipids inside the cell by using 180, 90, 60, 36, 30, 15, 9 projections, respectively. We observed that the quality of the reconstructed image was well maintained until the projection number was reduced to 30. A little distortion could still be seen when the number of the projection was decreased to 15, with little differences observed in the fine structures at the edge of the cell. Obvious differences were observed in morphology until the number of projections was reduced to 9. The SSIM and RMSE values as a function of the projection number were also calculated and shown in Fig. 9, and these results also indicated that a minimum projection number of 30 was required to remain a relatively good reconstruction fidelity. The SSIM value between the images which are reconstructed by the FBP and the PVDM based SART algorithms was 0.9551, and the RMSE between them was 0.3978. By combining the analysis results shown in Fig. 8 and Fig. 9, we could obtain an acceptable reconstructed image when the number of projections was set at around 20. These cell-based experimental results demonstrate the applicability potential of the proposed PVDM based SART algorithm for fast imaging of a single adipose cell. Currently, the image reconstruction scheme used in SRPT is based on the FBP algorithm that requires 180 projection images to rebuild a volumetric image [37]. Thus, the algorithm improved the imaging speed by a factor of 9–12 without sacrificing discernible imaging details.