## Abstract

Mesoscopic fluorescence molecular tomography (MFMT) is a novel imaging technique that aims at obtaining the 3-D distribution of molecular probes inside biological tissues at depths of a few millimeters. To achieve high resolution, around 100-150μm scale in turbid samples, dense spatial sampling strategies are required. However, a large number of optodes leads to sizable forward and inverse problems that can be challenging to compute efficiently. In this work, we propose a two-step data reduction strategy to accelerate the inverse problem and improve robustness. First, data selection is performed via signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR) criteria. Then principal component analysis (PCA) is applied to further reduce the size of the sensitivity matrix. We perform numerical simulations and phantom experiments to validate the effectiveness of the proposed strategy. In both *in silico* and *in vitro* cases, we are able to significantly improve the quality of MFMT reconstructions while reducing the computation times by close to a factor of two.

© 2017 Optical Society of America

## 1. Introduction

Mesoscopic Fluorescence Molecular Tomography (MFMT), also known as Fluorescence Laminar Optical Tomography (FLOT), is an emerging optical imaging modality operating between conventional microscopic and macroscopic imaging scales [1–3] for investigations in scattering samples. By collecting fluorescence signals exiting at multiple distances away from the illumination spot, MFMT can discriminate the depth location of the biomarkers of interest [4-5]. Thanks to its higher resolution performance compared to Diffuse Optical Tomography (DOT), MFMT has shown its utility for *in vivo* imaging applications [6–8] as well as *in vitro* applications such as evaluating engineered tissues [9]. However, MFMT is an inverse problem-based modality in which imaging performance is highly dependent on the accuracy of the forward model, the spatial sampling acquired, and the inverse solvers selected. Resolution of the technique can be improved by using advanced solver such as sparsity-enhancing methodologies [10], but it is still critically dependent on the number of optodes (source-detector pairs) acquired. Typically, MFMT is achieved via raster scanning of a laser beam over the sample and detection of light on discrete detectors (PMT or APD line array) in a de-scanned configuration [1,5,11]. Recently, we have developed a new system in which the spatial sampling is greatly increased by using an Electron Multiplying Charged Coupled Device (EMCCD) camera [12]. Such a system potentially leads to measurements in upwards of 10^{8-9} source-detector pairs. This amount of data cannot be directly employed in the inverse problem due to the sheer size of the computations involved. Hence, there is the need to develop post-processing methodologies that identify the most useful subset of measurements for accurate and computationally efficient reconstructions.

In this work, we propose a hybrid approach to reduce redundant and noisy detector readings and minimize their effect on the reconstruction results. The computation time to solve the inverse problem is reduced accordingly, without compromising the resolution and robustness of MFMT. First we investigate the effect of noisy data on reconstruction performance. Thresholds for both Signal-to-Noise Ratio (SNR) and Contrast-to-Noise Ratio (CNR) are established and detectors below the thresholds are discarded from further use. Second we explore Principal Component Analysis (PCA) to further extract the most useful and independent components of the sensitivity matrix as well as the measurements. In addition, we exploit the sparsity of fluorophores and apply a ${L}_{1}$ -norm regularization algorithm to solve the ill-posed inverse problem after the data reduction procedures.

In Section 2, we first give a brief introduction of our second generation MFMT system and relevant techniques, *i.e.*, building the forward model and solving the inverse problem. The two-step data reduction strategy and the metrics for evaluating reconstruction results are then described in detail in Section 3. To ascertain the performance of the proposed method, we carry out both numerical simulations and phantom experiments and compare reconstruction results before and after data reduction in Section 4. Finally, we discuss the relevance of the results and potential directions for future work.

## 2. High density MFMT

We have demonstrated the utility of MFMT *in vivo* and for tissue engineering applications, but our first generation imaging system was limited to an 8 detector array in a de-scanned acquisition mode [5,8–10]. Recently, we have upgraded the detection channel of this system by incorporating an EMCCD camera that can yield up to 512 × 512 measurements in parallel. In this section, we briefly report on the design and main features of our MFMT imaging platform as well as on the inverse problem formulation.

