Accelerating nonlinear reconstruction in laminar optical tomography by use of recursive SVD inversion

Mengyu Jia; Jingying Jiang; Wenjuan Ma; Chenxi Li; Shuang Wang; Huijuan Zhao; Feng Gao

doi:10.1364/BOE.8.004275

1. Introduction

Functional mesoscopic optical tomography, such as the laminar optical tomography (LOT) or the quantitative optoacoustic mesoscopy, bridges the gap of spatial-resolution and penetration-depth between microscopic and macroscopic imaging, and thus is emerging as a powerful clinical tool to detect chromophores (hemoglobin, melanin, water, etc.) and fluorochromes (endogenous fluorophores and exogenous agents) in superficial tissues [1,2]. In addition, miniaturized mesoscopic instrument has the potential of performing endoscopic imaging [3,4].

As a model-based biophotonic imaging technique, LOT is an adaptation of the epi-illumination diffuse optical tomography (DOT) principle to mesoscopic-scale applications [1]. Image reconstruction in LOT is inherently an ill-posed nonlinear inverse problem. Although the Newton-type iteration has been widely developed as a canonical method for DOT nonlinear reconstruction, and greatly facilitated by the adoption of the diffusion approximation (DA) model, its implementation in the high-resolution LOT is computationally difficult since a dense sampling and a physically-rigorous photon-transport model are desirable in this scenario [5–8]. Take the raster-scanning LOT for instance, to achieve 200-µm spatial resolution, at least 50 × 50 scanning spots (or source sites) per square centimeter are requisite, which normally leads to 2500 forward calculations in each iteration stage and pseudo-inversion of a 2500-row sensitivity matrix [1]. As a compromise, LOT generally employs a one-step linear reconstruction scheme, where the forward modeling is performed merely for a specific scanning spot with the assumed “homogeneous background”, and obtained for the rest scanning spots via translational transforms [1,4,6]. However, the linear reconstruction has been demonstrated to only accommodate small-perturbation scenarios with an accurate enough initial guess of the background [9–11].

In order to efficiently perform the nonlinear reconstruction for improved quantitation and resolution, efforts have been paid on enhancing both the forward and inverse calculations in DOT and can be naturally applied to LOT. Parallel implementation of the photon-migration modeling is typical strategy for the forward acceleration, where a combined multicore CPU and multistream GPU has been used [12–14]. For the scenarios depending on the delicate transport theory, the scaling scheme has been considered to be skillful technique that significantly speeds-up the repeated forward predictions through substituting the original re-modeling process with calibration operations, based on the established scaling relations between the perturbed and unperturbed results [15, 16]. By taking advantage of the scaling strategy, the forward model in the traditional one-step linear reconstruction for LOT has evolved into a recursive computation form for the nonlinear scheme to accommodate scenarios with “optically heterogeneous background” [17].

More works have been also devoted to improving the large-scale DOT inverse calculation. A gradient-optimization based solution has been proposed to circumvent the low-efficient inversion of the incurred large-scale and ill-conditioned sensitivity matrix (i.e., the linear inversion).The desired unknowns are obtained by using an iterative line minimization, where the searching direction is set to be the gradient of the objective function [18, 19]. However, it generally takes multiple iterations to reach the minimization. In contrast to the optimization strategy, the inverse problem is (piecewisely) linearized in the computationally robust perturbation method, i.e., the Newton-Raphson scheme, and the parameters updating is achieved through linear inversion [1, 5, 6]. Many approaches have been developed to reach a computationally and physically efficient inversion. One most commonly-utilized strategy is to reduce the model order. For instance, the adaptive meshing scheme refines meshes only in regions of interest to ensure necessary resolution while keeping the reconstruction computationally tractable [20, 21]. In view of the great sparsity of sensitivity matrix, truncated measurement sensitivity was proposed to perform selective generation only for those elements larger than a predefined threshold-value of sensitivity [22, 23]. Another typical method that has been demonstrated in the high-density DOT is development of the Schwarz-type domain decomposition (DD), to partition the whole domain into several overlapped subdomains and perform both the forward and inverse calculations over the subdomains in an alternating way [24]. The scheme can be in nature significantly accelerated with a combined multi-CPU and GPU parallelization [14]. With the downscaled sensitivity matrix, the linear inversion can be more efficiently achieved in a whole-matrix-fashioned way, e.g., the Tikhonov regularization or the truncated Singular Value Decomposition (SVD).

For the ill-posed linear inversion in the Newton-Raphson schemed DOT, SVD has been considered a powerful numerical tool, since the threshold of efficient singular value, as a significant reference value of regularization parameter, can be practically determined according to the systematic signal-to-noise ratio (SNR) [25, 26]. However, calculation of a complete SVD for both the singular values and vectors is time-consuming, and typically leads to a computational complexity of O(4MN² + 22M³) and O(4MN² + 8NM² + 9M³) for the Golub-Reinsch-SVD and R-SVD of a M × N (M≤N) matrix, respectively [27, 28]. As such, Yalavarthy et al. proposed a linear iterative reconstruction without iteratively updating the sensitivity matrix and its SVD terms, and demonstrated its conditional feasibility and superiority to the linear reconstruction [11].

In this work, a practical nonlinear LOT scheme for the absorption reconstruction is proposed within the framework of the Newton-Raphson linear iteration, where the original set of the governing equations is commonly linearized with a Taylor series expansion and iteratively solved based on the SVD-based linear inversion. To mitigate the requirement for SVD at each iteration, an explicitly recursive SVD version is derived from the scaling nature of the sensitivity computation, sharply reducing the overall computation complexity of the SVD-based linear inversion to O(NM² + MN²) ( $M ≪ N$ ). The proposed reconstruction scheme has been validated numerically and experimentally, and shown to be equivalent to the traditional nonlinear one.

2. Methods

In the rest of this paper, the element of a matrix or vector $ζ$ is referred to as $ζ_{i j}$ or $ζ_{i}$ , respectively; $ζ^{T}$ denotes the transpose of matrix $ζ$ ; $R^{M \times N}$ represents the set of M-by-N real matrices.

2.1 SVD-based framework for nonlinear inversion

To mitigate the ill-conditioned nature of the inverse problem, the linearized inversion for updating the absorption coefficient at k-th (k>0) iteration is typically regularized with the Tikhonov strategy, leading to the following updating equation [1,5,6]

