1. Introduction and rationale

Modelling light propagation in absorbing and scattering media (e.g. biological tissues) finds important applications in biomedical optics, notably in diffuse optical tomography (DOT), photoacoustic tomography, and photodynamic therapy to name a few. Time-domain (TD) measurements provide for the information richest data for imaging [1–3] since the shape of the measured TD profiles directly depends on the medium’s optical properties, namely absorption and scattering. This is the motivation for considering TD light propagation here. Modelling propagation in the time domain significantly increases the computational burden compared to the continuous-wave (CW) or frequency-domain (FD) cases. It is thus important to develop efficient numerical schemes for TD modelling in terms of algorithmic simplicity, memory usage/access, and computing time. This is the general objective pursued here. This is of major importance in various applications, such as in DOT imaging of small animals or of human tissues relying on model-based iterative image reconstruction (MOBIIR) paradigms [4], whereby the light propagation model needs to be solved at each iteration of an algorithm for retrieving maps of the medium’s absorption and scattering coefficients.

The optical diffusion equation (DE) is ubiquitous in biomedical optics for modelling light propagation [5]. The DE is a partial differential equation (PDE) derived from the radiative transfer equation (RTE) using the so-called diffusion approximation [1], which assumes that the light field propagating in the medium is isotropic at each point. The DE, having no angular dependence, is simpler to solve and allows reducing the computation time, but at the cost of lesser accuracy than the RTE [6,7]. Nevertheless, the DE has proved to be a very useful model, and is adopted here for its simplicity.

Different numerical techniques exist for solving PDEs: the finite difference method (FDM), the finite element method (FEM), the finite volume method (FVM), and the boundary element method (BEM). The FEM is generally favored as it has been developed to naturally handle curved boundaries [5,8] through the use of unstructured meshes (the FEM can also be used on meshes that correspond to structured grids, see for instance [9], but that is not the typical use of the FEM). Performing 3D volume unstructured meshing, as in the FEM can be a computationally very intensive task, often carried out using third party software owing to the inherent complexity, whereas generating a structured grid such as required by the FDM is comparatively much simpler and more accessible. Further, setting up the numerical propagation equations on an unstructured mesh incurs memory and computational overhead since node connectivity information must be explicitly provided and processed [10]. The same applies to the FVM relying on unstructured meshes. The FVM can also be used in conjuction with structured meshes, but it then essentially amounts to the FDM. The BEM is also well-suited for arbitrary shapes, but it assumes that the medium is piecewise homogeneous, which is generally not the case in biological tissues, the media of focus here.

All light propagation problems in biomedical optics, and in particular in DOT, are volumetric and require to compute the light field inside a tissue volume and at its boundary. The geometries dealt with in biomedical optics are relatively simple, even if 3D and with curved boundaries. It thus makes sense to use the FDM, as it is a much simpler method to implement than the FEM and FVM using unstructured meshes, and which can handle curved geometries with the blocking-off technique as will be seen. For intricate geometries, such as the whole structure of an aircraft and its surface in aerodynamics, the FEM or FVM should be used as these methods were developed for such highly complex situations. However, for biomedical optics problems, the FDM is sufficient, and is efficient for curved boundaries; the FEM or FVM with the complexity of unstructured meshes are not absolutely needed.

The main contribution of the present work is to arrive at a complete numerical FDM-based model of light propagation in the time domain in 3D absorbing and scattering media with curved geometries, taking into account realistic refractive index mismatch boundary conditions. Such an FDM-based model has not yet been developed for the TD case to the extent considered here. One work in the past [11] has presented a DE-based TD light propagation numerical model in absorbing and scattering media (this was part of a paper on a DOT image reconstruction algorithm based on TD data). That previous work, however, used Dirichlet zero boundary conditions, which are not realistic, and which do not require computing surface normals. Also, for the 3D case, a very simple cube geometry (hence flat boundaries in 3D) was considered. That past work shall not be underestimated, but the current paper goes beyond as regards light propagation modelling.

The FDM is very easy to understand and implement, and can be made to run very fast, which is a strong point, notably for image reconstruction. The reason for this is that the FDM uses a 3D regular structured Cartesian grid for spatial discretization. Hence the equations are simple to set up and have a simple structure which allows exploiting the connectivity between neighboring grid points. By resorting to a structured grid, FDM solvers are also computationally very efficient, since the system of equations obtained with the FDM gives rise to matrices with a simple banded highly sparse structure. Such structure can be exploited by highly efficient solvers [12]. The FDM thus has highly attractive features, especially for image reconstruction.