#### 2.1 Data acquisition

The optical diagram of our 2nd generation MFMT system is shown in Fig. 1. Briefly, the optical path starts with introduction of an excitation light (Ex) into the system through a linear polarizer (P). A Polarizing Beam Splitter (PBS) reflects ~90% of the S-polarized light onto a Galvanometer Mirror pair (GM). The GM controls the scanning area and dwell time for each excitation point through a Scan Lens (SL). The backscattered light is collected by the same SL and filtered by the PBS. The PBS allows ~90% transmission of P-polarized emission and minimizes the specular reflection of S-polarized excitation. After the PBS, the backscattered light is further filtered by another polarizer (A – to additionally reduce specular reflection) and then spectrally filtered using the appropriate interference filter (F). As this de-scanned configuration collects light exiting the tissue, 0.6-1.8mm away from the illumination spot, the signals acquired exhibit a large dynamical range. To mitigate this drawback, we introduce a custom-made reflection block (RB) right before relaying the signal onto the camera, effectively blocking the light originating at the same location as the illumination spot. RB ensures an adequate dynamical range and SNR for distal detectors. The light is then collected by an EMCCD (iXon^{EM+} DU-897 back illuminated, Andor Technology).

The relative positions of source and detector locations as imaged on the specimen by our system are shown in Fig. 2, where the blue square represents the source, green squares represent the detectors, and the blue dots represent the scanning positions on the surface of the specimen. This depiction renders an acquisition mode with 2 × 2 binning, leading to a total of 256 × 256 detectors covering a FOV of ~6.4 × 6.4 mm^{2}. Of those detectors, we selected an appropriate source-detector separation (~0.6 mm) for 49 detectors in a square grid formation (central detector is occluded due to RB). The source and detectors are moving together as indicated by the blue arrow and the source is always in the center of the detector array. The software-controlled raster scanning step size and dwell times are typically set to 100μm and 20ms, respectively. The typical FOV for illumination is 3.1 × 3.1 mm^{2}, leading to a total of 961 scanning points [12].

#### 2.2 The forward model

An accurate forward model, describing the photon propagation inside biological tissues and matching experimental configurations with collected measurements, is critical for every inverse-based tomography technique [13]. Since MFMT operates in the mesoscopic regime, the Monte Carlo (MC) method is used to build the forward model. More specifically, we employ a graphics processing unit (GPU)-based MC code to calculate continuous-wave Green’s functions at excitation and emission wavelengths ${G}^{x}$ and ${G}^{m}$ [14], and compute the weight function *W* with an adjoint formulation for efficiency [15]:

The 3-D distribution of the fluorophore’s effective quantum yield $\eta \left(r\right)$ can be revealed by solving the integral equation below:

#### 2.3 The inverse problem

MFMT is by nature an ill-posed inverse problem due to the diffuse nature of the light collected and the limited projection angle offered by the reflectance geometry. As a result, regularization methods have to be applied to obtain robust and accurate reconstructions. Herein, instead of the commonly used ${L}_{2}$-norm regularization, also known as Tikhonov regularization, we take advantage of the sparsity of fluorophore distribution and adopt an ${L}_{1}$-norm regularization scheme [17, 18], expressed as the following optimization problem:

where $\Vert {x}_{1}\Vert $ is the ${L}_{1}$ norm of the reconstructed fluorophore distribution vector, and $\lambda $ is the regularization parameter, usually determined through an $L$-curve. The minimization problem is convex and then solved with an iterative shrinkage-thresholding algorithm [19].## 3. Data reduction strategies for MFMT

The use of an EMCCD camera enables collecting a massive amount of spatial measurements within a short period, leading to potentially increased quality of MFMT reconstruction. However, such a data set is characterized by high redundancy due to scattering and potentially high variance in SNR among different detectors due to large dynamical range. Hence, it is expected that data reduction methodologies aiming at extracting the most informative and robust measurements would lead to better reconstruction accuracy. Herein, we investigate a post-processing methodology that first selects a sub-set of data based on quality measurement criteria and second, uses PCA for further data reduction.

