Mesh-based Monte Carlo method in time-domain widefield fluorescence molecular tomography

Jin Chen; Qianqian Fang; Xavier Intes

doi:10.1117/1.JBO.17.10.106009

1 October 2012 Mesh-based Monte Carlo method in time-domain widefield fluorescence molecular tomography

Jin Chen, Qianqian Fang, Xavier Intes

Author Affiliations +

Journal of Biomedical Optics, Vol. 17, Issue 10, 106009 (October 2012). https://doi.org/10.1117/1.JBO.17.10.106009

Abstract

We evaluated the potential of mesh-based Monte Carlo (MC) method for widefield time-gated fluorescence molecular tomography, aiming to improve accuracy in both shape discretization and photon transport modeling in preclinical settings. An optimized software platform was developed utilizing multithreading and distributed parallel computing to achieve efficient calculation. We validated the proposed algorithm and software by both simulations and in vivo studies. The results establish that the optimized mesh-based Monte Carlo (mMC) method is a computationally efficient solution for optical tomography studies in terms of both calculation time and memory utilization. The open source code, as part of a new release of mMC, is publicly available at http://mcx.sourceforge.net/mmc/.

1. Introduction

Optical imaging has been the subject of intense investigation over the past decades, largely due to the fact that light at visible and near-infrared wavelengths is a nonionizing, noninvasive probe with numerous applications in medicine.¹^,² In fluorescence molecular tomography (FMT), an emerging optical imaging technique for preclinical applications, one can solve for three-dimensional (3-D) in vivo fluorescent marker distribution maps using only noncontact surface optical measurements and tomographic image reconstructions.³^,⁴ This makes FMT a suitable method for small animal research where fluorophores are designed to label the drugs of interest, and enables the tracking of their delivery process. A successful tomographic reconstruction in FMT typically requires the forward modeling of photon transport inside optically complex tissue. The diffusion equation, an approximation to the more general radiative transport equation (RTE), is usually employed as the forward model for simplicity. However, this approximation becomes invalid when modeling low-scattering and highly absorbing tissue, resulting in potentially inaccurate quantification in small animal imaging where specimens exhibit wide variation in optical properties with possible presence of void regions (lung or bladder, for instance).⁵ Moreover, the diffusion equation cannot accurately represent the propagation of short light pulses in tissue for the portion of photons arriving at the surface early,⁶^,⁷ which has been shown as an effective datatype to improve the resolution in reconstruction when no a priori information is available.⁸^–¹⁰

The Monte Carlo (MC) method is an accurate light propagation modeling approach in dealing with general media such as those in small animals and early photons.¹¹^,¹² This method approximates the RTE solution via random sampling of large numbers of photons, thus it is often used as the gold standard for modeling time-resolved light propagation for all optical properties encountered in these scenarios. Nevertheless, MC-based techniques suffer from low computational efficiency (hours or days of computation)¹³ because numerous photon simulation trials are required to obtain satisfactory statistics. When time-domain data are considered, the photons are split into small time-bins, which leads to an increase in the number of photons to be simulated to achieve reasonable statistics per time window, and make calculations even more time-consuming. Only recently, our group developed a computationally efficient MC-based reconstruction technique for FMT utilizing time-gated data sets.¹⁴ This voxel-based Monte Carlo (vMC) method allows for computation of functional and fluorescent Jacobians for whole-body tomography within a couple of hours. However, difficulty arises when applying the vMC algorithm, as the commonly employed cubic shape voxels in a regular 3-D grid cannot accurately simulate the curved boundaries. One remedy is to raise the voxel density, but doing so drastically increases the memory and computational burden. This burden is further amplified when using widefield illumination strategies,⁹ where photons interact with a much larger boundary area compared to the conventional punctual light source scheme. Therefore, an MC algorithm that has flexibility in representing the arbitrary domain shape is highly desired.

In this work, we present an optimized software platform to solve the inverse problem in FMT with spatially extended arbitrary sources. The novelty of this approach lies in the accurate boundary representation and rapid calculation when employing the mesh-based Monte Carlo (mMC) method.¹⁵ Based on a currently available mMC code,¹⁶ an extended hybrid parallel version was implemented to model arbitrary illumination patterns. A comparison was performed between this method and our previously developed vMC method in terms of quantitative accuracy and computational efficiency for whole-body tomographic reconstruction in time-domain. This method was then validated and evaluated by small animal full-body simulations and experiments in tomographic settings.

2. Method

2.1.

Widefield Mesh-Based Monte Carlo Method

