## Abstract

Dynamic fluorescence molecular tomography (FMT) is an attractive imaging technique for three-dimensionally resolving the metabolic process of fluorescent biomarkers in small animal. When combined with compartmental modeling, dynamic FMT can be used to obtain parametric images which can provide quantitative pharmacokinetic information for drug development and metabolic research. However, the computational burden of dynamic FMT is extremely huge due to its large data sets arising from the long measurement process and the densely sampling device. In this work, we propose to accelerate the reconstruction process of dynamic FMT based on principal component analysis (PCA). Taking advantage of the compression property of PCA, the dimension of the sub weight matrix used for solving the inverse problem is reduced by retaining only a few principal components which can retain most of the effective information of the sub weight matrix. Therefore, the reconstruction process of dynamic FMT can be accelerated by solving the smaller scale inverse problem. Numerical simulation and mouse experiment are performed to validate the performance of the proposed method. Results show that the proposed method can greatly accelerate the reconstruction of parametric images in dynamic FMT almost without degradation in image quality.

© 2015 Optical Society of America

## 1. Introduction

Dynamic fluorescence molecular tomography (FMT) is a promising imaging technique that allows dynamically and three-dimensionally investigating the metabolic process of fluorescent biomarkers within small animals *in vivo* [1–7]. Meanwhile, compartmental modeling is a well-known method for pharmacokinetic analysis, and the pharmacokinetic parameters defined in the compartmental model can be utilized to describe the uptake and excretion of substances such as drugs, radiotracers and fluorophores in the body [6]. The boundary measurements of dynamic FMT can be transformed into three-dimensional (3-D) images of pharmacokinetic parameters through compartmental modeling [4–7], as has been done in positron emission tomography (PET) [8–10] and magnetic resonance imaging (MRI) [11–13]. Images of pharmacokinetic parameters are also known as parametric images, which can provide valuable physiological information for tumor detection [14–16], therapy assessment [17], and organ function evaluation [5–7, 18].

In order to obtain high-quality parametric images in dynamic FMT problem, we have proposed a novel full-direct method by applying regularization on the parametric images [7]. Compared with conventional indirect methods which reconstruct a sequence of static FMT images first and then estimate parametric images from the concentration-time curve of each voxel in a second step [5], the proposed direct method can reconstruct parametric images directly from boundary measurements by combining dynamic FMT reconstruction and compartmental modeling into one step. Therefore, the proposed direct method can make full use of temporal correlations of boundary measurements to improve the quality of parametric images [7]. Besides, the proposed direct method incorporates structural priors obtained from an X-ray computed tomography (XCT) system based on Laplace regularization to mitigate the inherent ill-posedness of FMT [19, 20], which can further improve the image quality.

However, the computational burden of the proposed direct method is extremely huge. This is because the boundary measurements used in the dynamic FMT problem are much larger than that in the static FMT problem. For dynamic FMT problem, the imaged small animal needs to be continuously rotated for multiple circles in a free-space full-angle FMT system [21], in order to monitor the metabolic process of fluorephores in the body. While the imaged small animal only needs to be rotated for one circle to acquire boundary measurements for the conventional static FMT problem. Additionally, the utilization of charge-coupled device (CCD) camera in the FMT system aggravates this problem. In each circle, multiple projections are acquired by CCD camera during the rotational process of the imaged small animal, and each projection can produce lots of measurement points.

To accelerate the static FMT inverse problem with huge amounts of boundary measurements, several compression approaches have been previously proposed based on Fourier transform [22], wavelet transform [23–25] and principal component analysis (PCA) [26, 27]. These approaches can obtain a smaller scale FMT inverse problem using data compression and dimensional reduction, thus the reconstruction process is accelerated. Despite of these successful applications, related work on how to accelerate the reconstruction process of parametric images in dynamic FMT has not been reported previously, although the acceleration for dynamic FMT inverse problem is more urgent than the static FMT inverse problem. One main reason is that the complexity of dynamic FMT inverse problem is greatly increased due to the incorporation of compartmental modeling, which makes it difficult to develop acceleration methods for the dynamic FMT inverse problem.

In this work, we present an acceleration method based on PCA for fast reconstructing parametric images in dynamic FMT problem. First, in the forward problem of dynamic FMT, the boundary measurements are acquired based on projection angles, and a sub weight matrix is assembled and used to formulate the matrix equation for directly mapping the parametric images to the grouped boundary measurements in each projection angle. The sub weight matrix is the same for all circles in the same projection angle. Second, in the inverse problem of dynamic FMT, PCA is used to analyze the sub weight matrix. The rows of the sub weight matrix have considerable correlations because they are from neighbor source-detector pairs in the same projection, thus the effective information of the sub weight matrix can be compressed into the first few principal components with the higher order components containing little useful information [27]. Finally, the dimension of sub weight matrix is reduced by discarding the less significant components and only retaining the first few larger components. Therefore, the reconstruction of parametric images in dynamic FMT problem can be accelerated by solving the smaller scale inverse problem. Numerical simulation and mouse experiment are carried out to test the performances of the proposed method. Results show that the proposed method can greatly accelerate the reconstruction process of dynamic FMT almost without degradation in reconstructed image quality.