#### 3.1 Data reduction based on SNR/CNR

First, we define the Signal-to-Noise Ratio (SNR) and Contrast-to-Noise Ratio (CNR) of one specific detector as:

Detectors with either low average SNR or low CNR readings are then discarded, and the corresponding rows in the sensitivity matrix are truncated as well. As a result, the overall quality of the measurements is improved and the size of the inverse problem is also reduced.

#### 3.2 Data reduction based on PCA

After removing measurements deemed poor based on $\text{SNR}$ and $\text{CNR}$ criteria, we address the data redundancy by using Principal Component Analysis (PCA). PCA has been proved to be an effective method of extracting useful information from large data sets in numerous fields, including Fluorescence Molecular Tomography (FMT) [20–22]. To filter out correlated measurements and further reduce the size of the inverse problem, we apply PCA on MFMT data sets following the methodology described below.

First, the covariance matrix of the original sensitivity matrix is calculated as $C=A{A}^{\text{T}}$, followed by an eigenvalue decomposition of the covariance matrix:

where**P**is the matrix of columns of eigenvectors of

**C**and

**Λ**is a diagonal matrix with the elements of eigenvalues ${\lambda}_{1}\ge {\lambda}_{2}\ge {\lambda}_{3}\ge \cdots $ of the covariance matrix. The inverse problem we want to solve changes from (3) to the following form:If we only retain the first $k$ principal components and leave out the less significant ones, the dimensions of the sensitivity matrix as well as the measurement vector are reduced and (8) becomes:where ${P}_{k}^{\text{T}}$ is a sub-matrix of ${P}^{\text{T}}$, containing only the first $k$ rows of ${P}^{\text{T}}$,

*i.e.*, the largest $k$ eigenvectors (out of $m$ in total) of the covariance matrix. The number of preserved principal components is typically determined by the cumulative percent of total variance ($\text{CPV}$), defined as:

#### 3.3 Evaluation of MFMT reconstruction results

Four metrics are adopted to quantify the difference between reconstruction results and the ground truth, namely, normalized sum squared difference ($\text{nSSD}$), normalized sum absolute difference ($\text{nSAD}$), normalized disparity ($nD$), and correlation ($R$) [23].

Suppose we have $N=X\times Y\times Z$ voxels in a cubic phantom, where $X,Y,Z$ are the number of voxels along the $x,y$ and $z$ axes. $A(i,j,k)$, $B(i,j,k)$ are the normalized numerical model and reconstructed results with values between 0 and 1 after all negative reconstructed values are set to 0 where $\text{nSSD}$ and $\text{nSAD}$ are defined as:

## 4. *In silico* and experimental validation

#### 4.1 Simulation settings

A numerical phantom was designed to evaluate the performance of the proposed algorithm and displayed in Fig. 3. The imaging domain had a surface area of 3.1 × 3.1mm^{2} with a depth of 3mm. This domain was discretized into 31 × 31 × 30 cubic voxels, leading to 100μm^{3} element of volumes. A vascular tree with main trunk and three groups of offshoots was simulated deep within the phantom (*z* = 1mm). The diameter of the trunk was set to 400μm whereas it was set to 200μm for the offshoots. The separation between two adjacent off shoots was set at 100μm (discretization level). The fluorophore concentration was assumed to be uniform in the vessel with effective quantum yield equal to 1 (for simplicity). The whole domain was assumed to be homogeneous at the excitation wavelength with absorption coefficient *μ _{a}* = 0.002mm

^{−1}, scattering coefficient

*μ*= 1 mm

_{s}'^{−1}, index of refraction

*n*= 1.34, and anisotropy factor

*g*= 0.81. These optical properties are derived from the collagen scaffold typically employed in our bio-printing application at collagen density of 6-9 mg/ml [5,9].

