= MMCLAB: MMC for MATLAB and GNU Octave =
*Author: Qianqian Fang <q.fang at neu.edu>
*License: GNU General Public License version 3 (GPLv3)
*Version: this package is part of Mesh-based Monte Carlo (MMC) 2.2, v2024.2
*URL: https://mcx.space/wiki/?Learn#mmc
<toc>
== # Introduction ==
MMC is a mesh-based Monte Carlo (MC) photon simulation software. It can
utilize a tetrahedral mesh to model a complex anatomical structure,
thus, has been shown to be more accurate and computationally efficient
than the conventional MC codes.
MMCLAB is the native MEX version of MMC for MATLAB and GNU Octave. By
converting the input and output files into convenient in-memory variables,
MMCLAB is very intuitive to use and straightforward to be integrated with
mesh generation and post-simulation analyses.
Because MMCLAB contains the same computational codes for multi-threading
photon simulation as in a MMC binary, running MMCLAB inside MATLAB is expected
to give similar speed as running a standalone MMC binary. If your CPU
supports, running mmclab with the 'sse' option can be 25% faster than
the standard mode.
== # Installation ==
Installation of MMCLAB is straightforward. You first download the
MMCLAB package and unzip it to a folder; then you add the folder path
into MATLAB's search path list. This can be done with the "addpath"
command in a working session; if you want to add this path permanently,
use the "pathtool" command, or edit your startup.m (~/.octaverc for
Octave).
After installation, please type "help mmclab" in MATLAB/Octave to print the
help information.
== # How to use MMCLAB in MATLAB/Octave ==
To learn the basic usage of MMCLAB, you can type
help mmclab
and enter in MATLAB/Octave to see the help information regarding how to use this
function. The help information is listed below. You can find the input/output
formats and examples. The input cfg structure has very similar field names as
the verbose command line options in MMC.
<pre>
#############################################################################%
Mesh-based Monte Carlo (MMC) - OpenCL %
Copyright (c) 2010-2024 Qianqian Fang <q.fang at neu.edu> %
https://mcx.space/#mmc & https://neurojson.io/ %
%
Computational Optics & Translational Imaging (COTI) Lab [http://fanglab.org]%
Department of Bioengineering, Northeastern University, Boston, MA, USA %
#############################################################################%
The MCX Project is funded by the NIH/NIGMS under grant R01-GM114365 %
#############################################################################%
Open-source codes and reusable scientific data are essential for research, %
MCX proudly developed human-readable JSON-based data formats for easy reuse.%
%
Please visit our free scientific data sharing portal at https://neurojson.io/%
and consider sharing your public datasets in standardized JSON/JData format %
#############################################################################%
$Rev::4ca899$v2024.2$Date::2023-09-11 17:37:07 -04$ by $Author::Qianqian Fang%
#############################################################################%
Format:
[fluence,detphoton,ncfg,seeds,traj]=mmclab(cfg);
or
fluence=mmclab(cfg);
newcfg=mmclab(cfg,'prep');
[fluence,detphoton,ncfg,seeds,traj]=mmclab(cfg, options);
Input:
cfg: a struct, or struct array. Each element in cfg defines
a set of parameters for a simulation.
option: (optional), options is a string, specifying additional options
option='preview': this plots the domain configuration using mcxpreview(cfg)
option='opencl': force using OpenCL (set cfg.gpuid=1 if not set)
instead of SSE on CPUs/GPUs that support OpenCL
cfg may contain the following fields:
== Required ==
*cfg.nphoton: the total number of photons to be simulated (integer)
*cfg.prop: an N by 4 array, each row specifies [mua, mus, g, n] in order.
the first row corresponds to medium type 0 which is
typically [0 0 1 1]. The second row is type 1, and so on.
*cfg.node: node array for the input tetrahedral mesh, 3 columns: (x,y,z)
*cfg.elem: element array for the input tetrahedral mesh, 4 columns
*cfg.elemprop: element property index for input tetrahedral mesh
*cfg.tstart: starting time of the simulation (in seconds)
*cfg.tstep: time-gate width of the simulation (in seconds)
*cfg.tend: ending time of the simulation (in second)
*cfg.srcpos: a 1 by 3 vector, the position of the source in mesh node length unit
*cfg.srcdir: if defined as [vx, vy, vy], it specifies the incident vector
if defined as [vx, vy, vy, focus], the first 3 elements define
the incident vector; focus controls the convergence or
divergence of the beam:
focus=0: collimated beam
focus<0: diverging beam from an imaginary src at c0-|focus|*[vx vy vz]
focus>0: converging beam, focusing to a point at c0+|focus|*[vx vy vz]
where c0 is the centroid of the source domain. Setting focus does
not impact pencil/isotropic/cone sources.
== MC simulation settings ==
cfg.seed: seed for the random number generator (integer)
if set to a uint8 array, the binary data in each column is used
to seed a photon (i.e. the "replay" mode), default value: 1648335518
cfg.isreflect: [1]-consider refractive index mismatch, 0-matched index
2 - total absorption on exterior surface
3 - prefect reflection (mirror) on exterior surface
cfg.isnormalized:[1]-normalize the output fluence to unitary source, 0-no reflection
cfg.isspecular: [1]-calculate specular reflection if source is outside
cfg.ismomentum: [0]-save momentum transfer for each detected photon
cfg.method: ray-tracing method, ["plucker"]:Plucker, "havel": Havel (SSE4),
"badouel": partial Badouel, "elem": branchless Badouel (SSE),
"grid": dual-grid MMC
cfg.mcmethod: 0 use MCX-styled MC method, 1 use MCML style MC
cfg.nout: [1.0] refractive index for medium type 0 (background)
cfg.minenergy: terminate photon when weight less than this level (float) [0.0]
cfg.roulettesize:[10] size of Russian roulette
cfg.steps: [dx, dy, dz], defines the DMMC grid voxel size,
must be isostropic, i.e. dx=dy=dz, only used when
cfg.method = 'grid', by default dx=dy=dz=1
cfg.unitinmm: defines the default length unit (to interpret mesh nodes, src/det positions
the default value is 1.0 (mm). For example, if the mesh node length unit is
in cm, one should set unitinmm to 10.
cfg.basisorder: [1]-linear basis, 0-piece-wise constant basis
== Source-detector parameters ==
cfg.detpos: an N by 4 array, each row specifying a detector: [x,y,z,radius]
cfg.maxdetphoton: maximum number of photons saved by the detectors [1000000]
cfg.srctype: source type, the parameters of the src are specified by cfg.srcparam{1,2}
'pencil' - default, pencil beam, no param needed
'isotropic' - isotropic source, no param needed
'cone' - uniform cone beam, srcparam1(1) is the half-angle in radian
'gaussian' - a gaussian beam, srcparam1(1) specifies the waist radius
(in default length unit); if one specifies a non-zero focal length
using cfg.srcdir, the gaussian beam can be converging to or
diverging from the waist center, which is located at srcpos+focus*srcdir;
optionally, one can specify the wavelength lambda (in cfg.unitinmm mm),
using srcparam1(2). This will rescale the Gaussian profile according
to w(z)=w0*sqrt(1-(z/z0)^2), where w0 is the waist radius, z is the
distance (in mm) to the waist center (focus), and z0 is the Rayleigh
range (in mm), and z0 is related to w0 by z0=w0^2*pi/lambda
'planar' - a 3D quadrilateral uniform planar source, with three corners specified
by srcpos, srcpos+srcparam1(1:3) and srcpos+srcparam2(1:3)
'pattern' - a 3D quadrilateral pattern illumination, same as above, except
srcparam1(4) and srcparam2(4) specify the pattern array x/y dimensions,
and srcpattern is a floating-point pattern array, with values between [0-1].
if cfg.srcnum>1, srcpattern must be a floating-point array with
a dimension of [srcnum srcparam1(4) srcparam2(4)]
Example: <demo_photon_sharing.m>
'fourier' - spatial frequency domain source, similar to 'planar', except
the integer parts of srcparam1(4) and srcparam2(4) represent
the x/y frequencies; the fraction part of srcparam1(4) multiplies
2*pi represents the phase shift (phi0); 1.0 minus the fraction part of
srcparam2(4) is the modulation depth (M). Put in equations:
S=0.5*[1+M*cos(2*pi*(fx*x+fy*y)+phi0)], (0<=x,y,M<=1)
'arcsine' - similar to isotropic, except the zenith angle is uniform
distribution, rather than a sine distribution.
'disk' - a uniform disk source pointing along srcdir; the radius is
set by srcparam1(1) (in default length unit)
'fourierx' - a general Fourier source, the parameters are
srcparam1: [v1x,v1y,v1z,|v2|], srcparam2: [kx,ky,phi0,M]
normalized vectors satisfy: srcdir cross v1=v2
the phase shift is phi0*2*pi
'fourierx2d' - a general 2D Fourier basis, parameters
srcparam1: [v1x,v1y,v1z,|v2|], srcparam2: [kx,ky,phix,phiy]
the phase shift is phi{x,y}*2*pi
'zgaussian' - an angular gaussian beam, srcparam1(1) specifies the variance in
the zenith angle
cfg.{srcparam1,srcparam2}: 1x4 vectors, see cfg.srctype for details
cfg.srcpattern: see cfg.srctype for details
cfg.srcnum: the number of source patterns that are
simultaneously simulated; only works for 'pattern'
source, see cfg.srctype='pattern' for details
Example <demo_photon_sharing.m>
cfg.replaydet: only works when cfg.outputtype is 'jacobian', 'wl', 'nscat', or 'wp' and cfg.seed is an array
0 replay all detectors and sum all Jacobians into one volume
a positive number: the index of the detector to replay and obtain Jacobians
cfg.voidtime: for wide-field sources, [1]-start timer at launch, 0-when entering
the first non-zero voxel
by default, mmc assumes the mesh and source position settings are all in mm unit.
if the mesh coordinates/source positions are not in mm unit, one needs to define
cfg.unitinmm (in mm) to specify the actual length unit.
== Optional mesh data ==
-cfg.facenb: element face neighbohood list (calculated by faceneighbors())
-cfg.evol: element volume (calculated by elemvolume() with iso2mesh)
-cfg.e0: the element ID enclosing the source, if not defined,
it will be calculated by tsearchn(node,elem,srcpos);
if cfg.e0 is set as one of the following characters,
mmclab will do an initial ray-tracing and move
srcpos to the first intersection to the surface:
'>': search along the forward (srcdir) direction
'<': search along the backward direction
'-': search both directions
== Output control ==
cfg.issaveexit: [0]-save the position (x,y,z) and (vx,vy,vz) for a detected photon
cfg.issaveref: [0]-save diffuse reflectance/transmittance on the exterior surfaces.
The output is stored as flux.dref in a 2D array of size [#Nf, #time_gate]
where #Nf is the number of triangles on the surface; #time_gate is the
number of total time gates. To plot the surface diffuse reflectance, the output
triangle surface mesh can be extracted by faces=faceneighbors(cfg.elem,'rowmajor');
where 'faceneighbors' can be found in the iso2mesh toolbox.
Example: see <demo_mmclab_basic.m>
cfg.issaveseed: [0]-save the RNG seed for a detected photon so one can replay
cfg.isatomic: [1]-use atomic operations for saving fluence, 0-no atomic operations
cfg.outputtype: 'flux' - output fluence-rate
'fluence' - fluence,
'energy' - energy deposit,
'jacobian' - mua Jacobian (replay mode)
'wl'- weighted path lengths to build mua Jacobian (replay mode)
'wp'- weighted scattering counts to build mus Jacobian (replay mode)
cfg.debuglevel: debug flag string, a subset of [MCBWDIOXATRPE], no space
cfg.debugphoton: print the photon movement debug info only for a specified photon ID
cfg.maxjumpdebug: [10000000|int] when trajectory is requested in the output,
use this parameter to set the maximum position stored. By default,
only the first 1e6 positions are stored.
fields marked with * are required; options in [] are the default values
fields marked with - are calculated if not given (can be faster if precomputed)
type: omit or 'omp' for multi-threading version; 'sse' for the SSE4 MMC,
the SSE4 version is about 25% faster, but requires newer CPUs;
if type='prep' with a single output, mmclab returns ncfg only.
Output:
fluence: a struct array, with a length equals to that of cfg.
For each element of fluence, fluence(i).data is a 2D array with
dimensions [size(cfg.node,1), total-time-gates] if cfg.basisorder=1,
or [size(cfg.elem,1), total-time-gates] if cfg.basisorder=0.
The content of the array is the normalized fluence-rate (or others
depending on cfg.outputtype) at each mesh node and time-gate.
If cfg.issaveref is set to 1, fluence(i).dref is not empty, and stores
the surface diffuse reflectance (normalized by default). The surface mesh
that the dref output is attached can be obtained by faces=faceneighbors(cfg.elem,'rowmajor');
detphoton: (optional) a struct array, with a length equals to that of cfg.
Starting from v2016.5, the detphoton contains the below subfields:
detphoton.detid: the ID(>0) of the detector that captures the photon
detphoton.nscat: cummulative scattering event counts in each medium
detphoton.ppath: cummulative path lengths in each medium (partial pathlength)
one need to multiply cfg.unitinmm with ppath to convert it to mm.
detphoton.mom: cummulative cos_theta for momentum transfer in each medium
detphoton.p or .v: exit position and direction, when cfg.issaveexit=1
detphoton.w0: photon initial weight at launch time
detphoton.prop: optical properties, a copy of cfg.prop
detphoton.data: a concatenated and transposed array in the order of
[detid nscat ppath mom p v w0]'
"data" is the is the only subfield in all mmclab before 2016.5
ncfg: (optional), if given, mmclab returns the preprocessed cfg structure,
including the calculated subfields (marked by "-"). This can be
used in the subsequent simulations to avoid repetitive preprocessing.
seeds: (optional), if give, mmclab returns the seeds, in the form of
a byte array (uint8) for each detected photon. The column number
of seed equals that of detphoton.
trajectory: (optional), if given, mmclab returns the trajectory data for
each simulated photon. The output has 6 rows, the meanings are
id: 1: index of the photon packet
pos: 2-4: x/y/z/ of each trajectory position
5: current photon packet weight
6: enclosing element's ID
By default, mcxlab only records the first 1e7 positions along all
simulated photons; change cfg.maxjumpdebug to define a different limit.
Example:
cfg.nphoton=1e5;
[cfg.node face cfg.elem]=meshabox([0 0 0],[60 60 30],6);
cfg.elemprop=ones(size(cfg.elem,1),1);
cfg.srcpos=[30 30 0];
cfg.srcdir=[0 0 1];
cfg.prop=[0 0 1 1;0.005 1 0 1.37];
cfg.tstart=0;
cfg.tend=5e-9;
cfg.tstep=5e-10;
cfg.debuglevel='TP';
% preprocessing to populate the missing fields to save computation
ncfg=mmclab(cfg,'prep');
cfgs(1)=ncfg; % when using struct array input, all fields must be defined
cfgs(2)=ncfg;
cfgs(1).isreflect=0;
cfgs(2).isreflect=1;
cfgs(2).detpos=[30 20 0 1;30 40 0 1;20 30 1 1;40 30 0 1];
% calculate the fluence and partial path lengths for the two configurations
[fluxs,detps]=mmclab(cfgs);
This function is part of Mesh-based Monte Carlo (MMC) URL: http://mcx.space/#mmc
License: GNU General Public License version 3, please read LICENSE.txt for details
</pre>
== # Examples ==
We provided several examples to demonstrate the basic usage of MMCLAB;
some of the examples also perform validations of MMC algorithm in both
homogeneous and heterogeneous domains. These examples are explained below:
==== demo_mmclab_basic.m ====
In this example, we show the most basic usage of MMCLAB. This includes
how to define the input configuration structure, launch MMC simulations
and plot the output data.
==== demo_example_validation.m ====
In this example, we validate MMCLAB with a homogeneous medium in a
cubic domain. This is the same as in mmc/examples/validation;
the details of the simulation is described in Fig.2 of [Fang2010].
==== demo_example_meshtest.m ====
In this example, we validate the MMCLAB solver with a heterogeneous
domain and the analytical solution from the diffusion model. The
domain is consisted of a 6x6x6 cm box with a 2cm diameter sphere
embedded at the center.
This test is identical to the simulations under mmc/examples/meshtest,
which are described by Fig. 3 in [Fang2010].
==== demo_example_onecube.m ====
In this example, we demonstrate how to use the debug flags
to print photon trajectories. The domain is a simple 10x10x10mm
cube.
This test is similar to the simulations under mmc/examples/onecube.
==== demo_example_replay.m ====
In this example, we run MMC using a homogeneous cubic domain,
same as in demo_example_validation.m. We save the seeds of the
detected photons, and then rerun these photons again (i.e.
the "replay" part). This feature allows one to identify features
that are only related to a source/detector pair.
==== demo_compare_mmc_mcx.m ====
In this example, we compare MMC and MCX using both simple and
widefield sources.
==== demo_sfdi_2layer.m ====
In this example, we simulate an spatial-frequency domain imaging source
using a 2-layer brain model.
==== demo_mmclab_slit.m ====
In this example, we show how to use a slit (line light source) in
an MMC simulation.
== # How to compile MMCLAB ==
To compile MMCLAB from source code, you need to make sure your
computer have the following requirements:
* your computer should have MATLAB C compiler (mex) and/or octave3.x-headers installed
* to compile the SSE versions, your computer should support SSE4 instructions
* the commands "mex" and "mkoctfile" should be in the search path ($PATH)
* you should have both gcc and g++ with version 4.4 or newer
* for windows, you should install Cygwin64 with mingw64-i686-gcc and make packages
* for windows, you need to create a mexopts.bat file as \
C:\Users\<username>\AppData\Roaming\MathWorks\MATLAB\RXXXX\mexopts.bat \
with content modified from mmc/src/mexopts_cygwin64_gcc.bat
To compile MMCLAB for MATLAB, you need to cd mmc/src directory, and type
make mex
or
make mexsse
from a shell window. You need to make sure your MATLAB is installed and
the command <tt>mex</tt> is included in your PATH environment variable. Similarly,
to compile MMCLAB for Octave, you type
make oct
or
make octsse
The command <tt>mkoctfile</tt> must be accessible from your command line
and it is provided in a package named "octave3.x-headers" in Ubuntu (3.x
can be 3.2 or 3.4 etc).
Compiling MMCLAB on Windows and Mac OS, you should use the following
command:
make ... -f makefile_sfmt ...
If you have multiple gcc/g++ installed on your system, and the default
gcc and g++ is older than 4.4 (use gcc -v to print version), you can
use the following command to specify the correct gcc command:
make ... CC=/path/to/gcc/4.4 CXX=/path/to/g++/4.4
== # Screenshots ==
Screenshot for using MMCLAB in MATLAB:
http://mcx.sourceforge.net/upload/matlab_mmclab.png
Screenshot for using MMCLAB in GNU Octave:
http://mcx.sourceforge.net/upload/octave_mmclab.png
== # Reference ==
[Fang2019] Qianqian Fang* and Shijie Yan, "GPU-accelerated mesh-based \
Monte Carlo photon transport simulations," J. of Biomedical Optics, in press, 2019. \
Preprint URL: https://www.biorxiv.org/content/10.1101/815977v1
[Brain2Mesh2020] Anh Phong Tran† , Shijie Yan† , Qianqian Fang*, (2020) "Improving \
model-based fNIRS analysis using mesh-based anatomical and light-transport models," \
Neurophotonics, 7(1), 015008, URL: https://doi.org/10.1117/1.NPh.7.1.015008
[Fang2010] Fang Q, "Mesh-based Monte Carlo method using fast ray-tracing \
in Plucker coordinates," Biomed. Opt. Express 1(1), 165-175 (2010)