The FDM is not well appreciated in the biomedical optics community, where curved boundaries are considered to be difficult to handle with the FDM because of the use of Cartesian grids. Such belief has been conveyed by some tutorial reviews, e.g. [13], and may have also been exacerbated by the fact that researchers in biomedical optics have been used to rely on FEM-based codes that have been made publicly available, such as NIFRAST [5] and TOAST++ [14], and that there have been numerous published research papers in which these codes have been used. This is unfortunate as the blocking-off method, less well known in the biomedical optics field, easily allows the FDM to handle curved boundaries [15]. Notably, it has been applied with success in computational fluid dynamics and heat transfer [16–19]. Note that other techniques have been developed for dealing with curved boundaries with the FDM, such as body-fitted grids [20] and block-structured grids [21]. Body-fitted grids rely on a change of coordinates for mapping the domain of interest (in the present case the biological tissue) onto a another domain with a simple shape amenable to a simple grid discretization. In practical cases involving irregular curved geometries (e.g. a small animal), such a mapping is hardly feasible, if not impossible [10]. As regards block-structured grids, which is a special case of adaptive mesh refinement, they require sophisticated automation and intricate programming to lay out the block topology in order to avoid user interaction [22].

We herein show how accessible and relatively easy it is to implement the FDM to solve a complex light propagation problem of importance in biomedical optics, and we detail the implementation steps. This should encourage biomedical optics researchers in using the FDM and developing FDM-based code for their specific needs. The FDM has proved its value in fluid dynamics, a highly demanding field, and we believe it is important to further popularize the FDM in the biomedical optics field, which the present work should contribute to. Hence, apart from its scientific contribution, this paper also serves a pedagogical purpose.

2. Overview

The TD-DE is a parabolic PDE, and the complete numerical model of light propagation in heteregenous absorbing and scattering media presented herein relies on the Crank-Nicolson scheme and blocking-off. The Crank-Nicolson scheme was chosen because it provides second-order accuracy in both time and space. As mentioned previously, blocking-off allows handling 3D media with curved boundaries. Blocking-off has been used to solve the DE or the RTE in the CW or frequency domains, which are essentially the same as regards the numerical solution, but this is the first time it is used to solve the DE in the time domain.

To handle the refractive index mismatch at the boundary of a tissue (typically an air-tissue interface), Robin-type (or mixed) boundary conditions (BCs) are considered. This requires the calculation of the normal vector at each boundary node. There exists different ways to do this for curved boundaries, notably by locally fitting an analytical expression to the boundary points to compute an expression for the normal vector. However, we introduce here a simpler way using of the Sobel operator of image processing, which is easy to implement, fast, and allows obtaining a smooth field of normal vectors along the boundary. This is the first time the Sobel operator is used to compute surface normals for BCs in the solution of PDEs with the FDM.

A computer code to implement the numerical light propagation model described herein has been developed in-house in the C++ programming language in view of its integration in a DOT MOBIIR code. Such MOBIIR code requires the modelling code to be directly incorporated into it; using black-box, already implemented, modelling codes is precluded. The C++ language was chosen as it allows optimization for fast code execution and efficient memory usage, two important criteria for MOBIIR.

The structure of the paper is as follows. First, the DE along with the Robin BCs used for TD light propagation modelling are presented in Sect. 3. This is followed in Sect. 4 with details on the numerical model (spatial discretization, Crank-Nicolson scheme, blocking-off, and boundary surface normal vectors computation with the Sobel operator). Sect. 5 discusses results obtained with the numerical model and their validation against gold-standard Monte Carlo simulations. Finally, a summary along with conclusions and an outlook of future work are given in Sect. 6.

3. Diffusion equation for time-domain light propagation in biological media

3.1 Diffusion equation

The DE is a low order approximation to the RTE [1], the latter being an integral- partial differential equation for the radiance that describes the strength and angular dependence of the light field at each point in space. The radiance depends on three spatial coordinates, two angular variables at each point in space, and time. The RTE is notoriously difficult to solve, both analytically and numerically owing to the angular dependence of the field [23], so-called the radiance. For this reason, the DE is often used instead, when applicable, namely in the case of high levels of scattering for which the field can be assumed isotropic at each point. The DE is obtained from the RTE via the diffusion approximation [1]; it is a parabolic PDE when time is considered. It is much simpler to solve than the RTE both analytically and numerically. In the DE, the light field is represented by the fluence rate, to be denoted $\psi$. The fluence rate has no angular dependance and is thus isotropic at each point. For time-domain light propagation, it depends on three spatial coordinates and time. Mathematically, the TD-DE is given by:

3.2 Boundary conditions

Appropriate boundary conditions (BCs) for the DE are necessary to adequately model light propagation in absorbing and scattering media. Realistic situations are considered here, in which a volume of biological tissue (e.g. the body of a mouse) is embedded in air. To account for the effect of refractive index mismatch at the boundary of a biological tissue, the so-called refractive index mismatch BCs are considered here. The mathematical derivation of these BCs can be found in [24]. These are Robin-type BCs, also called mixed or type 3 BCs, given by

Here $\mathbf {r'}$ is a point on the boundary, $\mathbf {\hat {n}}$ is the unit outward normal vector at $\mathbf {r'}$, $b_T$ represents a boundary source at $\mathbf {r'}$ (if any), and

Here, $\theta _c$ is the critical angle that needs to be accounted for since the refractive index of the diffusing medium is higher than that outside (air). There will thus be total internal reflection for directions of incidence with angles greater than the critical angle given by

Equations (3) and 4 are calculated numerically with the trapezoidal integration rule.

3.3 Boundary measurements

Related to the BCs is the model of measurements made at the boundary of the medium (detector readings). The measurable quantity at a position $\mathbf {r'}$ on the boundary used in imaging is the outward normal current (also called the exitance) given by [25]

Owing to the Robin BCs, and since in imaging there are no sources at points where measurements are made ($b_T(\mathbf {r'}, t) = 0$), the exitance is proportional to the fluence rate and is given more simply by

Since normalized measurements will be considered (notably in the image reconstruction algorithm to be developed in our laboratory from the work herein), the prefactor multiplying the fluence rate in the previous equation will be dropped. Hence, the fluence rate can directly be taken as the detector reading in the sequel (in any case, for validating the model later on, it does not matter whether the prefactor is accounted for or not since it is constant).

4. Light propagation numerical model

4.1 Finite difference spatial discretization scheme

To solve the DE, finite difference approximations are used to evaluate the derivatives. Considering the medium to be heterogeneous, the diffusion coefficient $D(\mathbf {r})$ is not constant. The following discretization centered scheme with half steps for evaluating the second order spatial derivatives containing the diffusion coefficient is resorted to [11,26]:

4.2 Crank-Nicolson method

Both explicit and implicit schemes have been discarded in this work. While the former is very efficient as it does not require solving large systems of equations, stability conditions impose the use of unreasonably small time steps [27]. As regards the implicit scheme, it is unconditionally stable, but is only first-order accurate in time [27]. The Crank-Nicolson method, is an average of the implicit and explicit methods. Its main advantage method is the stability for any value of $\Delta t$, as for an implicit scheme, and the fact that it is second order accurate in both time and space [26,27]. Note that although the Crank-Nicolson scheme is unconditionally stable, it may give rise to decaying oscillatory behavior in certain cases, depending on the value of the diffusion coefficient $D$, the size of the spatial discretization (represented by $\Delta x$), and the size of the time step $\Delta t$ (for details, see [28]). We have, however, never observed such oscillations in our results, and the Crank-Nicolson method is the one chosen here and used in the sequel. The finite difference scheme for the DE is then given by

Put in matrix form, these equations can be written as

4.3 Blocking-off method

In solving PDEs with the FDM, a regular Cartesian grid of nodes is used for discretizing the object. For boundaries that are aligned with the grid axes (rectangular objects), BCs are very easy to account for, and these are the situations in which the FDM is mostly used. It is often conveyed that handling curved irregular boundaries is difficult within the FDM framework [13]. This is not so, as this can be done with the blocking-off method [16–19]. Blocking-off first consists in dividing the grid nodes into three categories: nodes at the boundary (boundary nodes), nodes inside the medium (interior nodes) and nodes outside the medium (exterior nodes). Boundary and interior nodes form the active region and exterior nodes form the inactive region. For each grid node, a value is stored in a data array according to the category to which the node belongs. These values were arbitrarily chosen here to be 0 for exterior nodes (inactive), 1 for boundary nodes (active), 2 for interior nodes (active) (Fig. 1).

This array allows creating the vector $\mathbf{\Psi }$ of unknowns in Eq. (13), based on the FDM regular grid, by using only active nodes. To set up the system of equations to be evolved in time and which results in the matrix equation given in Eq. (13), the discretized DE (Eq. (11)) is used for interior nodes, whereas the discretized BC equation (Eq. (12)) is used for boundary nodes.

 

Fig. 1. Example of the application of the blocking-off method for a 2D object.

Download Full Size | PDF

The array of nodes with their respective values (0, 1, or 2) for blocking-off is obtained from a 3D binary image of the subject to be imaged. Building this binary image starts from a voxelized 3D image of the subject, for instance obtained from X-ray CT or MRI (in the case of Digimouse considered later in this paper, an X-ray CT image is used). The voxels in such images are generally easy to segment into voxels being part of the subject (subject voxels) or that are exterior to it (exterior voxels), since exterior voxels are generally very dark in such images. Then, determining whether a voxel’s value is lower or higher than a threshold allows categorizing it as an exterior or a subject voxel. This leads to a binary image, with, say, exterior voxels being set to a value of 0 and subject voxels being set to a value of 1. Subject voxels are further categorized as interior or boundary voxels. To determine whether a subject voxel at discrete position $(i, j, k)$ is an interior voxel, 6-connectivity can be resorted to. In 6-connectivity, voxels that are considered to be neighbors of voxel $(i, j, k)$ are those at positions $(i-1, j, k)$, $(i+1, j, k)$, $(i, j-1, k)$, $(i, j+1, k)$, $(i, j, k-1)$, and $(i, j, k+1)$. If all neighboring voxels are subject voxels, then voxel $(i, j, k)$ is an interior voxel, and labeled with a value of 1, otherwise it is a boundary voxel, which is then labeled with a value of 1. In the binary image, each voxel is associated to a node, with the node being considered to be at the center of the voxel. Hence, the nodes are also categorized as interior, exterior, or boundary (see Fig. 1).

4.4 Boundary surface normal vectors computation - Sobel operator

To determine the normal vector at each boundary node as required by the BC equation (Eq. (12)), the Sobel operator borrowed from image processing [29] can be resorted to and applied to the 3D binary image of the subject (see previous section). The Sobel operator efficiently computes an estimate of the gradient, which points into the direction of the normal. To compute the normal vector components at a given boundary node, the three 3D Sobel $3 \times 3 \times 3$ arrays (see Fig. 2) are convolved with the local data array of the binary image comprising the node considered and its 26 nearest neighbors. The convolution with the first Sobel array allows computing the component of the gradient along the $x$ direction at that node, and similarly for the second and third Sobel arrays for the $y$ and $z$ directions respectively. The gradient vector thus estimated is then normalized to a length of 1 as required by the boundary condition equation (Eqs. (2) and 12). This is repeated at all boundary nodes. To avoid any confusion, note that although an $xyz$ Cartesian grid is used, the normal vectors are not restricted with this technique to point along the $x$, $y$, or $z$ directions, but rather can be in any direction in space (the one direction which is computed at a given node being that of the normal at that node). This is illustrated in Fig. 3 for a 3D mouse, where it can be seen that a consistent set of smoothly varying normal vectors is obtained over the boundary surface.

 

Fig. 2. The $3 \times 3 \times 3$ array (tensor) of 3D Sobel operator $S_d$ along a given direction $d = x, y,$ or $z$, and its different layers along this direction.

Download Full Size | PDF

 

Fig. 3. Normal vectors (red needles) on the boundary surface of a mouse (blue).

Download Full Size | PDF

In past work [15], surface normals were computed with Matlab’s function surfnorm [30]. This function requires the prior generation of $X$, $Y$, and $Z$ 3D array data structures of surface coordinates (boundary surface here). Then at each $(x_b, y_b, z_b)$ point on the surface it performs a bicubic fit to the surface using that point and its neighbors on the surface to obtain a local analytical representation of the surface around $(x_b, y_b, z_b)$. This representation allows derivatives to be computed analytically to obtain two linearly independent tangent vectors at $(x_b, y_b, z_b),$ from which the normal vector can be computed using the vector cross product. Irrespective of the fact that it is implemented in Matlab, this approach is expensive as it requires generating the 3D surface array data structures, and then at each point of the surface contained in these arrays, a bicubic interpolation must be performed. With the Sobel operator, it is only needed to convolve each of the $x$, $y$, and $z$ Sobel 3D matrices at each point $(x_b, y_b, z_b)$ of the boundary and from that directly come the $x$, $y$, and $z$ components of the normal at $(x_b, y_b, z_b)$. There are no fits to be performed, and no specific data structures to generate; the approach is direct.

5. Results

The TD-DE numerical model described above was implemented in a computer code using the C++ programming language. The numerical model requires a solver for a 7-diagonal banded matrix system of equations. LU decomposition implemented in the sparseLU solver of the EIGEN library was used for this purpose [31].

5.1 Comparison with Monte Carlo

To validate the numerical model, the fluence rate obtained therewith was compared with that computed with a Monte Carlo simulation carried out with the MMCLAB software package that uses unstructured finite element meshes [32,33]. The exact same nodes were used for the numerical model as for the MMCLAB simulations, only that MMCLAB generated a tetrahedral mesh from the nodes.

Two situations were considered: a cylinder of radius 10 mm and height 60 mm (Fig. 4) and a complete 3D numerical mouse (Digimouse model [34]) (Fig. 5). The optical coefficients for the homogeneous cylinder are $\mu _a$ = 0.136 cm$^{-1}$, $\mu _s$ = 86 cm$^{-1}$, $g$ = 0.9 and $n_{in}$ = 1.37, which correspond to the optical coefficients of skin [34]. The optical coefficients for the mouse are those typical of real mouse tissues and are provided as part of the Digimouse model.

 

Fig. 4. Positions of the source S and of detectors D1 to D6 within the slice shown in red of the cylinder.

Download Full Size | PDF

 

Fig. 5. Positions of the source S and of detectors D1 to D6 within the slice shown in red of the mouse.

Download Full Size | PDF

For all simulations, $\Delta t$ = 10 ps, with 100 time steps, and a step of 0.6 mm is used in all three spatial directions. The sources used are isotropic point sources located as is customary at a distance of one reduced scattering mean free path (1/$\mu _s'$) of the boundary inside the medium [24,35] and on the same slice as the detectors, which is the configuration corresponding to that of the scanner developed in our group [36].

For each situation, the DE numerical model was solved, Monte Carlo simulations were performed, and the absolute difference between the fluence rate obtained from the DE numerical model and the Monte Carlo simulations was computed. Results are shown in Fig. 6 for the cylinder and Fig. 7 for the numerical mouse. The fluence rates obtained with the DE model and with the Monte Carlo simulations are both normalized to the same value at the time of the peak of the DE model. These results show that the DE numerical model developed herein agrees well with Monte Carlo code except at early times, where it is known that the DE is less accurate [37]. The comparison of the curves and of the absolute value of the difference between the fluence rate of the model and that of the Monte Carlo code allows to conclude that the model reproduces with acceptable error the Monte Carlo data. More quantitatively, apart from early times, the worst-case discrepancies between the numerical model and the Monte Carlo was computed to be around 6%.

 

Fig. 6. Time-domain profiles of the fluence rate for the DE numerical model (blue) and a Monte Carlo simulation (red) along with the absolute difference (green) between the DE and Monte Carlo for various detector positions on the surface of a homogeneous cylinder (Fig. 4).

Download Full Size | PDF

 

Fig. 7. Time-domain profiles of the fluence rate for the DE numerical model (blue) and a Monte Carlo simulation (red) along with the absolute difference (green) between the DE and Monte Carlo for various detector positions on the surface of Digimouse (Fig. 4).

Download Full Size | PDF

 

Fig. 8. Positions of slices A (axial) and B (coronal) in the numerical mouse.

Download Full Size | PDF

The maps of the fluence rate for the whole numerical mouse obtained respectively with the DE and MMCLAB at a given time (and this for different times) were also compared for different slices, and the amplitude of the absolute difference between the two was calculated. Figures 9 and 10 show such results for a time $t$ = 0.4 ns. Both agree extremely well with negligible absolute error ($\approx 10^{-5}$).

 

Fig. 9. Fluence rate map at $t$ = 0.4 ns in slice A for DE and MMCLAB and absolute difference (Fig. 8).

Download Full Size | PDF

 

Fig. 10. Fluence rate map at $t$ = 0.4 ns in slice B for DE and MMCLAB and absolute difference (Fig. 8).

Download Full Size | PDF

As an indication of the computing time, it takes 21 s with the algorithm ran on an Intel Core i7-4790K processor to model the light propagation in the complete 3D numerical mouse over 1000 ps with a time step of 10 ps when only one thread is used, and this for a very fine grid (0.6 mm spatial resolution in all directions). This is very fast considering the fine grid used, that this time includes setting-up the set of equations, and the fact that full time-domain simulations are carried out.

5.2 Influence of absorbing inclusion

As a final experiment, a situation showing the influence of an absorbing inclusion on the fluence rate is considered. The purpose of this is not to provide a quantitative assessment, but rather to illustrate the behavior of the fluence rate and see that it qualitatively corresponds to what is expected with the numerical code developed. The medium considered is the same cylinder as in the previous section with the same parameters, except that an absorbing inclusion with $\mu _a$ = 1.36 cm$^{-1}$ was added (Fig. 11). Maps of the fluence rate for the cylinder with and without the inclusion are depicted for two different slices in Figs. 12 and 13. As expected, it is clearly seen that the fluence rate is more attenuated at the inclusion’s site, while the propagation has a spherical shape in the homogeneous case.

 

Fig. 11. Positions of slices A and B within the cylinder with the absorbing inclusion depicted in white.

Download Full Size | PDF

 

Fig. 12. Fluence rate map at $t$ = 0.4 ns in slice A for DE with and without an absorbing inclusion (Fig. 11).

Download Full Size | PDF

 

Fig. 13. Fluence rate map at t = 0.4 ns in slice B for DE with and without absorbing inclusion (Fig. 11).

Download Full Size | PDF

6. Summary, conclusions, and outlook

This paper presented an efficient finite difference method for solving numerically the time-domain optical diffusion equation for heterogeneous media with curved boundaries, taking into account the refractive index mismatch at the boundary. The method relies on blocking-off to enable the finite difference scheme to deal with irregular curved boundaries. To account for realistic refractive index mismatch boundary conditions, which require the surface normal at each boundary point to be computed, the 3D Sobel operator is resorted to for the first time in the numerical solution of partial differential equations. This operator provides a simple and efficient way to estimate surface normals in a smooth way. A comparison with gold standard Monte Carlo simulations confirmed the accuracy of the numerical model to within the limitations of the diffusion equation. This is, to be best of the authors’ knowledge, the first time blocking-off is applied to solve this equation in the time domain using realistic boundary conditions for media with curved boundaries.

The present work focused on the development of an efficient numerical model in view of its use as the forward model in a DOT model-based iterative image reconstruction algorithm that will be presented in upcoming work. Such an algorithm requires several runs of forward modelling owing to its iterative nature, and because a plurality of illumination positions and associated measurement projections need to be considered at each iteration of the algorithm. It was strived here to resort to efficient methods in all aspects (finite difference method, blocking-off, Sobel operator) in view of performing image reconstruction in reasonable time and to make such an algorithm tractable for 3D imaging with time-domain data. It is again stressed that here the light propagation modelling equation was solved on regular grids containing equidistant points. The finite difference method combined with blocking-off makes possible the use of such structured grids for modelling irregular curved geometries. Structured grid calculations are highly efficient since data storage and access to neighboring grid points of regular or structured grids is direct. Neighboring grid points of regular grids are simply found in a 3D data array by adding or subtracting 1 from grid point indices. Irregular grids, however, require storage of cell-to-cell pointers in a one-dimensional data array, which provides a connectivity list between neighboring cells, leading to more storage and additional code instructions. Another appealing feature of the finite difference method with its use of structured grids is that it leads to simple banded matrices whose structure can be exploited by highly efficient numerical solvers for linear systems of equations [12,31], whereas this is more difficult for the sparse, but nevertheless denser, matrices obtained when unstructured meshes are used [13], such as in the finite element method and variants of the finite volume method. The approach presented here is efficient, and can be made to run fast. This is significant, in particular in the context of the use of numerical light propagation models within model-based iterative image reconstruction paradigms. Further, results show that the approach developed here is accurate. It showed good quantitative agreement with gold-standard Monte Carlo simulations. Discrepancies of around 6% were observed, which is excellent. Such discrepancies may be caused by the approximate nature of the diffusion equation, whereas Monte Carlo simulates the more precise radiative transfer equation.

It is hoped that the present work, will contribute to further popularize the use of finite difference methods in biomedical optics. It should now be clear that any type of boundaries can be handled rather easily with the finite difference method when used in conjuction with blocking-off and the Sobel operator. As a final note, the numerical approach developed here for solving the diffusion equation can obviously be adapted with no conceptual changes to more sophisticated models, such as the $SP_N$ approximation to the radiative transfer equation (RTE) [23,38,39], the RTE itself, or other PDEs of mathematical physics.