To generate *in silico* measurements, we replicated the imaging system configuration as used in real experiments: 961 illumination scanning positions and 48 detectors at each source location. This configuration yields a Jacobian with size of 46,128 by 28,830 elements. For computational ease and since the study is comparative by nature, we generated the measurements by multiplying the sensitivity matrix with the bio-printed vasculature model. Though, this approach leads to measurements with excellent SNR whereas experimental data can yield SNRs as low as 2. Hence, we added a spatially variant Gaussian noise that replicates the SNR levels seen experimentally. In simulation, we generated the noisy measurements by adding the pure measurements with the noise produced from a normal distribution with mean 0 and standard deviation determined by SNR. For the synthetic measurements we also investigate the dynamic range of the specific noise level in terms of SNR and CNR (see Fig. 3 (d) and (e)). Note that we are not able to calculate the SNR value of one source-detector pair for real experimental data, so the SNR value of one detector (961 source positions) is calculated as in Eq. (5) and used for filtering instead. Then reconstructions are performed as described above. For each of the following simulations or experiment cases, the optimal regularization parameter for the ${L}_{1}$-norm reconstruction algorithm is chosen from an $L$-curve analysis as well as visual inspection.

All computations carried to generate the Jacobians and synthetic measurements were performed under less than 5 minutes in a desktop computer (CPU Intel® Core i7-3820 Quad-Core 3.60 GHz 10MB Intel Smart Cache LGA2011, GPU Quadro K620). Although the reconstructions were performed on the same computer mentioned above, computational time is highly depends on the scale of measurements, Jacobians, and the algorithms adopted as exemplified in the following sections.

#### 4.2 In silico results under different data reduction methods

### 4.2.1 Data reduction via random sampling

To objectively assess the efficiency of the proposed data reduction, we first established the performance of data reduction via random depopulation of the sensitivity matrix. To ensure that the image-derived metrics are not affected by the stochastic nature of the random approach, ten trials were performed for each set level of data reduction reported. A summary of the evaluation metrics as well as computational times versus level of data reduction is provided in Fig. 4(a). Figures 4(b) and 4(c) also provide a visual validation for the effect of random sampling for two specific levels of data reduction.

Overall, although random sampling reduces the computation time, there are no obvious improvements on evaluation metrics when a certain data reduction level is implemented. Especially, as the number of measurements used in reconstruction is reduced, the offshoots cannot be resolved and the integrity of the main trunk is compromised (Fig. 4. (b)).

### 4.2.2 Data reduction based on SNR

MFMT is an ill-conditioned inverse problem and hence, the quality of the reconstructions is very sensitive to noise. Such sensitivity is typically observed as large artefacts on the surface and/or degradation of the resolution, leading to reconstructions with extremely poor fidelity. Hence, beyond implementing regularization techniques, one can implement a pre-processing data selection methodology that aims at discarding any detector readings with low SNR before proceeding to the reconstruction. An example of SNR distribution of each of the 48 detectors over 961 sources is provided in Fig. 5(a). From this distribution, a low SNR threshold can be defined and implemented such that any detector reading below this threshold is not considered in the inverse problem. For different values of this threshold, the number of measurements left and computation times are shown in Fig. 5(b) while the evaluation metrics of corresponding reconstruction results are shown in Fig. 5(c). The reconstructions associated with three specific threshold values are provided in Figs. 5(d-f). From these simulation results, we can see that removing detectors with lower SNR improves the quality of reconstruction and outperforms the random data reduction approach. However, if the threshold value is set too high, the reconstruction quality is compromised, as visually illustrated in Figs. 5(d-f). Hence, an optimal SNR threshold has to be determined based on the SNR level of the acquired data set as well as the amount of data points available.

### 4.2.3 Data reduction based on CNR

