Readme and User Manual for MCmatlab

What is MCmatlab?

MCmatlab is a Monte Carlo simulation for modeling light propagation in a 3D voxel space. Fluorescence can optionally be simulated after simulation of the excitation light. Included is also a finite element simulation for temperature increase and heat diffusion in the same voxel space.

Primarily targeted for tissue optics, but can be used in any environment with turbid media in which the wave nature of light (interference phenomena) is negligible and a ray-tracing model is appropriate.

You can find the latest version of MCmatlab on https://github.com/ankrh/MCmatlab. You can also use MCmatlab directly in MATLAB online by clicking the button below. If you publish results obtained with this software, we would be thankful if you cited its accompanying article: doi:10.1117/1.JBO.23.12.121622

License

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/.

How do I get set up?

Requirements

• Windows 7, macOS 10.12 (Sierra) or newer or Linux
• MATLAB R2018a or newer (only the intel-versions of MATLAB are supported)
• (For GPU accelerated computation) A Windows PC with a CUDA-enabled graphics card and the MATLAB Parallel Computing Toolbox

Helper files

All the helper functions needed for running MCmatlab are located in the folder (MATLAB package) "+MCmatlab", and its parent folder therefore has to be on your MATLAB path when trying to run any MCmatlab model files. You will not need to modify any of the contained files. We recommend keeping the model files in a folder together with the "+MCmatlab" folder, such that you don't need to manually add the parent folder of "+MCmatlab" to your MATLAB path.

How do I use MCmatlab?

Compilation

The folders include all the executables necessary, so you don't need to compile anything. If, however, you want to change the routine in either the MCmatlab.c source code or the finiteElementHeatPropagator.c source code (both located in the folder "+MCmatlab/src"), you will need to recompile the respective mex-files. Check out those two source-files on how to do so. The source code is written in such a way that it can be compiled as either C or C++ code using either GCC, MSVC, clang or as CUDA source code using NVCC.

Model files

In MCmatlab, you set up your model in a single m-file called "Model File". You can find some examples of Model Files to get you started in the parent folder of "+MCmatlab". A complete list of parameters that can be set is provided below.

Model object

A "Model Object" stores all the parameters for a single Monte Carlo simulation: a single geometry and a single type of lightsource, together with all the results from the Monte Carlo (and optional heat transfer) simulation. The Model Object is an instance of the class MCmatlab.model, i.e., it's of the class "model" in the MATLAB-package "MCmatlab". One Model File can, if desired, generate multiple Model Objects with a different geometry and/or lightsource each.

Model Objects can be saved to and loaded from a standard binary MATLAB MAT-file like any other MATLAB variable. When loading a Model Object from a MATLAB MAT-file, the parent folder of "+MCmatlab" needs to be on your MATLAB path: prior to loading the Model Object, either navigate to the parent folder of "+MCmatlab" in the "Current Folder" view of MATLAB, or use the command "addpath" to add the parent folder of "+MCmatlab" to your path.

MATLAB function definitions

MCmatlab makes extensive use of user-defined functions, passed into the simulation in the form of function handles. For more information, see https://www.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html. Inside other functions, nested functions can be defined anywhere, but in the main model file script, MATLAB requires that all functions are defined at the end of the m file.

Media properties

The optical properties of all media are defined in a function at the end of the model file. The examples contain some media you can modify or copy-paste to your own model file, or you define your own media from scratch in your own models. All media properties functions must start with this line: mediaProperties = MCmatlab.mediumProperties; in order to initialize the mediaProperties variable as an array of mediumProperties objects. Make sure that each media in a single model file has its distinct "j"-number, as this is the number you will refer to when building your model in the geometry function (see below).

Each medium can be specified with up to 13 properties. Many of these may be specified either as a scalar or as function handles that can take wavelength, fluence rate (FR), temperature (T) and fractional damage (FD) as input:

• name - A char array containing the name of the medium, used for the visualization.
• mua - The absorption coefficient in units of 1/cm. May be (a) scalar, (b) a function handle taking wavelength as input or (c) a function handle taking wavelength, FR, T and FD as inputs.
• mus - The scattering coefficient in units of 1/cm. May be (a) scalar, (b) a function handle taking wavelength as input or (c) a function handle taking wavelength, FR, T and FD as inputs.
• g - The scattering anisotropy (the mean cosine to the scattering angle), to be used in the Henyey-Greenstein phase function. It has no units. May be (a) scalar, (b) a function handle taking wavelength as input or (c) a function handle taking wavelength, FR, T and FD as inputs.
• customPhaseFunc - An alternative to specifying g. If specified, this must be a function handle that takes wavelength and polar scattering angle as inputs and yields the scattering phase function as output. This function does not need to be normalized.
• n - The refractive index. May be (a) scalar or (b) a function handle taking wavelength as input.
• QY - The quantum yield of fluorescence, that is, photons of fluorescence emitted relative to excitation photons absorbed. It has no units. May be (a) scalar, (b) a function handle taking wavelength as input or (c) a function handle taking wavelength, FR, T and FD as inputs.
• ES - The emission spectrum of the fluorescence emission, if any. Doesn't need to be normalized. May be (a) scalar, (b) a function handle taking wavelength as input or (c) a function handle taking wavelength, FR, T and FD as inputs.
• VHC - The volumetric heat capacity of the medium in units of J/(cm^3 K). May be (a) scalar or (b) a function handle taking T and FD as inputs.
• TC - The thermal conductivity of the medium in units of W/(cm K). May be (a) scalar or (b) a function handle taking T and FD as inputs.
• E - The Arrhenius activation energy in units of J/mol. Must be scalar.
• A - The Arrhenius pre-exponential factor in units of 1/s. Must be scalar.
• nBins - The number of bins (sub-media) to split the medium into in case of media that have properties that depend on fluence rate (FR), temperature (T) or fractional damage (FD). Must be scalar.

For Monte Carlo simulations, name, mua, mus and either g or customPhaseFunc need be specified for all media.

For fluorescence Monte Carlo simulations, the QY and ES properties must additionally be specified for those media that fluoresce.

For heat simulations, VHC and TC must be specified for all media.

For heat simulations with calculation of tissue damage, E and A must be specified for those media that can be damaged.

Building the model

You build each model in a separate m-file. Each model requires the first two and optionally more of the following steps:

1. Build geometry

• Modify or add to "mediaPropertiesFunc" at the end of the model file to include the definitions of the media you're interested in. You may optionally include thermal and/or fluorescence properties.
• In your model file, specify the side lengths Lx, Ly, Lz of your model cuboid and the resolutions Nx, Ny, Nz.
• The "GeomFunc" you define at the end of the model file simply returns a 3D array containing the media definition for each voxel in the model cuboid. The media are referred to by a number, corresponding to the "j" defined in "mediaPropertiesFunc" at the end of the model file.
• You may optionally import your geometry from STL files as shown in example 18.
• After this step, when calling "plot(model,'G')", you will be shown a figure with the geometry you defined.

2. Calculate light distribution

• The section "%% Monte Carlo simulation" in the model file contains the definitions for the Monte Carlo simulation.
• After this step, when calling "plot(model,'MC')", you will be shown various figures; one with an overview of the optical properties, one with the absorbed light in the cuboid, one with the normalized fluence rate (NFR, roughly equivalent to normalized intensity or irradiance), and one showing the fluence rate at the cuboid boundaries (taking into account only the exiting photon packages, not the incident beam). The last figure will only be shown if some of your cuboid boundaries are "escaping" boundaries, i.e., boundary type 1 (all boundaries are escaping boundaries) or 2 (only the top is an escaping boundary).
• (Optional) If you set "model.MC.nExamplePaths" to an integ