## 1.

## Introduction

Fluorescence molecular tomography (FMT) is molecular imaging method that uses the fluorescence signal at the surface of imaging objects to inversely calculate the distribution of the fluorescent yield^{1}^{,}^{2} in 3D, such that the location and concentration of the fluorochromes can be determined. FMT has been adopted as a reliable technology for the research of cardiovascular,^{3} pulmonary, oncology, inflammation, and infectious disease.^{4}5.^{–}^{6}

The complex geometry of small animals has a tremendous influence on the accuracy of FMT reconstruction, and this complex structural information cannot be acquired by FMT alone. A compression of the small animals’ bodies^{7} or a match solution of optical properties^{8} is needed in single-modality FMT because the geometry of small animals can be considered as a simple slab. However, these methods ignore the inner boundaries of small animals; therefore, the reconstruction accuracy is limited in the case of small animals with tissue of complex geometry.

In order to increase the accuracy of FMT reconstruction for small animals, structural information from other modalities, such as microcomputed tomography (micro-CT), ultrasound and magnetic resonance imaging (MRI), is introduced into the FMT reconstruction. This information not only can provide the accurate boundary conditions for forward calculations but also can provide an accurate constraint condition to reduce the ill-posedness of reconstruction and increase the convergence speed and accuracy in the inverse problem. For ultrasound, acoustic attenuation may influence the quantification and image-fidelity, and MRI/optics system demands the optical system to be fiber coupling or low-performance magnetic resonance-compatible cameras to be used, which allows few sources and detectors. Combined with CT, optical system can be engineered onto the rotating gantry, because of which, a high performance charge-coupled device (CCD) can be applied.^{9} Therefore, dual-modality micro-CT/FMT has become the mainstream in FMT technology. The dual-modality micro-CT/FMT system can provide both structural and functional information for small animals in biological research.^{10}11.12.13.14.^{–}^{15}

There are still some factors limiting the reconstruction accuracy of the dual-modality micro-CT/FMT system. Firstly, the diversity of the optical properties of small animals indicates that there is not only high-scattering tissue but also low-scattering tissue, for example, rabbit muscle, colon adenocarcinomas,^{16} the cerebral spinal fluid layer in the brain, and cysts in the human breast.^{17} Therefore, those method based on diffusion approximation with the finite element method may lead to error^{18}^{,}^{19} in the reconstruction for small animals with arbitrary optical properties (including both low-scattering and high-scattering tissue). The reconstruction method based on the radiative transfer equation (RTE) was proposed to satisfy both the low-scattering and high-scattering cases. Yet, RTE is time consuming to attain a high degree of accuracy.^{18} On the other hand, the Monte Carlo (MC) method based on statistics was proposed to depict the propagation process of photons in tissue, which is as accurate as RTE,^{20}21.^{–}^{22} and it is suitable for small animals with arbitrary optical properties. However, the huge time consumption of MC greatly limits this method’s application.^{23} Therefore, a graphics processing unit (GPU) is introduced to MC to accelerate its calculation speed.^{17}^{,}^{24}^{,}^{25} Different methods based on MC accelerated by GPU have been introduced in the FMT reconstruction.^{26}27.^{–}^{28}

Secondly, since the resolution of micro-CT is approximately tens or hundreds of microns, for example the resolution of our micro-CT is 75 *μ*m, while the resolution of FMT is approximately 1 mm,^{29} an efficient reconstruction method that can correctly incorporate the high-resolution structural information of micro-CT into a low-resolution FMT reconstruction is still required. Due to the mismatch of the resolution between micro-CT and FMT, it is hard to choose the reconstruction pixels size of FMT. When we use the pixel size of FMT as that of micro-CT, the memory and calculation time will be quite huge, which is far beyond the calculation capability one can access in most case. Therefore, voxel binning is needed for images form micro-CT to match the pixel size of FMT that is about 1 mm. After voxel binning, the low-resolution data set of the structural information that served as prior information will be incorporated into the regularization method of FMT.^{30}^{,}^{31} However, voxel binning will inevitably lead to a loss of structural information, especially in the irregular boundary, which will cause inaccuracy in the boundary definition, and ultimately will affect the accuracy of the reconstruction.

In this paper, we propose a micro-CT guided non-equal voxel MC (NVMC) and a modified Laplacian regularization method for a dual-modality micro-CT/FMT system. The NVMC is based on Monte Carlo eXtreme (MCX).^{17} With the guidance of micro-CT, the voxel at the irregular boundary or region of interest is precisely merged, while the voxel at the same region is coarsely merged. Through modifying the MC and Laplacian regularization method, which was proposed based on an equal-space voxel and made improvements in quantitation and linearity of response,^{32} the structural information of micro-CT with non-equal voxels can be accurately incorporated into the FMT reconstruction.

First, a detailed description and theoretical derivation of this method will be shown; then, through a simulation experiment, the validity of micro-CT guided NVMC in the forward problem will be proven, and the advantage of NVMC will be shown in the case of a curved surface boundary. A phantom experiment will further validate the NVMC and modified Laplacian regularization method in the reconstruction. Finally, small animal experiments will prove the validity of our method.