It is imperative not to cull the data too drastically based on SNR. It remains that some detector readings, while having good SNR, contribute little to revealing the fluorophore distribution due to the lack of image contrast. Therefore, they can also be considered as redundant data and could be removed. Similar to the SNR-based case, we calculated the CNR values of 48 detectors from all 961 source positions, and plot the CNR distribution, numbers of measurements left and time costs in Figs. 6(a-b). The results under different CNR levels, as shown in Figs. 6(c-f), also follow the same trend as for the SNR-based case: the reconstructions improve as detectors with low CNR are discarded, but then worsen with increasing the CNR cutoff value.

### 4.2.4 Data reduction via PCA

To investigate the performance of PCA-derived data reduction on MFMT reconstructions, we retained the largest principal components *k*, from 5,000 to 40,000 out of 46,128. The cumulative percent of total variance is plotted in Fig. 7(a), while the computation times and values of evaluation metrics are shown in Fig. 7(b). As shown, the reconstruction quality improves significantly when *k* rises from 5,000 to 10,000, but then remains relatively stable when *k* varies from 10,000 to 35,000, and then degrades as *k* continues to increase to 40,000 until all components are used. These results are in accordance with the expectation that the majority of the information in the Jacobian can be represented with a much smaller dimension – in our case >80% of CPV explained with <33% of the largest components. Also because the measurements are noisy and highly redundant, introducing more principal components is not expected to bring additional useful information, but only amplify the noise effect instead, as visually validated in Figs. 7(c-g).

To further demonstrate the benefit of PCA, we provide a detailed comparison in Table 1 between when no PCA is applied and when only the 35,000 largest components are retained. Despite the decrease in Jacobian size and time cost when applying PCA, all four evaluation metrics are higher compared to the original result. PCA is thus proven to be effective in both accelerating the inverse problem and improving the quality of reconstruction.

#### 4.3 Summary of quantification results for data reduction methods on simulation data

Here we give a table to shortly summarize the quantification results for the above mentioned simulations. As shown in the Table 2, although reconstruction time is basically affected by the size of sensitivity matrix, it highly depends on the property of sensitivity matrix, such as sparse, correlation of components. In addition, the preset thresholds of SNR/CNR or CPV have a huge impact on the size of sensitivity matrix and reconstruction quality, as shown in Fig. 4 to Fig. 7. Note that we set different thresholds of SNR/CNR or CPV in order to achieve high quality reconstruction in different cases.

#### 4.4 Performance of data reduction on experimental data

Last we applied the proposed two-step data reduction algorithm on an experimental data set, acquired from a collagen sample with optical properties *μ _{a}* = 0.002 mm

^{−1},

*μ*= 1 mm

_{s}'^{−1},

*n*= 1.34, and

*g*= 0.81. The sample size was 3.1 × 3.1 × 3.0mm

^{3}, with four polystyrene fluorophore beads (GFP 488/509, Cospheric) placed 1.7~1.8$\text{mm}$ beneath the sample surface following standard protocols. The voxel size of 100 × 100 × 100μm

^{3}is also employed for the experimental phantom. Transversal micro-MRI slices taken across the x-y and y-z planes are shown in Figs. 8(a-b). As in the simulation configuration, the source scans over 31 × 31 positions and measurements of 48 detectors are collected at each source position in less than 20ms.

The reconstruction with the original sensitivity matrix and measurements was initially performed. Then, for the first step of data reduction, we discarded detector readings with low SNR and low CNR sequentially and carried out a second reconstruction based on the simulation results. According to the simulation results the thresholds of SNR and CNR were optimally set to 2 and 6.5, leading to the removal of 6 detectors (5,766 measurements) and 3 detectors (2,883 measurements), respectively (see Fig. 8 (g-h)). Lastly, for the second step of data reduction, we applied PCA to the trimmed Jacobian from the last step. The CPV value was set to 95%, corresponding to 28,123 principal components, and a third reconstruction with the data-reduced Jacobian and measurement vector was conducted. All three inverse problems were solved with our *L*_{1} reconstruction algorithm as described in section 2.3 and the optimal regularization parameters were determined through L-curves.

