Getting Started with MMC

1. Mesh generation
2. Define optical properties
3. Prepare the input file
4. Run the simulation

The best place to get started with MMC is the mmc/example/onecube folder. In this folder, we provide a toy-problem for you to play with.

1. Mesh generation

The first step to prepare for an MMC simulation is to generate a tetrahedral mesh. In this case, the matlab script "createmesh.m" does everything for you. Please look into the script to understand the key elements for this process.

sessionid='onecube';

addpath('../../matlab/');
[node,elem]=genT6mesh(0:1,0:1,0:1);
elem=sortrows(elem);
elem(:,1:4)=meshreorient(node,elem(:,1:4));

node=node*10;

srcpos=[2,8,0];
savemmcmesh(sessionid,node,elem,[]);
eid=tsearchn(node,elem(:,1:4),srcpos)

The first line defines a string called a "session id", this will be used in the later simulation to identify the input and output file names. The second line adds the folder to a matlab toolbox provided in the mmc package. This folder includes a simple mesh generator for uniform cubic domain (or grid), such as genT5mesh and genT6mesh. The next line calls the mesh script genT6mesh to produce the tetrahedral mesh: the output "node" stores the node coordinates and elem stores the tetrahedron node indices (starting from 1). You can plot the mesh with tetramesh(elem,node) in matlab.

In the next line, we sort the elements. This is not necessary. We include this line only because Octave and Matlab produce different element orders. The next line is important. Here we call the meshreorient function to test all elements and make sure the nodes are oriented consistently. This function is defined in the iso2mesh toolbox, you have to install it first. The reorientation is done by calculating the volumes of the element and make sure they are all positive.

Then we multiply node by 10 to scale the cubic domain. The next command defines the coordinate of the source. In this case, it is {x=2,y=8,z=0}. The savemmcmesh line is very important. It produces all required mesh files for an mmc simulation. This not only includes the node and element files, but also the volumes of each element (velem file) and the neighboring element indices (facenb file). In fact, inside savemmcmesh, we also call meshreorient to double-check the orientation. The last line do a search to identify the initial element ID (an integer starting from 1) that enclosing the source. This number will be printed in the command window and you will need it in preparing the input file.

This example show you how to make very simple mesh. To produce more realistic complex mesh, you may have to learn how to use iso2mesh or other mesh generators. In any case, you will be able to get a list of node coordinates (output as node), and a list of tetrahedron element indices (output as elem). The rest of the preparation is identical to the above example.

2. Define optical properties

If your mesh contains multiple tissue types or sub-domains, you want to mark them with different labels. We can add the 5-th column to elem to do this (if the property is associated element-wise). If you call iso2mesh to mesh a multi-domain/segmented volume, the output elem array carries this label automatically. Function savemmcmesh will recognize this and automatically put the labels in the right place.

If you elem only contains 4 columns, then you have only 1 tissue type; if your elem contains 5 columns, then the tissue type is max(elem(:,5)). For each tissue type, you need to define the optical properties, and name it as prop_{sessionid}.dat. The following is the content of prop_onecube.dat:

 1 1
 1 0.005 1.0101010101 0.01 1.0

The first row is a header. The first number 1 is the format version. The second 1 indicates how many tissue types are defined in this file. In this case, we have only 1. The following row contains 5 numbers, the first integer is the tissue index, the second number is the absorption coefficient (mua in 1/mm); the third is the scattering coeff (mus in 1/mm); the 4-th is the anisotropy (g); the last is the refractive index (n).

You can find another example at mmc/examples/meshtest/prop_sph1.dat :

 1 2
 1 0.002 1.0 0.01 1.37
 2 0.050 5.0 0.9  1.37

In this case you can see that we have two tissue types, and the they have different mua, mus and g values.

3. Prepare the input file

Input file defines most of other parameters that are not associated a mesh. The content of onecube.inp is shown below