The mMC method utilizes fast ray-tracing techniques to accelerate calculation and is capable of simulating point illumination on complex geometry.¹⁵ Here the widefield illumination was developed based on the version 0.8 release of mMC as provided at Ref. 16. In order to simulate a realistic widefield source with spatial intensity variance, the illumination patterns were represented as a two-dimensional (2-D)-grid ( $1 {mm}^{2} pixel$ ) images with an intensity value assigned at each grid element [Fig. 1(a)]. The initial positions of the total photons were uniformly distributed over the illumination image by using a uniform random number generator. Any photon falling into one pixel had the pixel’s intensity value assigned as its initial packet weight, allowing for arbitrary illumination settings [Fig. 1(b)]. A simulated photon was then projected along the $z$ -axis and entered the mesh at an intersection point on the mesh surface. To identify the surface injection point, we applied a ray-tracing step to obtain its 3-D coordinates by testing all surface triangles using a Havel-Herout ray-tracing algorithm.¹⁷ In case multiple intersection points were identified [Fig. 1(b) and 1(c)], the intersection with the shortest distance to the initial position [Fig. 1(b)] was selected as the injection point on the surface. After a photon’s injection point was determined, its propagation in the media was simulated following the same rules as in classical MC techniques, until it exited the surface or the total time of interest was reached.

Fig. 1

The illustration for generating the widefield illumination using a digital mouse phantom. (a) arbitrary illumination pattern; (b) photons projected to the surface; (c) surface region of interest; a: initial position of a photon; b: the actual injection point of this photon; c: a potential injection point of this photon (discarded since the distance to the initial position is longer than b).

Unfortunately, a complete surface triangle test of this additional ray-tracing step imposed a significant computational overhead, particularly when the surface mesh was dense. For example, the calculation time for the ray-triangle intersection test could be $50 \times$ more than that for the photon propagation on a mouse model. In order to accelerate the computation time for the widefield method to match the punctual illumination method, a few algorithmic optimizations were implemented, particularly considering the fact that the $x - y$ coordinate of the injection position is dependent on the photon’s initial position, while the $z$ coordinate is dependent on the intersection of the photon and the surface. First, we confined the initial positions to the region of interest (ROI) and the surface triangles in the ROI were picked out [Fig. 1(c), the dark gray surface]. Second, using the projected triangles on the $x - y$ plane, a subset of triangles was selected and stored in memory. For any launched photon, a culling technique¹⁸ was then implemented to confine ray-triangle tests only to the triangles with their bounding boxes containing the photon’s $x - y$ coordinate. These operations reduced redundant calculations thus allowed for acceleration in computational time.

2.2.

Inverse Model for Reconstruction

The 3-D distribution of a fluorophore’s effective quantum yield $η (r)$ can be obtained by solving an integral equation relating the fluorescence signals at time $t$ and $η (r)$ :

Eq. (1)

U_{F} (r_{s}, r_{d}, t) = \int_{Ω} W (r_{s}, r_{d}, r, t) η (r) d r,

where

U_{F} (r_{s}, r_{d}, t)

is the fluorescence detected by a detector located at

r_{d}

at

t

resulting from an excitation at the source

r_{s}

at

t_{0} = 0

, the integration domain

Ω

is defined as the entire imaging volume, and

W (r_{s}, r_{d}, r, t)

, referred to as the weight function or Jacobian, describes how sensitive a change in

η (r)

will result in a change in

U_{F} (r_{s}, r_{d}, t)

. There are a few MC-based methods to calculate

W

, and based on the comparison between these different methods, we selected the forward-adjoint MC method to generate time-gated Jacobians. It allowed for computing Jacobian in an efficient manner,¹⁹ thanks to the smaller dataset acquired with widefield illumination strategies.⁹^,²⁰ In this method,

W

is computed by convolving the Green’s functions (the light propagation for impulse sources) and the lifetime decay¹⁴:

Eq. (2)

