Simulation of optical light in objects with high heterogeneity, such as biological tissue, has been a long-standing challenge in the field of biomedical engineering. Many of the recent groundbreaking bioengineering technologies utilize some form of biomarker, many of which emit optical light, such as fluorescence, as their signal. Thus simulating the nature of the light traveling through biological tissue has very significant implications for the future of biomedical engineering research.
For many years, the main limitation to light (photon) simulation has been the highly variable nature of photon propagation in biological tissue (scattering), which makes the determination of the ray-boundary intersection prohibitively expensive from a computational perspective. The “TIM-OS” algorithm we proposed is a unique algorithm for light propagation in highly heterogeneous domain. In this algorithm, the photon propagation is guided by a tetrahedral mesh. Since the photon is highly scattered, the average step size of a photon is much smaller than the size of a tetrahedron, thus only one tetrahedron requires analysis per step. And, each tetrahedron has only four triangular surfaces, so it is a fairly straightforward geometry to compute the photon surface interaction.
A simplified TIM-OS version for MCML style input can be download here with source code.
To run the examples, please run:
The above runall.sh script is for machines with 32 cores. You may need to adjust the parameters to fit your machine.
Note: We also observed problems under certain Intel Compiler. Some versions of Intel C/C++ compiler have problems for pthreads. If you face this problem, please reduce the optimization level from -O3 to -O1.
timos -f finite_element_filename -p optical_filename -s source_filename -o output_filename -m [s|i|is|si|t] -T time_step num_time_step -t num_thread -r start_rand_id
TIM-OS needs at least four parameters to work:
- A finite element file contains the tetrahedral mesh for the specimen.
- A file contains optical parameters for each region of the tetrahedral mesh.
- A file defines the optical source.
- And a file to store the output.
The program supports time domain simulation:
-T time_step num_time_step
To enable time domain simulation, we need to provide the above parameters. The time_step is in nano second (10^-9 sec). num_time_step should be an integer number greater than 1. Note that time domain simulation is memory intensive. The main program needs about (#Tetrahedron + # Surface Triangle) * num_time_step * 8 to store simulation result. For the above demo, TIM-OS need about 20MB to store one time step. 100 time steps will need about 2GB memory. So, choose the number of time step wisely. To disable time domain simulation, simply do not supply -T parameter
The program supports several output formats:
s: output surface fluence
i: output internal fluence
si|is: output surface and internal fluences
The program supports multi-thread:
If the host system supports multi-core, TIM-OS can launch as many threads as the number of cores. It is also possible to launch more threads than the number of cores. The performance of multi-threads is depending on the system’s thread scheduling algorithm. We observed that for a dual Xeon 5520 machine (8 cores and 16 hyper-threads), launch 32 threads can reach maximum performance.
It is possible to split a large problem into several small jobs. Each job can be ran on a different machine. In this case, it is better to use different random stream for each thread across all the machines. The random number generator used in this program can have at most 1024 independent streams. Since each thread has a different stream, the number of machines can be used for a problem are: 1024/(num_thread_each_machine). For example, if each machine will use 16 threads, the first machine will have a start_rand_id = 1 and the second machine will have a start_rand_id = 17.
In the package, there are several samples can be used to test and help you understand the usage of TIM-OS. Look in
where you can find the following files:
- cube_5med.mesh: defines tetrahedral mesh for a cubic with five different regions.
- cube_5med.opt: defines optical parameters.
- cube_5med.source: defines source.
Format for tetrahedral mesh file.
The following data is from cube_5med.mesh. The comments are added later. TIM-OS does not allow comments in this version.
48000 | number of tetrahedra in this mesh
-5 -5 -5 | follows 9261 lines data
-5 -4.5 -5 | each line is a node’s position (X, Y, Z)
-5 -4 -5
-5 -3.5 -5
-5 -3 -5
-5 -2.5 -5
-5 -2 -5
-5 -1.5 -5
5 5 5
1 2 23 443 1 | follows 48000 lines of data
1 443 23 464 1 | each line defines a tetrahedron:
1 442 443 464 1 | node_idx1 node_idx2 node_idx3 node_idx4 region_idx
1 22 23 464 1 | node index starts from 1
1 22 442 464 1 | region index starts from 1
6020 6021 6042 6462 2 | tetrahedron belongs to region 2
Format for optical parameter file.
The following data is from cube_5med.opt.
| 1: assign the optical parameters to the regions of the phantom | 2: assign the optical parameters to the tetrahedra of the phantom | In this setting, each tetrahedron can have different optical parameters | Since there are five regions in the model file, we need to define five sets of optical parameters. | Each line contains four values: 0.05 20 0.9 1.3 | mu_a: absorption coefficient 0.20 20 0.8 1.2 | mu_s: scattering coefficient 0.10 10 0.9 1.4 | g: anisotropy factor 0.20 20 0.9 1.5 | n: refractive index | 1: the environment has an uniform refractive index | 2: the environment will match the refractive index of the boundary tetrahedron which means no reflection on the boundary. 1.0 | If the above line is 1, then this is the refractive index of the environment …
Format for source file.
The following data is from cube_5med.source
1 0.3 1.31 1.333 10000000 | the second line define a isotropic point source
Current version of TIM-OS accepts three types of sources: isotropic point, isotropic region, pencil beam
isotropic point source (type id 1)
To define a point source, we need first specify it is a isotropic point source (type id 1). Then, we need the position of the source (x y z) and how many photons will be simulated for this point source. Definitely, the position of the source should be inside the mesh.
Source Type | Position(X,Y,Z) | NumPhoton
isotropic region source (type id 2)
To define a region source, we need the source type (type id 2), the tetrahedron it in (the index of the tetrahedron), and the number of photons:
Source Type | Tetrahedron Index | NumPhoton
pencil beam source (type id 11)
To define a pencil beam source, we need the source type (type id 11), the tetrahedron index it hits on the surface, the position (x y z) it hits, the direction (ux uy uz) of the source, and the number of photons:
11 4537 0 0 0 0 0 1 100000000 Source Type | Tetrahedron Index | Position(X,Y,Z) | Direction (UX,UY,UZ) | NumPhoton
The support of pencil beam source is very limited in this version. The pencil beam should be normal to the surface triangle it hits on. TIM-OS will not check whether it is normal or not.
Now it is time to run TIM-OS:
The program will print some progress information:
Optical filename: cube_5med.opt
Fem data filename: cube_5med.mesh
Source filename: cube_5med.source
Begin read optical parameters.
Optical parameter per region
The environment has a uniform refractive index.
End read optical parameters.
Begin read finite element file.
End read finite element file.
End read finite element file.
Begin Preprocessing the finite element mesh.
End Preprocessing the finite element mesh.
Begin read source file.
Begin check source
End check source
End read source file.
Num Photon: 10000000
Time domain setting. Step size: 0.1 ns, Num Step: 2
Output filename: result_cube_5med.dat
Output format: surface and internal fluence
Number of threads: 32
Start random stream: 1
Absorbed Fraction: 0.750473
total running time: 125sec
total CPU time: 1999.49 sec
Write output file
The result will be saved in test.dat:
% Fem data filename: cube_5med.mesh
% Source filename: cube_5med.source
% Num Photon: 10000000
% Time domain setting. Step size: 0.1 ns, Num Step: 2
% Number of threads: 32
% Start random stream: 1
1 4800 2
| The first number “1” indicates that the follow section is for surface fluence
| The second number “4800” indicates that there are 4800 surface triangles
| The third number “2” indicates that the there are 2 time steps in the simulation 1 2 23 1.25000000e-01 0.00000000e+00 4.93664765e+01
| The first three numbers are the node index for the triangle
| The fourth number is the size of the triangle
| The rest numbers are the surface fluence 1 2 443 1.25000000e-01 2.90289403e+00 2.07782283e+01
1 22 23 1.25000000e-01 2.91742664e+00 4.77912430e+01
1 22 442 1.25000000e-01 3.03633752e+00 4.47211477e+01
1 442 443 1.25000000e-01 4.96746220e+00 4.13068283e+01
2 3 24 1.25000000e-01 5.33395510e+00 7.47500342e+01
2 48000 2
| The first number “2” indicates that the follow section is for volumetric fluence
| The second number “48000” indicates that there are 48000 tetrahedra
| The third number “2” indicates that the there are 2 time steps in the simulation1 2 23 443 2.08333333e-02 1.61099729e+01 2.10588505e+02
| The first four numbers are the node index for the tetrahedron
| The fifth number is the size of the tetrahedron
| The rest numbers are the volumetric fluence
1 23 443 464 2.08333333e-02 1.50813033e+01 2.80943125e+02
1 442 443 464 2.08333333e-02 1.94574351e+01 2.72147497e+02
1 22 23 464 2.08333333e-02 8.61062892e+00 2.38070462e+02
1 22 442 464 2.08333333e-02 1.07722744e+01 2.61379380e+02
22 442 463 464 2.08333333e-02 2.60714773e+01 3.14004137e+02
- H. Shen and G. Wang. Tetrahedron-based inhomogeneous Monte-Carlo optical simulator. Phys. Med. Biol. 55:947-962, 2010.
Featured Article in Physics in Medicine & Biology 2010.
- Y. Lu, B. Zhu, H. Shen, J.C. Rasmussen, G. Wang, E.M. Sevick-Muraca. Parallel Adaptive Finite Element Simplified Spherical Harmonics Approximations Solver for Frequency-domain Fluorescence Molecular Tomography. Physics in Medicine and Biology, 55:4624-4645 (2010).
Featured Article in Physics in Medicine & Biology 2010.
Highlighted in Physics in Medicine & Biology 2010.
- W.C. Vogt, H. Shen, G. Wang, and C.G. Rylander. Parametric study of tissue optical clearing by localized mechanical compression using combined finite element and Monte Carlo simulation. Journal of Innovative Optical Health Sciences, 3(3):203-211, (2010).
- H. Shen and G. Wang. A study on tetrahedron-based inhomogeneous Monte-Carlo optical simulation. Biomed. Opt. Express. 2(1):44-57 (2010).
- H. Shen and G. Wang. Reply to “Comment on ‘A study on tetrahedron-based inhomogeneous Monte-Carlo optical simulation'”. Biomedical Optics Express. Accepted.
A full version shows why tetrahedron-based fluence quantification scheme is better than node-based fluence quantification scheme can be download here.
For general information about optical Monte Carlo simulation, please refer to the original “MCML” paper by Dr. Lihong Wang.