- 1. validation - Validation of MMCM in a homogeneous cubic domain
- 2. meshtest - Validation of MMCM in a heterogeneous domain - a sphere inside a cube
- 3. onecube - Visual debugging of photon migration in a simple mesh
- 4. dcs - Validation of MMCM simulation of DCS auto-correlation measurements
- 5. statnoise: Fluence variation vs photon numbers (Experimental)
- 6. colin27 - Simulation with the Colin27 Human Brain Atlas
- 7. planar - Validation of uniform planar source
- 8. sphere - Validation of uniform planar source use a spherical domain
- 9. sfdi2layer - Validation of SFDI source
- 10. Reference:

The examples contained in each subfolders are explained below:

examples/

|-validation - homogeneous cubic domain 60x60x60, i.e. Fig. 2 in Fang2010 |-meshtest - 10mm sphere embedded in the cubic domain, Fig. 3 in Fang2010 |-mcxsph - scripts to run MCX for the above two cases |-onecube - debugging flag show case 1, simple mesh without reflection |-reftest - debugging flag show case 2, simple mesh with reflection |-rngtest - random number generator test |-dcs - validation of momentum transfer in DCS simulation |-colin27 - simulation with the Colin27 Human Brain Atlas |-statnoise - test statistical noise against launched photon numbers |-misctest - other misc. tests |-planar - validation of uniform planar source in a simple cube |-sphere - validation of uniform planar source use a spherical domain |-sfdi2layer - validation SFDI illumination on a two-layer brain model |-regression/exitangle - regression testing of photon exit angles

A complex human brain atlas mesh, i.e. Fig. 4 in [Fang2010], can be downloaded separately at the following URL:

http://mcx.sourceforge.net/cgi-bin/index.cgi?MMC/CollinsAtlasMesh

In this example, we validate the MMCM algorithm using a homogeneous cubic domain. This simulation generates the results shown in Fig. 2 in the paper.

The cubic domain has a dimension of 60x60x60 mm with optical properties mua=0.001, mus=1, n=1.0 and g=0.01. The analytical solution can be computed by the cwdiffusion (for CW solution) and tddiffusion (for time-domain solutions) functions in the package of Monte Carlo eXtreme (MCX) under the mcx/utils directory.

1. In order to run this example, you need to first download and install the svn version of iso2mesh from

http://iso2mesh.sourceforge.net/cgi-bin/index.cgi?Download

the password is anonymous_user

2. Start matlab, run createmesh to create all mesh files. Matlab will print a value for variable eid. This value indicates the initial element ID in the mesh that encloses the source. You need to manually set this number as the second integer on the 7-th line in cube.inp input file.

3. Run the simulation bash script run_test.sh, this will run 30000000 photons with the specified mesh.

4. To produce Fig. 2, you need to run MCX simulations. You have to download MCX package from

http://mcx.sourceforge.net/cgi-bin/index.cgi?Download

using the password "anonymous_user". You can also download the pre-compiled binary package from http://mcx.sf.net/cgi-bin/index.cgi?Download

5. Assuming you have compiled and installed MCX, you can run MCX simulation by

5.1 go to mmc/examples/mcxsph run createmcxbin from matlab

5.2 open benchbox.sh, edit the thread number (-t) and thread block size (-T) based on the compute capability of your card. For 8800GT/9800GT, the -T number can not exceed 128; for 280/285/295, -T can not be more than 256; for 470, -T can not be more than 576

5.3 run benchbox.sh to generate the MCX output box.mc2

6. When all simulations are down, you start matlab again, and run plotcuberes to generate the plots.

In this example, we validate MMCM algorithm using a sphere object in a homogeneous cubic background. This simulation generates the results shown in Fig. 3 in the paper.

The cubic domain has a dimension of 60x60x60 mm with optical properties mua=0.002, mus=1, n=1.37 and g=0.01. A sphere is centered at [30 30 30]mm with a radius of 10mm. The optical properties of the sphere is mua=0.05, mus=5, n=1.37 and g=0.9. The analytical solution, approximated by a sphere inside a infinite slab, can be computed by the sphdiffusionslab function in mmc/matlab/ directory (this has already been done, the result for plane y=30 is saved in the sphdiffsemiinf.mat file).