## 2.

## Method

In this section, a detailed description of micro-CT guided NVMC for the reconstruction of dual-modality micro-CT/FMT will be shown. This method can be listed as follows: first, with the guidance of micro-CT, we make a non-equal voxel binning, which makes precise voxel binning at the irregular boundary or region of interest but makes coarse voxel binning for voxels belonging to same region; then, the Green functions of FMT are calculated by NVMC, and with these Green functions, the Jacobian matrix is formed; last, the distribution of the fluorescence yield is reconstructed with a modified Laplacian regularization method for non-equal space voxels.

## 2.1.

### Micro-CT Guided Non-Equal Discretization of Voxels

The resolution of micro-CT in our dual-modality system is 75 *μ*m, and the resolution of FMT in our dual-modality system is 3 mm.^{13} To reduce the dataset and the time consumption in the reconstruction, the original data from micro-CT with voxel size of $75\times 75\times 75$ *μ*m^{3} are coarsely merged in to slices with a pixel pitch of $1\times 1\times 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$. Then, according to the thresholding image segmentation method, the coarsely merged slices are segmented into ${N}_{s}$ regions. Each value of the pixel in the same region is assigned to ${N}_{i}$, where ${N}_{i}\in [1,{N}_{S}]$. Figure 1(a) shows the coarsely merged slice with the voxel size of $1\times 1\times 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$.

To improve the accuracy and reduce the amount of voxels, the original micro-CT slices are used to guide the refinement of the coarsely merged slices, and two datasets of different sizes are prepared at the same time. The canny-edge detection method^{33} is applied in the original micro-CT slices to find the edge between different regions, such as bone, lung and muscle. Coarse voxels at the edge between different regions are refined into tiny voxels with a minimum edge length of $L$ ($L<1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$). The refined method is also applied in the region of interest. The index of the voxel in the coarsely merged slice that contains the refined voxels is saved in a file for the MC calculation. The regions that contain the refined voxel are also recorded in a file. Figure 1(b) shows the refined voxel.

## 2.2.

### Monte Carlo Simulation for Non-Equal Voxels

After the micro-CT guided non-equal discretization of the voxels, the value of the optical properties for each voxel must be assigned. Because the relationship between each voxel and its region index (${N}_{i}$) was determined in Sec. 2.1, we only assign the optical properties to each region. With all the structural information mentioned in Sec. 2.1 and the optical properties of each region, the propagation of photons in the tissue can be calculated. However, the classical MC is proposed for equal-space voxels.^{17}^{,}^{20}^{,}^{22} Therefore, we modified the classical MC method to be suitable for the case of non-equal voxels, and GPU is used to speed up the MC calculation. The differences between the NVMC and the classical MC based on equal voxels have been listed below, and Fig. 2 shows the flowchart of NVMC.

(1) Modification of the photon substeps in voxel. Because the dataset of non-equal space voxels contains voxels with different sizes, we assign the value of photon substeps in a voxel to $L$, which is the minimum edge length of the minimum refined voxel mentioned in Sec. 2.1. This modification can ensure the photons travel inside the minimum voxel without stepping across it. Therefore, the calculation accuracy in the minimum voxel can be guaranteed. This method is used in the papers of Fang

^{17}and Boas,^{21}but only in equal voxel MC method, and we’ll analyse the validity of the modification in the follow simulation experiment.(2) Modification of the storage space of raw probability. The storage space of raw probability for equal voxels can be sequentially accessed according to the position of the voxel. For non-equal voxels, this is not possible because the relationship between the position of the voxel and the index of the voxel is not linear. Therefore, two kinds of storage space of raw probability are created: one is used to store the raw probability of the coarse voxel, while the other is used to store the raw probability of the refined voxel.

(3) Modification of the normalization of raw probability. Because the volume of the voxels is different in the dataset of non-equal voxels, the calculated raw probability must be normalized by the volume of each voxel.

## 2.3.

### Modified Laplacian Regularization for Non-Equal Voxels

Calculation of the Jacobian matrix is the main task of FMT reconstruction; it can be expressed as the equation shown in Eq. (1), which is derived from Ref. 26. We use the NVMC to calculate the Green function in Eq. (1)^{28} and obtain the Jacobian matrix.

## (1)

$$J({R}_{d},{R}_{s})=\frac{\partial [{\varphi}_{f}({R}_{d},{R}_{s})]}{\partial x}=\int {g}^{\lambda em}(r,{R}_{d}){g}^{\lambda ex}(r,{R}_{s})\mathrm{d}r,$$The Laplacian regularization method is used to reconstruct the distribution of the fluorescence yield; the equation is shown in Eq. (2):

where $J$ is the Jacobian matrix shown in Eq. (1), and $y$ is the fluorescence signal detected by the experiment; $x$ is the reconstructed fluorescence yield $x=\eta {\mu}_{af}$; ${\mu}_{af}$ is the absorption coefficient of the fluorochromes; $\eta $ is the quantum efficiency of the fluorochromes; $\lambda $ is the relaxation factor determined by the L-curve method;^{34}the conjugate gradient method is used to solve Eq. (2).