This paper is organized as follows. Section 2 describes the proposed acceleration method for dynamic FMT as well as the evaluation metrics. Section 3 shows the results of numerical simulation and mouse experiment. Finally, Section 4 discusses the results and draws the conclusion.

## 2. Methods

#### 2.1. Diffusion model

In the near infrared spectral window, photon propagation in biological tissues can be modeled using the diffusion equation [28]. For the static FMT problem, a coupled diffusion equation is used to describe the propagation of excitation and emission lights. The coupled diffusion equation in continuous wave domain is given by [20]

*x*and

*m*denote the excitation and emission wavelengths, respectively; $S(r)$ is the excitation source; ${\Phi}_{x,m}(r)$ represents the optical field in the medium $\Omega $; ${\mu}_{ax,m}(r)$ represents the absorption coefficient, and ${D}_{x,m}(r)$ represents the diffusion coefficient; $n(r)$ denotes the fluorophore concentration, and $\eta $ is a constant accounting for the absorption and quantum yield of the fluorophores.

Here, the Robin boundary condition for the coupled diffusion equation is used [29]

In the diffusion equation, a collimated point source can be modeled as an isotropic source $S(r)\text{\hspace{0.17em}}=\delta (r-{r}_{s})$, where $\left({\mathrm{min}}^{-1}\right)$ is the point located one transport mean free path into the medium from the illumination spot. Similarly, a uniform line-shaped excitation source can be modeled as [21]

The fluorescence measurement detected at a point ${r}_{d}$ due to a line-shaped excitation source $L{S}_{s}$ can be obtained based on the first order Born approximation [30]

*s*to an arbitrary position $r$ inside the medium at the excitation wavelength; the Green’s function ${G}_{m}({r}_{d},r)$ describes the light propagation from the position $r$ inside the medium to the detector position ${r}_{d}$ at the emission wavelength; $\Theta $ is a calibration factor which accounts for the unknown gain of the system.

#### 2.2. Dynamic forward problem

In the dynamic FMT problem, the concentration of fluorophores varies with time continuously. The time-varying concentration $n(r,t)$ is usually caused by the metabolism of fluorophores in the body, and can provide valuable physiological information for biological studies, disease diagnosis and drug development [7]. Recently, indocyanine green (ICG) has been widely used in the field of optical imaging because of its low toxicity and approval for human use [4, 5, 14, 15]. Based on compartmental modeling, the metabolic processes of ICG in the organs and tissues can be described by a bi-exponential model [5, 14]

where the parameters*A*,

*B*, $\alpha $ and $\beta $ are the pharmacokinetic parameters describing the metabolic process of ICG [14]. When the gain of the FMT system is unknown, the pharmacokinetic parameters $A$ and $B$ have an arbitrary unit (a.u.) [5]. $\alpha $ and $\beta $ (${\mathrm{min}}^{-1}$) are the uptake and excretion rates of fluorophores, which have physiological significance for quantitative evaluation of organ function [31].

In our previous work [7], ICG has been used to investigate the metabolic features of organs in the mouse. Each organ in the body plays a different role in the metabolization of fluorophores, and thus has a distinctive temporal pattern of fluorophores uptake and release when fluorophores is injected into the body [32, 33]. Therefore, parameters *A*, *B*, $\alpha $ and $\beta $ can be used to quantitatively evaluate the functions of organs and tissues, such as liver function [31]. The images of pharmacokinetic parameters are the parametric images to be reconstructed.

When the imaged domain is discretized into $N$ voxels, and ${r}_{j}$ ($j=1,\cdots ,N$) is employed to denote the spatial locations of voxels, the parametric images in the discrete domain can be defined as

where each parametric image ${x}_{u}$ ($u=1,2,3,4$) is given by ${x}_{u}\text{=}{\left[{x}_{u}({r}_{1}),\cdots ,{x}_{u}({r}_{N})\right]}^{T}$.The boundary measurements of the dynamic FMT problem are obtained based on a hybrid FMT/XCT system as shown in Fig. 1 [34]. The imaged small animal is fixed on a rotation stage. Around the rotation stage, a line-shaped source and a CCD camera are positioned to generate and acquire the fluorescence data sets respectively, while the XCT system is used to acquire the anatomical information of the small animal. For dynamic FMT problem, the small animal is continuously rotated for multiple circles in order to monitor the metabolic process of the fluorophores. Figure 2 shows the schematic diagram of the measurement scheme. Suppose that *S* projections are acquired by the CCD camera in each circle and the small animal is rotated for *L* circles continuously. For projection *s* ($s=1,\cdots ,S$), ${M}_{s}$ measurement points are acquired. And for each circle, $M={\displaystyle {\sum}_{s=1}^{S}{M}_{s}}$measurement points are collected. Let ${t}_{k}$ ($k=1,\cdots ,K$) denote the time at which the *k-*th projection is acquired, where $K\text{=}LS$ is the total number of projections. Note that the subscript *k* can also be defined as $k=(l-1)S+s$($s=1,\cdots ,S;\text{\hspace{0.17em}}\text{\hspace{0.17em}}l=1,\cdots ,L$) (see Fig. 2). Then, the time-varying concentration at time ${t}_{k}$ can be defined in the discrete domain as