To validate MMCM, we generate an FE-mesh for the sphere and the cube. Three meshes are tested: sph1: a coarse FE mesh with only 10,000 nodes, sph2: a uniform dense FE mesh with 60000 nodes, and sph3: a mesh with higher density around the sphere surface and near the source. The case sph2 and sph3 correspond to "MMCM Mesh 1" and "MMCM Mesh 2" in the paper (Table 1), respectively.

1. First, you need to create the 3 meshes to run this example. You need to first download and install the svn version of iso2mesh from

http://iso2mesh.sourceforge.net/cgi-bin/index.cgi?Download

the password is anonymous_user

2. Start matlab, run createmesh to create all mesh files

3. Open initial_elem.txt, and type the integer in each line to replace the second integer at the 7-th line of sph{1,2,3}.inp (this sets the initial element ID, because the mesh generator may generate a different mesh on your system, you have to update this manually)

4. run the simulation bash script run_test.sh

5. when all 3 simulations are complete, you start matlab again, and run plotmmcsph to generate the plots.

In this example, we run photon migration in a toy-problem to show the basic steps and behaviors of the simulation code. The mesh used in this case is a simple cube, splitted into 6 tetrahedra. The edge length of the cube is 10mm. A total of 20 photons are simulated starting from a source position [2,8,0] along an incident vector [0 0 1]. We run the simulation and print out the moving trajectories of the photon, and the exit and accumulation sites. Using a matlab script, one can visualize the photon moving from the output of the simulation.

1. The mesh files for this example has already been generated. if you want to regenerate the mesh files, please open matlab or octave and run "createmesh"

2. Run the simulation by typing

./run_test.sh

3. Start matlab, run "plotmmcdebug" to see the photon trajectory plot

- Author: Stefan Carp (carp <at> nmr.mgh.harvard.edu or carp <at> alum.mit.edu)
- Date: October 18, 2011

Monte Carlo simulations of light propagation are generally used to estimate fluence distributions, or the intensity of diffusely reflected light from a turbid medium. By storing the history of scattering angles on a per-photon basis (momentum transfer) in addition to the photon pathlenghts it is possible to derive the expected electric field auto-correlation on the surface of the medium for a given distribution of scatterer dynamics in the medium.

Diffuse correlation spectroscopy (DCS) as Diffusing Wave Spectroscopy (DWS) is often called in the biomedical arena, is a non-invasive method for measuring tissue blood flow based on the extensions of the dynamic light scattering (DLS) to the multiple scattering regime.

DCS measurements consist of measuring diffusely reflected light using photon counting photo detectors, and computing the temporal intensity auto-correlation curve g2(tau) for typical correlation delays from 100 ns to 1 s. These curves can be fitted using analytical expressions derived from the correlation diffusion equation to estimate scatterer motion parameters. These expression give the field auto- correlation function g1, related to the measured g2 by the Siegert relationship g2=1+beta*g1^2, where beta is a factor that depends on the detection geometry and light polarization (usually lost in tissue due to multiple scattering).

Here we use an MMCM simulation of photon propagation in a thick slab (considered semi-infinite) for a symmetrical pattern of 4 2mm diameter detectors 15 mm away from a source. We turn on momentum transfer recording (the "-m" parameter in the command line), and compute the field auto-correlation curve g1 assuming Brownian diffusion like scatterer motion (shown to approximate biological tissue measurements well). We also simulate g2 by assuming a beta of 0.4 (typical for room light experimental conditions). We then fit both simulated curves using the analytical model and attempt to recover the Brownian diffusion coefficient. The 2nd example where a simulated intensity auto- correlation is fit represents the way experimental DCS data would be processed to obtain a blood flow index.

1. First, you need to create the mesh to run this example. If you have not done it already, install iso2mesh from Sourceforge:

http://iso2mesh.sourceforge.net/cgi-bin/index.cgi?Download

2. Start Matlab, run createmesh and make a note of the value of the element ID that encloses the source (reported by the createmesh script). Edit "dcs.inp" to set line 7 to this number. If you receive errors, make sure iso2mesh is in your Matlab/Octave path

3. Run the simulation shell script "run_test.sh" from the command line. It will run 20 million photons.

