1. Introduction

Monte Carlo method for photon transport (hereafter referred to as ’Monte Carlo’) is used for simulating how photons are scattered and absorbed within a medium. It has been used to simulate light propagation in tissue for a variety of biomedical applications, e. g. near-infrared spectroscopy, diffuse optical imaging, photoacoustic tomography and light-based therapies [1,2]. A key feature of the method, as the name Monte Carlo implies, is utilisation of random numbers in propagating photons in a stochastic manner. Due to its reliability, Monte Carlo has been regarded as the ’gold standard’ for simulating light transport in biological tissue. It is generally used as the reference approach for deterministic methods which are typically analytical or numerical solutions of the radiative transfer equation or its approximations such as the diffusion approximation.

Numerous implementations of the Monte Carlo algorithm exist [1,2]. One of the first widely utilised Monte Carlo method in biomedical optical imaging was the photon packet method, introduced for biomedical optics by Prahl et al. [3]. In photon packet method, the efficiency of a Monte Carlo program is increased by propagating many photons, i.e. a photon packet, along each pathway. Instead of considering the absorption of a single photon as a stochastic event, it is treated by reducing the weight (size) of the photon packet. Various implementations of the photon packet method followed.

In the early implementations, photon fluence was computed in geometries with a predefined shape such as layered spheres [4] and multi-layered slab geometries [5]. Later, implementations in other predefined shapes [6,7] and in more general complex geometries followed [8–10]. Typically in complex geometries, Monte Carlo implementations have been implemented using two-dimensional (2D) pixel or three-dimensional (3D) voxel discretisations [8–15]. It has been shown that approximating boundaries in voxel basis may lead to large errors in the solution of light transport simulation [16]. Therefore, Monte Carlo implementations in unstructured meshes consisting of triangles or tetrahedra have been constructed for better and more accurate solutions in complex heterogeneous geometries [17–24]. In addition to the discretisation, Monte Carlo implementations differ also depending on how they sample scattering length and reduce the weight of the photon packet. For more information on the comparison of different Monte Carlo implementations see e. g. articles [1,25] where it has also been shown that these different implementations are statistically equivalent.

During the past decades, the huge development of computers has enabled parallised and Graphics Processing Unit (GPU) advanced Monte Carlo implementations that can be used to obtain solutions to complex three-dimensional simulation problems in a time-frame from few seconds to a few minutes [12,26–30]. Thus, unlike at the early 1990s, when the available computational resources limited stochastic simulation methods such as Monte Carlo as a reference approach which took hours or even weeks to run, the current Monte Carlo implementations can be easily utilised for various problems. This has made them a tempting alternative for other type of problems as well. Monte Carlo simulations have been used, for example, in studying photoacoustics [31,32], fluorescence tomography [33,34], time resolved diffuse spectroscopy [35,36], statistical models for photon transport [37], biophotonics of corals [38], dose optimization for photodynamic therapy [39] and laser treatment of port-wine stains [40]. Furthermore, the first studies in which Monte Carlo is utilised in the solution of the inverse problem of an optical imaging problem have been implemented [15,41–46].

Although various implementations of Monte Carlo algorithm exists, there is only a limited number of free software available, see e. g. [5,12,14,17,22,29]. In this paper, our Monte Carlo software ValoMC is described. The geometry is defined using unstructured meshes, i. e. triangles in 2D and tetrahedra in 3D. The software consists of the main simulation code, written in C++, and auxiliary MATLAB (The MathWorks Inc., Natick, MA) functions. These functions serve as an easy to use interface for anyone who wants to run Monte Carlo simulations for their problem but are not comfortable with technical programming. The simulation code is based on our earlier Monte Carlo software that has previously been utilised in simulating light propagation for example in [20,23,47–51]. A summary of other Monte Carlo simulation software that have MATLAB interface can be found in Ref. [52].

The rest of the paper is organised as follows. First, the theory is shortly reviewed in Section 2. Then, program description and usage guidelines for the MATLAB interface are given in Section 3, followed by usage examples and validation studies in Section 4. Software and hardware specifications and availability are discussed in Section 5. A summary is given in Section 6.

2. Theory and computational methods