δ μ_{a}^{(k)} = {[{(J^{(k)})}^{T} J^{(k)} + α^{(k)} I]}^{- 1} {(J^{(k)})}^{T} b^{(k)},

where

μ_{a}^{(k)} \in R^{N}

denotes the absorption perturbations at the N meshing voxels;

J^{(k)} \in R^{M \times N}

is the sensitivity matrix; the regularization parameter

α^{(k)}

controls the influence of the model mismatch, noise and systematic errors on the solutions of

μ_{a}^{(k)}

;

b^{(k)} \in R^{M}

indicates the residuals between the boundary measurements

Γ

and the model predictions

F^{(k)}

, i.e.,

b^{(k)} = Γ - F^{(k)}

.

M ≪ N

is generally achievable in LOT reconstruction. To perform a reliable regularization, we transform Eq. (1) into the following form of the SVD inversion

δ μ_{a}^{(k)} = V^{(k)} \frac{Λ^{(k)}}{{(Λ^{(k)})}^{2} + α^{(k)} I} {(U^{(k)})}^{T} b^{(k)},

where

U^{(k)} \in R^{M \times M}

and

V^{(k)} \in R^{N \times N}

are orthonormal matrices containing the left and right singular vectors of

J^{(k)}

, respectively, and

Λ^{(k)} \in R^{M \times N}

is a diagonal matrix with the non-negative diagonal elements being the singular values of

J^{(k)}

, i.e.,

J^{(k)} = U^{(k)} Λ^{(k)} {(V^{(k)})}^{T}

. By such an expression, it is now clear to set

α^{(k)}

equal to the lower-bound of the efficient singular values, which is practically determined as per the systematic SNR.

On account of the high complexity of computing SVD, we propose herein to achieve $U^{(k)}$ , $Λ^{(k)}$ and $V^{(k)}$ with an explicitly recursive scheme, respectively. Without loss of generality, the recursive expression is formed as

Θ^{(k + 1)} = δ Θ^{(k)} + Θ^{(k)},

where

Θ \in {U, V, Λ}

. The small enough perturbation

δ μ_{a, n}^{(k)}

holds the following approximation

δ Θ_{i, j}^{(k)} \approx \sum_{n = 1}^{M} \sum_{m = 1}^{N} \sum_{l = 1}^{N} \frac{\partial Θ_{i, j}^{(k)}}{\partial J_{m, l}^{(k)}} \frac{\partial J_{m, l}^{(k)}}{\partial μ_{a, n}^{(k)}} δ μ_{a, n}^{(k)} .

2.2 Recursive SVD scheme for sensitivity matrix

In view of the photon-migration scaling principle, the forward prediction for the voxelized absorption perturbations $δ μ_{a, n}^{(k)}$ can be iteratively updated as formulated below

F_{m}^{(k + 1)} = F_{m}^{(k)} \prod_{n = 1}^{N} e^{- s_{m, n}^{(k)} δ μ_{a, n}^{(k)}},

where

s_{m, n}^{(k)}

(

0 < m \leq M

) denotes the averaged pathlength of the detected photons inside the nth voxel under the mth source-detector configuration [15,16]. For the proof-of-concept survey,

F_{m}^{(0)}

is obtained through solving the delicate radiative transport equation (RTE) with the discrete-ordinates-method (DOM). Combing Eq. (5) and the definition formula of the measurement sensitivity, we have (Appendix A)

J_{m, n}^{(k + 1)} = J_{m, n}^{(k)} e^{- s_{m, n}^{(k)} δ μ_{a, n}^{(k)}} .

Note that $\partial J_{m, l}^{(k)} / \partial μ_{a, n}^{(k)} = 0$ in case of $l \neq n$ . Inserting Eq. (6) into Eq. (4), we have

δ Θ_{i, j}^{(k)} \approx \sum_{m = 1}^{M} \sum_{n = 1}^{N} \frac{\partial Θ_{i, j}^{(k)}}{\partial J_{m, n}^{(k)}} δ J_{m, n}^{(k)},

where

δ J_{m, n}^{(k)} = J_{m, n}^{(k)} (e^{- s_{m, n}^{(k)} δ μ_{a, n}^{(k)}} - 1)

, and

\partial Θ_{i, j}^{(k)} / \partial J_{m, n}^{(k)}

is detailed to the following explicit expressions [27,28]