^{35}Laplacian regularization matrix $L$ is determined by the structural information provided by micro-CT, and $L$ can relieve the ill-posedness of the reconstruction and increase the convergence speed. The classical Laplacian matrix is derived from a finite difference approximation on an ‘$N$’ number of equal spaces for a Laplace equation.

^{30}It can be expressed as Eq. (3):

## (3)

$${\nabla}^{2}x(r){h}^{2}\asymp {x}_{1}+{x}_{2}+\text{\hspace{0.17em}}\dots -N{x}_{N/2}+\dots +{x}_{N}=0,$$The classical Laplacian regularization method is proposed based on the equal-space voxel, and it may not be suitable for the case of non-equal voxels. Therefore, we modified the classical Laplacian regularization method. If the region contains a voxel with the same size, Eq. (3) is correct, but if a region contains voxels with different sizes, a new discretion of the bigger voxel must be made. Because the bigger voxel is divided into several minimum voxels, the space of the region is equal again, but the volume of the equal-space voxel is ${V}_{\mathrm{min}}$, and total number of voxels is ${N}_{\mathrm{min}}$, where, ${N}_{\mathrm{min}}=\sum _{i=1}^{N}\frac{{V}_{i}}{{V}_{\mathrm{min}}}$ and the step size is ${h}_{\mathrm{min}}$, such that Eq. (3) can be expressed as follows:

## (4)

$${\nabla}^{2}x(r){h}_{\mathrm{min}}^{2}\asymp {x}_{1}+{x}_{2}+\text{\hspace{0.17em}}\dots -{N}_{\mathrm{min}}{x}_{{N}_{\mathrm{min}}/2}+\dots +{x}_{{N}_{\mathrm{min}}}=0.$$## (5)

$${x}_{Vi}\approx \frac{\sum _{j=1}^{{V}_{i}/{V}_{\mathrm{min}}}{x}_{j}}{{V}_{i}/{V}_{\mathrm{min}}}.$$Substituting Eq. (5) into Eq. (4) and dividing by $-{N}_{\mathrm{min}}$ on both sides of Eq. (4), Eq. (6) can be derived:

## (6)

$$\frac{-{x}_{1}}{{N}_{\mathrm{min}}}+\frac{-({V}_{i}/{V}_{\mathrm{min}}){x}_{Vi}}{{N}_{\mathrm{min}}}+\text{\hspace{0.17em}}\dots +{x}_{{N}_{\mathrm{min}}/2}+\dots \text{\hspace{0.17em}}\frac{-({V}_{j}/{V}_{\mathrm{min}}){x}_{Vj}}{{N}_{\mathrm{min}}}+\frac{-{x}_{{N}_{\mathrm{min}}}}{{N}_{\mathrm{min}}}=0.$$Therefore, the Laplacian regularization matrix can be expressed as follows:

## (7)