Four parameters are used to represent the optical properties of medium [53]: absorption coefficient $\mu _a$, scattering coefficient $\mu _s$, scattering anisotropy in the Henyey-Greenstein phase function $g$ and refractive index $n$, and are all function of the position within the medium. Absorption $\mu _a$ and scattering $\mu _s$ are defined so that the probability for the photon to be absorbed or scattered within an infinitesimal length ${\textrm {d}}x$ along the propagation direction are $\mu _a {\textrm {d}}x$ and $\mu _s {\textrm {d}}x$, respectively, and assumed to be independent. Further, the Henyey-Greenstein phase function represents the directionality of the scattering [54]. In contrast to e. g. Rayleigh and Mie phase functions, it is not derived from a physical theory. Rather it simply serves as mathematical representation of the scattering medium. For example, when $g = 0$ the scattering is isotropic and when it is close to 1, the photon scattering angles become small. Its numerical value can be determined using experimental data or fits to theoretical considerations. Henyey-Greenstein is the most commonly used phase function for biological tissue but other functions are used as well [55]. The refractive index determines the amount of light that is reflected or refracted on tissue interface. In this work, the spatial dependence of all the parameters is described using a triangular or tetrahedral mesh, and considering them as piecewise constant values.

2.1 Geometry description

A coordinate matrix $r$, a topology matrix $H$ and a boundary topology matrix $B_H$ represent a triangular or a tetrahedral mesh. Examples of such 2D and 3D meshes are shown in Fig. 1. Although we consider here very simple shapes for illustration purposes, general complicated shapes can be constructed and utilised in ValoMC. Let us consider a simple example geometry shown in Fig. 1(b). To set up a quadrilateral region in 2D that consists of two triangles, the following coordinate matrix can be used

 

Fig. 1. Illustration of the matrices $r$ (coordinates), $H$ (element topology) and $B_H$ (boundary element topology). (a) Shown is a 2D circle constructed from triangles. (b) The geometry for the simple example discussed in the text. (c) 3D half-sphere constructed from tetrahedrons. (d) A tetrahedron $H_1 = (\,1\enspace 3\enspace 5\enspace 2\,)$ indicated with coloured region and a boundary element ${B_H}_1 = (\,1\enspace 3\enspace \,5)$ indicated with a dashed line.

Download Full Size | PDF

A neighbourhood matrix $H_{N}$ is used to make finding the neighbouring elements and boundary elements for a given element more efficient. This matrix is only used internally and the description is provided here for completeness. Along the same lines, each row in $H_{N}$ refers to a single element and each element to a neighbouring element.

Negative indices are used to distinguish boundary elements from regular elements. The neighbourhood matrix for the current example is

The generalization to a cube in 3D is straightforward, as visualised in Fig. 1. In 3D, the fundamental element is a tetrahedron. The four corners of the tetrahedron are similarly defined using a four-columned matrix $H$ while the boundary elements are triangles. They can be set up using a three column $B_H$ matrix which has the same structure as $H$ in 2D.

2.2 Algorithm

The flowchart of the algorithm is shown in Fig. 2. Our implementation follows Ref. [3] with additions due to the mesh based geometry and difference in path weighting. Random numbers are generated using the Mersenne-Twister algorithm [56].

 

Fig. 2. Flowchart of the algorithm.

Download Full Size | PDF

The algorithm starts with the creation of a photon packet (Fig. 2(a)). A random location from those boundary elements that act as a light source is selected. Further, the propagation direction (angle) is drawn from a probability density distribution that can be defined for each boundary element independently (see Section 3). The photon packet is initialised with a given weight $w=1$. This number describes the relative number of photons in the packet.

After the photon packet is initiated, the program enters in the main loop that ends in the termination of the packet. In the loop, a random scattering length $s$, is drawn so that the drawn lengths follow an exponential probability density distribution

The scattering length $s$ is exhausted by moving the packet at a straight line either to the interface of the next element or, if the distance to the interface is greater than the remaining length, to the end of the distance (Fig. 2(c,d)). While the packet moves within the medium, it is affected by the optical coefficients in each element on its path. If the scattering coefficient changes along the path, the remaining length is adjusted to account for the change by scaling the remaining length by $\mu _{s,\textrm {curr}} / \mu _{s,\textrm {next}}$ where $\mu _{s,\textrm {next}}$ refers to the scattering coefficient of the element the packet is entering. At the interface of two elements, the photon packet may undergo Fresnel reflection or be transmitted to the neighbouring element with a modified propagation direction due to a change in the refractive index (Fig. 2(e)). Throughout the propagation, the weight of the photon packet is reduced to