{\begin{cases} \frac{\partial Λ_{i}^{(k)}}{\partial J_{m, n}^{(k)}} = u_{m, i}^{(k)} v_{n, i}^{(k)} \\ \frac{\partial U^{(k)}}{\partial J_{m, n}^{(k)}} = U^{(k)} ϕ^{(k)} (m, n) \\ \frac{\partial V^{(k)}}{\partial J_{m, n}^{(k)}} = - V^{(k)} φ^{(k)} (m, n) \end{cases},

with

ϕ^{(k)}

and

φ^{(k)}

being the solution to the following 2 × 2 matrix equation

[\begin{matrix} Λ_{j}^{(k)} & Λ_{i}^{(k)} \\ Λ_{i}^{(k)} & Λ_{j}^{(k)} \end{matrix}] [\begin{matrix} ϕ_{i, j}^{(k)} (m, n) \\ φ_{i, j}^{(k)} (m, n) \end{matrix}] = [\begin{matrix} u_{m, i}^{(k)} v_{n, j}^{(k)} \\ - u_{m, j}^{(k)} v_{n, i}^{(k)} \end{matrix}] .

It is noteworthy that, Eq. (9) becomes ill-conditioned as

Λ_{i}^{(k)}

is close to

Λ_{j}^{(k)}

, and the Tikhonov strategy is used to regularize the linear inversion. Substituting Eqs. (7) and (8) into Eq. (3) produces the following matrix-fashioned recursive expressions

{\begin{cases} Λ^{(k + 1)} = d i a g (Ω^{(k)}) + Λ^{(k)} \\ U^{(k + 1)} = U^{(k)} [ξ^{(k)} \otimes Ω^{(k)} - η^{(k)} \otimes {(Ω^{(k)})}^{T} + I] \\ V^{(k + 1)} = V^{(k)} [ξ^{(k)} \otimes {(Ω^{(k)})}^{T} - η^{(k)} \otimes Ω^{(k)} + I] \end{cases},

where the operator “⊗” denotes the Hadamard/Schur product,

Ω^{(k)} = {(U^{(k)})}^{T} δ J^{(k)} V^{(k)},

and

[\begin{matrix} ξ_{i, j}^{(k)} & η_{i, j}^{(k)} \end{matrix}] = \frac{[\begin{matrix} {(Λ_{i}^{(k)})}^{2} + {(Λ_{j}^{(k)})}^{2} + {α^{'}}^{(k)} & - 2 Λ_{i}^{(k)} Λ_{j}^{(k)} \end{matrix}]}{{[{(Λ_{i}^{(k)})}^{2} + {(Λ_{j}^{(k)})}^{2} + {α^{'}}^{(k)}]}^{2} - 4 {(Λ_{i}^{(k)} Λ_{j}^{(k)})}^{2}} [\begin{matrix} Λ_{i}^{(k)} & Λ_{j}^{(k)} \\ Λ_{j}^{(k)} & Λ_{i}^{(k)} \end{matrix}] (1 - Δ_{i j}),

with Δ_ij being the Kronecker delta function, and

α'^{(k)}

the Tikhonov regularization parameter. The most expensive step in the recursive calculation is computing Eqs. (10) and (11), making the leading order of updating

δ μ_{a}^{(k)}

as O(2M²N + MN²).

3. Simulative investigations

The performance of the proposed method is evaluated using the simulated data for a semi-infinite tissue model at mesoscopic scale. The model contains a stick-like heterogeneity target, with only absorption contrast, as shown in Fig. 1. The cross-section area of the heterogeneity is 1 × 1mm². The background optical properties of the model is μ_a0 = 0.01 mm⁻¹, μ_s0 = 2.5 mm⁻¹, and g = 0.85. 100 × 100 raster-scanning spots of 0.15-mm diameter are uniformly distributed on the model surface covering a field of view (FOV) of 8.3 × 8.3mm², i.e., the scanning interval is 0.083 mm. There are totally 8 detecting sites corresponding to each scanning spot, each with an increment of ~0.24 mm in the source-detector separation (SDS). To achieve an imaging depth of ~1.7 mm and also to eliminate the residual specular reflections in practice [1], the maximum and minimum SDSs are set to ~2.19 mm and ~0.25 mm, respectively.

Fig. 1 Tissue model for generating simulated data.

Download Full Size | PDF

The simulated data to the raw measurement are generated by the MC method. To be realistic, a Gaussian-type noise is added to obtain a conventionally achievable signal-to-noise ratio (SNR) of 30 dB in the LOT measurements [4]. In reconstruction procedure, 100 × 100 raw measurements for a fixed SDS are decimated into 20 × 20 with a low-pass filter. To achieve a mesoscopic resolution of ~200 μm for the reconstruction, the computation domain of 8.3 × 8.3 × 3mm³ is discretized into 60 × 60 × 43 cubic voxels. With this meshing scheme, no evident numerical error would be observed according to our previous practice. The above settings for the downsampling ratio and discretization lead to a large-scale updating equation [i.e., Eq. (1)] with a 3200 × 154800 coefficient/sensitivity matrix. The regularization parameters in Eq. (1) are set as $α = 0.1, 0.05, a n d 0.01$ for target depths being 0.4mm, 0.9mm, and 1.2mm, respectively [1]. For a fair comparison, the iteration number of the Newton-Raphson optimization in all the following reconstructions is set to 50. As an example, only the reconstructed Y-Z and X-Y sectional images are illustrated. It is noteworthy that, the Y-Z images is an averaged result from x = 0mm to 8.3mm to achieve numerical stabilization. All reconstructions were performed on a workstation with a 3.60 GHz CPU (Intel Core i7-3820) and a GPU card (Nvidia GeForce GTX 970).

To evaluate the efficiency and accuracy of the proposed method, we compare the reconstructions using the one-step linear (‘linear’ for short), linear iterative, proposed, and traditional methods, in terms of computational time, quantitativeness ratio (QR), and full width at half-maximum along X-, Y-, or Z- axis (referred to as FWHM_ax with the subscript ‘ax’ being ‘x’, ‘y’, or ‘z’). To quantify the performance of position reconstruction, we introduce a measure, defined as $ε_{a x} = | γ_{a x} - γ_{a x, r e a l} | / γ_{a x, r e a l}$ , with $γ_{a x}$ being the distribution centroid along the referred axis of the reconstructed μ_a image. With the definition, $ε_{a x} = 0$ indicates no position error, which, however, is hard to achieve in DOT/LOT due to the inherent skin-effect [4, 29].

3.1 Comparisons for model with different target-absorption-contrast

The model in Fig. 1 is employed to compare the four reconstructions for different target-absorption-contrast (TAC). The simulated data are generated for TACs of 2, 3 and 5, i.e., μ_a = 0.02, 0.03 and 0.05 mm⁻¹, respectively, and a fixed target depth of 0.9mm. Figure 2 shows the reconstructed μ_a images and the corresponding profiles along the Z-axis (Z-profile) at y = 4.165mm (indicated by the white dashed lines). Table 1 presents the FWHM_z, QR, and $ε_{z}$ of the reconstructed target calculated from the images. Furthermore, the computation time per iteration taken to evaluate Eq. (1) is also calculated, which is averaged for all the Newton-Raphson iterations in one reconstruction scheme. The computation time will be fixed in all the following experiments, as long as the reconstruction-associated parameters, e.g., measurement configuration and discretization scheme, are unchanged.

Fig. 2 Reconstructed images and the corresponding Z-profiles along the white dashed lines with the four schemes, for (a) TAC = 2, (b) TAC = 3, and (c) TAC = 5, respectively, with a fixed target depth of 0.9mm. The black dotted boxes indicate the original targets.

Download Full Size | PDF

Table 1. Comparisons among the reconstructions for three TACs

View Table | View all tables in this article

It is observed that, the linear reconstruction show the poorest performance, e.g., lower QR, larger FWHM_z, and more obvious skin-effect. For the linear reconstruction achieved by solving Eq. (1) with k = 1, $δ μ_{a}$ depends on the regularization parameter, which generally has to be changed with the imaging depth [1,4,30]. As a result, the position-related measures and QR are hard to be simultaneously optimized through the one-step procedure. For the linear iterative scheme, although it spends less computation time than the proposed method dose, it has lower QR, larger FWHM_z and $ε_{z}$ . In general, the linear iterative reconstruction leads to a larger solution error in cases where the perturbation is large, e.g., an initial guess is not close to the real value [9–11]. However, non-linear approaches allow larger perturbations to be adequately modeled by imposing iterative updates on searching direction, i.e., the derivative $J_{m, n}^{(k)}$ [9]. In comparison to the traditional non-linear one, the proposed nonlinear reconstruction shows a similar performance, while taking merely 26.7% computation time. Seeing from the reconstructed Y-Z images, the bottom edge of the target is smoothed and blurred more severely than the top, which is a common problem inherent in LOT or DOT and can be physically ascribed to the increased number of detected photon scatterings with deepening target depth [1,4,30].

3.2 Comparisons for model with different target depth

The imaging depth of the four methods is also compared. Figures 3(a) to 3(c) show the reconstructed μ_a images for target depths of 0.4 mm, 0.9 mm, and 1.2 mm, respectively, as well as the corresponding Z-profile at y = 4.165 mm, with a fixed TAC of 5, i.e., μ_a = 0.05 mm⁻¹. Table 2 shows the corresponding QR, FWHM_z and $ε_{z}$ of the reconstructions calculated from the images.

Fig. 3 Reconstructed images and the corresponding Z-profiles along the white dashed lines with the four schemes, for the target depths of (a) 0.4mm, (b) 0.9mm, and (c) 1.2mm, respectively. The black dotted boxes indicate the original targets.

Download Full Size | PDF

Table 2. Comparisons among the reconstructions for different target depths

View Table | View all tables in this article

We can see that, for the case of superficial target, e.g., depth = 0.4mm, all the four schemes show accurate position-associated reconstructions. With increase of target depth, image quality reduces with rising FWHM_z and $ε_{z}$ . Furthermore, artifacts can be found near the surface, as indicated with white arrowheads in Fig. 3(c), which might be attributed to the twofold aspects below. On one hand, slight differences between the raster-scanning positions in the meshing for forward calculation and those in the simulated data may produce serious artifacts [30]. On the other hand, with increase of target depth, the regularization parameter extremely reduces, imposing numerical instability on solving Eq. (1) [30]. For a typical LOT, the highest measurement sensitivity along Z-axis locates at several mean-free-pathlength below incident sites due to the Gaussian shape of snapshotting incident photon-density along Z-axis [31]. Consequently, the artifacts shown in Y-Z cross-sections emerge at a small depth.

Due to the decreasing measurement sensitivity, the QR reduces in theory with the target depth increasing. However, several “exceptions” can be observed in Table 2, since the above principle makes sense merely for well-recovered target position, i.e., a reasonably small ε_z herein. In this case, a higher QR with none-zero ε_z should be worse depending on ε_z. Consequently, it might be more sensible to evaluate the reconstruction quality firstly in terms of ε_z rather than QR.

Since LOT is originally designed to detect blood hemoglobin in epithelial tissue, visible excitation is employed to achieve a high target-absorption-contrast (TAC), further reducing the maximum imaging depth (typical ~2 mm) [1,4]. To improve target depth localization, depth-compensation scheme can be incorporated into the reconstruction process through modifying the sensitivity matrix to optimally counterbalance the decay nature of light propagation in turbid medium [32].

3.3 Comparisons of noise robustness

The noise robustness of the four methods is evaluated using the same model with target depth and absorption contrast being 0.9mm and 5, respectively. To intuitively represent the influence of the measurement noise to the reconstructions, X-Z cross-sections at y = 4.165mm (crossing the target center) are presented in Fig. 4 instead of the averaged Y-Z ones, for the SNR of 10, 20 and 30dB, respectively. The last row in Fig. 4 shows the corresponding Z-profiles along the white dashed lines at y = 4.165 mm in the above images.

Fig. 4 Reconstructed slices along y = 4.165 mm and the corresponding Z-profiles along the white dashed lines, with the four schemes for a fixed target depth of 0.9 mm, when SNRs of the simulated data are (a) 10 dB, (b) 20 dB, and (c) 30 dB, respectively. The black dotted lines indicate the edges of original targets.

Download Full Size | PDF

As shown in Fig. 4, the recovered absorption distribution in the target region becomes unreasonable with the SNR deceasing due to the statistical instability. For noise-robust quantitation, the mean QR, FWHM_z and $ε_{z}$ along X-axis are further calculated for the Y-Z cross-sections (not shown), and averaged for 10 statistically independent reconstructions in Table 3. It can be found that, the QR has reduction rate of 26.38%, 28.62%, 27.67% and 27.67% for the linear, linear iterative, proposed and traditional methods, respectively, and correspondingly, the FWHM_z and $ε_{z}$ slightly increase for all the four methods, as the SNR goes down from 30 dB to 10 dB. Nevertheless, the target quantitativeness, size and location recovered with the proposed method are in good agreement with those recovered with the traditional one, indicating that both the algorithms have nearly the same noise robustness and reconstruction fidelity for SNR≥10 dB. Additional simulations to assessing the noise robustness have also been made for models embedded with 0.4 mm and 1.2 mm inclusions to support a SNR threshold of 10 dB.

Table 3. Comparisons among the reconstructions for noise robustness

View Table | View all tables in this article

For all the above simulative investigations, the iterative reconstructions (i.e., the linear iterative and the two nonlinear approaches) extremely outperforms the one-step linear one. In contrast to nonlinear reconstructions, the linear iterative method using a fixed searching direction to achieve economic computation with sacrificed accuracy. The proposed approach represents an extremely similar performance with the traditional nonlinear method, while considerably cuts down the computation time.

4. Experimental validations

In this work, experimental measurements are provided by a condensed dip-angle LOT (abbreviated as cdaLOT) system as illustrated in Fig. 5 [33]. The system uses a confocal-microscopy-based setup while remains the off-axis back-scattered light measurements to allow noncontact imaging with a nearly 200-μm spatial resolution over depth of ~2.5 mm. A focused laser beam is densely raster- or line- scanned over the sample surface, and the de-scanned lights emerging from the sample at an increasing distance from 0 to ~3 mm relative to the focus are detected by a linear-array photomultiplier (PMT) or CCD to form the raw measurements. Comparing to the traditional LOT setup, the cdaLOT substitutes the achromatic doublet or microscope objective lens with a bi-telecentric lens (MVTC23200, Thorlabs) to provide a nearly constant incident angle assumed by general forward model during raster-scanning. The polarized and collimated light source is provided by a 632.8 nm He-Ne laser. The SDSs of the eight channels are configured as same as those in simulative investigations. Note that the SDS increment is equal to the proportion between the channel interval of the linear-array PMT and the optical magnification of system [1,33]. For all the experiments below, the measurement geometry and source-detector configuration are the same as that in Fig. 1. For all the following reconstructions, the mesh- and iteration- regarding parameters are the same as those in the simulation situation.

Fig. 5 Schematic of the cdaLOT system.

Download Full Size | PDF

4.1 Phantom experiment

A semi-infinite phantom for LOT is fabricated from epoxy resin, with which titanium oxide particles (TiO₂) and dye (QS675B) are mixed as the scatter and absorber, respectively. After careful calibration, the optical properties at 632.8nm is set to μ_a0 = 0.01mm⁻¹, μ_s = 3mm⁻¹, and g = 0.9 to mimic those of human cervical tissues [34]. Vessel-mimic targets are parallelly embedded along X-axis at depths of 0.4mm and 0.9mm, with fixed cross-section area of 1 × 1mm² and fixed TAC of 5, i.e., $μ_{a} / μ_{a 0} = 5$ . It is noteworthy that, the center-to-center separation of the two targets is large enough to guarantee only one target is to be imaged in a single experiment. To simulate photon-migration in semi-infinite medium, computational domain adopted in reconstruction is extended along X-Y plane with a distance of $\sim 3 / {μ^{'}}_{t}$ where ${μ^{'}}_{t}$ is the reduced mean-free-path [35]. Figure 6(a) shows the reconstructed absorption images using the four schemes from left to right columns, respectively, with upper and lower rows showing the results for target depths of 0.4mm and 1mm, respectively. The corresponding Y- and Z- profiles along the white dashed lines in the images are plotted in the left and right columns, respectively, in Figs. 6(b) and 6(c).

Fig. 6 Phantom experimental results. (a) Reconstructed $μ_{a}$ images for target depths of 0.4mm (upper row) and 1mm (lower row), with the four schemes from left to right columns, respectively. (b)-(c) Corresponding profiles of the reconstructed $μ_{a}$ images along the white dashed lines marked in the above images for target depths of 0.4mm and 1mm, respectively.

Download Full Size | PDF

It can be seen from the upper row of Fig. 6(a) that, the target size and position for 0.4-mm scenario can be accurately reconstructed with the four schemes. Among them, the one-step linear scheme renders the lowest QR as expected. Deviations of the reconstructed target shape might be caused by the fabrication error of the phantom, since it is hard to guarantee the optical uniformity for such a slim heterogeneity in the epoxy-resin-based medium. With increase of the target depth, the measurement SNR reduces, as well as the measurement sensitivity. As shown in Fig. 6 and Table 4, the 1-mm reconstructions show inferior performance, e.g., more noise can be found in the target region and background, and the edge of the target reconstructed by the linear method is smoothed and diffused seriously in X-Y view. These are common problems in imaging technique based on reflecting measurement [4,31]. Nevertheless, the nonlinear methods outperform the linear (iterative) methods. For the above reconstruction errors, the biased estimations to the optical properties of the background and the instability of the system also contribute to them.

Table 4. Comparisons among the reconstructions for different target depths

View Table | View all tables in this article

4.2 In-vivo mouse experiment

It has been reported for tumor development that, in addition to oncogene-induced cell proliferation, another attractive candidate for a cellular change necessary is angiogenesis— the induction of neovascularization [36]. All entity tumors including pancreatic are highly vascularized, and the activation of neovascularization appears to be occurring prior to overt tumor formation in a subset of the preneoplastic nodules. Since the hemoglobin in blood has a strong absorption effect in visible wavelength, the tumor at an early stage might be detected via functional optical imaging technique with a mesoscopic-scale resolution, e.g., LOT.

A suspension of 5 × 10⁵ PanCO2 murine pancreatic tumor cells was subcutaneously injected into the abdominal region of a 3-week-old BALB/C mouse, and the in-vivo imaging with LOT was performed 2 weeks later. All animal procedures were reviewed and approved by the subcommittee on research animal care at Tianjin Medical University Cancer Institute & Hospital, where these experiments were performed.

The photograph of the tumor-bearing mouse is shown in Fig. 7(a), where the tumorous region in the murine abdomen is given a close-up description boxed with white dotted line in the right subplot. Superficial information can be found in the photo, such as skin texture, ecchymosis-suspected regions (region #1), and sub-surface capillaries (region #2). The subcutaneous haemorrhage might be caused by soft tissue strain. Figure 7(b) shows the averaged cdaLOT raw images over 100 frames for eight SDSs. Pixels of LOT raw image are composed of the raw measurements collecting from a fixed channel in PMT. Actually, the raw image is the integral of X-Y slices from z = 0 to a SDS-corresponding value. To objectively evaluate the reconstructions, anatomy experiment was conducted afterward near the tumor region. Anatomy experiments show that, the pancreatic tumor has a central depth of ~0.7mm and an area of ~4mm². As shown in Fig. 7(b), in addition to the superficial information, the curve-shaped pancreatic tumor can be found in those of large SDS, indicated by white arrowheads. With the increase of SDS, the shadow of the tumor (indicated by yellow arrowhead) can be seen and the separation between them increases correspondingly. This is caused by twice strong attenuations to the incident and emerged light, respectively [1,4]. The lateral direction between the twice mutations is orthogonal to that of scanning, and the lateral displacement is equal to SDS. Although these raw images are qualitatively useful, surface features as well as the slight lateral displacement of each raw LOT depth image can complicate interpretation of the true depth of the observed features [1]. Therefore they need to be further processed via the image reconstruction/deconvolution procedures.

Fig. 7 Imagings of the murine cancerous region: (a) photograph and (b) cdaLOT raw images.

Download Full Size | PDF

Given that the murine skin is two-layer structured, i.e., 100μm-thickness epidermis with an optical property set of {μ_a0 = 0.5mm⁻¹, μ_s0 = 52.26mm⁻¹, g = 0.9} and semi-infinite dermis with {μ_a0 = 0.5mm⁻¹, μ_s0 = 33.41mm⁻¹, g = 0.9} [34,37]. In the forward modeling, target surface to be imaged within such a limited FOV in LOT is herein treated to be flat. The reconstructed images and corresponding Z-profiles using the four schemes are presented in Fig. 8. In Fig. 8(a), the reconstructed X-Y slices at z = 0.2mm, 0.5mm, 0.71mm, and 0.9mm are shown from the first to fourth columns, respectively, and Z-Y slices across the tumor region are shown in the last column. Figure 8(b) shows the corresponding Z-profiles along the white dashed lines marked in the Y-Z images, and the red dashed line indicates z-position of the tumor absorption centroid. It is seen from Fig. 8(a) that, the tumors can be clearly found (as encircled by black dashed lines) in the X-Y and Z-Y slices at z = 0.71mm and x = 4.1mm with the linear iterative, proposed and traditional reconstructions. The tumor size and location recovered by the three iterative schemes are nearly the same and in agreement with the result of anatomy experiment. However, the tumor imaged with the one-step linear scheme has an extremely low absorption contrast and blurred edge. Similar with the previous experimental results, linear iterative one provides a lower QR comparing to the proposed and traditional approaches. To further flatten the target surface, the mice abdominal region is slightly pressed with a high-transparence coverslip as shown in Fig. 7(a), leading to a compressed tumor in Z-Y view.

Fig. 8 In-vivo experimental results. (a) Reconstructed images sliced at z = 0.2mm, z = 0.5mm, z = 0.71mm, z = 0.9mm, and x = 4.1mm, from left to right columns, respectively, with the four reconstruction schemes. (b) Corresponding Z-profiles along the white dashed lines marked in the above Y-Z slices. The red dashed line indicates z-position of the tumor centroid.

Download Full Size | PDF

5. Discussions and conclusions

In this paper, a practical nonlinear absorption reconstruction based on the recursive SVD inversion is proposed for LOT, where the recursions for calculating SVD components are derived based on an established scaled expression of sensitivity matrix. The computation complexity of solving the SVD-based linear inversion reduces from typical O(4NM² + 22N³) or O(4NM² + 8MN² + 9N³) to O(2NM² + MN²), with $M ≪ N$ being general occurrence in DOT/LOT. Time cost of iteratively updating $δ μ_{a}^{(k)}$ using the proposed (T_prop) and traditional schemes (T_trad) are compared in Fig. 9, for different size of sensitivity matrix J. On average, $T_{p r o p} / T_{t r a d}$ is merely 26.27%. The feasibility and availability of the proposed reconstruction strategy is borne out through numerical and experimental investigations with an established cdaLOT system. In contrast, linear iterative method shortens computation time at the cost of image quality, and might never be a reasonable choice in practice.

Fig. 9 Computation time per iteration VS. size of (J).

Download Full Size | PDF

To further enhance the proposed approach, rank degeneracy principle can be adopted to reduce the dimension of (singular) matrices in calculations [22, 23, 38]. In special, for the SVD-based linear inversion, this approach reduces to the truncated SVD (TVSD) regularization [22, 23, 38]. Furthermore, most prevalent accelerations, e.g., model-order reduction and parallel computing, etc., can be added. Last but not least, to accommodate the simultaneous perturbations of absorption and scattering, the scattering-dependent term needs to be added in Eq. (5) [17]. The voxelized scaling factor of averaged collision-number can be estimated as averaged photon-pathlength within the corresponding voxel divided by mean-free-path.

The traditional reconstruction can be accelerated with advanced SVD algorithms, e.g., the widely used randomized technique [39–41]. For the proof-of-concept survey, improved reconstructions using two randomized SVD methods, referred to as rSVD [40] and fSVD [41], respectively, are herein investigated to compare the proposed and traditional reconstructions, as shown in Fig. 10. The number of the retained singular values in rSVD and fSVD are set to 400 and 500, respectively, to keep almost the same computation time as the proposed method. It can be seen in Table 5 that, the traditional reconstructions with rSVD and fSVD have lower QR and larger $ε_{z}$ . Actually, for these two approaches, iterative modification almost fails due to absence of the significant components in singular values and vectors of the original sensitivity matrix, since the randomized technique is just to seek a nearly optimal low-rank approximation to achieve high computational efficiency. As a result, enhancement by such approximation-based methods depends on accuracy required, and might not be applicable for the mesoscopic reconstruction. In contrast, for the proposed methodology, the acceleration is based on the physically reasonable scaling principle of the photon-migration model, and therefore the derived recursive formulation is mathematically rigorous.

Fig. 10 (a) Reconstructed $μ_{a}$ images along y = 4.165mm with the proposed, traditional, and improved approaches using the rSVD and fSVD, respectively, for TAC = 2, SNR = 10dB, and a fixed target depth of 0.9 mm. (b) Corresponding Z-profiles along the white dashed lines marked in the above X-Z slices.

Download Full Size | PDF

Table 5. Comparisons among the reconstructions for different target depths. The numbers bracketed in Table 5 are the corresponding 1-iteration results.

View Table | View all tables in this article

Disclosures

The authors declare that there are no conflicts of interest related to this article.

Appendix A derivation of scaled sensitivity

Inserting Eq. (5) into the following limit form of the sensitivity definition

J_{m, n}^{(k)} = \lim_{δ μ_{a, n}^{(k)} \to 0} \frac{F_{m}^{(k + 1)} - F_{m}^{(k)}}{δ μ_{a, n}^{(k)}},

we obtain (N = 1)

J_{m, n}^{(k)} = \lim_{δ μ_{a, n}^{(k)} \to 0} \frac{F_{m}^{(k)} e^{- s_{m, n}^{(k)} δ μ_{a, n}^{(k)}} - F_{m}^{(k)}}{δ μ_{a, n}^{(k)}} = - F_{m}^{(k)} s_{m, n}^{(k)} .

Inserting Eq. (5) into the following derivative form of the sensitivity definition

J_{m, n}^{(k + 1)} = {\frac{\partial F_{m}^{(k + 1)}}{\partial μ_{a, n}} |}_{μ_{a, n}^{(k)}},

we have

J_{m, n}^{(k + 1)} = - s_{m, n}^{(k)} F_{m}^{(k)} e^{- s_{m, n} δ μ_{a, n}^{(k)}}

Combining Eq. (A2) and Eq. (A4) leads to

J_{m, n}^{(k + 1)} = J_{m, n}^{(k)} e^{- s_{m, n}^{(k)} δ μ_{a, n}^{(k)}} .

Funding

National Natural Science Foundation of China (81271618, 81371602, 61475115, 61475116, 81401453, 61575140, 81571723, 81671728), the 111 Project (B07014), and Tianjin Municipal Government of China (14JCQNJC14400, 15JCZDJC31800, 15JCQNJC14500, 16JCZDJC31200, 17JCZDJC32700).

References and links

1. E. M. Hillman, D. A. Boas, A. M. Dale, and A. K. Dunn, “Laminar optical tomography: demonstration of millimeter-scale depth-resolved imaging in turbid media,” Opt. Lett. 29(14), 1650–1652 (2004). [CrossRef] [PubMed]

2. M. Schwarz, A. Buehler, J. Aguirre, and V. Ntziachristos, “Three-dimensional multispectral optoacoustic mesoscopy reveals melanin and blood oxygenation in human skin in vivo,” J. Biophotonics 9(1-2), 55–60 (2016). [CrossRef] [PubMed]

3. A. Taruttis and V. Ntziachristos, “Advances in real-time multispectral optoacoustic imaging and its applications,” Nat. Photonics 9(4), 219–227 (2015). [CrossRef]

4. D. A. Boas, C. Pitris, and N. Ramanujam, Handbook of Biomedical Optics (Taylor and Francis, 2010).

5. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25(12), 123010 (2009). [CrossRef]

6. A. Dunn and D. Boas, “Transport-based image reconstruction in turbid media with small source-detector separations,” Opt. Lett. 25(24), 1777–1779 (2000). [CrossRef] [PubMed]

7. I. Seo, C. K. Hayakawa, and V. Venugopalan, “Radiative transport in the delta-P1 approximation for semi-infinite turbid media,” Med. Phys. 35(2), 681–693 (2008). [CrossRef] [PubMed]

8. L. Yao, Y. Sun, and H. Jiang, “Quantitative photoacoustic tomography based on the radiative transfer equation,” Opt. Lett. 34(12), 1765–1767 (2009). [CrossRef] [PubMed]

9. E. M. Hillman, H. Dehghani, J. C. Hebden, S. R. Arridge, M. Schweiger, and D. T. Delpy, “Differential imaging in heterogeneous media: limitations of linearization assumptions in optical tomography,” Proc. SPIE 4250, 327–338 (2001). [CrossRef]

10. A. P. Gibson, J. C. Hebden, J. Riley, N. Everdell, M. Schweiger, S. R. Arridge, and D. T. Delpy, “Linear and nonlinear reconstruction for optical tomography of phantoms with nonscattering regions,” Appl. Opt. 44(19), 3925–3936 (2005). [CrossRef] [PubMed]

11. S. Gupta, P. K. Yalavarthy, D. Roy, D. Piao, and R. M. Vasu, “Singular value decomposition based computationally efficient algorithm for rapid dynamic near-infrared diffuse optical tomography,” Med. Phys. 36(12), 5559–5567 (2009). [CrossRef] [PubMed]

12. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef] [PubMed]

13. C. Gong, J. Liu, L. Chi, H. Huang, J. Fang, and Z. Gong, “GPU accelerated simulations of 3D deterministic particle transport using discrete ordinates method,” J. Comput. Phys. 230(15), 6010–6022 (2011). [CrossRef]

14. X. Yi, X. Wang, W. Chen, W. Wan, H. Zhao, and F. Gao, “Full domain-decomposition scheme for diffuse optical tomography of large-sized tissues with a combined CPU and GPU parallelization,” Appl. Opt. 53(13), 2754–2765 (2014). [CrossRef] [PubMed]

15. Q. Liu and N. Ramanujam, “Scaling method for fast Monte Carlo simulation of diffuse reflectance spectra from multilayered turbid media,” J. Opt. Soc. Am. A 24(4), 1011–1025 (2007). [CrossRef] [PubMed]

16. C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26(17), 1335–1337 (2001). [CrossRef] [PubMed]

17. M. Jia, S. Cui, X. Chen, M. Liu, X. Zhou, H. Zhao, and F. Gao, “Image reconstruction method for laminar optical tomography with only a single Monte-Carlo simulation,” Chin. Opt. Lett. 12(3), 031702 (2014). [CrossRef]

18. A. H. Hielscher, A. D. Klose, and K. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imaging 18(3), 262–271 (1999). [CrossRef] [PubMed]

19. S. Arridge and M. Schweiger, “A gradient-based optimisation scheme foroptical tomography,” Opt. Express 2(6), 213–226 (1998). [CrossRef] [PubMed]

20. X. Gu, Y. Xu, and H. Jiang, “Mesh-based enhancement schemes in diffuse optical tomography,” Med. Phys. 30(5), 861–869 (2003). [CrossRef] [PubMed]

21. A. Joshi, W. Bangerth, and E. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004). [CrossRef] [PubMed]