By combining Eq. (4) and Eq. (8), a mapping from the pharmacokinetic parameters to the boundary measurements is obtained in the discrete domain, as

*s*, and its entries are defined aswhere $\Delta V$ is the volume of each individual voxel and arises from the discretization of the volume integral in Eq. (4). Similarly, $W={\left\{{W}_{s}\right\}}_{s=1}^{S}$ of size $M\times N$ is defined as the whole weight matrix for each circle. Note that ${W}_{s}$ in projection

*s*is the same for all

*L*circles, and $W$ is also the same for all

*L*circles.

Let $f(X)$ denote the forward model which maps the parametric image matrix to the boundary measurements. Then the boundary measurements predicted by the forward model can be obtained by assembling the values calculated in Eq. (9), as follows

*s*at time ${t}_{k}$ can be written as the following matrix form

Similarly, the boundary measurements obtained from CCD can be defined as

*L*circles.

#### 2.3. Acceleration of inverse problem

In order to fully utilize structural priors and temporal correlations of boundary measurements to improve the quality of parametric images, we have proposed a full-direct method under new regularization concept in a previous work [7]. The proposed method directly applies the regularization matrix on the parametric images, and the objective function is given by [7]

Although the proposed method has been proven to be effective in improving the reconstruction quality, the computational burden of the inverse problem is extremely large. This is because huge amounts of boundary measurements are collected by the CCD camera for the dynamic FMT problem, in which the imaged small animal has to be continuously investigated for a long time in order to monitor the metabolic process of fluorophores in the body. Therefore, acceleration for the dynamic FMT problem is necessary. In this work, we use PCA to reduce the dimension of the weight matrix and thus achieve acceleration for the reconstruction of parametric images in the dynamic FMT problem.

PCA is a common statistical processing method for compressing high-dimensional data into a lower-dimensional form by choosing only the highest variance components of the data set. Using orthogonal projection, PCA linearly transforms an original set of variables into another set of uncorrelated variables, while the first few uncorrelated variables can contain most of the information in the original set of variables. From the correlation point of view, the row vectors of the sub weight matrix ${W}_{s}$ are highly correlative since they correspond to the neighboring detectors with the same excitation source [27]. By retaining only the first few principle components with higher variance, PCA can effectively discard the correlative information in the weight matrix and thus reduce the dimension of the weight matrix.

Since the sub weight matrix ${W}_{s}$ with a size of ${M}_{s}\times N$ is constructed based on the ${M}_{s}$ detectors in projection *s*, in the view of PCA, the rows of ${W}_{s}$ correspond to the ${M}_{s}$variables, and the columns of ${W}_{s}$ correspond to the *N* observations. Then, the covariance matrix $C={W}_{s}{W}_{s}^{T}$ of the sub weight matrix ${W}_{s}^{T}$ can be diagonalized to find eigenvalues (variances) and eigenvectors (components), as follows

With the matrix of eigenvectors, the principal components ${\overline{W}}_{s}$ of the sub weight matrix ${W}_{s}$ can be obtain by

Besides, for projection*s*at time ${t}_{k}$ ($k=(l-1)S+s;\text{\hspace{0.17em}}\text{\hspace{0.17em}}l=1,\cdots ,L$), by projecting the boundary measurements predicted by the forward model onto the eigenvectors, the corresponding score vector is obtained as followsSimilarly, the score vector of the boundary measurements obtained from CCD is given byTherefore, with the implementation of PCA, Eq. (13) can be modified as

By taking advantage of the compression property of PCA, we can retain the first $\theta $ largest principal components and leave out the rest less significant principal components of ${\overline{W}}_{s}$ and $\overline{f}(X,{t}_{k})$, thus a dimension-reduced matrix equation can be given as follows

where the size of ${\overline{W}}_{s}^{\theta}$ is $\theta \times N$. In this work, the cumulative percent of variance ($CPV$) is used to determine the number of retained principal components $\theta $ [27]*CPV*reaches to a preset threshold (e.g., 0.90). Note that the values of $\theta $ could be different in different projections although the preset threshold of

*CPV*is the same for all projections. Thus, after retaining only the first $\theta $ largest principal components in Eq. (23), the boundary measurements of all circles defined in Eqs. (12) and (15) are compressed into

*CPV*will lead to a higher compression degree.

With the dimensional reduction of boundary measurements, the objective function (17) can be modified as

#### 2.4. Evaluation metrics