After the photon packet has propagated the distance $s$, a new scattering event occurs. In the scattering event, the packet obtains a new direction that is drawn so that the probability of the direction follows Henyey-Greenstein phase function. This is achieved by using the inversion method [57], i.e. uniformly distributed random numbers in the inverse of the Henyey-Greenstein phase function. The current implementation draws the scattering angle in 2D and the cosine of the scattering angle in 3D (similarly as in Ref. [3]). Then, a new scattering length is drawn, and the main loop of the simulation continues.

The photon packet is terminated if the packet exits the computation domain. In addition, the packet can be terminated in a procedure known as ’survival roulette’ [3] before drawing a new scattering length. The survival roulette is performed in order to avoid the computation of excessively long propagation paths with a low influence on the computation results. Once the photon packet weight has decreased below a threshold value (by default $10^{-3}$), it is either terminated with a probability of $1 - 1/m$ (by default $m=10$) or its weight is promoted to $mw$ with probability $1/m$. Thereby, the total weight of all photon packets is well conserved despite the terminations.

The main outcome of the simulation is the fluence distribution in the elements of the computation domain and the exitance at the boundary elements. As discussed previously, an initial weight of one is assigned to each packet, and throughout the propagation, the weight of the packet is reduced. Conversely, a weight of $w (1-\exp (-\mu _a \Delta s))$ is added to the element through which the packet propagates. When the packet reaches a boundary element and exits the computation domain (i.e. does not undergo Fresnel reflection at the boundary, c.f. Fig. 2(f)), the remaining weight is added to the boundary element. Since a photon packet has an initial weight of one, the total weight that is deposited in the computational domain (including its boundary) is equal to the number of photon packets $N$. The piecewise constant values obtained in this way represent the weight $W_i$ absorbed into the element $i$ and the weight $W_j$ escaped through boundary element $j$, per unit time. At the end of the simulation, the fluence $\Phi _i$ at element $i$ is obtained from

2.3 Frequency domain calculation

Frequency domain Monte Carlo (i.e. intensity modulated light sources) are implemented by using a complex weight for the photon packet [2]. This approach is motivated by noticing that the radiative transfer equation in the frequency domain can be interpreted as time-independent one but with a complex absorption coefficient $\mu _a + \frac {\mathrm {i} \omega }{c}$ where $\omega$ is the angular frequency of the intensity modulation, $c$ is the speed of light in the medium and $\mathrm {i}$ is the imaginary unit. Therefore, given an initial weight $w$ and a distance $\Delta s$ propagated within an element, the new weight is computed according to

3. Program description

The program consists of two parts: the simulation code (written in C++) and the MATLAB interface. While the C++ program can be used independently without MATLAB, this section focuses on describing the MATLAB interface part that is intended to facilitate the initialisation of the simulations and the visualisation of the results. Documentation for the C++ code can be found in ”cpp/README”.

A summary of the currently implemented features of the software is given Table 1. It should be noted ValoMC does not provide a fixed pipeline from meshing to plotting the results. It focuses only on the simulation part. Usage examples on how to construct a custom pipeline are provided. ValoMC stores the mesh and the results in a format that is directly compatible with the plotting tools in MATLAB. Functions to import meshes and export results to specialized software are, in addition, provided. As an output, the computation returns the exitance at each boundary element and the fluence at each element in the computation domain. Therefore, in the current version, detectors (or light sources) with a more complicated spatial or angular dependence can be implemented by making the mesh finer or extending the mesh. Spatial sensitivity profiles of detectors can be implemented by weighting the boundary solution from a finer mesh in post-processing. Angular sensitivity can be implemented by creating vacuum extensions ($\mu _s = 0$, $\mu _a = 0$) to the mesh at the detector locations. Since photons propagate in straight lines in vacuum, the extension can be adjusted so that the boundary elements at the end of the extension collect photons at desired angles only.

Table 1. A summary of currently implemented features in ValoMC.