W (r_{s}, r_{d}, r, t) = \int_{0}^{t} e^{- (t - t^{'}) / τ} d t^{'} \int_{0}^{t^{'}} G^{x} (r_{s}, r, t^{'} - t^{''}) G^{m} (r, r_{d}, t^{''}) d t^{''},

where

G_{x}

and

G_{m}

are the time-dependent background Green’s functions calculated by mMC simulations at excitation and emission wavelength, respectively,²¹ and

τ

is the lifetime of the fluorophore. The system of linear equations representing the measurements detected at different positions and time can then be solved to obtain the fluorophore distribution.

2.3.

Computational Settings and Efficiency Evaluation

For forward simulations, we incorporated the widefield-related computation in the mMC software platform. In addition to the previously implemented multithreading and single-instruction multiple-data (SIMD) parallelism, we further enhanced the modeling efficiency by integrating Message Passing Interface (MPI)-based parallel computing technique with the code. This allows one to run mMC in a distributed memory system, such as a high-performance cluster. To utilize the maximum computational power of platforms with both multithreads and multinodes, an adaptive hybrid parallelization method using the communication protocol MPI and OpenMP APIs was implemented. The MC approach is highly parallelizable, that is, the large number of photons can be broken up into smaller sets and calculated independently. In terms of distribution, the MPI protocol had photons equally divided and distributed to all nodes, while the OpenMP protocols dynamically set the assignment of photons for the threads on a node. The seed for the random number generator was assigned based on the thread and node ID. A reduce operation was performed to sum up the result to a single node (master node) after all nodes finished computing.

In this particular work, the computational efficiency of the hybrid parallel code was tested using single/dual core CPUs (BlueGene/Opteron on CCNI, RPI). To estimate the scalability of multithreading and multi-CPU calculation, quantitative metrics were computed to measure the performance of the parallel programs. We define the speedup as the ratio of the execution time on a single processor (the sequential version) to that on a multinode cluster. This ratio is defined as

Eq. (3)

S (n) = \frac{t_{s}}{t_{p} (n)},

where

n

is the number of threads/processes used in the parallel simulation,

t_{s}

is the execution time on a single-thread processor, and

t_{p}

is the execution time on a parallel computer.

S (n)

therefore describes the scalability of the system as the number of threads/processes increases. The efficiency was computed as

Eq. (4)

E (n) = \frac{t_{s}}{n \times t_{p} (n)},

which describes the fraction of the time that is being used by the processors for the computation.

2.4.

Accuracy Comparison

Due to a wider illuminated surface, the simulation result of widefield illumination can be more sensitive to errors of boundary modeling than conventional point illumination used in optical tomography. To quantitatively assess the performance of mMC and vMC dealing with different types of illumination on complex models, a simulation using high-resolution voxel-based model was first set up as the reference. The light propagation profiles using both a mesh-based model and a resampled low resolution voxel-based model with the same number of tetrahedral elements/valid voxels (voxels that are not air) were then computed.

To mimic a preclinical tomographic reconstruction, the high-resolution Digimouse model ( $0.1 {mm}^{3}$ ) generated from CT data²² was employed to create the mesh model and low-resolution voxelized model (Fig. 2). The software package Iso2mesh²³ was used to convert the volumetric model to a mesh model containing 4075 node, 22,379 tetrahedron elements. The corresponding low-resolution voxel-based model had a comparable number of 22,433 valid voxels resulting in $1 {mm}^{3}$ voxels, which is the typical voxel size in preclinical applications. Notice here the mesh model was homogeneously tesselated with no denser sampling around the boundary or internal organs to assure an impartial comparison with the homogeneously voxelized model. The illumination pattern covered a 5- to 7-cm area with full coverage of the transverse plane. A cross-sectional plot is shown in Fig. 2(b) for mesh- and voxel-based model, respectively. In all simulations above, the numbers of photons were kept at $10^{8}$ according to Ref. 19. The optical properties were set to $μ_{a} = 0.3 {cm}^{- 1}$ , $μ_{s}^{'} = 15 {cm}^{- 1}$ , $g = 0.9$ , $n = 1.37$ , which are typical values for mouse tissue in the NIR spectral region in the literature.²⁴ The time profiles were recorded with a 300 ps gate width and 20 ps time shift between the gates to replicate the RPI imaging platform for optimal experimental performance.²⁵

Fig. 2

(a) The high-resolution mouse atlas employed in the accuracy comparison. (b) The selected slice in (a) for boundary visualization overlapped with the mesh model and the low resolution voxel model. The red line shows the boundary of the high-resolution CT image, while the black and blue lines represent that of the mesh- and low-resolution voxel model, respectively.

The time-gated profiles of light propagation using the mesh- and high-resolution models were scaled down to the size of the low-resolution profile, then the percentage error of mesh- and low-resolution models in relation to the high-resolution model was calculated to assess the accuracy quantitatively. Further evaluation of the methods were conducted by calculating the overall root mean square difference (RMS difference) between the results using the two methods and high-resolution standard:

Eq. (5)

RMS (t) = \sqrt{\frac{\sum_{i = 1}^{N} {[p (t) - p_{0} (t)]}^{2}}{\sum_{i = 1}^{N} p_{0} {(t)}^{2}}},

where

N

is the number of voxels in the low-resolution voxel model,

t

is the time,

p_{0}

is the high resolution result, and

p

is that calculated by low-resolution mMC or vMC.

2.5.

Tomographic Reconstruction Settings

2.5.1.

Simulation reconstructions

The above mentioned models were utilized for a full tomographic reconstruction study. The high-resolution mouse model was employed to generate fluorescence measurements with the pre-segmented kidneys labeled with simulated fluorescent markers (effective quantum yield 0.1; $τ = 0.8 ns$ ). 64 sliding half-space illumination patterns and 64 point detectors were simulated within the abdomen section (x: 50 to 70 mm, y: 8.5 to 28.5 mm, 400 pixels for each pattern) using 512 single-thread nodes on BlueGene. The detectors were evenly spanned over the region of interest with a 2.857 mm separation [Fig. 3(a)]. A 25% rising gate (early gate) was selected and the Born normalization²⁶ of the corresponding time point spread functions (TPSFs) was employed as the datatype for reconstruction herein. The conjugate gradient method was used to solve the inverse problem and the iteration was stopped when the relative error was less than 0.02.

Fig. 3

(a) Mesh-based mouse model with kidneys depicted in red. The patch of black and red shows one sample scanning pattern out of a total of 64. The black areas indicate where the photon source is off, while the red area shows where it is on. The black circles at the bottom represent the detectors. (b) The segmented boundary of CT-scanned images from the in vivo experiment overlaid with the half-space illumination patterns and normalized Born measurements (Video 1, MPEG, 1.61 MB) [URL: http://dx.doi.org/10.1117/1.JBO.17.10.106009.1]. The red box and the black circles show the entire illuminated area and the detector positions, respectively.

2.5.2.

Experimental validation

The performance of this method was also evaluated with experimental data for FMT. The data were acquired from an in vivo experiment using a time-resolved widefield tomography platform (the details of the system can be found in Ref. 25). The system employed a tunable femtosecond laser as the source and a time-gated ICCD camera as the detector. A pico-projector module was used for source pattern generation, allowing for rapid acquisition of spatially and temporally dense measurements. The experiment was performed on a freshly euthanized mouse with a 13-mm-long tube with 1 mm diameter placed in the thoracic cavity [Fig. 4(a)].²⁷ The tube contained 14 pmol of IRDye 800CW R(LiCOR) in 1 μL ethanol. 64 sliding bar-shaped illumination patterns were projected over a $33 \times 18 mm$ area (594 pixels to represent each pattern) on the abdomen with 97 point detectors ( $0.5 {mm}^{2}$ areas) measuring the transmitted excitation and emission photons. The detectors were empirically selected from the high-resolution CCD camera images over the area where fluorescence signals were detected, with denser detectors at the highly fluorescing region. We acquired a profile covering 4.6 ns with 40 ps offset between two successive gates. The animal was then subsequently scanned by Micro CT (Scanoco Viva CT40) and the CT images were registered with the optical platform to locate the fluorescent tube.

Fig. 4

(a) 3-D reconstruction showing the 50% isovolume of the reconstructed object (blue) and the fluorescent vessel (red). Axial and frontal slices at $z = 6.5 mm$ and $y = 21.5 mm$ for (b) node-based, (c) element-based, and (d) voxel-based reconstructions, respectively.

The surface geometry of the region of interest ( $40 \times 31 mm$ ) was extracted to generate a homogeneously tesselated model with 92,713 elements and 15,581 nodes for the mesh-based weight matrix calculation. The entire model was assigned the average background properties of the small animal acquired by spectroscopic fitting of the excitation measurements ( $μ_{a} = 0.3 {cm}^{- 1}$ and $μ_{s}^{'} = 25 {cm}^{- 1}$ ). The results of the forward simulations were stored in both node-based format (photon package weights are accumulated at each node) and element-based format (the result is saved for each element) to calculate the weight matrix. The number of nodes was only 16.8% of that of elements, thus the node-based format had advantages in both storing space and speed for data processing. The reconstructions employing weight matrices using the two formats were compared, to assess the possible loss of accuracy of using node-based format. The early-gate at 25% of the maximum value was selected in reconstruction, similar to the synthetic study. For comparison, a reconstruction employing the same measurement dataset on a voxel-based model with 92,831 non-air isotropic $0.5 {mm}^{3}$ voxels was performed.

3. Results

3.1.

Computational Efficiency

Our implementation allows the generation of spatially complex illumination patterns over arbitrary oriented surfaces to accurately model noncontact widefield illumination strategies. After the optimizations mentioned in Sec. 2.1, the overall computational overhead of simulating extended time-resolved sources in this complex scenario was only marginally larger (5% to 10%) over a point source configuration, while providing great flexibility to define complex sources.

Figure 5(a) shows the speedup of mMC calculation while using different threading settings with $10^{7}$ and $10^{8}$ photons. For simulations using a larger number of photons, the speedup curve has a marked improvement toward the ideal curve, implying higher efficiency ( $E (16) = 93.2 %$ for $10^{7}$ photons and $E (16) = 99.3 %$ for $10^{8}$ photons with single-thread computation setting). The difference between the three threading methods is minimal ( $< 3 %$ ) when the same number of photons is employed, with pure threading being the closest to ideal. However, the current generation of multicore hardware has a limitation on total number of threads and thus limits the speedup potential of threading-only technique.

Fig. 5

(a) The speedup using widefield mMC with $10^{7}$ and $18^{8}$ photons on single-thread nodes, four-thread nodes and an eight-thread node; and (b) the time cost for mMC and vMC on single-tread nodes.

Figure 5(b) depicts the time cost for mMC and vMC using different number of processes under single-thread parallel settings while having similar numbers of elements/valid voxels. Under all available processes settings (128 to 4096), mMC remains roughly four to five times faster than vMC, while for large number of processes (4096) the program approaches its scaling bottleneck.

3.2.

Accuracy of Time-Domain Jacobians

Figure 6(a) and 6(b) shows the forward photon propagation profiles at $y = 60 mm$ for the full field and point illumination, respectively. As expected, mMC fits the boundaries much better compared to vMC for widefield case, while both mMC and vMC matches well to the standard for point source case. Figure 6(c) presents the TPSFs for these two illumination settings at the same detector indicated in Fig. 6(a) and 6(b). Both the results using vMC and mMC methods fit the standard well for the point illumination, while the TPSFs for widefield illumination result in increased error for vMC compared to mMC, especially around the rising part (early gates) of the TPSF. A comparison of the different datatype extracted from the time-resolved measurements at these two detectors are listed in Table 1, with the continuous wave (CW) measurement being a percentage error and the other parameters as the absolute errors in time. The error in the first moment of TPSF (CW measurements) reflects the integrated intensity difference and that in the second and third moments (mean flight time and variance) represent a change in position and shape of the TPSF. Overall, the mMC results are more accurate than the vMC ones, while the inaccuracy at the side detector (detector 1) is greater than that at the center detector (detector 2), due to the fact that the side detector is closer to the edge, thus more affected by the mismatched boundary area. Comparing the widefield and point illumination, the error is greater in the widefield case as expected, as a result of an illumination pattern probing larger boundary area where mismatch boundaries exist.

Fig. 6

Forward simulation comparison for (a) widefield and (b) point source illumination. mMC: blue, vMC: gray; high-resolution standard (std): solid-color-filled contours (background). (c) The normalized time-resolved detector readings (black crosses and circles in (a) and (b). The solid and dashed lines represent the TPSFs in the point source and widefield simulations, respectively. In the bottom figure the dashed and solid lines are overlapped due to similar detector responses for simulations with the point source and the widefield source.

Table 1

Error table comparing the low-resolution vMC and mMC with the standard. The CW measurement error is in percentage, the others are in absolute time.

		CW(%)		Tmax (ps)		Mean flight time (ps)		Variance (ps)
		Det 1	Det 2	Det 1	Det 2	Det 1	Det 2	Det 1	Det 2
Widefield	mMC	$- 1.29$	$- 0.81$	40	20	20.81	15.94	$- 3.93$	$- 2.19$
Widefield	vMC	$- 3.35$	1.99	80	0	67.83	18.81	$- 8.51$	$- 3.00$
Point	mMC	$- 1.62$	$- 0.38$	0	20	$- 6.96$	20.93	1.55	$- 3.43$
Point	vMC	2.48	0.73	0	20	21.12	18.81	$- 3.30$	$- 3.00$

The contours of continuous wave Jacobians using both the mMC and vMC methods are shown in Fig. 7(a) and 7(b) for widefield and point illumination, respectively. In Fig. 7(c) and 7(d), the Jacobian contour comparison at the 25% rising gate is displayed. In both cases, contours using the mMC and vMC methods match the reference contour well for the point-point SD pair. However, the contour using the mMC method fits the reference much better than that using the vMC method in a widefield-point SD setting, especially around the boundaries where the voxel model experiences increased mismatch to the high-resolution model. An average RMS difference of 10.8% and 9.4% is achieved in the point-point SD pair Jacobians using vMC and mMC, respectively, while the RMS difference is 23.5% and 13.4% for the widefield-point SD pair. The values may be compared to a baseline RMS difference of 2.4% and 3.8%, obtained using two high-resolution simulations with $10^{9}$ photons and different random seeds. The Jacobian accuracy comparison in time-domain is shown in Fig. 7(c). For the point-point SD pair, the RMS differences using the two methods do not experience significant change with respect to the changes in gate position. However, for the widefield case, the mMC has a considerably lower RMS difference when compared to the vMC for early gates.

Fig. 7

The contour comparison of Jacobians using mMC and vMC for (a) and (c) widefield and (b) and (d) point illumination in (a) and (b) continuous wave and (c) and (d) time gated modes. The contours filled with solid colors are from the high-resolution simulation. (e) The RMS comparison in time-domain for four cases and all time gates.

3.3.

Simulation Reconstruction Performance

The reconstruction using mesh-based MC methods is shown in Fig. 8. The kidneys are accurately reconstructed in regard to their original discretization. The maximum reconstructed value for the mesh-based MC is 91.0% of the expected value. The centroids of the two kidneys have position error of 1.02 mm and 0.78 mm, respectively. The maximum reconstructed value for the mMC method is 3.4% more accurate toward the expected value compared to the vMC method. The run time for creating the time-domain mesh-based Jacobian was 3.47 h for $64 \times 64$ SD pairs (in total 4096) on 512 nodes.

Fig. 8

The simulation reconstruction result. (a) The 50% isovolume of the reconstruction. (b) 5 slices from $x = 50 mm$ to $x = 70 mm$ with a separation of 5 mm.

3.4.

Experimental Comparison

The 50% isovolume of the reconstructed effective quantum yield using the mMC node-based results is shown in Fig. 4(a). This isovolume has a maximum length of 11.3 mm and mean diameter of 2.7 mm. Moreover, sample slices of cross-sections overlaid on the corresponding CT slices demonstrates the accurate localization of the inclusion. The element-based (92,713 elements) and voxel-based (15,581 nodes) reconstructions are also shown for comparison in Fig. 4(b) and 4(c), respectively. The node-based and element-based weight matrices result in very close reconstructions in both quantification ( $< 2 %$ difference in maximum reconstructed value) and position of the inclusion. The centroid position error of the 50% isovolume using vMC is 2.83 mm while that using mMC is 0.56 mm.

4. Discussion

In this study, we successfully implemented the widefield illumination for calculating the propagation of light through arbitrary 3-D media in mesh-based models. After implementing an optimized ray-test step for determining the injecting point for each photon, computational time for widefield illumination simulations becomes only marginally longer (5% to 10%) than point-source simulations. This is a remarkable improvement compared to the 50 times increase when no optimization was employed in the widefield illumination.

Moreover, the parallelization of the code is a significant advancement of the computation capacity. The aim of developing a mixed mode MPI/OpenMP code was to utilize the full computation availability of single symmetric multiprocessors (SMPs) or a cluster of SMPs, and achieve a performance improvement over a pure MPI code. As expected, when the number of threads/processes is the same, the most efficient parallelization implementation is using pure multithreading (OpenMP). Although the number of threads on a single SMP is restricted at current time (usually less than eight), larger systems of single SMP could become available in the foreseeable future with 64-bit memory addressing. The efficiency comparison of hybrid code with the pure MPI implementation reveals that a performance improvement is achieved: the overall computational time has reduced and the scaling of the code with increasing thread number is better. Hence a combination of MPI and OpenMP parallelization paradigms provides more efficient parallelization strategy. In addition, the hybrid code makes full use of the memory on each node due to the shared memory mechanism, allowing for the possibility of running larger-scale programs. However, the improvement in time is quite small since the pure MPI parallelization algorithm for forward MC method is optimized and scales very well already, with no significant impairment due to load imbalance or memory limitation. Overall, the three parallelization implementations have excellent scalability with almost linear acceleration when a parallelized platform is available. Thus when more threads/processes are utilized, the faster the simulation will be.

The widefield implementation enables the ability to perform simulations with spatial variation in illumination intensity, which is essential for the emerging pattern light tomography applications.⁹^,²⁸^–³⁰ Compared to the conventional punctual illumination, the widefield light illumination results in more errors on the boundaries where the model mismatch exists due to the larger illumination area, in both the forward propagation and Jacobians. When the time-gated information is considered, the RMS difference comparison represents a significant improvement in accuracy for mMC at early gates, yet with an acceleration of $\sim 5 x$ compared to vMC. This is of importance as the early gates are not modeled accurately by diffusion equation while the incentive of using early gates to increase spatial resolution in FMT is rising. However, this accuracy improvement in Jacobian does not necessarily relate to the reconstruction performance (13.2% discrepancy in the RMS differences for Jacobian at the selected time gate, 3.4% difference in the maximum reconstructed value comparing the mMC and vMC reconstruction). This can be explained by the fact that the disparity mainly exists around the boundaries while the expected fluorescent objects are closer to the center.

In the current version of the code, the time-of-flight between the illumination plane and the object surface was ignored, however this can be easily corrected with a postprocessing. In addition, according to the light propagation velocity in air, the maximum error in time for the typical specimen thickness (around 2 cm) in preclinical applications is $\sim 3 ps$ , which is negligible compared to the minimum time interval of 20 ps for the time-resolved platform. Another potential difficulty for in vivo assessment may reside in the availability of an anatomical atlas, as high-resolution imaging modalities are required to obtain accurate surfaces for mesh-based reconstructions.

As in a few examples, this code can be used for investigating the fluorophore distribution with time-gated data type while serving as the forward solver for MRI-guided FMT. In this work, the mesh model for the forward and reconstruction processes was confined in homogeneous tissue and uniformly tessellated; however, it is noteworthy that it can also present in a highly heterogeneous medium with arbitrary internal/external boundaries while retaining a small number of nodes and elements, thanks to its flexibility. The representation of mesh-based models for high-resolution 3-D anatomic images can be further optimized (e.g., reducing the nodes while keeping the geometry unchanged) to benefit the accuracy of the reconstructions.

The rationale of similar results using a node- and element-based weight matrix is that these methods render essentially the same spatial profiles of photon propagation, though the number of nodes for both cases are much less than that of tetrahedron elements (15,581 nodes compared to 92,713 elements). Furthermore, in the ill-posed inverse problem, the effective information mainly relies on the measurements but not the image-space spatial distribution when the number of measurements (dependent on the number of SD pairs and time gates selected) is much less than that of the unknowns. Casting nodes as the unknowns in the inverse problem leads to a much smaller unsolved linear system,³¹ which is superior as the memory constraint is usually a burden in reconstructions, especially when the time-resolved data are employed to provide extra information. Due to this advantage in memory utilization, employing nodes instead of elements also allows for the use of denser source-detector pairs, which is crucial for improving resolution of the reconstructed objects.

In addition, the formulation for diffuse optical tomography (DOT) using the MC method is readily available for the adjoint method, thus the mMC method can be easily applied to DOT to quantitatively determine the functional parameters (blood volume and hemoglobin oxygenation) related to the optical properties as well.³² As FMT/DOT is becoming frequently used in combination with high-resolution imaging methods such as MRI, CT, x-ray, and ultrasound,³³ mMC can serve as an effective tool to incorporate high-resolution structural information into FMT/DOT with limited spatial resolution for quantitative functional and molecular studies.

5. Conclusion

This work focuses on assessing the feasibility and evaluating the computational performance of the widefield mMC algorithm for time-resolved FMT, and the improvements on the quality of models with complex boundaries. We demonstrate that the mMC method is a computationally efficient solution for tomographic use when modeling of general media is required (about five times faster than vMC). Simulations on a complex digital mouse phantom establish that both mMC and vMC match well with the high-resolution vMC output in temporal and spatial profiles for punctual SD settings; however, for widefield illumination, the accuracy for mMC is improved compared to that for vMC. Moreover, with the recent progress of MC simulation using GPU,³⁴^,³⁵ it is expected that a full reconstruction can be finished with a typical personal computer in the time frame comparable to the classic DE-based models. The current version of this code is available for download as MMC v1.0 at http://mcx.sourceforge.net/mmc/, and more features may be added, such as curved illumination patterns using a surface mesh.

Acknowledgments

This work was supported by National Institutes of Health (NIH) grant R21 EB013421 and by the National Science Foundation (NSF) grant CBET-1149407. Qianqian Fang wishes to thank the funding support from the Bill & Melinda Gates Foundation (#OPP1035992).

References

1.

V. Ntziachristoset al., “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol., 23 (3), 313 –320 (2005). http://dx.doi.org/10.1038/nbt1074 NABIF9 1087-0156 Google Scholar

2.

A. GibsonJ. HebdenS. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol., 50 (4), R1 –R43 (2005). http://dx.doi.org/10.1088/0031-9155/50/4/R01 PHMBA7 0031-9155 Google Scholar

3.

J. R. Mansfield, “Distinguished photons: a review of in vivo spectral fluorescence imaging in small animals,” Curr. Pharm. Biotechnol., 11 (6), 628 –638 (2010). http://dx.doi.org/10.2174/138920110792246474 CPBUBP 1389-2010 Google Scholar

4.

A. H. Hielscher, “Optical tomographic imaging of small animals,” Curr. Opin. Biotech., 16 (1), 79 –88 (2005). http://dx.doi.org/10.1016/j.copbio.2005.01.002 CUOBE3 0958-1669 Google Scholar

5.

K. M. YooF. LiuR. R. Alfano, “When does the diffusion approximation fail to describe photon transport in random media?,” Phys. Rev. Lett., 65 (22), 2210 –2211 (1990). http://dx.doi.org/10.1103/PhysRevLett.65.2210 PRLTAO 0031-9007 Google Scholar

6.

M. Xuet al., “Photon migration in turbid media using a cumulant approximation to radiative transfer,” Phys. Rev. E, 65 (6), 066609 (2002). http://dx.doi.org/10.1103/PhysRevE.65.066609 PLEEE8 1063-651X Google Scholar

7.

M. SakamiK. MitraT. Vo-Dinh, “Analysis of short-pulse laser photon transport through tissues for optical tomography,” Opt. Lett., 27 (5), 336 –338 (2002). http://dx.doi.org/10.1364/OL.27.000336 OPLEDP 0146-9592 Google Scholar

8.

M. J. Niedreet al., “Early photon tomography allows fluorescence detection of lung carcinomas and disease progression in mice in vivo,” Proc. Natl. Acad. Sci., 105 (49), 19126 –19131 (2008). http://dx.doi.org/10.1073/pnas.0804798105 PNASA6 0027-8424 Google Scholar

9.

J. Chenet al., “Time-resolved diffuse optical tomography with patterned light illumination and detection,” Opt. Lett., 35 (13), 125112 (2010). http://dx.doi.org/10.1364/OL.35.002121 OPLEDP 0146-9592 Google Scholar

10.

J. Wuet al., “Fluorescence tomographic imaging in turbid media using early-arriving photons and laplace transforms,” Proc. Natl. Acad. Sci., 94 (16), 8783 –8788 (1997). http://dx.doi.org/10.1073/pnas.94.16.8783 PNASA6 0027-8424 Google Scholar

11.

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

12.

G. Maet al., “Comparison of simplified Monte Carlo simulation and diffusion approximation for the fluorescence signal from phantoms with typical mouse tissue optical properties,” Appl. Opt., 46 (10), 1686 –1692 (2007). http://dx.doi.org/10.1364/AO.46.001686 APOPAI 0003-6935 Google Scholar

13.

D. 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). http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-10-3-159 OPEXFF 1094-4087 Google Scholar

14.

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

15.

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

16.

Q. Fang, “Mesh-based Monte Carlo (MMC),” (2012) http://mcx.sourceforge.net/mmc/ Google Scholar

17.

J. HavelA. Herout, “Yet faster ray-triangle intersection (using sse4),” IEEE Trans. Vis. Comput. Graphics, 16 (3), 424 –428 (2010). http://dx.doi.org/10.1109/TVCG.2009.73 TVCG 1077-2626 Google Scholar

18.

I. Wald, “Realtime ray tracing and interactive global illumination,” Saarland University, (2004). Google Scholar

19.

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

20.

S. Belangeret al., “Realtime diffuse optical tomography based on structured illumination,” J. Biomed. Opt., 15 (1), 016006 (2010). http://dx.doi.org/10.1117/1.3290818 JBOPFO 1083-3668 Google Scholar

21.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl., 15 (2), 41 –93 (1999). http://dx.doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar

22.

B. Dogdaset al., “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol., 52 (3), 577 –587 (2007). http://dx.doi.org/10.1088/0031-9155/52/3/003 PHMBA7 0031-9155 Google Scholar

23.

Q. FangD. Boas, “Tetrahedral mesh generation from volumetric binary and gray-scale images,” in Proc. IEEE Int. Symp. Biomedical Imaging 2009, 1142 –1145 (2009). Google Scholar

24.

G. AlexandrakisF. R. RannouA. F. Chatsiioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol., 50 (17), 4225 –4241 (2005). http://dx.doi.org/10.1088/0031-9155/50/17/021 PHMBA7 0031-9155 Google Scholar

25.

V. VenugopalJ. ChenX. Intes, “Development of an optical imaging platform for functional imaging of small animals using widefield excitation,” Biomed. Opt. Express, 1 (1), 143 –156 (2010). http://dx.doi.org/10.1364/BOE.1.000143 BOEICL 2156-7085 Google Scholar

26.

V. NtziachristosR. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett., 26 (12), 893 –895 (2001). http://dx.doi.org/10.1364/OL.26.000893 OPLEDP 0146-9592 Google Scholar

27.

V. Venugopalet al., “Full-field time-resolved fluorescence tomography of small animals,” Opt. Lett., 35 (19), 3189 –3191 (2010). http://dx.doi.org/10.1364/OL.35.003189 OPLEDP 0146-9592 Google Scholar

28.

A. Joshiet al., “Fully adaptive FEM based fluorescence optical tomography from time-dependent measurements with area illumination and detection,” Med. Phys., 33 (5), 1299 –1310 (2006). http://dx.doi.org/10.1118/1.2190330 MPHYA6 0094-2405 Google Scholar

29.

J. Duttaet al., “Illumination pattern optimization for fluorescence tomography: theory and simulation studies,” Phys. Med. Biol., 55 (10), 2961 –2982 (2010). http://dx.doi.org/10.1088/0031-9155/55/10/011 PHMBA7 0031-9155 Google Scholar

30.

C. D’Andreaet al., “Fast 3D optical reconstruction in turbid media using spatially modulated light,” Biomed. Opt. Express, 1 (2), 471 –481 (2010). http://dx.doi.org/10.1364/BOE.1.000471 BOEICL 2156-7085 Google Scholar

31.

P. K. YalavarthyD. R. LynchB. W. Pogue, “Implementation of a computationally effcient least-squares algorithm for highly under-determined three-dimensional diffuse optical tomography problems,” Med. Phys., 35 (5), 1683 –1697 (2008). http://dx.doi.org/10.1118/1.2889778 MPHYA6 0094-2405 Google Scholar

32.

J. ChenX. Intes, “Time gated perturbation Monte Carlo for whole body functional imaging in small animals,” Opt. Express, 17 (22), 19566 –19579 (2009). http://dx.doi.org/10.1364/OE.17.019566 OPEXFF 1094-4087 Google Scholar

33.

F. AzarX. Intes, “Translational multimodality optical imaging,” Introduction to Clinical Optical Imaging, 1 –13 Artech House, Norwood, MA (2008). Google Scholar

34.

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

35.

Q. FangA. B. 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

Citation Download Citation

Jin Chen, Qianqian Fang, and Xavier Intes "Mesh-based Monte Carlo method in time-domain widefield fluorescence molecular tomography," Journal of Biomedical Optics 17(10), 106009 (1 October 2012). https://doi.org/10.1117/1.JBO.17.10.106009

Published: 1 October 2012

Access the abstract

JOURNAL ARTICLE
8 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 63 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

Tomography

Sensors

Computer simulations

Photon transport

Picosecond phenomena

1.

Introduction

2.

Method

2.1.

Widefield Mesh-Based Monte Carlo Method

Fig. 1

2.2.

Inverse Model for Reconstruction

Eq. (1)

Eq. (2)

2.3.

Computational Settings and Efficiency Evaluation

Eq. (3)

Eq. (4)

2.4.

Accuracy Comparison

Fig. 2

Eq. (5)

2.5.

Tomographic Reconstruction Settings

2.5.1.

Simulation reconstructions

Fig. 3

2.5.2.

Experimental validation

Fig. 4

3.

Results

3.1.

Computational Efficiency

Fig. 5

3.2.

Accuracy of Time-Domain Jacobians

Fig. 6

Table 1

Fig. 7

3.3.

Simulation Reconstruction Performance

Fig. 8

3.4.

Experimental Comparison

4.

Discussion

5.

Conclusion

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years