In order to quantitatively evaluate the reconstruction quality, the normalized root-mean-square error (NRMSE) was used as an index

where ${X}_{PCA}$ is the 3-D parametric images reconstructed by method with PCA acceleration; ${X}_{org}$ is the 3-D parametric images reconstructed by method without PCA acceleration. The NRMSE was used to quantitatively evaluate the differences between the parametric images obtained with the two methods.Besides, the time cost of the reconstruction process was calculated and used to quantitatively evaluate the computational burden of the two methods.

## 3. Results

Numerical simulation and *in vivo* mouse experiment were carried out to evaluate the performances of the proposed method with PCA acceleration. The programs were performed on a personal computer with Intel(R) Core (TM) i7-2600 CPU (3.4 GHz) and 16 GB memory. The software environment was Matlab R2011a (MathWorks, Natick, MA, USA). All reconstruction results were 3-D parametric images with a voxel size of 1 mm$\times $1 mm$\times $1 mm. Based on the structural priors, by assigning the corresponding optical properties to the segmented organs [6], heterogeneous FMT forward models were constructed and used in both the numerical simulation and mouse experiment.

#### 3.1. Numerical simulation

### 3.1.1. Simulation setups

Numerical simulation was first implemented to validate the performance of the proposed method with PCA acceleration. In the numerical simulation, the Digimouse atlas [35] was employed to construct a 3-D simulation model. As shown in Fig. 4(a), the simulation model includes three different organs: liver, lungs and kidneys. In the liver, an ischemia-reperfusion injury (IRI) region (3 mm$\times $3 mm$\times $3 mm) with different metabolic features was simulated [31]. Figure 4(b) shows the ICG concentration curves which mimic the metabolic processes of ICG in different organs and other tissues. The ICG curves are obtained according to Eq. (5), and the corresponding pharmacokinetic parameters are shown in Table 1 [32, 33].

Figure 5 shows the segmented XCT results of the numerical simulation. Figure 5(a) shows the segmented XCT slice across the center of the IRI region. This representative slice of the liver is indicated by the red line shown in Fig. 4(a). It was segmented into the liver and other tissues, and the IRI region cannot be seen in the XCT images. Figures 5(b) and 5(c) show the segmented XCT slices indicated by the cyan and blue lines shown in Fig. 4(a), and they are the representative slices of the lungs and kidneys, respectively. The whole simulation model was segmented to the liver, lungs, kidneys and other tissues. The segmentation results were used as the structural priors to construct the Laplacian-type matrix [19].

The dynamic FMT boundary measurements of the simulation model were generated based on the scheme shown in Fig. 2. A line-shaped excitation source with the same length as the simulation model was used. The simulation model was continousely rotated for 60 circles (*L* = 60) within 60 min. For each circle, 24 projections (*S* = 24) were acquired with an angular increment of 15°. The ICG concentrations in different organs and tissues varied in each projection following the curves shown in Fig. 4(b). Then, $K\text{=}LS\text{=}1440$ projections in total were acquired. The simulation model was discretized into 8743 voxels (*N* = 8743) with a voxel size of $1\text{\hspace{0.17em}}\text{mm}\times 1\text{\hspace{0.17em}}\text{mm}\times 1\text{\hspace{0.17em}}\text{mm}$. For each circle, 9443 measurement points (*M* = 9443) were collected from the 24 projections. Thus, *P = LM* = 566580 measurement points (i.e., the boundary measurements) in total were acquired for the 60 circles.

### 3.1.2. Simulation results

The proposed method was implemented to process the boundary measurements and obtain the parametric images. *CPV* = 0.90 was used to determine the number of retained large principal components ($\theta $) of each sub weight matrix ${\overline{W}}_{s}$. The reconstruction results are shown in Figs. 6–8.

Figure 6 shows the parametric images corresponding to the slice shown in Fig. 5(a). Figures 6(a)–6(d) are the true parametric images obtained according to Table 1. The IRI region in the liver has different pharmacokinetic parameters and thus can be distinguished from the liver in the parametric images. Figures 6(e)–6(l) show the parametric images reconstructed by the direct methods without and with PCA acceleration, respectively. There are almost no visual differences in the results reconstructed from the two methods. Figures 7 and 8 show the parametric images corresponding to the slices shown in Figs. 5(b) and 5(c), respectively. Similarly to Fig. 6, the two direct methods without and with PCA acceleration obtain almost the same reconstruction results.

Table 2 shows the sizes of the original and reduced weight matrices, as well as the time cost of the two methods for reconstructing the parametric images. As shown in the table, the time cost in the calculation of forward problem is the same, since this step has no difference for the two methods. The time cost in performing PCA is only available for the method with PCA acceleration, but this step takes very little time. In the step of solving the inverse problem, the time cost of the method with PCA acceleration is much lower than the method without PCA acceleration, because the reduced weight matrix has a much smaller scale than the original weight matrix. As a result, the total time cost of the proposed method is only about 1/8 of the total time cost by using the original weight matrix. These results demonstrate that the proposed method can effectively accelerate the reconstruction process.