$${L}_{ij}=\{\begin{array}{lll}1& & i\text{\hspace{0.17em}}\mathrm{and}\text{\hspace{0.17em}}j\text{\hspace{0.17em}}\mathrm{is}\text{\hspace{0.17em}}\mathrm{the}\text{\hspace{0.17em}}\mathrm{same}\text{\hspace{0.17em}}\mathrm{point}\\ -\frac{{V}_{i}/{V}_{\mathrm{min}}}{{N}_{\mathrm{min}}}& & i\text{\hspace{0.17em}}\mathrm{and}\text{\hspace{0.17em}}j\text{\hspace{0.17em}}\mathrm{is}\text{\hspace{0.17em}}\mathrm{in}\text{\hspace{0.17em}}\mathrm{the}\text{\hspace{0.17em}}\mathrm{same}\text{\hspace{0.17em}}\mathrm{region}\\ 0& & i\text{\hspace{0.17em}}\mathrm{and}\text{\hspace{0.17em}}j\text{\hspace{0.17em}}\mathrm{is}\text{\hspace{0.17em}}\mathrm{not}\text{\hspace{0.17em}}\mathrm{in}\text{\hspace{0.17em}}\mathrm{the}\text{\hspace{0.17em}}\mathrm{same}\text{\hspace{0.17em}}\mathrm{region}\end{array}.$$We also analysed the error term when the classical Laplacian regularization was applied in the non-equal voxel case.

The term $-N{x}_{{N}_{\mathrm{min}}/2}$ is added into both sides of Eq. (4). Then, by substituting Eq. (5) into Eq. (4) and dividing both sides of Eq. (4) by $N$ (the total number of voxels), Eq. (8) can be derived:

## (8)

$$\frac{-{x}_{1}}{N}+\frac{-{x}_{Vi}}{N}+\text{\hspace{0.17em}}\dots +{x}_{{N}_{\mathrm{min}}/2}+\dots \frac{-{x}_{Vj}}{N}+\frac{-{x}_{{N}_{\mathrm{min}}}}{N}=\phantom{\rule{0ex}{0ex}}\frac{-({N}_{\mathrm{min}}-N){x}_{{N}_{\mathrm{min}}/2}}{N}+\frac{({V}_{i}/{V}_{\mathrm{min}}-1){x}_{Vi}}{N}+\dots \text{\hspace{0.17em}}\phantom{\rule{0ex}{0ex}}+\frac{({V}_{j}/{V}_{\mathrm{min}}-1){x}_{Vj}}{N}.$$The left side of Eq. (8) is the Laplacian regularization matrix proposed in the paper.^{30} Comparing Eq. (8) with Eq. (3), we can see the right side of Eq. (8) is not equal to zero. This result means constant noise caused by the right side of Eq. (8) will be added in each iteration of the conjugate gradient method, which will ultimately affect the accuracy of the reconstruction.

## 3.

## Validity of Micro-CT Guided NVMC in the Forward Problem

In this section, through simulation experiments, the validity of micro-CT guided NVMC in the forward problem will be proven. The advantage of this method for a curved surface boundary is compared with the classical MC based on equal-space voxels. Furthermore, besides the standard coarse voxel size $1\times 1\times 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$, a fine equal voxel size ($0.5\times 0.5\times 0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$) is also used.

## 3.1.

### Validity in Regular Plane Boundary

The experimental conditions of this simulation study are listed below, which is similar to Refs. 16, 20, and 34. For a clear comparison, two kinds of media are placed in a $60\times 60\times 60\text{-}{\mathrm{mm}}^{3}$ cubic domain. The source is placed in the position of [30, 30, 0] mm with the incident direction along the $z$ direction. The optical coefficients of medium 1 are as follows: absorption coefficient: ${\mu}_{\alpha}=0.005\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$; scattering coefficient: ${\mu}_{s}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\phantom{\rule{0ex}{0ex}}{\mathrm{mm}}^{-1}$; refractive index: $n=1.37$; and anisotropy: $g=0.01$. The optical coefficients of medium 2 are as follows: absorption coefficient: ${\mu}_{\alpha}=0.002\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$; scattering coefficient: ${\mu}_{s}=\phantom{\rule{0ex}{0ex}}5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$; refractive index: $n=1.37$; and anisotropy: $g=0.01$.

The cubic domain is firstly discretized into $60\times 60\times 60$ voxels. Suppose with the structural information of micro-CT, we know medium 1 is at the domain $(x\in [20,40],y\in \phantom{\rule{0ex}{0ex}}[20,40],z\in [0,10]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm})$, which is the region of interest, and the boundaries between medium 1 and 2 are plane boundaries. The voxels in medium 1 are refined into $0.5\times 0.5\times 0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$. A quantity of $1.3\times {10}^{7}$ photons is simulated. MCX is used for coarse equal voxel simulation, and NVMC is used for non-equal voxel simulation and fine equal voxel simulation, for which the medium 1 and medium 2 are both refined into $0.5\times 0.5\times 0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$. The contour plot of the fluence (with 10-db spacing) in the middle slice is shown in Fig. 3.

In the case of a regular plane boundary, the loss of high-resolution structural information caused by voxel binning is rare for both equal-space voxels and non-equal space voxels. The results are shown in Fig. 3(a) to 3(d). For that reason, the accuracy of both MCX with coarse equal-voxel and NVMC is almost the same. On the other hand, these results prove the validity of NVMC, although it is as accurate as equal-voxel in the plane surface boundary.

## 3.2.

### Validity in the Curved Surface Boundary

In this section, we will prove that the loss of structural information caused by the voxel binning in the curved surface boundary will affect the accuracy of MC in the forward problem. Comparing NVMC with the result of coarse and fine equal-space voxels MC, the advantage of NVMC for the curved surface boundary will be shown. The experimental conditions have been listed below:

Suppose with the structural information of micro-CT, we know that a sphere is located in a $60\times 60\times 60\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ cubic domain, the centre of the sphere is at [30, 30, 30] mm, and the radial of the sphere is 10 mm. The source is placed at position [30, 30, 0] mm with the incident direction along the $z$-axis. The optical coefficients of the cubic domain are as follows: absorption coefficient, ${\mu}_{\alpha}=0.002\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$; scattering coefficient, ${\mu}_{s}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$; refractive index, $n=1.37$; and anisotropy, $g=0$. The optical coefficients of the spherical domain are as follows: absorption coefficient, ${\mu}_{\alpha}=0.05\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$; scattering coefficient, ${\mu}_{s}=5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$; refractive index, $n=1.37$; and anisotropy, $g=0$. The cubic domain is first discretised into voxels of $1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$. Then, voxels at the domain $(x\in [15,45],\phantom{\rule{0ex}{0ex}}y\in [15,45],z\in [15,45])\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, which is the region of interest, are refined into $0.25\times 0.25\times 0.25\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ voxels. A quantity of $3\times {10}^{7}$ photons is simulated. The contour plots of the continuous wave (CW) fluence (with 10-db spacing) in the plane $\mathrm{z}=30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ are shown in Fig. 4. Besides the three methods used in Sec. 3.1, analytical solution from the diffusion model is also added.^{36}^{,}^{37} As the boundary of a cube for MC is finite, the analytical solution for a sphere in a semi-infinite medium is only plotted near the source and inside the sphere.