22. M. E. Eames, B. W. Pogue, P. K. Yalavarthy, and H. Dehghani, “An efficient Jacobian reduction method for diffuse optical image reconstruction,” Opt. Express 15(24), 15908–15919 (2007). [CrossRef] [PubMed]

23. X. Wu, A. T. Eggebrecht, S. L. Ferradal, J. P. Culver, and H. Dehghani, “Fast and efficient image reconstruction for high density diffuse optical imaging of the human brain,” Biomed. Opt. Express 6(11), 4567–4584 (2015). [CrossRef] [PubMed]

24. A. Quarteroni and A. Valli, Domain Decomposition Methods for Partial Differential Equations (Oxford Science, 1999).

25. G. H. Golub and C. F. Van Loan, Matrix Computation (Johns Hopkins Press, Baltimore, 1983).

26. L. Trefethen and D. Bau, Numerical Linear Algebra (Siam, 1997).

27. T. Papadopoulo and M.I. A. Lourakis, “Estimating the Jacobian of the Singular Value Decomposition: Theory and Applications,” Research Report 3961, INRIA Sophia-Antipolis, June 2000.

28. G. W. Stewart and J.-G. Sun, Matrix Perturbation Theory (Academic, 1990).

29. Y. Yamada and S. Okawa, “Diffuse optical tomography: present status and its future,” Opt. Rev. 21(3), 185–205 (2014). [CrossRef]