View Table | View all tables in this article

In MATLAB, the simulations are initiated with the command

solution = ValoMC (vmcmesh, vmcmedium, vmcboundary, vmcoptions)

The MATLAB interface consists of auxiliary functions that help to set up the input structures. These structures are documented in Table 2 while a more detailed description with additional features can be found from the homepage. The fields of vmcmesh are the same matrices as described in Section 2.1, vmcmedium contains the optical coefficients for each element and vmcboundary the boundary conditions for each boundary element. A typical modelling task consists of creating a suitable mesh and setting up the optical coefficients of each element in vmcmedium and light sources or other boundary conditions of each boundary element in vmcboundary. The fields of vmcoptions can be used to set a number of global options, for example, to adjust the total number of photon packets launched ($10^6$ by default) or to enable a frequency domain calculation. If the options structure is not provided, default values will be used.

Table 2. List of the mandatory fields in the input and output of the main function of ValoMC. More detailed options are provided in the homepage and documentation.

View Table | View all tables in this article

Simulations can also be run outside MATLAB by using the functions exportValoMC and importValoMC. These functions write an input file for the C++ program or read its output for analysis in MATLAB, respectively, and enable computations, for example, on a computing cluster without MATLAB.

3.1 Creating a mesh, a medium, and a boundary

ValoMC comes with a few convenient functions to create simple meshes for prototyping purposes, such as createRectangularMesh. For creating more complicated meshes, multiple dedicated freeware tools are available, such as iso2mesh [58] and NetGen [59]. It is possible to import NetGen ’.vol’ ASCII files using importNetGenMesh. The importer returns the region and boundary definitions in addition to the mesh. Commercial software are available [60] that produce meshes in a format easily convertible for ValoMC (c.f. Table 2) within MATLAB. It should be noted that meshing tools often introduce multiply defined boundaries between different regions of the mesh (i.e. boundaries between connected elements) as a part of the boundary definition. Currently, the overlapping parts must be removed before using them as $B_H$. Alternatively, the function createBH can be used to create a $B_H$ matrix for a given $H$ matrix and overwrite the boundary given by the meshing tool. Furthermore, ValoMC has been made compatible to mesh formats used in some biomedical optics simulation software such as Toast++ [61] and k-Wave [62].

Function createGridMesh is provided for working with pixel or voxel format data. It creates a structured mesh from a grid so that the mesh pixels or voxels can be filled by two triangles or six tetrahedrons, respectively. When the mesh has been created with createGridMesh, the fields of vmcmedium can be introduced as 2D or 3D arrays. Furthermore, the resulting fluence field is returned as a 2D/3D array in solution.grid_fluence.

Once the mesh has been created, functions createMedium and createBoundary can be used to create an accompanying vmcmedium (used to set material parameters) and vmcboundary (used to set boundary conditions) structure, respectively. They take vmcmesh as an input. If createMedium is given a vmcmedium as the second input parameter, it will repeat the arrays so that they match in size with the mesh. A convenient way to initialise all elements to uniform values is to give a structure with scalar optical coefficients as an input. createMedium will set the refractive index at the boundary so that there is no mismatch at the boundary (see Section 3.3).

3.2 Finding boundaries and regions

All material parameters and boundary conditions are set individually for all elements and boundary elements, respectively. However, it is often convenient to group elements for setting material parameters and boundary conditions using a single instruction. This is done by creating vectors that hold indices to $B_H$ (boundary conditions) or to $H$ (material parameters). These matrices were introduced in Fig. 1(b). If such vectors are not provided by the meshing tool, functions findBoundaries and findElements can be used. These functions take the mesh as an input, a keyword such as ’circle’ and optional arguments. A list of keywords is provided in the code documentation and by invoking ’help findElements’. For example, in 2D all elements within a circle of $r < 10 \, \textrm {mm}$ with a centre at the origin $(0 \, \, 0)$ can be retrieved by invoking

circle = findElements (vmcmesh, ’circle’, [0 0], 10);

The return value, circle, is a vector that holds indices to vmcmesh.H. These indices could be now used to set e.g. absorption coefficient for the elements within the circle by

vmcmedium.absorption_coefficient (circle) = 0.02;