The 3D reconstruction results for the above three cases are shown in Table 3 and Figs. 8(d-f). Compared to the original case, the artifacts are remarkably suppressed after data points with low SNR/CNR measurements are removed. However the accurate locations and shapes of four beads are still not satisfactory compared to micro-MRI. After applying the PCA data reduction approach, the artifacts are completely eliminated while the locations and depths of beads match perfectly with those in the micro-MRI images, shown in Figs. 8(a-c). Hence, these results demonstrate that the sequential combination of SNR/CNR with PCA data reduction approaches leads to faster as well as improved reconstruction performances.

## 5. Conclusions

In conclusion, we have presented a comprehensive data reduction method on an MFMT data set, where measurements with low SNR/CNR are first removed and PCA is then applied for dimension reduction. With the proposed method, the noise and redundancy in MFMT raw data are minimized while the most useful information revealing fluorophore distribution is retained. We have tested the performance of the data reduction algorithm on data sets from numerical simulations as well as real experiments. In both cases, the performances of the reconstruction were significantly improved in terms of spatial accuracy as well as computational efficiency. In the case of the experimental study, locations and shape of the imaged objects were retrieved with high fidelity (resolution <200 µm) even though they were deeply embedded (*z* = 1.8 mm) in a highly scattering medium (*μ _{s}'* = 1 mm

^{−1}). We expect that this methodology will support our efforts in monitoring longitudinally bio-printed tissues [12] and

*in vivo*tumor xenografts [8]. Additionally, beyond MFMT and/or CW data sets, we will investigate the utility of this data processing methodology for applications relying on temporal and/or spectral data sets such as those acquired in our novel hyperspectral time-resolved imager [24, 25].

## Funding

National Institutes of Health (NIH) (R01-EB019443 and R01 BRG-CA207725); National Science Foundation (NSF) (Career Award CBET 1149407); Natural Science Foundation of Shandong Province (ZR2014FL031); National Natural Science Foundation of China (61472227); Science and Technology Project of Yantai (2015ZH059).

## 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. **S. Yuan, Q. Li, J. Jiang, A. Cable, and Y. Chen, “Three-dimensional coregistered optical coherence tomography and line-scanning fluorescence laminar optical tomography,” Opt. Lett. **34**(11), 1615–1617 (2009). [CrossRef] [PubMed]

**3. **M. S. Ozturk, C.-W. Chen, R. Ji, L. Zhao, B.-N. B. Nguyen, J. P. Fisher, Y. Chen, and X. Intes, “Mesoscopic Fluorescence Molecular Tomography for Evaluating Engineered Tissues,” Ann. Biomed. Eng. **44**(3), 667–679 (2016). [CrossRef] [PubMed]

**4. **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]

**5. **L. Zhao, V. K. Lee, S.-S. Yoo, G. Dai, and X. Intes, “The integration of 3-D cell printing and mesoscopic fluorescence molecular tomography of vascular constructs within thick hydrogel scaffolds,” Biomaterials **33**(21), 5325–5332 (2012). [CrossRef] [PubMed]

**6. **S. Björn, K. H. Englmeier, V. Ntziachristos, and R. Schulz, “Reconstruction of fluorescence distribution hidden in biological tissue using mesoscopic epifluorescence tomography,” J. Biomed. Opt. **16**(4), 046005 (2011). [CrossRef] [PubMed]

**7. **Q. Tang, V. Tsytsarev, A. Frank, Y. Wu, C. W. Chen, R. S. Erzurumlu, and Y. Chen, “In Vivo Mesoscopic Voltage-Sensitive Dye Imaging of Brain Activation,” Sci. Rep. **6**(1), 25269 (2016). [CrossRef] [PubMed]

**8. **M. S. Ozturk, D. Rohrbach, U. Sunar, and X. Intes, “Mesoscopic fluorescence tomography of a photosensitizer (HPPH) 3D biodistribution in skin cancer,” Acad. Radiol. **21**(2), 271–280 (2014). [CrossRef] [PubMed]

