Micro-computed tomography-guided, non-equal voxel Monte Carlo method for reconstruction of fluorescence molecular tomography

Guotao Quan; Kan Wang; Xiaoquan Yang; Yong Deng; Qingming Luo; Hui Gong

doi:10.1117/1.JBO.17.8.086006

8 August 2012 Micro-computed tomography-guided, non-equal voxel Monte Carlo method for reconstruction of fluorescence molecular tomography

Guotao Quan, Kan Wang, Xiaoquan Yang, Yong Deng, Qingming Luo, Hui Gong

Author Affiliations +

Journal of Biomedical Optics, Vol. 17, Issue 8, 086006 (August 2012). https://doi.org/10.1117/1.JBO.17.8.086006

Abstract

The study of dual-modality technology which combines microcomputed tomography (micro-CT) and fluorescence molecular tomography (FMT) has become one of the main focuses in FMT. However, because of the diversity of the optical properties and irregular geometry for small animals, a reconstruction method that can effectively utilize the high-resolution structural information of micro-CT for tissue with arbitrary optical properties is still one of the most challenging problems in FMT. We develop a micro-CT-guided non-equal voxel Monte Carlo method for FMT reconstruction. With the guidance of micro-CT, precise voxel binning can be conducted on the irregular boundary or region of interest. A modified Laplacian regularization method is also proposed to accurately reconstruct the distribution of the fluorescent yield for non-equal space voxels. Simulations and phantom experiments show that this method not only effectively reduces the loss of high-resolution structural information of micro-CT in irregular boundaries and increases the accuracy of the FMT algorithm in both forward and inverse problems, but the method also has a small Jacobian matrix and a short reconstruction time. At last, we performed small animal imaging to validate our method.

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¹^,² 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,³ pulmonary, oncology, inflammation, and infectious disease.⁴^–⁶

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⁷ or a match solution of optical properties⁸ 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.⁹ 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.¹⁰^–¹⁵

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,¹⁶ the cerebral spinal fluid layer in the brain, and cysts in the human breast.¹⁷ Therefore, those method based on diffusion approximation with the finite element method may lead to error¹⁸^,¹⁹ 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.¹⁸ 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,²⁰^–²² and it is suitable for small animals with arbitrary optical properties. However, the huge time consumption of MC greatly limits this method’s application.²³ Therefore, a graphics processing unit (GPU) is introduced to MC to accelerate its calculation speed.¹⁷^,²⁴^,²⁵ Different methods based on MC accelerated by GPU have been introduced in the FMT reconstruction.²⁶^–²⁸

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,²⁹ 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.³⁰^,³¹ 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).¹⁷ 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,³² 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.¹³ 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³ are coarsely merged in to slices with a pixel pitch of $1 \times 1 \times 1 {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 {mm}^{3}$ .

Fig. 1

One slice of the discretization of non-equal voxels. (a) Coarse voxel binning. (b) Non-equal voxel binning. The shadow area in both (a) and (b) is region one; outside of the shadow area is region two; the curved line of the shadow area is the boundary; and the boundary crosses over three voxels in (a).

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³³ 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 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.¹⁷^,²⁰^,²² 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.

Fig. 2

Flowchart of non-equal voxel Monte Carlo (NVMC) method; shadow box indicates the different between NVMC and MCX.

(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¹⁷ and Boas,²¹ 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)²⁸ and obtain the Jacobian matrix.

Eq. (1)

J (R_{d}, R_{s}) = \frac{\partial [ϕ_{f} (R_{d}, R_{s})]}{\partial x} = \int g^{λ e m} (r, R_{d}) g^{λ e x} (r, R_{s}) d r,

where

J (R_{d}, R_{s})

is the Jacobian matrix with source

R_{s}

and detector

R_{d}

;

ϕ_{f} (R_{d}, R_{s})

is the flux with source

R_{s}

and detector

R_{d}

;

g^{λ e m} (r, R_{d})

is the Green function with wavelength

λ_{e m}

, and

g^{λ e x} (r, R_{s})

is the Green function with wavelength

λ_{e x}

.

The Laplacian regularization method is used to reconstruct the distribution of the fluorescence yield; the equation is shown in Eq. (2):

Eq. (2)

\min ({‖ J x - y ‖}^{2} + λ ‖ L x ‖),

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 = η μ_{a f}

;

μ_{a f}

is the absorption coefficient of the fluorochromes;

η

is the quantum efficiency of the fluorochromes;

λ

is the relaxation factor determined by the L-curve method;³⁴ the conjugate gradient method is used to solve Eq. (2).³⁵ 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.³⁰ It can be expressed as Eq. (3):

Eq. (3)

\nabla^{2} x (r) h^{2} ≍ x_{1} + x_{2} + \dots - N x_{N / 2} + \dots + x_{N} = 0,

where

N

is the number of equal spaces, and

h

is the step size.

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_{\min}$ , and total number of voxels is $N_{\min}$ , where, $N_{\min} = \sum_{i = 1}^{N} \frac{V_{i}}{V_{\min}}$ and the step size is $h_{\min}$ , such that Eq. (3) can be expressed as follows:

Eq. (4)

\nabla^{2} x (r) h_{\min}^{2} ≍ x_{1} + x_{2} + \dots - N_{\min} x_{N_{\min} / 2} + \dots + x_{N_{\min}} = 0 .

When the volume of the voxel is small, the value of

x

for a bigger voxel with a volume of

V_{i}

can be expressed by the value of the minimum voxel. The relationship between them has been shown in Eq. (5).

Eq. (5)

x_{V i} \approx \frac{\sum_{j = 1}^{V_{i} / V_{\min}} x_{j}}{V_{i} / V_{\min}} .

Substituting Eq. (5) into Eq. (4) and dividing by $- N_{\min}$ on both sides of Eq. (4), Eq. (6) can be derived:

Eq. (6)

\frac{- x_{1}}{N_{\min}} + \frac{- (V_{i} / V_{\min}) x_{V i}}{N_{\min}} + \dots + x_{N_{\min} / 2} + \dots \frac{- (V_{j} / V_{\min}) x_{V j}}{N_{\min}} + \frac{- x_{N_{\min}}}{N_{\min}} = 0.

Therefore, the Laplacian regularization matrix can be expressed as follows:

Eq. (7)

L_{i j} = {\begin{cases} 1 & i and j is the same point \\ - \frac{V_{i} / V_{\min}}{N_{\min}} & i and j is in the same region \\ 0 & i and j is not in the same region \end{cases} .

We also analysed the error term when the classical Laplacian regularization was applied in the non-equal voxel case.

The term $- N x_{N_{\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:

Eq. (8)

\frac{- x_{1}}{N} + \frac{- x_{V i}}{N} + \dots + x_{N_{\min} / 2} + \dots \frac{- x_{V j}}{N} + \frac{- x_{N_{\min}}}{N} = \frac{- (N_{\min} - N) x_{N_{\min} / 2}}{N} + \frac{(V_{i} / V_{\min} - 1) x_{V i}}{N} + \dots + \frac{(V_{j} / V_{\min} - 1) x_{V j}}{N} .

The left side of Eq. (8) is the Laplacian regularization matrix proposed in the paper.³⁰ 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 {mm}^{3}$ , a fine equal voxel size ( $0.5 \times 0.5 \times 0.5 {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 - {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: $μ_{α} = 0.005 {mm}^{- 1}$ ; scattering coefficient: $μ_{s} = 1 {mm}^{- 1}$ ; refractive index: $n = 1.37$ ; and anisotropy: $g = 0.01$ . The optical coefficients of medium 2 are as follows: absorption coefficient: $μ_{α} = 0.002 {mm}^{- 1}$ ; scattering coefficient: $μ_{s} = 5 {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 [20, 40], z \in [0, 10] 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 {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 {mm}^{3}$ . The contour plot of the fluence (with 10-db spacing) in the middle slice is shown in Fig. 3.

Fig. 3

Validity of non-equal voxel Monte Carlo (NVMC) method in the plane boundary. (a)–(c) Detail with an enlarged scale of the rectangular region in (d)–(f). (d)–(f) Coarse equal, none-equal and fine equal voxel discretization of the voxels in the middle slice of the cube. (g) Contour plots of the fluence (with 10-db spacing) in the middle slice with coarse equal voxel, non-equal voxel, fine equal voxel.

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 {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, $μ_{α} = 0.002 {mm}^{- 1}$ ; scattering coefficient, $μ_{s} = 1 {mm}^{- 1}$ ; refractive index, $n = 1.37$ ; and anisotropy, $g = 0$ . The optical coefficients of the spherical domain are as follows: absorption coefficient, $μ_{α} = 0.05 {mm}^{- 1}$ ; scattering coefficient, $μ_{s} = 5 {mm}^{- 1}$ ; refractive index, $n = 1.37$ ; and anisotropy, $g = 0$ . The cubic domain is first discretised into voxels of $1 {mm}^{3}$ . Then, voxels at the domain $(x \in [15, 45], y \in [15, 45], z \in [15, 45]) mm$ , which is the region of interest, are refined into $0.25 \times 0.25 \times 0.25 {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 $z = 30 mm$ are shown in Fig. 4. Besides the three methods used in Sec. 3.1, analytical solution from the diffusion model is also added.³⁶^,³⁷ 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.

Fig. 4

Accuracy of non-equal voxel Monte Carlo (NVMC) method in a curved surface boundary. (a)–(c) Detail with the enlarged scale of the rectangular region in (d)–(f). (c) Coarse equal, none-equal and fine equal voxel discretization of the voxels in the middle slice of the cube. (g) Contour plots of the fluence (with 10-db spacing) in the middle slice with coarse equal voxel Monte Carlo (MC), non-equal voxel MC, fine equal voxel MC, and diffusion mode with analytical solution.

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 mm$ , $y \in [20, 40] mm$ , $z = 30 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³⁸ is placed in a rectangular glass box ( $17 \times 44 \times 60 {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 {mm}^{- 1}$ ; reduced scattering coefficient, $1.8 {mm}^{- 1}$ . These optical properties refer to the optical coefficients for the skin of rabbits.¹⁶

The dual-modality imaging system combined with micro-CT and FMT¹³^,³⁹ 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³ were merged into slices with a $1 \times 1 \times 1 - {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 {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;²⁸ the whole process took only 13 min on the GPU cluster. The experimental data are treated with normalized Born ratio.⁴⁰

The reconstruction region ( $17 \times 44 \times 20 {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 $λ = 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).

Fig. 5

Reconstruction of fluorescence concentration for non-equal voxels with classical and modified Laplacian regularization methods. (a) Reconstructed concentration changes with the iteration number. (b) Reconstructed middle slice at the 29th iteration with the modified Laplacian regularization method. (c) Reconstructed middle slice at the 29th iteration with the classical Laplacian regularization method.

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 {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 {mm}^{3}$ ) and fine equal-space voxel ( $0.5 \times 0.5 \times 0.5 {mm}^{3}$ ) reconstruction based on MC²⁸ was compared with our NVMC approach. The voxels for NVMC were non-equal voxels ( $0.5 \times 0.5 \times 0.5 {mm}^{3}$ for the region of glass tube and $1 \times 1 \times 1 {mm}^{3}$ for the other regions).

The hardware for reconstruction can be found in our paper.²⁸ 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 $λ = 0.001$ , which was determined by the $L$ -curve method, and the conjugate gradient method was stopped at the 18th iteration.

Fig. 6

Reconstructed slices. (a) Is the reconstructed slice from micro-CT. (b) Shows the reconstructed slice with Monte Carlo (MC) for coarse equal voxels. (c) Represents the reconstructed slice with MC for fine equal voxels. Panel d is the reconstructed slice with non-equal voxel Monte Carlo (NVMC) method.

Fig. 7

Reconstructed concentration of DiR fluorescence dye with respect to the real concentration.

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³ were merged into voxels with sizes of $1 \times 1 \times 1 {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 \times 44, 660$	4.98	10
Fine equal voxel	119,680	$119, 680 \times 44, 660$	39.8	30
Non-equal voxel	24,380	$24, 380 \times 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.

Fig. 8

Imaging result of the mouse. (a) Volume rendering image of the mouse from micro-CT. (b) Micro-CT reconstruction slice of slice 1, whose position is indicated in (a) with white line. (c) Fluorescence molecular tomography (FMT) reconstruction result based on the coarse equal-voxel MC ( $1 \times 1 \times 1 {mm}^{3}$ ) according to (b). (d) FMT reconstruction result based NVMC according to (b). (e) Micro-CT reconstruction slice of slice 2, whose position is indicated in (a) with black line. (f) FMT reconstruction result based the coarse equal-voxel MC ( $1 \times 1 \times 1 {mm}^{3}$ ) according to (e). (g) FMT reconstruction result based NVMC according to (e). The FMT reconstruction values are in the same range from 0.1 to 1. The black and white circles in (c), (d), (f), and (g) are the outlines of glass tube and the mouse, respectively.

A dual-modality FMT/micro-CT system was used to image the mouse on the rotating stage. Meanwhile, combined with the co-calibration method³⁹ 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 - {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 {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
Slice1	$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

1.

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

2.

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). http://dx.doi.org/10.1364/BOE.2.000871 BOEICL 2156-7085 Google Scholar

3.

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

4.

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

5.

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

6.

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

7.

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

8.

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). http://dx.doi.org/10.1109/TMI.2007.904662 ITMID4 0278-0062 Google Scholar

9.

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

10.

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

11.

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

12.

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

13.

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

14.

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). http://dx.doi.org/10.1109/TMI.2011.2112669 ITMID4 0278-0062 Google Scholar

15.

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

16.

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). http://dx.doi.org/10.1088/0031-9155/42/11/017 PHMBA7 0031-9155 Google Scholar

17.

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). http://dx.doi.org/10.1364/OE.17.020178 OPEXFF 1094-4087 Google Scholar

18.

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

19.

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). http://dx.doi.org/10.1109/JSTQE.2009.2038112 IJSQEN 1077-260X Google Scholar

20.

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). http://dx.doi.org/10.1016/0169-2607(95)01640-F CMPBEK 0169-2607 Google Scholar

21.

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). OPEXFF 1094-4087 Google Scholar

22.

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/S1793545810000927 Google Scholar

23.

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

24.

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). http://dx.doi.org/10.1117/1.3041496 JBOPFO 1083-3668 Google Scholar

25.

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

26.

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

27.

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

28.

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). http://dx.doi.org/10.1117/1.3544548 JBOPFO 1083-3668 Google Scholar

29.

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

30.

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

31.

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

32.

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

33.

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

34.

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). http://dx.doi.org/10.1016/S0377-0427(00)00414-3 JCAMDI 0377-0427 Google Scholar

35.

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). http://dx.doi.org/10.1109/TAP.2007.901993 IETPAK 0018-926X Google Scholar

36.

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). http://dx.doi.org/10.1073/pnas.91.11.4887 PNASA6 0027-8424 Google Scholar

37.

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

38.

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

39.

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). http://dx.doi.org/10.1118/1.3658727 MPHYA6 0094-2405 Google Scholar

40.

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). http://dx.doi.org/10.1109/TMI.2005.857213 ITMID4 0278-0062 Google Scholar

Citation Download Citation

Guotao Quan, Kan Wang, Xiaoquan Yang, Yong Deng, Qingming Luo, and Hui Gong "Micro-computed tomography-guided, non-equal voxel Monte Carlo method for reconstruction of fluorescence molecular tomography," Journal of Biomedical Optics 17(8), 086006 (8 August 2012). https://doi.org/10.1117/1.JBO.17.8.086006

Published: 8 August 2012

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 5 scholarly publications and 1 patent.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Monte Carlo methods

Luminescence

Optical properties

Glasses

Fluorescence tomography

Tomography

Tissue optics

1.

Introduction

2.

Method

2.1.

Micro-CT Guided Non-Equal Discretization of Voxels

Fig. 1

2.2.

Monte Carlo Simulation for Non-Equal Voxels

Fig. 2

2.3.

Modified Laplacian Regularization for Non-Equal Voxels

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

3.

Validity of Micro-CT Guided NVMC in the Forward Problem

3.1.

Validity in Regular Plane Boundary

Fig. 3

3.2.

Validity in the Curved Surface Boundary

Fig. 4

4.

Validity of the Modified Laplacian Regularization Method in the Inverse Problem

Fig. 5

5.

Experimental Results of Micro-CT Guided NVMC for FMT Reconstruction

5.1.

Phantom Results

Fig. 6

Fig. 7

Table 1

5.2.

Small Animal Experimental Results

Fig. 8

Table 2

6.

Conclusion

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years