### 3.1.3. Influence of CPV

In order to test the influence of the retained principal component number for each sub weight matrix on the reconstructed parametric image quality and the computation time, we investigated the reconstructed results using different *CPVs*. The NRMSE was calculated to quantitatively evaluate the differences between the parametric images reconstructed by the two methods without and with PCA acceleration.

Figure 9 shows the reconstructed parametric images corresponding to the slice shown in Fig. 5(a) when the *CPV* is set as 0.98, 0.95, 0.90, 0.80, 0.70 and 0.60, respectively. It can be seen that the reconstruction results generally have no visual difference between the parametric images reconstructed by the two methods without and with PCA acceleration when *CPV* is above 0.80. Nevertheless, when *CPV* is set as 0.70 and 0.60, the IRI region in the parametric images begins to become distorted. It means that part of valuable information has been discarded when leaving out less significant principal components of the sub weight matrix.

To quantify the above observations, Fig. 10 shows the corresponding NRMSE and computation time as a function of *CPV*, respectively. It can be seen that a larger *CPV* decreases the NRMSE and improves the reconstruction quality of parametric images at the cost of more computation time. A compromise method is to choose an appropriate *CPV* which can achieve a balance between image quality and computation time. It can be observed in Fig. 10 that the reduction of computation time becomes less significant when *CPV* is less than 0.90. Thus, the *CPV* of 0.90 is suggested as a compromise between image quality and computation time.

#### 3.2. Mouse experiment

### 3.2.1. Experiment setups

Mouse experiment was conducted under the protocol approved by the Institutional Animal Care and Use Committee of Tsinghua University. Briefly, a healthy BALB/c nude mouse (about 8 weeks) was fixed on the rotation stage of the hybrid FMT/XCT system [35] and anesthetized during the experiment. A bolus of ICG ($\text{0}\text{.1}\text{mL}$ of $\text{50}\text{\hspace{0.17em}}\text{\mu g/mL}$) was injected through the tail vein. Then the mouse was continuously rotated for 70 circles (*L* = 70) to monitor the metabolic process of ICG in the body. For each circle, 24 projections (*S* = 24) were acquired with an angular increment of 15°. In total, $K\text{=}LS\text{=}1680$ projections were acquired. A line-shaped excitation source with a length of 4 cm provided by a Xenon lamp was used as the excitation source. Fluorescence images were acquired with a 775$\pm $6 nm excitation filter and an 840$\pm $6 nm emission filter. All projections were recorded by a 512$\times $512 pixel, −70°C cooled CCD camera (Andor, Belfast, Northern Ireland, U.K.). The exposure time of CCD was set to 1 s, and the CCD binning was set to 2$\times $2. The rotational speed of the mouse was 8°/s, and the total rotation time of full 360° was 45 s. As a result, the FMT imaging time for each circle was about 69 s, and the total FMT imaging time for all circles was about 80 min. After collecting the fluoresence datesets, a hepatobiliary contrast agent for XCT imaging, Fenestra LC (Advanced Research Technologies, Montreal, CA) was injected at a dose of 15 mL/kg body weight through the tail vein. One hour after the injection, XCT data sets were collected to obtain the anatomical information of the mouse. Finally, a steel anchor point was plastered on the surface of the mouse to provide the relative height information for data set registration between the XCT and FMT systems.

For dynamic FMT reconstruction, the investigated region of the mouse was discretized into 5662 voxels (*N* = 5662) with a voxel size of $1\text{\hspace{0.17em}}\text{mm}\times 1\text{\hspace{0.17em}}\text{mm}\times 1\text{\hspace{0.17em}}\text{mm}$. For each circle, 5058 measurement points (*M* = 5058) were collected from the 24 projections. Thus, *P = LM* = 354060 measurement points (i.e., the boundary measurements) in total were acquired for the 70 circles.

Figure 11 shows the XCT results of the mouse. Figure 11(a) is a representative coronal XCT slice, and shows the region selected for 3-D reconstruction. Figures 11(b)–11(d) are the transversal XCT slices indicated by the green, cyan and blue dashed lines in Fig. 11(a), respectively. All transversal XCT slices were manually segmented to the liver, lungs, kidneys and other tissues. Figures 11(e)–11(g) are the segmentation results for Figs. 11(b)–11(d), respectively. A heterogeneous FMT forward model was constructed by assigning the corresponding optical parameters to the segmented organs [5]. The segmentation results were used as the structural priors to construct the Laplacian-type matrix.

### 3.2.2. Experiment results

Similarly to the numerical simulation, the proposed method was implemented to obtain parametric images of the mouse experiment. *CPV* = 0.90 was used to determine the number of retained large principal components ($\theta $) of each sub weight matrix ${\overline{W}}_{s}$.