4. Run "dcs_example.m" in Matlab to produce the test-figures and verify the reported fitting error is less than 5%.

Boas D.A., "Diffuse Photon Probes of Structural and Dynamical Properties of Turbid Media: Theory and Biomedical Applications", Ph.D. Thesis, Univ. of Pennsylvania, 1996

1. First, you need to create the mesh files to run this example. You need to first download and install iso2mesh version 1.0 from

http://iso2mesh.sourceforge.net/cgi-bin/index.cgi?Download

or, you can download the git version of iso2mesh from

https://github.com/fangq/iso2mesh

2. Start matlab, run createmesh to create all mesh files. Matlab will print a value for variable eid. This value indicates the initial element ID in the mesh that encloses the source. You need to open vartest.json with a text editor and replace the number after "InitElem".

3. Run the simulation bash script run_test.sh, this will run the simulation for 10 repetitions and save fluence maps for 1E4, 1E5, 1E6 and 1E7 photons.

4. When all simulations complete, you start matlab again, and run plotstdhist to generate a plot between signal standard deviation and photon number.

- This test is experimental. The interpretations to the results need

In this example, we test MMC with a complex human brain atlas: the Colin27 atlas. This example is explained in details in the original MMC paper [Fang, BOE 1(1),2010]. Running this example allows you to reproduce Fig. 4 in the paper.

This example requires the Colin27 atlas mesh. It can be downloaded at

http://mcx.sourceforge.net/cgi-bin/index.cgi?MMC/Colin27AtlasMesh

In the original MMC paper, the V1 atlas mesh was used. In this example, we configure to use the V2 atlas mesh. We believe the difference caused by the mesh versions is negligible.

1. First, you need to download the Colin27 atlas mesh from

http://mcx.sourceforge.net/cgi-bin/index.cgi?MMC/Colin27AtlasMesh

after download, you should extract file "MMC_Collins_Atlas_Mesh_Version_2L.mat" and put it under this folder

2. Start matlab, run createmesh to generate all the mesh files.

3. run the simulation bash script run_test.sh

4. The time-resolved fluence/flux maps will be saved in a file named brain.dat. You can then load it to make plots. You will need the qmeshcut function from iso2mesh toolbox. Please refer to examples/meshtest/plotmmcsph.m script for more details.

In this example, we test a simple widefield source -- uniform planar on a digimouse mesh model. Three corners of the 3D quadrilateral are specified by srcpos, srcpos+srcparam1(0:2), srcpos+srcparam2(0:2). The original mesh is re-tessellated after adding source nodes in addsource.m . The digimouse mesh model can be downloaded at http://mcx.sourceforge.net/cgi-bin/index.cgi?MMC/DigimouseMesh

In this example, we test a simple widefield source -- uniform planar over a sphere model. Three corners of the 3D quadrilateral are specified by srcpos, srcpos+srcparam1(0:2), srcpos+srcparam2(0:2). The original mesh is re-tessellated after adding source nodes in addsource_sphere.m .

In this example, we test a spatial frequency domain (SFDI) source -- on a two-layer cubic model.

Three corners of the 3D quadrilateral are specified by srcpos, srcpos+srcparam1(0:2), srcpos+srcparam2(0:2), k-numbers of the spatial frequency along two directions are specified by srcparam1(3) and srcparam2(3).

The original mesh is saved as "**_media.dat' and the retessellated mesh is saved as "**_sfdi.dat".

The optical properties mimic that of the skull and CSF of a human brain model, and can be found at

http://mcx.sourceforge.net/cgi-bin/index.cgi?MMC/Colin27AtlasMesh#Tissue_optical_properties.

Run "createmesh.m" first to generate mesh (before/after retessellation), then run "run_test.sh" to perform the simulation. Finally, run "plot_result.m" to see the continuous wave (CW) fluence distribution inside the model.

[Fang2010] Qianqian Fang, "Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates," Biomed. Opt. Express 1, 165-175 (2010) [Yao2016] Ruoyang Yao, Xavier Intes, Qianqian Fang, "Generalized mesh-based Monte Carlo for wide-field illumination and detection via mesh retessellation," Biomed. Optics Express, 7(1), 171-184 (2016)