30. B. Yuan, S. A. Burgess, A. Iranmahboob, M. B. Bouchard, N. Lehrer, C. Bordier, and E. M. C. Hillman, “A system for high-resolution depth-resolved optical imaging of fluorescence and absorption contrast,” Rev. Sci. Instrum. 80(4), 043706 (2009). [CrossRef] [PubMed]

31. M. Jia, H. Zhao, J. Li, L. Liu, L. Zhang, J. Jiang, and F. Gao, “Coupling between radiative transport and diffusion approximation for enhanced near-field photon-migration modeling based on transient photon kinetics,” J. Biomed. Opt. 21(5), 050501 (2016). [CrossRef] [PubMed]

32. H. Niu, Z. J. Lin, F. Tian, S. Dhamne, and H. Liu, “Comprehensive investigation of three-dimensional diffuse optical tomography with depth compensation algorithm,” J. Biomed. Opt. 15(4), 046005 (2010). [CrossRef] [PubMed]

33. H. Zhao, S. Wang, M. Jia, X. Chen, J. Qi, J. Tian, W. Ma, J. Li, and F. Gao, “A modified laminar optical tomography system with small dip-angle and the initial validation,” Proc. SPIE 9700, 97001B (2016). [CrossRef]

34. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef] [PubMed]