Figures 12, 13, and 14 show the parametric images of the mouse experiment obtained with the methods without and with PCA acceleration, respectively. Since the pharmacokinetic parameters $A$ and $B$ only have an arbitrary unit without knowing the gain of the FMT system, the two parameters are normalized by the mean value of $A$ in the liver [7]. Figure 12 shows the parametric images corresponding to the representative slice of the liver shown in Fig. 11(e). Similarly to the numerical simulation, it can be observed in Fig. 12 that the reconstruction results obtained with the two methods almost have no visual difference. Figure 12 shows the parametric images corresponding to the representative slice of the lungs shown in Fig. 11(f), while Fig. 14 shows the parametric images corresponding to the representative slice of the kidneys shown in Fig. 11(g). Similarly to Fig. 12, the two methods obtain almost the same reconstruction results by visual assessment.

Table 3 shows the sizes of the original and reduced weight matrices, as well as the time cost of the two methods for reconstructing the parametric images. As shown in the table, because of using the reduced weight matrix to solve the inverse problem, the total time cost of the method with PCA acceleration is only about 1/5 of the total time cost by using the method without PCA acceleration. Similarly to the numerical simulation, the proposed method can effectively accelerate the reconstruction process in the mouse experiment almost without degrading the qualities of the reconstructed parametric images.

### 3.2.3. Influence of CPV

In the mouse experiment, we also test the influence of *CPV* on the reconstructed parametric image quality and the computation time. The NRMSE was used to quantitatively evaluate the differences between the parametric images reconstructed by the two methods without and with PCA acceleration.

Figure 15 shows the reconstructed parametric images corresponding to the slice shown in Fig. 11(e) when the *CPV* is set as 0.98, 0.95, 0.90, 0.80, 0.70 and 0.60, respectively. The reconstructed parametric images almost have no visual difference between the results acquired by the two methods without and with PCA acceleration when *CPV* is higher than 0.80. Nevertheless, when *CPV* is set as 0.70 and 0.60, the results begin to have a little distortion, especially for parametric images $B$ and $\beta $. This is in accordance with the results in numerical simulation.

To quantify the above results, Fig. 16 shows the corresponding NRMSE and computation time as a function of *CPV*, respectively. It can be seen that the NRMSE increases and the computation time decreases when *CPV* is reduced. The reduction of computation time becomes less significant when *CPV* is less than 0.90, while the parametric images obtained with *CPV* = 0.90 have a relatively small NRMSE compared with results reconstructed by method without PCA acceleration. Thus, the *CPV* of 0.90 can also achieve an appropriate compromise between image quality and computation time in the mouse experiment.

## 4. Discussion and conclusion

Dynamic FMT can provide pharmacokinetic information of fluorescent biomarkers in small animals with parametric images, and thus is important for metabolic research and drug development. However, the computational burden of dynamic FMT is extremely huge due to the large data sets arising from the long measurement process and the densely sampling device. The main aim of this paper is to accelerate the reconstruction process of dynamic FMT based on the data compression property of PCA.

It has been shown in the numerical simulation and mouse experiment that, when *CPV* = 0.90, parametric images reconstructed by the proposed method with PCA acceleration have almost the same image quality as those reconstructed by the method without PCA acceleration (see Figs. 6–8 and Figs. 12–14). Nevertheless, the reconstruction process of dynamic FMT has been greatly accelerated by the proposed method, as shown in Tables 2 and 3. Besides, the influence of compression degree on the image quality and computation time has also been studied. It can be observed in Figs. 10 and 16 that image quality and computation time are simultaneously reduced when the compression degree of the sub weight matrix is increased (i.e., lower *CPV*). By choosing an appropriate compression degree, the computation time can be greatly reduced with little loss in the image quality, because PCA can compress most of the effective information of the sub weight matrix into the first few principal components.

Although the proposed acceleration direct method for parametric image reconstruction is based on regularization framework, the acceleration strategy of the proposed method can also be extended to other direct methods [3, 4, 36]. In addition, the reconstruction process of dynamic FMT may be further accelerated by applying compression approaches in parametric image domain [25].

In conclusion, we have presented a novel method based on PCA for accelerating the reconstruction process of parametric images in dynamic FMT. By taking advantage of the compression property of PCA, the dimension of the sub weight matrix is reduced and thus the reconstruction process of dynamic FMT is accelerated. The results of numerical simulation and mouse experiment demonstrate that the proposed method can accelerate the dynamic FMT reconstruction process almost without quality degradation of parametric images.

## Acknowledgments

This work is supported by the National Basic Research Program of China (973) under Grant No. 2011CB707701; the National Natural Science Foundation of China under Grant No. 81227901, 81271617, 61322101, 61361160418, 61401246; the National Major Scientific Instrument and Equipment Development Project under Grant No. 2011YQ030114; the Fundamental Research Funds for the Central Universities under Grant No. 2015JBM026, 2013JBZ014; the China Postdoctoral Science Foundation under Grant No. 2015M570032, 2014M550073. F. Liu is supported in part by the Postdoctoral Fellowship of Tsinghua-Peking Center for Life Sciences.