**9. **M. S. Ozturk, V. K. Lee, L. Zhao, G. Dai, and X. Intes, “Mesoscopic fluorescence molecular tomography of reporter genes in bioprinted thick tissue,” J. Biomed. Opt. **18**(10), 100501 (2013). [CrossRef] [PubMed]

**10. **F. Yang, M. S. Ozturk, L. Zhao, W. Cong, G. Wang, and X. Intes, “High-Resolution Mesoscopic Fluorescence Molecular Tomography Based on Compressive Sensing,” IEEE Trans. Biomed. Eng. **62**(1), 248–255 (2015). [CrossRef] [PubMed]

**11. **N. Ouakli, E. Guevara, S. Dubeau, E. Beaumont, and F. Lesage, “Laminar optical tomography of the hemodynamic response in the lumbar spinal cord of rats,” Opt. Express **18**(10), 10068–10077 (2010). [CrossRef] [PubMed]

**12. **M. S. Ozturk, X. Intes, and V. K. Lee, “Longitudinal Volumetric Assessment of Glioblastoma Brain Tumor in 3D Bio-Printed Environment by Mesoscopic Fluorescence Molecular Tomography,” in *Optics and the Brain*, (Optical Society of America, 2016), JM3A. 46.

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

**14. **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]

**15. **J. Chen and X. Intes, “Comparison of Monte Carlo methods for fluorescence molecular tomography-computational efficiency,” Med. Phys. **38**(10), 5788–5798 (2011). [CrossRef] [PubMed]

**16. **Q. Fang, “Monte Carlo eXtreme (MCX),” http://mcx.space/.

**17. **L. Zhao, H. Yang, W. Cong, G. Wang, and X. Intes, “Lp regularization for early gate fluorescence molecular tomography,” Opt. Lett. **39**(14), 4156–4159 (2014). [CrossRef] [PubMed]

**18. **R. Yao, Q. Pian, and X. Intes, “Wide-field fluorescence molecular tomography with compressive sensing based preconditioning,” Biomed. Opt. Express **6**(12), 4887–4898 (2015). [CrossRef] [PubMed]

**19. **A. Beck and M. Teboulle, “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems,” SIAM J. Imaging Sci. **2**(1), 183–202 (2009). [CrossRef]

**20. **P. Mohajerani and V. Ntziachristos, “Compression of Born ratio for fluorescence molecular tomography/x-ray computed tomography hybrid imaging: methodology and in vivo validation,” Opt. Lett. **38**(13), 2324–2326 (2013). [CrossRef] [PubMed]

**21. **X. Cao, X. Wang, B. Zhang, F. Liu, J. Luo, and J. Bai, “Accelerated image reconstruction in fluorescence molecular tomography using dimension reduction,” Biomed. Opt. Express **4**(1), 1–14 (2013). [CrossRef] [PubMed]

**22. **J. Zhang, J. Shi, S. Zuo, F. Liu, J. Bai, and J. Luo, “Fast reconstruction in fluorescence molecular tomography using data compression of intra- and inter-projections,” Chin. Opt. Lett. **13**(7), 071002 (2015). [CrossRef]

**23. **J. S. Silva, J. Cancela, and L. Teixeira, “Fast volumetric registration method for tumor follow-up in pulmonary CT exams,” J. Appl. Clin. Med. Phys. **12**(2), 3450 (2011). [CrossRef] [PubMed]

**24. **Q. Pian, R. Yao, N. Sinsuebphon, and X. Intes, “Compressive Hyperspectral Time-resolved Wide-Field Fluorescence Lifetime Imaging,” Nat. Photonics . **11**411–417 (2017).

**25. **Q. Pian, R. Yao, L. Zhao, and X. Intes, “Hyperspectral time-resolved wide-field fluorescence molecular tomography based on structured light and single-pixel detection,” Opt. Lett. **40**(3), 431–434 (2015). [CrossRef] [PubMed]