35. H. Fujii, S. Okawa, Y. Yamada, and Y. Hoshi, “Hybrid model of light propagation in random media based on the time-dependent radiative transfer and diffusion equations,” J. Quant. Spectrosc. Radiat. Transf. 147, 145–154 (2014). [CrossRef]

36. J. Kandel, E. Bossy-Wetzel, F. Radvanyi, M. Klagsbrun, J. Folkman, and D. Hanahan, “Neovascularization is associated with a switch to the export of bFGF in the multistep development of fibrosarcoma,” Cell 66(6), 1095–1104 (1991). [CrossRef] [PubMed]

37. T. J. Muldoon, S. A. Burgess, B. R. Chen, D. Ratner, and E. M. Hillman, “Analysis of skin lesions using laminar optical tomography,” Biomed. Opt. Express 3(7), 1701–1712 (2012). [CrossRef] [PubMed]

38. M. E. Wall, A. Rechtsteinner, and L. M. Rocha, “Singular value decomposition and principal component analysis,” in A Practical Approach to Microarray Data Analysis, D. P. Berrar, W. Dubitzky, and M. Granzow ed. (Kluwer, Norwell, 2003).

39. C. Musco and C. Musco, “Randomized block krylov methods for stronger and faster approximate singular value decomposition,” in Proceedings of Neural Information Processing Systems 28 (NIPS, 2015), pp. 1396–1404.