Note that circle can be a single number as well to set the properties of an individual element. Similarly, findBoundaries returns indices to vmcmesh.BH, i.e. boundary elements. Using command

lightsource1 = findBoundaries (vmcmesh, ’direction’,…[0 0], [5 0], 1); will return boundary elements that fall within a rectangle defined by a line from (0 0) to (5 0) that intesects it in the middle and has a width of 1.

3.3 Light sources and other boundary conditions

Light emitting boundary elements in vmcboundary are defined by a cell array lightsource. Each element of the array is a string that determines the type of the light source. For example, vmcboundary.lightsource (lightsource1) = {’direct’}; will set a light source with a type ’direct’ on each boundary element whose index is contained in vector lightsource1. All photons in ’direct’ light sources are launched in the same direction (by default, in the direction of the normal).

Light sources can be directed using the field lightsource_direction. This is a $N_b \times 2$ / $N_b \times 3$ size matrix where each row contains a direction vector ($N_b$ is the number of boundary elements). The direction can be given either as an absolute direction or relative to the normal. The light source types affect the distribution of initial propagation angles with respect to the direction vector. For example, in a {’cosinic’} light source initial angles are random but follow a cosine distribution and in a {’gaussian’}, a Gaussian.

By default, the initial position of the photon packet is selected from all light emitting boundary elements with a uniform probability. This means that all boundary elements that contain a light source have an equal irradiance. Alternatively, the relative strength of each light source element can be determined using the field lightsource_irradiance. This field adjusts the probability for a particular boundary element to be selected. For example, if the value of the field on a single boundary element is set to $2$ and the rest of the elements are $1$, the element will emit twice as many photons per unit length than the rest of the elements.

The refractive index at the boundary of the computation domain (i.e. the refractive index of the material surrounding the computation domain) can be set for each boundary element individually using the field exterior_refractive_index. Only the photon packets that are exiting the computation domain can be reflected at the boundary. In other words, light sources are considered to lie inside the computation domain (a convention that varies among different software).

3.4 Simulation output and visualisation

The solution of a simulation (fluence at each element and exitance at each boundary element) is returned in solution.element_fluence and solution.boundary_fluence. Visualisation of the solution arrays is straightforward with MATLAB built-in functions such as patch, trimesh and tetramesh. Various usage examples are provided with the code. Visualisation of large tetrahedral meshes can be inefficient using tetramesh. For this purpose, the function exportX3D is provided. It converts tetrahedral meshes into triangular meshes so that each four triangular face of tetrahedron will have the same value. The resulting X3D file can be opened and post-processed for example in MeshLab [63]. The triangular representation is also convenient for making even cuts to the mesh for visualisation purposes, e. g. using the Geometry Processing Toolbox [64] for MATLAB.

4. Simulations

The simulation code has been validated against analytical and numerical solutions of the radiative transfer equation and other Monte Carlo software. Below, three example studies are described. We have used

4.1 Validation of the 2D code

The results of the simulation code were compared against analytical solution of the radiative transfer equation in a 2D semi-infinite layered media. The analytical solution is presented in the article by A. Liemert and A. Kienle [65]. The geometry of the problem and the optical parameters are illustrated in the inset of Fig. 3 (see also Figs. 1 and 5 of Ref. [65]). The slab consisted of two layers with different optical coefficients: the first layer was $5 \, \textrm {mm}$ thick while the second was infinite. The optical properties of the medium were $\mu _a = 0.08 \, \textrm {mm}^{-1}$, $\mu _s = 5 \, \textrm {mm}^{-1}$, g = 0.9 within the first layer and $\mu _a = 0.01 \, \textrm {mm}^{-1}$, $\mu _s = 12 \, \textrm {mm}^{-1}$, g = 0.9 within the second layer. A beam with a spatially Gaussian profile was incident on the slab with an incident angle $\theta _0$ which had the values $0$, $45$ and $70$ degrees.

 

Fig. 3. Comparison of ValoMC against the analytical solution of the radiative transfer equation in a semi-infinte layered medium. The unit of exitance ($1/\textrm {mm}$) is used for consistency with Ref. [65]. The data for the analytical solution was provided by A. Liemert and A. Kienle [65].

Download Full Size | PDF