## References and links

**1. **S. Patwardhan, S. Bloch, S. Achilefu, and J. Culver, “Time-dependent whole-body fluorescence tomography of probe bio-distributions in mice,” Opt. Express **13**(7), 2564–2577 (2005). [CrossRef] [PubMed]

**2. **X. Liu, F. Liu, Y. Zhang, and J. Bai, “Unmixing dynamic fluorescence diffuse optical tomography images with independent component analysis,” IEEE Trans. Med. Imaging **30**(9), 1591–1604 (2011). [CrossRef] [PubMed]

**3. **A. B. Milstein, K. J. Webb, and C. A. Bouman, “Estimation of kinetic model parameters in fluorescence optical diffusion tomography,” J. Opt. Soc. Am. A **22**(7), 1357–1368 (2005). [CrossRef] [PubMed]

**4. **B. Alacam and B. Yazici, “Direct reconstruction of pharmacokinetic-rate images of optical fluorophores from NIR measurements,” IEEE Trans. Med. Imaging **28**(9), 1337–1353 (2009). [CrossRef] [PubMed]

**5. **G. Zhang, F. Liu, B. Zhang, Y. He, J. Luo, and J. Bai, “Imaging of pharmacokinetic rates of indocyanine green in mouse liver with a hybrid fluorescence molecular tomography/x-ray computed tomography system,” J. Biomed. Opt. **18**(4), 040505 (2013). [CrossRef] [PubMed]

**6. **G. Zhang, F. Liu, H. Pu, W. He, J. Luo, and J. Bai, “A direct method with structural priors for imaging pharmacokinetic parameters in dynamic fluorescence molecular tomography,” IEEE Trans. Biomed. Eng. **61**(3), 986–990 (2014). [CrossRef] [PubMed]

**7. **G. Zhang, H. Pu, W. He, F. Liu, J. Luo, and J. Bai, “Full-direct method for imaging pharmacokinetic parameters in dynamic fluorescence molecular tomography,” Appl. Phys. Lett. **106**(8), 081110 (2015). [CrossRef]

**8. **M. E. Kamasak, C. A. Bouman, E. D. Morris, and K. Sauer, “Direct reconstruction of kinetic parameter images from dynamic PET data,” IEEE Trans. Med. Imaging **24**(5), 636–650 (2005). [CrossRef] [PubMed]

**9. **G. Wang and J. Qi, “Generalized algorithms for direct reconstruction of parametric images from dynamic PET data,” IEEE Trans. Med. Imaging **28**(11), 1717–1726 (2009). [CrossRef] [PubMed]

**10. **J. Tang, H. Kuwabara, D. F. Wong, and A. Rahmim, “Direct 4D reconstruction of parametric images incorporating anato-functional joint entropy,” Phys. Med. Biol. **55**(15), 4261–4272 (2010). [CrossRef] [PubMed]

**11. **C. C. Martin, R. F. Williams, J. H. Gao, L. D. Nickerson, J. Xiong, and P. T. Fox, “The pharmacokinetics of hyperpolarized xenon: Implications for cerebral MRI,” J. Magn. Reson. Imaging **7**(5), 848–854 (1997). [CrossRef] [PubMed]

**12. **B. M. Kelm, B. H. Menze, O. Nix, C. M. Zechmann, and F. A. Hamprecht, “Estimating kinetic parameter maps from dynamic contrast-enhanced MRI using spatial prior knowledge,” IEEE Trans. Med. Imaging **28**(10), 1534–1547 (2009). [CrossRef] [PubMed]

**13. **L. Chen, P. L. Choyke, T. H. Chan, C. Y. Chi, G. Wang, and Y. Wang, “Tissue-specific compartmental analysis for dynamic contrast-enhanced MR imaging of complex tumors,” IEEE Trans. Med. Imaging **30**(12), 2044–2058 (2011). [CrossRef] [PubMed]

**14. **M. Gurfinkel, A. B. Thompson, W. Ralston, T. L. Troy, A. L. Moore, T. A. Moore, J. D. Gust, D. Tatman, J. S. Reynolds, B. Muggenburg, K. Nikula, R. Pandey, R. H. Mayer, D. J. Hawrysz, and E. M. Sevick-Muraca, “Pharmacokinetics of ICG and HPPH-car for the detection of normal and tumor tissue using fluorescence, near-infrared reflectance imaging: A case study,” Photochem. Photobiol. **72**(1), 94–102 (2000). [CrossRef] [PubMed]

**15. **D. J. Cuccia, F. Bevilacqua, A. J. Durkin, S. Merritt, B. J. Tromberg, G. Gulsen, H. Yu, J. Wang, and O. Nalcioglu, “In vivo quantification of optical contrast agent dynamics in rat tumors by use of diffuse optical spectroscopy with magnetic resonance imaging coregistration,” Appl. Opt. **42**(16), 2940–2950 (2003). [PubMed]