40. A. Liutkus, “Randomized SVD,” MATLAB Central File Exchange, 2014.

41. V. Vijayan, “Fast SVD and PCA,” MATLAB Central File Exchange.

Method	Time (s)	TAC = 2			TAC = 3			TAC = 5
Method	Time (s)	QR (%)	FWHM_z (mm)	ε_z (%)	QR (%)	FWHM_z (mm)	ε_z (%)	QR (%)	FWHM_z (mm)	ε_z (%)
Linear	\	54.54	1.78	19.90	43.66	1.59	19.90	36.22	1.49	19.90
Linear iterative	38.4	86.60	1.50	4.39	69.72	1.38	4.39	58.49	1.32	4.39
Proposed	84.8	104.42	1.24	1.81	84.07	1.17	4.39	70.48	1.13	4.39
Traditional	317.3	104.41	1.24	1.81	84.05	1.17	4.39	70.46	1.13	4.39

Method	Depth = 0.4mm			Depth = 0.9mm			Depth = 1.2mm
Method	QR (%)	FWHM_z (mm)	ε_z (%)	QR (%)	FWHM_z (mm)	ε_z (%)	QR (%)	FWHM_z (mm)	ε_z (%)
Linear	39.89	1.18	3.10	36.22	1.49	19.90	33.70	1.54	14.97
Linear iterative	46.49	1.15	3.10	58.49	1.32	4.39	59.86	1.26	8.43
Proposed	50.74	1.18	0.78	70.48	1.13	4.39	66.39	1.18	8.43
Traditional	50.74	1.18	0.78	70.46	1.13	4.39	66.38	1.18	8.43