The incident power per unit length through a line that is perpedicular to $\theta _0$ is $S(y_{\perp })$ where S is a Gaussian function (width $\sigma = 0.4\, \textrm {mm}$) and $y_{\perp }$ is a position along that line. Therefore, we computed the strength of the light source (given by the field lightsource_irradiance) at $y$ as $S(\cos (\theta _0) y)) \cos (\theta _0)$, where $y$ is a position of the boundary element at the edge (x=0) of the domain. This expression follows from the geometry and preserves the beam profile in the perpendicular direction when tilted. Note that the light source type was {’direct’} and not {’gaussian’} since all photons have the same initial propagation direction.

Since infinite domains are not supported in ValoMC, the finite size error was made small by creating a suitably large computation domain with thickness of $55 \,\textrm {mm}$ and width of $200 \, \textrm {mm}$. Element width of the rectangular mesh (createGridMesh) was $0.1 \, \textrm {mm}$ in y-direction and $1 \, \textrm {mm}$ in x-direction. $10^8$ photon packets were launched into the medium. The relative errors of the Monte Carlo exitance computed according to Eqn. (10) in the domain $y = -5 \cdots 5\;\textrm {mm}$ were $E (0^{\circ })\approx 0.46\;\%$, $E (45^{\circ })\approx 0.35\;\%$, $E (70^{\circ })\approx 0.17\;\%$. The resulting exitance, given by the field boundary_exitance is shown in Fig. 3. The simulations demonstrate an excellent agreement with the analytical solution.

4.2 Validation of the 3D code

The performance of ValoMC was tested in a 3D cubic simulation domain and compared against MMC (mesh-based Monte Carlo) software by Q. Fang et al. [19]. The geometry of the calculation is illustrated in Fig. 4(a) (see also Fig. 3 of Ref. [19]). The domain consisted of a cube as a background (edge length $60 \, \textrm {mm}$, $\mu _a = 0.005 \, \textrm {mm}^{-1}$, $\mu _s = 1 \, \textrm {mm}^{-1}$, $g = 0.01$, $n = 1.37$) with a sphere inside of the cube (radius $10 \, \textrm {mm}$, $\mu _a = 0.05 \, \textrm {mm}^{-1}$, $\mu _s = 5 \, \textrm {mm}^{-1}$, $g = 0.9$, $n = 1.37$). A pencil beam was incident from the middle of the bottom face. For compatibility of the two results, all MMC simulations were run on a piecewise constant basis and performed using a large time window and step ($0.5 \, \textrm {s}$). The mesh had 11853 nodes and 70749 elements and it can be found as a part of the distribution of MMC.

 

Fig. 4. Comparison of ValoMC against Mesh-based Monte Carlo (MMC) [19]. (a) Fluence distribution when the refractive index of the sphere was set to $n = 1.0$ calculated with ValoMC. (b) Difference (relative error) between fluence obtained using ValoMC and MMC with different number of photon packets calculated using Eqn. (10). (c) Total energy absorption rate in the sphere as a function of its refractive index calculated with ValoMC and MMC.

Download Full Size | PDF

The results were in a good agreement with each other, showing a maximal difference of $1.16 \times 10^{-4} \, \textrm {Wmm}^{-2}$ in fluence after $10^8$ photon packets were launched into the domain. Differences in the results calculated according to Eqn. (10) over the entire computation domain were $1.57~\%$, $0.16~\%$ and $0.05~\%$ after $10^5$, $10^7$ and $10^8$ packets launched in to the medium, respectively. A plot of the difference is shown in Fig. 4(b). A comparison of the computation speeds is given in Section 5.2.

To validate the computation of internal reflections due to a refractive index mismatch, we varied the refractive index of the sphere. The refractive index of the background medium and at the boundary of the computation domain were kept unchanged. The total energy absorption rate into the sphere was computed with both codes using $10^7$ photon packets for each refractive index and is shown in Fig. 4(c). Changing the refractive index did not affect the difference between the results (0.16 $\pm$ 0.1 % for all computations).

4.3 Frequency domain calculation

The use of intensity modulated light sources was simulated following an example given in the book by L.V. Wang [66]. Two light sources, $180$ degrees out of phase, were placed adjacent to each other (see Fig. 5) and the resulting amplitude was measured on a line further away. Due to the symmetry of the problem, the photon fluence from the two light sources should be at an opposite phase on a line in half-way between the light sources, which cancels the amplitude on that line.

 