**16. **M. Hsing, Y. Lin, and G. Gulsen, “Pharmacokinetic analysis for tumor characterization using MR-guided dynamic contrast enhanced diffuse optical tomography,” Biomed. Opt. and 3D Imag. BTu2A.3. (2012).

**17. **M. Choi, K. Choi, S. W. Ryu, J. Lee, and C. Choi, “Dynamic fluorescence imaging for multiparametric measurement of tumor vasculature,” J. Biomed. Opt. **16**(4), 046008 (2011). [CrossRef] [PubMed]

**18. **C. B. Amoozegar, T. Wang, M. B. Bouchard, A. F. H. McCaslin, W. S. Blaner, R. M. Levenson, and E. M. C. Hillman, “Dynamic contrast-enhanced optical imaging of in vivo organ function,” J. Biomed. Opt. **17**(9), 096003 (2012). [CrossRef] [PubMed]

**19. **R. B. Schulz, A. Ale, A. Sarantopoulos, M. Freyer, E. Soehngen, M. Zientkowska, and V. Ntziachristos, “Hybrid system for simultaneous fluorescence and x-ray computed tomography,” IEEE Trans. Med. Imaging **29**(2), 465–473 (2010). [CrossRef] [PubMed]

**20. **G. Zhang, X. Cao, B. Zhang, F. Liu, J. Luo, and J. Bai, “MAP estimation with structural priors for fluorescence molecular tomography,” Phys. Med. Biol. **58**(2), 351–372 (2012). [CrossRef] [PubMed]

**21. **F. Liu, X. Liu, D. Wang, B. Zhang, and J. Bai, “A parallel excitation based fluorescence molecular tomography system for whole-body simultaneous imaging of small animals,” Ann. Biomed. Eng. **38**(11), 3440–3448 (2010). [CrossRef] [PubMed]

**22. **J. Ripoll, “Hybrid Fourier-real space method for diffuse optical tomography,” Opt. Lett. **35**(5), 688–690 (2010). [CrossRef] [PubMed]

**23. **T. J. Rudge, V. Y. Soloviev, and S. R. Arridge, “Fast image reconstruction in fluorescence optical tomography using data compression,” Opt. Lett. **35**(5), 763–765 (2010). [CrossRef] [PubMed]

**24. **N. Ducros, C. D’andrea, G. Valentini, T. Rudge, S. Arridge, and A. Bassi, “Full-wavelet approach for fluorescence diffuse optical tomography with structured illumination,” Opt. Lett. **35**(21), 3676–3678 (2010). [CrossRef] [PubMed]

**25. **T. Correia, T. Rudge, M. Koch, V. Ntziachristos, and S. Arridge, “Wavelet-based data and solution compression for efficient image reconstruction in fluorescence diffuse optical tomography,” J. Biomed. Opt. **18**(8), 086008 (2013). [CrossRef] [PubMed]

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

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

**28. **A. Ishimaru, *Wave Propagation and Scattering in Random Media* (Academic, 1978).

**29. **M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. **22**(11), 1779–1792 (1995). [CrossRef] [PubMed]

**30. **R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissues with noncontact measurements,” IEEE Trans. Med. Imaging **23**(4), 492–500 (2004). [CrossRef] [PubMed]

**31. **H. Shinohara, A. Tanaka, T. Kitai, N. Yanabu, T. Inomoto, S. Satoh, E. Hatano, Y. Yamaoka, and K. Hirao, “Direct measurement of hepatic indocyanine green clearance with near-infrared spectroscopy: separate evaluation of uptake and removal,” Hepatology **23**(1), 137–144 (1996). [CrossRef] [PubMed]

**32. **E. M. C. Hillman and A. Moore, “All-optical anatomical co-registration for molecular imaging of small animals using dynamic contrast,” Nat. Photonics **1**(9), 526–530 (2007). [CrossRef] [PubMed]

**33. **L. Chen, T. H. Chan, P. L. Choyke, E. M. C. Hillman, C. Y. Chi, Z. M. Bhujwalla, G. Wang, S. S. Wang, Z. Szabo, and Y. Wang, “CAM-CM: A signal deconvolution tool for in vivo dynamic contrast-enhanced imaging of complex tissues,” Bioinformatics **27**(18), 2607–2609 (2011). [PubMed]

**34. **X. Guo, X. Liu, X. Wang, F. Tian, F. Liu, B. Zhang, G. Hu, and J. Bai, “A combined fluorescence and microcomputed tomography system for small animal imaging,” IEEE Trans. Biomed. Eng. **57**(12), 2876–2883 (2010). [CrossRef] [PubMed]

**35. **B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: A 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. **52**(3), 577–587 (2007). [CrossRef] [PubMed]

**36. **G. Zhang, H. Pu, W. He, F. Liu, J. Luo, and J. Bai, “Bayesian framework based direct reconstruction of fluorescence parametric images,” IEEE Trans. Med. Imaging. in press.