100                   # total photon
17182818            # RNG seed, negative to generate
2 8 0.0        # source position (mm)
0. 0. 1.                # initial directional vector
0.e+00 5.e-09 5e-10  # time-gates(s): start, end, step
onecube              # volume ('uchar' format)
3                    # index of the element that enclosing the source (starting from 1)
4	1            # detector number and radius (mm)
30.0	20.0	1.0  # detector 1 position (mm)
30.0	40.0	1.0  # ...
20.0	30.0	1.0
40.0	30.0	1.0

For each row, the meanings of the numbers are explain in the comments. You must use the same "sessionid", the string you defined for mesh generation, in the 6th row, and the value of "eid" computed in the meshing script as the 7th row. Although we define the detector parameters in the input file, the current version of mmc simply ignore them. We will add support to detectors in the future.

4. Run the simulation

If you have succeeded for all the above steps, now it is time to run mmc to perform the simulation. The mmc binary accepts various command line options. To list all the options, you simply type mmc without any parameters. Here is the output:

usage: mmc <param1> <param2> ...
where possible parameters include (the first item in [] is the default value)
 -i 	       (--interactive) interactive mode
 -s sessionid  (--session)     a string to label all output file names
 -f config     (--input)       read config from a file
 -n [0.|float] (--photon)      total photon number
 -b [0|1]      (--reflect)     1 do reflection at int&ext boundaries, 0 no ref.
 -e [0.|float] (--minenergy)   minimum energy level to trigger Russian roulette
 -U [1|0]      (--normalize)   1 to normalize the fluence to unitary,0 save raw
 -d [1|0]      (--savedet)     1 to save photon info at detectors,0 not to save
 -S [1|0]      (--save2pt)     1 to save the fluence field, 0 do not save
 -C [1|0]      (--basisorder)  1 piece-wise-linear basis for fluence,0 constant
 -u [1.|float] (--unitinmm)    define the length unit in mm for the mesh
 -h            (--help)        print this message
 -l            (--log)         print messages to a log file instead
 -D [0|int]    (--debug)       print debug information (you can use an integer
  or                           or a string by combining the following flags)
 -D [''|MCBWDIOXATRP]          1 M  photon movement info
                               2 C  print ray-polygon testing details
                               4 B  print Bary-centric coordinates
                               8 W  print photon weight changes
                              16 D  print distances
                              32 I  entering a triangle
                              64 O  exiting a triangle
                             128 X  hitting an edge
                             256 A  accumulating weights to the mesh
                             512 T  timing information
                            1024 R  debugging reflection
                            2048 P  show progress bar
      add the numbers together to print multiple items, or one can use a string
example:
       mmc -n 1000000 -f input.inp -s test -b 0 -D TP

For any simulation, the required parameter is the "-f" option. It must be followed by the input file you prepared before. In this case, if you run

 mmc -f onecube.inp
it will read in mesh files labeled by "onecube" (as specified in the line 6 of onecube.inp), then launch 100 photons (as specified in the line 1 of onecube.inp) and produce an output file named "onecube.inp.dat".

You can use other parameters to overwrite the settings in the input file. For example, if you use "-n 10000" in the command line, it will launch 10000 photons instead; if you add "-s test", it will save the output file as "test.dat". When you simulate a large number of photons and want to have some idea on the percentage of completion, you can add "-D P" in the command line to print a progress bar. If you use "-D T", mmc will print the timing information after each simulation step, so you can benchmark for speed. Multiple debug flags can be used together, for example, you can do "-D TP" or "-D PT" to combine the above two functions.

In this particularly example, we have prepared a shell script, run_test.sh, to run mmc with "-D M" and "-D MA" to print out the moving trajectories of the simulated photons. The results can then be visualized by running the "plotmmcdebug.m" script in matlab or octave. You can see the tracks of photons lifting from the injecting source position and scatter and move around inside the cubic domain. The red dots mark the locations of weight accumulating (energy deposit). All photons terminate either at the exterior boundary or it reaches the maximum simulation time window (defined in the input file).

Powered by Habitat