Fig. 5. Frequency domain calculation in 2D. Coloured image shows the phase of the photon fluence. The blue line shows the amplitude of the exitance.

Download Full Size | PDF

Since the options related to the intensity modulation affect all light sources in ValoMC, two simulation runs were performed to demonstrate the effect. To activate a sinusoidally modulated light source, the field options.frequency was set to $200 \, \textrm {MHz}$ in both simulations. Using options.phase0, the initial phase of the first source was set to $\pi /2$ in the first simulation and the initial phase of the second source was set to $-\pi /2$ in the second simulation. The total fluence and exitance were formed by summing the individual results from the two simulations. The results are shown in Fig. 5.

5. Software and hardware specifications and availibilty

The simulation code is written in C++ and parallelized using OpenMP. A 64-bit platform is required to run the software. ValoMC has been tested on the platforms listed in Table 3. The software is available for download from the homepage https://inverselight.github.io/ValoMC/ and the source code from GitHub https://github.com/InverseLight/ValoMC. The project is distributed under the MIT-license, excluding the external components. LICENCE file is provided for detailed information. Several examples demonstrating the various features of the code are included.

Table 3. List of the environments ValoMC has been tested on.

View Table | View all tables in this article

5.1 Installation

Installation instructions are provided in the homepage and in the README.md file. Currently, binaries are not distributed and the program must be compiled before use. We have distributed the program to a small number of test users and found that most problems in installation arose due to MATLAB not supporting the version of the compiler that was installed. These problems are avoided e.g. by using CMake [67] to compile the program and changing the default compiler (instructions provided in the README files).

5.2 Performance

On our Linux test system (Ubuntu GCC 8.1.0, Intel(R) Core(TM) i7-6500U processor, using 4 threads), the computation times of the simulations described in Section 4.2 were somewhat slower than with MMC but similar in magnitude. ValoMC had the computation rate of about 55 000 packets/second, whereas for MMC it was 99 000 packets/second under comparable conditions.

In the current implementation, ValoMC sacrifices some of the optimisations (e.g. by using double precision numbers and standard math functions, whereas MMC uses single precision and custom math functions) for numerical accuracy. Furthermore, ValoMC (unlike MMC) stores by default the boundary solution. On the other hand, MMC performs a time domain calculation. The compiler switches were those reported optimal for MMC in its user manual. All simulations were performed with the stand-alone executable and not in MATLAB. Thus, although while we have made efforts to make the comparison as fair as possible, some differences in the computations remain. The results suggest no significant gain in accuracy by using double precision numbers and standard math functions.

Furthermore, we found that on our Windows test system, the ValoMC binary compiled with Visual Studio 2017 compiler (optimized for speed in the settings) performed considerably better than the GCC 8 compiler without losing accuracy ($\sim 1.25 \times$), while we were unable to compile MMC with Visual Studio to perform the comparison also with that. For users wanting to achieve the maximal performance, we therefore recommend trying out multiple compilers. This is straightforward with CMake.

6 Discussion and conclusions

We have developed an open source Monte Carlo software for studying photon transport. While similar programs exist and are readily available, the code packages may require substantial learning efforts. To address this issue, we have also developed a simple MATLAB interface to the simulation code and provide multiple examples demonstrating the features of the software. Functions for facilitating the use of the software with certain biomedical optics simulation software such as Toast++ [61] and k-Wave [62] and meshing software such as NetGen [59] are provided. The simulation code and the interface are suitable for fast prototyping as well as setting up large-scale problems. Near-term feature extension plans include time domain computations, internal light sources (for modelling e.g. fluorescence), coordinate based detectors and optimizations. In particular, we consider that there is still considerable optimization potential left in e.g. improving cache performance, utilising the vector instructions on CPU and math libraries, while GPU implementations are also in consideration.

With this project, we aim to bring high-performance Monte Carlo modelling of light transport accessible to a wide range of researchers, teachers and students, including those not previously familiar with the method or lacking background in low-level programming languages. We vision ValoMC as a lightweight, efficient and flexible tool that is developed according to the principles of open source to suit increasingly diverse uses in biomedical and other applications.