Because the geometry of a boundary is mostly irregular, constituted by different curved surface boundaries in small animals, this simulation experiment is practiced to prove the advantage of NVMC for a curved surface boundary. From Fig. 4(a) and 4(b), we can find the boundary definition for coarse equal-space voxels is obviously inaccurate. This inaccuracy is caused by the loss of structural information during voxel binning, which will directly affect the calculation accuracy in the forward problem and ultimately affect the reconstruction accuracy. Whereas, with the micro-CT guided non-equal voxel discretization, the boundary definition is more accurate. The influence data along the source direction and inside sphere domain ($x=30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, $y\in [20,40]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, $z=30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$) is gotten and then compared with the fine equal voxel MC, the mean square error (MSE) of NVMC and coarse equal voxel MC are calculated. The coarse equal voxel MC and NVMC are 0.5442 and 0.3297, respectively. Therefore, the accuracy of the NVMC in the forward problem is higher than for MC based on coarse equal-space voxels. This result can be seen in Fig. 4.

## 4.

## Validity of the Modified Laplacian Regularization Method in the Inverse Problem

In this section, we will prove the validity of the modified Laplacian regularization method for non-equal voxels. The experimental conditions are listed below:

One glass tube (2.2 mm in diameter) with 0.8-*μ*mol/L DiR fluorescence dye^{38} is placed in a rectangular glass box ($17\times 44\times 60\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$, 2-mm thick glass) with a mixed solution of Intralipid and Indian ink.

The major optical coefficients of the phantom are as the follows: absorption coefficient, $0.07\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$; reduced scattering coefficient, $1.8\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. These optical properties refer to the optical coefficients for the skin of rabbits.^{16}

The dual-modality imaging system combined with micro-CT and FMT^{13}^{,}^{39} was used to obtain the anatomic information and fluorescence images. The micro-CT subsystem consists of an x-ray tube (Ultrabright, Oxford Instruments, USA) and a flat panel x-ray detector (PaxScan2520V, Varian MedicalSystems, USA). The FMT subsystem contains a diode laser (BWF1-750, B&W TEK, USA), a laser scanner and a CCD camera (DU-897, Andor, UK). An F-theta lens after the scanner is used to focus the laser. Emission filter is located before the CCD camera. A rotation stage (ADRS-100, Aerotech, USA) which is used to hold the imaging object is fixed at the center of the system. A thresholding method was applied on the original micro-CT slice. Three regions were segmented (synthetic glass box, glass tube with fluorescence dye, mixed solution of Intralipid and Indian ink); the segmented information was used to form the L matrix shown in Eq. (7). Through pixel combination, the original micro-CT images with a pixel size of $75\times 75\times 75$ *μ*m^{3} were merged into slices with a $1\times 1\times 1\text{-}{\mathrm{mm}}^{3}$ pixel interval; then, the voxel in the region of the glass tube was refined to a smaller voxel ($0.5\times 0.5\times 0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$). All of the non-equal voxel information was stored in a volume file for the calculation of NVMC. The calculation of the Green function with NVMC was processed on a GPU cluster;^{28} the whole process took only 13 min on the GPU cluster. The experimental data are treated with normalized Born ratio.^{40}

The reconstruction region ($17\times 44\times 20\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$) was smaller than the real size of the phantom in the $z$-axis because the fluorescence signal only appeared in that area. In these experiments, there were 154 sources and 290 detectors; therefore, 444 Green functions were calculated. Therefore, the amount of data from the experiment ($154\times 290=44660$) was larger than the number of reconstruction fluorescent yields, which was equal to the number of voxels.

With the same experimental data, two different Laplacian regularization methods were compared with a relaxation factor of $\lambda =0.001$ decided by the $L$-curve method. The variation trend of the reconstructed concentration of this fluorescence tube was observed along different iterations.

From Fig. 5, with the increment of the iteration, we can see that the reconstructed concentration with the classical Laplacian regularization method increases first, but the increment trend is restrained after the tenth iteration. Finally, the reconstructed concentration deviates far from the real value. However, the result of the modified Laplacian regularization method steps closely to the real value. This result confirms the theoretical analysis of Eq. (8).

## 5.

## Experimental Results of Micro-CT Guided NVMC for FMT Reconstruction

In this section, through a phantom experiment, we will show the advantage of micro-CT guided NVMC for a dual-modality micro-CT/FMT system. In addition, small animal study will further validate this method with irregular boundaries and complex distributions of optical properties.

## 5.1.

### Phantom Results

The experimental conditions are the same as that mentioned in Sec. 4. The difference is that four glass tubes (2.2 mm in diameter) with different concentrations of DiR fluorescence dye (from 1.8 to 1.2 *μ*mol/L) were placed in a rectangular glass box ($17\times 44\times 60\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$, 2-mm thick glass). To investigate the advantage of NVMC for the FMT reconstruction, coarse equal-space voxel ($1\times 1\times 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$) and fine equal-space voxel ($0.5\times 0.5\times 0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$) reconstruction based on MC^{28} was compared with our NVMC approach. The voxels for NVMC were non-equal voxels ($0.5\times 0.5\times 0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ for the region of glass tube and $1\times 1\times 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ for the other regions).