Method	SNR = 10dB			SNR = 20dB			SNR = 30dB
Method	QR (%)	FWHM_z (mm)	ε_z (%)	QR (%)	FWHM_z (mm)	ε_z (%)	QR (%)	FWHM_z (mm)	ε_z (%)
Linear	25.73	1.54	14.05	32.88	1.48	13.11	36.17	1.46	12.14
Linear iterative	38.47	1.34	3.40	49.29	1.32	1.81	53.90	1.22	1.81
Proposed	52.33	1.14	1.81	65.05	1.13	1.81	72.35	1.06	1.81
Traditional	52.32	1.14	1.81	65.03	1.13	1.81	72.34	1.06	1.81

Method	Depth = 0.4mm				Depth = 1mm
Method	QR (%)	FWHM_y (mm)	FWHM_z (mm)	ε_z (%)	QR (%)	FWHM_y (mm)	FWHM_z (mm)	ε_z (%)
Linear	13.29	1.54	1.45	9.30	24.14	1.50	1.65	16.28
Linear iterative	58.00	1.11	1.24	9.30	46.69	1.33	1.19	11.63
Proposed	60.15	1.08	1.29	2.33	59.68	1.33	1.07	6.97
Traditional	60.16	1.08	1.29	2.33	57.89	1.33	1.09	6.98

Method	QR (%)	FWHM_y (mm)	FWHM_z (mm)	ε_z (%)
Trad + rSVD	33.78 (31.58)	1.16 (1.19)	1.24 (1.24)	27.65 (27.65)
Trad + fSVD	32.45 (29.90)	1.16 (1.20)	1.46 (1.42)	12.14 (17.31)
Proposed	52.18 (30.83)	1.02 (1.15)	1.10 (1.40)	6.98 (17.31)
Traditional	52.18 (30.83)	1.02 (1.15)	1.10 (1.40)	6.98 (17.31)

Accelerating nonlinear reconstruction in laminar optical tomography by use of recursive SVD inversion

Abstract

1. Introduction

2. Methods

2.1 SVD-based framework for nonlinear inversion

2.2 Recursive SVD scheme for sensitivity matrix

3. Simulative investigations

3.1 Comparisons for model with different target-absorption-contrast

3.2 Comparisons for model with different target depth

3.3 Comparisons of noise robustness

4. Experimental validations

4.1 Phantom experiment

4.2 In-vivo mouse experiment

5. Discussions and conclusions

Disclosures

Appendix A derivation of scaled sensitivity

Funding

References and links

Cited By

Figures (10)

Tables (5)

Equations (17)

Biomedical Optics Express