The hardware for reconstruction can be found in our paper.^{28} The reconstructed slice and the concentration of the fluorochromes are shown in Figs. 6 and 7, and the classical Laplacian regularization method was applied in the reconstruction based on MC for equal-space voxels. The relaxation factor was $\lambda =0.001$, which was determined by the $L$-curve method, and the conjugate gradient method was stopped at the 18th iteration.

According to Fig. 6(b) to 6(d), the average value of the reconstructed fluorescent yield of the four fluorochromes was calculated. The reconstructed concentrations of four fluorochromes can be calculated according to the linear relationship between the reconstructed fluorescent yield and the concentration of fluorochromes. The results are shown in Fig. 7.

The influence of the structural information loss on the inverse problem can be found in Fig. 6. Because of the coarse equal-voxel binning, the original voxels of micro-CT with sizes of $75\times 75\times 75$ *μ*m^{3} were merged into voxels with sizes of $1\times 1\times 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$. Since the diameter of the fluorescence tube was only 2 mm, part of the structural information of the circular boundary was lost. Therefore, the circular region of this fluorescence tube was reconstructed into a rectangle. Although we can create fine equal-voxel binning for the original voxels of micro-CT to reduce the loss of structural information and increase the accuracy of the reconstruction (reconstructed result has been shown in Fig. 6(c), from Table 1, we can see that the accuracy of fine equal-voxel binning is at the cost of high occupancy of the memory and long-time consumption for the reconstruction. To ensure the accuracy of the reconstruction and reduce the consumption of time of the reconstruction, we propose a micro-CT-guided NVMC, which can generate fine voxel binning in an irregular boundary or region of interest, while creating coarse binning in homogenous regions. Therefore, the loss of structural information in an irregular boundary can be reduced, and the accuracy of calculation in the boundary can be increased, while the size of the Jacobian matrix increases only slightly, which can reduce the time consumption of the reconstruction.

## Table 1

Comparison of three kinds of voxel configurations for FMT reconstruction.

Number of Voxels | Size of Jacobian matrix | Memory for Jacobian (GB)* | Reconstruction time (min) | |
---|---|---|---|---|

Coarse equal voxel | 14,960 | 14,960×44,660 | 4.98 | 10 |

Fine equal voxel | 119,680 | 119,680×44,660 | 39.8 | 30 |

Non-equal voxel | 24,380 | 24,380×44,660 | 8.11 | 13 |

## *

double float

## 5.2.

### Small Animal Experimental Results

A 44-g mouse was anaesthetized by injecting excessive chloral hydrate in the abdominal cavity. A glass tube (2 mm in diameter) with 1 DiR fluorescence dye was inserted into the thorax from the neck with a little angle as shown in Fig. 8(a), and the depth from the surface was about 3 mm.

A dual-modality FMT/micro-CT system was used to image the mouse on the rotating stage. Meanwhile, combined with the co-calibration method^{39} for FMT/micro-CT system, the data was acquired from multiple angle projections. The reconstructed micro-CT slices were segmented into four regions (air, soft tissue, bone, DiR fluorescence dye). The $L$$L$ matrix shown in Eq. (7) was calculated according to this segmented information. Through pixel combination, the original micro-CT slices with a 75 *μ*m pixel interval were merged into slices with a 1-mm pixel interval. Because the glass tube was in the subcutaneous tissue, the voxel of the soft tissue near the glass tube was refined into a $0.25\times 0.25\times 1\text{-}{\mathrm{mm}}^{3}$ voxel. All of this information was stored in a volume file for the NVMC simulation.

As shown in Fig. 8(c) and 8(f) and Table 2, for the coarse equal-voxel MC ($1\times 1\times 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$) based reconstruction method, the shape and position of reconstructed fluorescence yields deviate their true value, which is related to the relative position of the yield and the grid [Fig. 1(a)] of image segmentation. But the NVMC based reconstruction method exactly reconstructs the outlines of the fluorescence yields and the positions. The fine equal-voxel MC method is also employed, but the reconstruction is failed because of the sharp increased number of voxels (several times as many as that of coarse voxels).

## Table 2

Coordinate value of tube center in micro-CT, reconstruction based coarse equal-voxel MC and NVMC.

Micro-CT | Coarse equal-voxel MC | NVMC | ||
---|---|---|---|---|

Slice1 | X/mm | 16.9 | 16.8 | 16.8 |

Y/mm | 18.4 | 18.5 | 18.5 | |

Slice 2 | X/mm | 15.0 | 14.6 | 15.0 |

Y/mm | 17.2 | 17.2 | 17.2 |

The *in vivo* small animal experiment demonstrates that micro-CT guided NVMC for a dual-modality micro-CT/FMT system is suitable for small animals with irregular boundaries and complex distributions of optical properties.

## 6.

## Conclusion

In conclusion, we proposed a micro-CT guided NVMC in the forward problem and a modified Laplacian regularization method in the inverse problem. Simulations and phantom experiments have shown that micro-CT guided NVMC can efficiently reduce the loss of structural information in the irregular boundary; therefore, the accuracy of the calculation can be increased in both the forward problem and inverse problem, while the size of the Jacobian matrix and time consumption are greatly reduced. The phantom experiment also demonstrated that the modified Laplacian regularization method correctly incorporates the structural information of non-equal space voxels into the reconstruction. Lastly, the small animal experiment showed the validity of the method for irregular boundaries and complex distributions of optical properties.

After the theoretical analysis and the investigation of the experiment, the validity of micro-CT guided NVMC for both the forward problem and inverse problem is demonstrated in Figs. 3 to 4 and 6 to 7, respectively. Since NVMC is proposed in the forward problem, in order to correctly incorporate the structural information of non-equal voxels into the reconstruction, we modified the classical Laplacian regularization method. The theoretical derivation was shown in Sec. 2.3. An inference was derived when the classical Laplacian regularisation method was applied in the case of non-equal space voxels, and an error term may be incorporated into the reconstruction, which finally restrained the reconstruction and affected the accuracy of the reconstruction.

## Acknowledgments

This work was supported by the National Key Technology R&D Program of China (Grant No. 2012BAI23B02), the Science Fund for Creative Research Group of China (Grant No. 61121004), the National Natural Science Foundation of China (Grant No. 61078072), and Specific International Scientific Cooperation (Grant No. 2010DFR30820).

## References

A. X. CongG. Wang, “A finite-element-based reconstrution method for 3D fluorescence tomography,” Opt. Express 13(24), 9847–9857 (2005).OPEXFF1094-4087http://dx.doi.org/10.1364/OPEX.13.009847Google Scholar

J. ChenV. VenugopalX. Intes, “Monte Carlo based method for fluorescence tomographic imaging with lifetime multiplexing using time gates,” Biomed. Opt. Express 2(4), 871–886 (2011).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.2.000871Google Scholar

F. K. Swirskiet al., “Identification of splenic reservoir monocytes and their deployment to inflammatory sites,” Science 325(5940), 612–616 (2009).SCIEAS0036-8075http://dx.doi.org/10.1126/science.1175202Google Scholar

V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed Eng. 8, 1–33 (2006).ARBEF71523-9829http://dx.doi.org/10.1146/annurev.bioeng.8.061505.095831Google Scholar

X. Montetet al., “Tomographic fluorescence imaging of tumor vascular volume in mice,” Radiol. 242(3), 751–758 (2007).RADLAX0033-8419http://dx.doi.org/10.1148/radiol.2423052065Google Scholar

G. Zacharakiset al., “Volumetric tomography of fluorescent proteins through small animals in vivo,” Proc. Natl. Acad. Sci. A 102(51), 18252–18257 (2005).PNASA60027-8424http://dx.doi.org/10.1073/pnas.0504628102Google Scholar

S. V. Patwardhanet al., “Time-dependent whole-body fluorescence tomography of probe bio-distributions in mice,” Opt. Express 13(7), 2564–2577 (2005).OPEXFF1094-4087http://dx.doi.org/10.1364/OPEX.13.002564Google Scholar

T. Lasseret al., “Surface reconstruction for free-space 360 degrees fluorescence molecular tomography and the effects of animal motion,” IEEE Trans. Med. Imag. 27(2), 188–194 (2008).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2007.904662Google Scholar

V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Meth. 7(8), 603–614 (2010).1548-7091http://dx.doi.org/10.1038/nmeth.1483Google Scholar

R. Holtet al., “A hybrid approach combining microCT and fluorescence tomography: imaging workflow and system of coordinate registration,” Proc. SPIE 7892, 789213 (2011).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.881449Google Scholar

D. Kepshireet al., “A microcomputed tomography guided fluorescence tomography system for small animal molecular imaging,” Rev. Sci. Instrum. 80(4), 043701 (2009).RSINAK0034-6748http://dx.doi.org/10.1063/1.3109903Google Scholar

Y. Linet al., “Quantitative fluorescence tomography using a combined tri-modality FT/DOT/XCT system,” Opt. Express 18(8), 7835–7850 (2010).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.18.007835Google Scholar

X. Yanget al., “Combined system of fluorescence diffuse optical tomography and microcomputed tomography for small animal imaging,” Rev. Sci. Instrum. 81 (5), 054304 (2010).RSINAK0034-6748http://dx.doi.org/10.1063/1.3422252Google Scholar

F. Stukeret al., “Hybrid small animal imaging system combining magnetic resonance imaging with fluorescence tomography using single photon avalanche diode detectors,” IEEE Trans. Med. Imag. 30(6), 1265–1273 (2011).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2011.2112669Google Scholar

R. B. Schulzet al., “Hybrid system for simultaneous fluorescence and x-ray computed tomography,” IEEE Trans. Med. Imag. 29(2), 465–473 (2010).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2009.2035310Google Scholar

J. F. Beeket al., “In vitro double-integrating-sphere optical properties of tissues between 630 and 1064 nm,” Phys. Med. Biomed. 42(11), 2255–2261 (1997).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/42/11/017Google Scholar

Q. FangD. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.17.020178Google Scholar

A. Joshiet al., “Radiative transport-based frequency-domain fluorescence tomography,” Phys. Med. Biol. 53(8), 2069–2088 (2008).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/53/8/005Google Scholar

H. K. KimJ. H. LeeA. H. Hielscher, “PDE-constrained fluorescence tomography with the frequency-domain equation of radiative transfer,” IEEE J. Sel. Top. Quant. 16(4), 793–803 (2010).IJSQEN1077-260Xhttp://dx.doi.org/10.1109/JSTQE.2009.2038112Google Scholar

L. WangS. L. JacquesL. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Biomed. 47(2), 131–146 (1995).CMPBEK0169-2607http://dx.doi.org/10.1016/0169-2607(95)01640-FGoogle Scholar

D. A. Boaset al., “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10(3), 159–170 (2002).OPEXFF1094-4087Google Scholar

T. LiH. GongQ. Luo, “MCVM: Monte Carlo modeling of photon migration in voxelized media,” J. Innov. Opt. Health Sci. 3(2), 1–12 (2010).http://dx.doi.org/10.1142/S1793545810000927Google Scholar

J. ChenX. Intes, “Comparison of Monte Carlo methods for fluorescence molecular tomography-computational efficiency,” Med. Phys. 38(10), 5788–5798 (2011).MPHYA60094-2405http://dx.doi.org/10.1118/1.3641827Google Scholar

E. AlerstamT. SvenssonS. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. 13(6), 060504 (2008).JBOPFO1083-3668http://dx.doi.org/10.1117/1.3041496Google Scholar

W. C. Y. Loet al., “GPU-accelerated Monte Carlo simulation for photodynamic therapy treatment planning,” Proc. SPIE 7373, 737313(2009).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.831944Google Scholar

A. T. N. Kumaret al., “A time domain fluorescence tomography system for small animal Imaging,” IEEE Trans. Med. Imag. 27(8), 1152–1163 (2008).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2008.918341Google Scholar

X. Zhanget al., “Development of a noncontact 3-D fluorescence tomography system for small animal in vivo imaging,” Proc. SPIE 7191, 71910D (2009).SPIECJ0361-0748http://dx.doi.org/10.1117/12.808199Google Scholar

G. Quanet al., “Monte Carlo-based fluorescence molecular tomography reconstruction method accelerated by a cluster of graphic processing units,” J. Biomed. Opt. 16(2), 026018 (2011).JBOPFO1083-3668http://dx.doi.org/10.1117/1.3544548Google Scholar

E. E. Graveset al., “A submillimeter resolution fluorescence molecular imaging system for small animal imaging,” Med. Phys. 30(5), 901–911 (2003).MPHYA60094-2405http://dx.doi.org/10.1118/1.1568977Google Scholar

P. K. Yalavarthyet al., “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express 15(13), 8043–8058 (2007).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.15.008043Google Scholar

D. Hydeet al., “Data specific spatially varying regularization for multimodal fluorescence molecular tomography,” IEEE Trans. Med. Imag. 29(2), 365–374 (2010).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2009.2031112Google Scholar

S. C. Daviset al., “Image-guided diffuse optical fluorescence tomography implemented with Laplacian-type regularization,” Opt. Express 15(7), 4066–4082 (2007).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.15.004066Google Scholar

J. Canny, “A computational approach to edge detection,” IEEE Trans. Pattern Anal. Mach. Intell. PAMI-8(6), 679–698 (1986).ITPIDJ0162-8828http://dx.doi.org/10.1109/TPAMI.1986.4767851Google Scholar

D. Calvettiet al., “Tikhonov regularization and the L-curve for large discrete ill-posed problems,” J. Comput. Appl. Math. 123(1–2), 423–446 (2000).JCAMDI0377-0427http://dx.doi.org/10.1016/S0377-0427(00)00414-3Google Scholar

T. Rubaeket al., “Nonlinear microwave imaging for breast-cancer screening using Gauss-Newton’s method and the CGLS inversion algorithm,” IEEE T. Antenn. Propag. 55(8), 2320–2331 (2007).IETPAK0018-926Xhttp://dx.doi.org/10.1109/TAP.2007.901993Google Scholar

D. A. Boaset al., “Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications,” Proc. Natl. Acad. Sci. U. S. A. 91(11), 4887–4891 (1994).PNASA60027-8424http://dx.doi.org/10.1073/pnas.91.11.4887Google Scholar

Q. Fang, “Mesh-based Monte Carlo method using fast raytracing in Plücker coordinates,” Bio. Opt. Express 1(1), 165–175 (2010).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.1.000165Google Scholar

Z. Zhanget al., “Biomimetic nanocarrier for direct cytosolic drug delivery,” Angew. Chem. Int. Ed. 48(48), 9171–9175 (2009).ACIEF51433-7851http://dx.doi.org/10.1002/anie.v48:48Google Scholar

J. W. Fuet al., “A generic, geometric cocalibration method for a combined system of fluorescence molecular tomography and microcomputed tomography with arbitrarily shaped objects,” Med. Phys. 38(12), 6561–6570 (2011).MPHYA60094-2405http://dx.doi.org/10.1118/1.3658727Google Scholar

A. SoubretJ. RipollV. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities: study of the normalized Born ratio,” IEEE Trans. Med. Imag. 24(10), 1377–1386 (2005).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2005.857213Google Scholar