New nonlocal forward model for diffuse optical tomography

Wenqi Lu; Wenqi Lu; Jinming Duan; Joshua Deepak Veesa; Iain B Styles; Iain B Styles

doi:10.1364/BOE.10.006227

1. Introduction

In diffuse optical tomography (DOT), near-infrared light (650-900 nm) is injected into an object through optical fibers placed on its surface. The light is injected through each fibre in turn and propagates through the object. The spatial distribution of light remitted from the object’s surface is measured for each source fibre, and this information is used to estimate the object’s internal optical properties by iteratively refining the optical properties of a forward model of light propagation in the object until the model predictions match the measured surface remittance. As such, the forward model of light propagation must be able to accurately model the main interactions (i.e. absorption and scattering) between light and the object so as to recover internal properties faithfully.

When its wave nature is neglected and light is interpreted as a stream of particles (photons), the main interactions between light and biological tissue are characterized as absorption and scattering and are modelled by the radiative transfer equation (RTE) which is generally accepted to accurately describe how light propagates in biological tissues [1,2]. Analytical solutions exist for the RTE only for simple geometries with nearly homogeneous interior structure [3,4]. Although a number of algorithms exist to find numerical solutions for more complex inhomogeneous domains [5,6], they are extremely computationally expensive, especially for large 3D volumes. Due to the complexity of the RTE and the limitations of existing algorithms, a range of stochastic or deterministic approximation schemes are frequently adopted to simplify the RTE. The Monte Carlo method is the most commonly used stochastic model [7,8]. It is however costly in computational time, because millions of photons need to be tracked to acquire meaningful statistics. Light propagation can also be modelled by a deterministic diffusion equation (DE) using the $P_1$ approximation of the RTE. It is based on the assumption that the radiance in an optical medium in which multiple scattering occurs is almost isotropic, and scattering in that medium is dominant over absorption. The modified Beer–Lambert law (MBLL) is also sometimes used to model thick tissues [9,10]. However Boas et al. [11] observed that standard MBLL analysis cannot accurately quantify relative changes in the concentration of chromophores. Bhatt et al. [12] proposed a generalized Beer–Lambert model to overcome this limitation and have applied this method widely to near-infrared spectroscopy (NIRS) studies.

In DOT, technically, such interactions can be accurately described by a diffusion equation (DE) which is derived from the radiative transfer equation (RTE) [1] under the assumption that the radiance in an optical medium is almost isotropic, and that the scattering interactions dominate over absorption [13]. Defining a computational domain $\Omega$ with boundary surface $\Gamma$ and internal domain $\Omega '$ (i.e. $\Omega = {\Omega '} \cup \Gamma$ and $\Omega ' \cap \Gamma = \emptyset$), the DE for a continuous wave (CW) imaging system is given as

(1)$$- \nabla \cdot\left( {\kappa \left( x \right)\nabla \Phi \left( x \right)} \right) + {\mu _a}\left( {x} \right)\Phi \left( x \right) = q_0 \left(x\right)\;\;\; \textrm{for}\;x \in \Omega.$$

$\Phi \left ( x \right )$ is the photon fluence rate as a function of position $x$. The diffusion coefficient $\kappa \left ( {x} \right )=1/(3(\mu _a\left ( {x} \right )+\mu _s^{\prime }\left ( {x} \right )))$, where $\mu _a$ and $\mu _s$ are the spatially varying absorption and scattering coefficients and $\mu _s^{\prime }\left ( {x} \right )=(1-g)\mu _s\left ( {x} \right )$ where $g$ is the anisotropy factor [14]. $q_0\left ( x \right )$ is the isotropic component of the source. $\nabla$ is the gradient operator and $\nabla \cdot (\cdot )$ denotes the differential divergence of a vector function (i.e. $\kappa \nabla \Phi$). This is usually solved under the Robin boundary condition (RBC) in which light that escapes the medium does not come back. The RBC is written as

(2)$$2A\widehat {\textrm{n}}\cdot\left( {\kappa \left( x \right)\nabla \Phi \left( x \right)} \right) + \Phi \left( x \right) = 0 \;\; \textrm{for}\;x \in \Gamma,$$

where $\widehat {\textrm{n}}$ denotes the outward unit normal on the boundary. $A$ is related to the relative refractive index mismatch between the medium and air and is derived from Fresnel’s law [14].

Mathematically, Eq. (1) is an elliptic partial differential equations, the differential operators (i.e. gradient or divergence) in which are defined using the classical vector calculus. A general approach to analytically solve the DE (with its RBC) is to apply the Green function, but analytical solutions are only known for homogeneous objects [15–17]. For more complex DOT geometries, the finite element method (FEM) [16] is commonly used to discretize the DE and its RBC. In this discretization, the computational domain $\Omega$ is divided into a series of elements (triangles in 2D, tetrahedra in 3D). The discretized geometry using FEM is normally termed as a FE mesh. However, FEM implementations can be difficult and time-consuming, especially when higher-order polynomial basis (shape) functions are used for non-linear interpolation between vertices of high-order elements [18]. In previous work [19], we introduced a graph representation to discretize a total-variation regularization term for the inverse problem in DOT. In this discretization, the object geometry is represented by an unstructured graph, defined by vertices, edges and weights. The graph was constructed by exploring neighborhood relationships between vertices.

In order to fully leverage the power of graph-based discretization, one must use the nonlocal vector calculus. In the classical local vector calculus, the differential operators are numerically evaluated using purely local information. In the nonlocal calculus, the operators include more pixel information in the domain. For example, in image processing, some well-known PDEs and variational techniques such as nonlocal image denoising [20,21] and inpainting [22,23] have explored the advantages of nonlocal vector calculus [21,24]. When applied to these problems, local operators include information from only neighbouring pixels whilst nonlocal methods include information from a wider area and are naturally formulated in a graph-based representation instead of in terms of the classical local difference operators.

In image processing, nonlocal methods are shown to have several advantages over local methods, including preservation of important image features such as texture and ability to handle unstructured geometries. It has also been observed that many PDE-based physical processes, minimizations and computational methods, such as CT image processing and reconstruction [25], can be generalized to be nonlocal. Therefore we expect that such a framework may be useful for the physical modelling in DOT.

As such, we propose a nonlocal diffusion equation (NDE) as a new forward model for DOT. The concept of differential operators under the nonlocal vector calculus [21–24] is used to formulate a new forward model that can accurately simulate light propagation in turbid media. The discretization for the NDE is performed using a graph-based numerical method (GNM). As a result, the proposed method naturally applies without modification to complex, unstructured DOT geometries in both two and three dimensions. The accuracy of the proposed model is compared against the conventional diffusion equation implemented by FEM and to the existing analytical solution on a homogeneous slab. We also compare the image reconstruction accuracy of different forward models on a 2D circular model and a 3D human head model. It should be noted that the diffusion equation is also used to model light propagation in imaging techniques such as diffuse correlation spectroscopy and near infrared spectroscopy, and our results can be applied to any technique that uses a diffusion-based model of light propagation.

2. Methodology

Our approach is based on reformulating the diffusion equation (Eq. (1)) in terms of nonlocal differential operators. We denote $\nabla _w(\cdot )$, $\mathrm {div}_w(\cdot )$ and ${\cal N}_w(\cdot )$ as the nonlocal gradient, the nonlocal divergence and the nonlocal normal derivative, respectively. Their definitions are given in Eq. (6), Eq. (7) and Eq. (8). We simply replace the differential operators in Eq. (1) with their nonlocal counterparts and solve the new NDE under the framework of nonlocal vector calculus:

(3)$$- \mathrm{div}_w \left( {\kappa \left( x \right)\nabla_w \Phi \left( x \right)} \right) + {\mu _a}\left( x \right) \Phi \left( x \right) = q_0 \left(x\right)\;\; \textrm{for}\;x \in \Omega '.$$

Similarly, we reformulate the RBC with the nonlocal normal derivative and the nonlocal gradient to give a nonlocal boundary condition (NBC):

(4)$$2A{\cal N}_w\left( {\kappa \left( x \right)\nabla_w \Phi \left( x \right)} \right) + \Phi \left( x \right) = 0 \;\; \textrm{for}\;x \in \Gamma.$$

We now formulate a graph-based numerical method to discretize the NDE with its NBC. Following established methods [21,24], we first discretize the computational domain $\Omega$ using a weighted graph $G = \left ( {V,E,w} \right )$, where $V = \left \{ {{V_k}} \right \}_{k = 1}^N$ denotes a finite set of $N$ vertices, and $E \in V \times V$ represents a finite set of weighted edges. Here $V=V_{\Omega '} \cup V_{\Gamma }$ with $V_{\Omega '}$ representing vertices in $\Omega '$ and $V_{\Gamma }$ vertices on boundary $\Gamma$. In this study, we assume that $G$ is an undirected simple graph (no multiple edges). Let $(i,j) \in E$ be an edge of $E$ that connects the vertices $i$ and $j$ in $V$. The weight $w_{ij}$ denotes the similarity between two vertices $i$ and $j$. The computation of this quantity is discussed later in this section. The nonlocal differential operators required by Eq. (3) and Eq. (4) on the graph $G$ are then defined as follows.

Definition (Nonlocal gradient). For a function $\Phi _i: V \to \mathbb {R}$ and a nonnegative and symmetric weight function $w_{ij}$: $V \times V \to \mathbb {R}$, the nonlocal partial derivative can be written as

(5)$${\partial_j}{\Phi _{i}} \triangleq \left( {{\Phi _j} - {\Phi _i}} \right)\sqrt {{w_{ij}}} :V \times V \to \mathbb{R}.$$

Therefore the nonlocal gradient ${\nabla _w}{\Phi _{i}}$ is defined as the vector of all partial derivatives:

(6)$${\nabla_w}{\Phi _{i,j}} \triangleq \left( {{\Phi _j} - {\Phi _i}} \right)\sqrt {{w_{ij}}} :V \times V \to \mathbb{R}.$$

Definition (Nonlocal divergence). Given a vector function $\boldsymbol {\nu }_i$: $V_{\Omega '} \to \mathbb {R}$ and a weight function $w_{ij}$: $V \times V \to \mathbb {R}$, the nonlocal divergence operator $\mathrm {div}_w$ acting on $\boldsymbol {\nu }_i$ is

(7)$$\mathrm{div}_w\, \boldsymbol{\nu} _i \triangleq \sum\limits_{j=1}^N {\left( {{\nu _{ij}} - {\nu _{ji}}} \right)\sqrt {{w_{ij}}} }: V_{\Omega'} \to \mathbb{R} ,$$

where ${\nu _{ij}}$ is the $j$’th element of $\boldsymbol {\nu } _{i}$.

Definition (Nonlocal normal derivative). Given a function $\boldsymbol {\nu }_i$: $V_{\Gamma } \to \mathbb {R}$ and a weight function $w_{ij}$: $V \times V \to \mathbb {R}$, the nonlocal normal operator acting on $\boldsymbol {\nu }_i$ is

(8)$${\cal N}_w \boldsymbol{\nu}_i \triangleq - \sum\limits_{j=1}^N{\left( {\nu _{ij} - {\nu _{ji}}} \right)\sqrt{w_{ij}}}: V_{\Gamma} \to \mathbb{R} .$$

Definition (Nonlocal Laplacian). Let $\Phi _i: V \to \mathbb {R}$ and $w_{ij}$: $V \times V \to \mathbb {R}$. The linear nonlocal Laplace operator acting on $\Phi _i$ is defined based on Eq. (6) and Eq. (7):

(9)$${\Delta_w}{\Phi _{i}} \triangleq \frac{1}{2}\mathrm{div}_w\,\left( {{\nabla}_w{\Phi _{i}}} \right) = \sum\limits_{j =1}^N {\left( {{\Phi _{j}} - {\Phi _{i}}} \right){w_{ij}}}: V \to \mathbb{R}.$$

The nonlocal normal derivative ${\cal N}_w$ in Eq. (8) is a nonlocal analogue of the normal derivative operator at the boundary encountered in the classical differential vector calculus (i.e. $\widehat {\textrm{n}}$ in Eq. (2)). Note that $\mathrm {div}_w$ in Eq. (7) and ${\cal N}_w$ in Eq. (8) have similar definitions but differ in their signs and the regions over which $\mathrm {div}_w\, \boldsymbol {\nu }_i$ and ${\cal N}_w \boldsymbol {\nu }_i$ are calculated. Also note that the mapping $\boldsymbol {\nu }_i \mapsto {\cal N}_w \boldsymbol {\nu }_i$ is scalar-valued which is analogous to the local differential divergence of a vector function in Eq. (1). Finally, with the definitions of $\mathrm {div}_w$ and ${\cal N}_w$, the nonlocal divergence theorem is $\int _{\Omega '} {\mathrm {div}_w\, \boldsymbol {\nu }dx} = \int _\Gamma {{\cal N}_w \boldsymbol {\nu }dx}$, which essentially relates the flow (i.e. flux) of a nonlocal vector field through a boundary/surface to the behaviour of the nonlocal vector field inside the boundary/surface.

It should be noticed from the nonlocal differential operator definitions (Eq. (6), Eq. (7), Eq. (8) and Eq. (9)) that, in a full non-local scheme, each vertex has connections with all the vertices in $V$ over $\Omega$ such that the constructed graph is fully connected. This can make the computational load extremely heavy and so approaches based on spectral graph theory [26,27] or nearest neighbors [28], are typically employed to partition the vertices in the computational domain into groups according to their similarities. For example, Bertozzi [27] used spectral approaches along with the Nyström extension method to efficiently calculate the eigendecomposition of a dense graph Laplacian. The second eigenvector of the graph Laplacian was used to initialize the partitioning so that the weights between vertices in different groups are small and the weights between vertices within the same group are large. In this paper, we build the graph by using the positions of the nodes and the connectivity between nodes in the finite element mesh as the vertices and edges in the graph to sparsify the graph for computational efficiency. We have learned from previous work [19] that the graph-based nonlocal inverse model with this sparse method can achieve accurate and stable reconstruction, regardless of the mesh resolution. Therefore for each vertex $i$, we consider only those vertices that are directly connected to the vertex $i$ for ${{\cal N}_i}$ (i.e. those vertices that share the same edge with $i$). With this structure and the nonlocal discrete differential operators, we can derive the following discretized versions of Eq. (3) and Eq. (4):

(10)$$\begin{aligned} \sum\limits_{j \in {{ {\cal{N}}}_i}} {\left( {{\kappa _i} + {\kappa _j}} \right)\left( {{\Phi _i} - {\Phi _j}} \right){w_{ij}}} + {{\mu_a}_i}{\Phi _i} = {q_{0i}}\;\;\;{\textrm{for}}\;i \in \Omega ' \\ 2A\sum\limits_{j \in {{{\cal{N}}}_i}} {\left( {{\kappa _i} + {\kappa _j}} \right)\left( {{\Phi _i} - {\Phi _j}} \right){w_{ij}}} + {\Phi _i} = 0\;\;\;\;\;{\textrm{for}}\;i \in \Gamma\\ \end{aligned}$$

The nonnegative and symmetric weight function $w_{ij}$ between two connected vertices $i$ and $j$ has many possible choices. In this work, we first obtain the similarity $w_{ij}$ by simply using the inverse of the Euclidean distance $d_{ij}$ between two nodes. Then we normalize the similarity using $w_{ij}/{\sum _{j \in {{{\cal {N}}}_i}}w_{ij}}$ to convert the similarities into probabilities and ensure that the probabilities sum to one.

We note that due to the nature of the graph representation, the implementation of Eq. (10) is identical for a 2D or 3D geometry. It should also be noted that increasing the number of vertices and edges will decrease the sparsity of the graph and increase the computational burden with no change in the implementation. Under these assumptions, Eq. (10) can be rewritten in matrix form as

(11)$${\cal{M}}\Phi = Q.$$

${\cal {M}}$ is a $N \times N$ sparse matrix and a symmetric, diagonally dominant and positive definite real-value matrix, whose entries are

{{\cal{M}}_{i,j}} = \left\{ \begin{array}{ll} \sum\limits_{j \in {{\cal {N}}_i}}^{} {\left( {{\kappa _i} + {\kappa _j}} \right){w_{ij}}} + {{\mu_a}_i} & {\textrm{if}}\; i=j \in \Omega '\\ \sum\limits_{j \in {{\cal {N}}_i}}^{} {\left( {{\kappa _i} + {\kappa _j}} \right){w_{ij}}} + \frac{1}{2A} & {\textrm{if}}\; i=j \in \Gamma \\ -{\left( {{\kappa _i} + {\kappa _j}} \right){w_{ij}}} & {\textrm{if}}\; i \neq j \; {\textrm{and}}\; {j \in {{\cal {N}}_i}} \\ \;\;\;\;\;\;\;\;\;\;\;{0} & {\textrm{otherwise}} \end{array} \right..

$Q$ is a $N \times N_s$ sparse matrix where $N_s$ is the number of sources and each column represents one distributed Gaussian source. The linear system (Eq. (11)) can be solved exactly by using a direct solver with Cholesky decomposition.

3. Experimental results

In this section, numerical experiments are conducted to quantitatively evaluate the performance of the proposed NDE method. The NDE method with the GNM implementation will be compared against the original DE with the FEM implementation. We evaluate the light propagation performance of the proposed method in a 3D homogeneous rectangular slab and then a heterogeneous two-layer rectangular slab where the analytical solutions are known, followed by two dimensional (2D) and three dimensional (3D) image reconstruction examples. All the experiments are performed using Matlab 2018b on a Windows 7 platform with an Intel Xeon CPU i7-6700 (3.40 GHz) and 64 GB memory.

3.1 Forward modelling on a 3D rectangular-slab model

To quantitatively compare our GNM method with classical FEM approaches, we model a rectangular-slab of size: 200$\times$100$\times$100 mm$^3$, as shown in Fig. 1. The mesh is composed of 442381 nodes corresponding to 2620541 tetrahedral elements, with the average nodal distance of 1.5 mm. For the forward model based on FEM, such a discrete structure can be directly employed for the finite element method. However, the forward model based on GNM requires only the vertices and edges of the mesh. It is well known that Monte Carlo method is the gold standard for modeling photon and electron transport in medium. However it is costly in computational time, because a large number of photons need to be simulated so as to acquire meaningful statistics. FEM based models have been shown to be capable of producing quantitatively correct boundary measurements when the mesh resolution is high [29]. Therefore in this work, we evaluate the light propagation performance of both methods (FEM and GNM) on a slab mesh in which the average nodal distance is as small as 1.5 mm. We then compare both methods with the well-known analytical solution. We conduct simulations using a CW source for which we can analytically calculate the photon flux measurement on the boundary (BF) as well as the fluence rate (FR) at each vertex.

Fig. 1. Rectangular-slab mesh with one source (red dot) and six detectors (green dots). The distance between the source and the six detectors varies from 15 mm to 40 mm, in 5 mm increments.

Download Full Size | PDF

3.1.1 3D homogeneous rectangular-slab model

For a 3D homogeneous rectangular-slab model, the optical parameters $\mu _a$ and $\mu _s^{\prime }$ in the slab were set to 0.01 mm$^{-1}$ and 1 mm$^{-1}$, respectively. The analytical solution of the BF has the form [15]:

(12)$$I\left( \rho \right) = \frac{1}{{4\pi }}\left[ {\frac{1}{{{\mu _a} + {\mu _{s}^\prime}}}}\left( {{\mu _{\mathrm{eff}}} + \frac{1}{{{r_1}}}} \right)\frac{{{e^{ - {\mu _{\mathrm{eff}}}{r_1}}}}}{{r_1}^2} + \frac{{3 + 4A}}{{3\left( {{\mu _a} + {\mu _{s}^\prime}} \right)}}\left( {{\mu _{\mathrm{eff}}} + \frac{1}{{{r_2}}}} \right)\frac{{{e^{ - {\mu _{\mathrm{eff}}}{r_2}}}}}{{r_2}^2} \right] ,$$

where $\rho$ represents the distance from the source, $A$ is the internal reflection parameter for the air-tissue interface, ${\mu _{\mathrm {eff}}}$ is the effective attenuation coefficient which is $\sqrt {3{\mu _a}\left ( {{\mu _a} + {\mu _{s}^\prime }} \right )}$, ${r_1} = \sqrt {1/( {\mu _a} + {\mu _{s}^\prime } )^2 + {\rho ^2}}$ and ${r_2} = \sqrt {{( 3 + 4A )^2}/{( {3( {\mu _a} + {\mu _{s}^\prime } )} )^2} + {\rho ^2}}$.

In Fig. 2(a), we plot the normalized photon flux at the boundary (NBF). We normalize the BF to remove any constant offset resulting from the use of different propagation models. It can be seen that the NBF from both forward models match the analytical solution. In order to observe the difference clearly, in Fig. 2(b), we plot the percentage of error between the analytical solution and the other two methods with regards to NBF. The percentage of error is calculated by, for each source-detector channel, dividing the absolute difference between each forward model and the analytical solution by the analytical solution. We average the percentage errors along the six source-detector pairs. The forward models based on FEM and GNM are both shown to reproduce the analytical solution to within 7% on average.

Fig. 2. The flux measurements on the boundary versus the source-detector distance. (a): NBF; (b): Percentage of error based on NBF.

Download Full Size | PDF

We then compare the FR calculated at the vertices inside of the medium. The analytical solution of the FR is [17]:

(13)$$\Phi \left( {r,z} \right) = \frac{P\mu_{\mathrm{eff}} ^2}{4\pi\mu_a}\left[ \left( {\frac{\exp \left\{ { - \mu_{\mathrm{eff}} {{\left[ {{{\left( z - z_0 \right)}^2} + r^2} \right]}^{1/2}}} \right\}}{{ - {\mu _{\mathrm{eff}}} {{\left[ \left( z - z_0 \right)^2 + r^2 \right]}^{1/2}}}}} \right) - \left( {\frac{{\exp \left\{ { - \mu_{\mathrm{eff}} {{\left[ \left( z + z_0 \right)^2 + {r^2} \right]}^{1/2}}} \right\}}}{{ - \mu_{\mathrm{eff}} {{\left[ {{{\left( z + z_0 \right)}^2} + r^2} \right]}^{1/2}}}}} \right) \right],$$

where $P$ is the source power. $z_0$ is the depth of the source which is $1/\mu '_s$. $z$ represents the depth under the surface which is $z=50$mm in our case. $r$ is the distance between a given vertex and the source on the X-Y plane. Note that $\sqrt {{( z - {z_0} )^2} + {r^2}}$ represents the distance between a given vertex and the source.

In Fig. 3, we compare the FR calculated using Eq. (13) and the FEM and GNM methods. We choose the vertical plane across the source-detector positions as the region of interest (ROI). For each method, in order to remove any constant offset resulting from the use of different propagation models, we rescaled FR onto the range $[0,1]$ by dividing the FR with the highest FR value in the ROI and name the rescaled FR as NFR. This is necessary because in FEM, point sources are distributed across the nodes belonging to the element in which the source is placed, whereas in GNM, the source is fully attached to the nearest node. The two methods can therefore have different initialization states for the same source. In Fig. 3(a)–(c), we plot the NFR at each vertex in the ROI calculated using the analytical method, and the FEM and GNM models, respectively. We also plot its logarithm in (d)–(f), corresponding to the NFR in (a)–(c) respectively. It can be observed that the light propagation in the medium modelled by the proposed forward model is comparable to the one modelled by the forward model based on FEM. In order to see the difference clearly, in Fig. 4, we plot the descending tendency of the NFR calculated by different propagation methods. Specifically, we plot the logarithm of NFR along the z axis starting from the source position. As can be seen, for all methods the fluence rate gradually drops as the light penetrates deeper. The descending tendency of the curves derived from the both forward methods are almost parallel to the one from the analytical solution. Therefore we can see that all the three models can generate the same NFR distribution.

Fig. 3. (a)–(c): NFR at each vertex in the ROI calculated using the analytical solution, forward models based on FEM and the one based on GNM, respectively; (d)–(f): logarithm of the NFRs, corresponding to (a)–(c).

Download Full Size | PDF

Fig. 4. Descending tendency of the NFR from the source to the medium along the z axis.

Download Full Size | PDF

3.1.2 3D heterogeneous rectangular-slab model

We further evaluate the light propagation performance of FEM and GNM on a heterogeneous model which has the same geometry as the previous example, but contains two optically different layers. We choose the thickness of the first layer to be relevant to one potential application of the layered model: DOT for measurements of cerebral oxygenation, where the thickness of the tissues above the brain is around 10 mm. Figure 5 gives the 2D section of the two-layered model in which the thickness of the first layer is 7.5mm. The optical properties ($\mu _a$, $\mu _s^{\prime }$) in the two layers were set to (0.015 mm$^{-1}$, 1.5 mm$^{-1}$) and (0.01 mm$^{-1}$, 1 mm$^{-1}$), respectively. The positions of sources and detectors are the same as that in Section 3.1.1. The analytical solution of the BF for two-layered medium is given in [30] and we use it as the ground truth to evaluate the accuracy of our proposed forward model. In Fig. 6, we first plot the normalized photon flux (NBF) at the boundary versus the source-detector separation and then give the percentage of error between the analytical solution and the other two methods. The forward models based on FEM and GNM are both shown to reproduce the analytical solution to within 8.5% on average. Similar to what we observed in the previous section, GNM has less error than FEM when the source-detector distance is short (15mm) while the error becomes greater when the source-detector distance is as long as 40mm. There is no significant difference in the percentage errors when the source-detector distance is between 20mm to 35mm. Therefore similar conclusions can be achieved from the experiments in Section 3.1.1 and 3.1.2.

Fig. 5. Scheme of the heterogeneous rectangular-slab model.

Download Full Size | PDF

Fig. 6. The flux measurements on the boundary versus the source-detector distance. (a): NBF; (b): Percentage of error based on NBF.

Download Full Size | PDF

3.1.3 Computational time

After evaluating the accuracy of the fluence rates and boundary measurements modelled by different forward models, in Fig. 7, we compare the computational efficiency of FEM and GNM forward models on the 3D homogeneous rectangular-slab model. We run each model on six meshes with different average nodal distance of 1.5, 2, 2.5, 3, 3.5 and 4 mm respectively. The mesh spatial resolution becomes lower when the nodal distance is larger. We run each forward modelling process ten times and record the mean and standard deviation of the CPU time consumed for computing one source-detector channel. For a fair comparison, we use a direct solver with Cholesky decomposition to solve the linear equation resulting from each forward model.

Fig. 7. CPU time (s) consumed at one source-detector channel using different forward models. ’A’ represents the FEM approach while ’B’ represents the GNM approach. Right figure is the zoomed-in plot of the area in the green dash line of the left figure.

Download Full Size | PDF

For all mesh resolutions, based on each source-detector channel, the CPU times required by the FEM model are larger than for GNM. When the mesh resolution is low (for example the case where the average nodal distance is 4 mm) the CPU time consumed by the FEM approach (0.11s) is 175% larger than the time required by the GNM approach (0.04s). When the mesh resolution is high (average nodal distance is 1.5 mm), the CPU time consumed by the FEM approach (14.6s) is only 14% longer than the GNM approach (12.7s). It should be noticed that the CPU time plotted in Fig. 7 is only for one source-detector channel. GNM is 14% faster than the FEM approach when average nodal distance is low (mesh resolution is high). For DOT which normally has thousands of source-detector channels in some clinical scenarios, this 14% improvement at each channel would make a significant difference in the total computational time. This result demonstrates the computational efficiency of the proposed forward model.

3.2 Image reconstruction using different forward models

We now consider the recovery of the optical properties at each vertex within the medium using both forward models. The image reconstruction process is implemented by iteratively refining the optical properties of the forward model until the forward model prediction matches the boundary measurements [14]. It can be implemented by solving the following minimization problem:

(14)$${\mu_a ^*} = \arg {\min\limits _{\mu_a} }\left\{\| {\Phi _{}^\textrm{M} - {\cal F}\left( \mu_a \right)\|_2^2 + \lambda {\cal R}\left( \mu_a \right)} \right\},$$

where $\Phi ^{\textrm {M}}$ represents the boundary measurements acquired from the optical detectors, ${\cal F}$ is the non-linear operator induced from the forward model, ${\cal R}$ is a general regularization term, and $\lambda$ is a weight that determines the extent to which regularization will be imposed on the solution $\mu _a^*$. In this paper, we adopt the popular quadratic Tikhonov-type regularization (${\cal R}(\mu _a)=\|\mu _a - \mu _{a,0} \|_2^2$) for all methods for fair comparison [14]. Four quantitative evaluation metrics are considered to evaluate the reconstruction results: the average contrast (AC) [31], peak signal-to-noise ratio (PSNR) [31], structural similarity index (SSIM) [32] and root mean square error (RMSE) [32]. If the reconstructed image is identical to the ground truth image, AC is equal to 1. For PSNR and SSIM, the recovered image has higher quality if higher PSNR or SSIM values are obtained. Lower RMSE represents better reconstruction results. Randomly generated Gaussian noise is added to the amplitude of the measurement vector to simulate real noise in a CW system. In order to reduce the randomness resulting from the randomly distributed Gaussian noise, we run each experiment ten times and record the average (mean) and standard deviation (SD) of the four evaluation metrics.

3.2.1 Image reconstruction on a homogeneous circle model

We consider a 2D homogeneous circular geometry containing one target activation region (Fig. 8(b)). The model has a radius of 43mm and is composed of 1785 nodes and 3418 linear triangle elements. Sixteen source-detector fibres are placed equidistant around the external boundary for data acquisition (Fig. 8(a)). When one fibre as a source is turned on, the rest are used as detectors, leading to 240 total boundary measurements. All sources were positioned one scattering distance within the outer boundary because the source is assumed to be spherically isotropic. The background absorption coefficient is set to 0.01 mm$^{-1}$. One 10mm radius target region is centred at (20mm, 0mm) with 0.03 mm$^{-1}$ absorption coefficient. The reduced scattering coefficient is set to be homogeneous throughout the whole computational domain with the value of 1 mm$^{-1}$. 1% normally distributed Gaussian noise was added to the amplitude of the measurement vector.

Fig. 8. (a): A typical circle mesh with sixteen co-located sources and detectors; (b): True distribution of $\mu _a$; (c): Images reconstruction of $\mu _a$ using the forward model based on FEM and GNM (from left to right column) on $0\%$ (top part) and $1\%$ (bottom part) noisy data.

Download Full Size | PDF

Figure 8(c) shows the reconstruction results using the forward model based on FEM (Eq. (1)) and GNM (Eq. (3)) on $0\%$ and $1\%$ noisy data respectively. By visual inspection, it is evident that for the same level of Gaussian noise, the image recovered using the GNM approach is similar to the one recovered using the FEM approach. Figure 9 gives the 1D cross section of the results recovered in Fig. 8 along the horizontal line across the centre of the target (20mm, 0mm). It can be seen that the curves resulting from different forward models have similar edge smoothing resulting from the Tikhonov regularization and slightly different peak values. This is consistent with our visual observation from the reconstructed images in Fig. 8.

Fig. 9. 1D cross sections of images recovered in Fig. 8 along the horizontal line across the centre of the target. Left to right column: $0\%$ and $1\%$ added Gaussian noise.

Download Full Size | PDF

In Table 1, the values of the metrics AC, PSNR, SSIM and RMSE are shown to qualitatively evaluate the results in Fig. 8. It can be observed that when the data is clean, GNM gives AC closer to 1, slightly higher PSNR and SSIM, and lower RMSE than FEM. For the noisy data, GNM achieves similar AC, PSNR, SSIM and RMSE values with the FEM approach. This experiment quantitatively validates the forward modelling capacity of our proposed model and the consistency between these two forward models.

Table 1. Evaluation metrics for the recovered results using FEM and GNM on data with 0% and 1% added noise.

View Table | View all tables in this article

3.2.2 Image reconstruction on a heterogeneous head model

We now evaluate both forward models on a physically realistic three dimensional heterogeneous head model. This head model is composed of three tissue layers which are scalp, skull and brain. The reconstruction mesh consists of 50721 nodes associated with 287547 tetrahedral elements, with the average element size 9.3mm$^{3}$. Each node is assigned to one of the three layers. Absorption coefficients assigned to each layer refer to an in vivo study [33] at 750nm.

A large rectangular imaging array with 36 sources and 37 detectors was placed over the back-head area (Fig. 10(a)), allowing use of multiple sets of overlapping measurements which can improve both the spatial resolution and quantitative accuracy [34]. The source-detector (SD) separation distances ranges from 1.3 to 4.8cm, leading to 590 overlapping, multi-distance measurements. One anomaly with 15mm radius is simulated in the brain (Fig. 10(b)). In order to simulate traumatic brain injury (TBI) cases where the cerebral tissue oxygen saturation ($\textrm{StO}_{2}$) is normally between 50% and 75% [35] (compared to 80% in healthy tissue), the absorption coefficient in the anomaly is calculated using Beer’s law [14] with 55% ${\textrm {StO}}_{2}$. In line with the current in vivo performance of the imaging system, 0.12$\%$, 0.15$\%$, 0.41$\%$ and 1.42$\%$ Gaussian random noise was added to first (13mm), second (30mm), third (40mm) and fourth (48mm) nearest neighbor measurements to provide realistic data [36]. Reconstructed absorption coefficients of the simulated anomaly using different models are displayed in the third to fourth column of Fig. 10. Corresponding 2D cross section is given in the second row. The visualization suggests that GNM can achieve better reconstruction performance with optical property values closer to the ground truth. This may be because the simulated anomaly is close to the outer surface and GNM gives lower errors when the source detector distance is low. We can also see that the results by both methods are smoothed and the volume sizes of the recovered anomaly are smaller than the ground truth. Evaluation metrics are given in Table 2. Even though there is slight visual difference between the two methods , no obvious difference between these two reconstruction models can be observed from the four evaluation metrics. These findings further quantitatively validate the consistency between these two forward models.

Fig. 10. (a): Three-dimensional head mesh and distribution of the rectangular imaging array with 36 sources (red dots) and 37 detectors (green dots); (b): Ground truth; (c) Reconstruction with the forward model based on FEM and GNM, respectively.

Download Full Size | PDF

Table 2. Evaluation metrics for $\mu _a$ on the recovered results shown in Figure 10.

View Table | View all tables in this article

4. Conclusion

We have proposed a new formulation of the forward model for DOT that is based on the concepts of differential operators under a nonlocal vector calculus. The discretization of the new forward model is performed using an efficient graph-based numerical method. Our proposed model is shown to be able to accurately model the light propagation in the medium and is quantitatively comparable with both analytical and FEM forward models. Compared with the conventional forward model based on FEM, our proposed model has the following two advantages: 1) according to the experiments in Section 3.1, our proposed model is shown to be more computationally efficient with an average speed improvement of 30% compare to the FEM forward model due to the simple graph-based discretization; 2) it allows identical implementation for geometries in different dimensions thanks to the nature of the graph representation. In addition, the proposed graph-based discretization method can also be applied to other imaging techniques which are modelled using a diffusion equation, where has potential to improve the computational efficiency and simplicity.

Funding

H2020 Marie Skłodowska-Curie Actions (675332).

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 Marie Sklodowska-Curie Innovative Training Networks (ITN-ETN) programme, under grant agreement no 675332, BitMap.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References

1. J. J. Duderstadt and W. R. Martin, Transport Theory (John Wiley & Sons Chichester, 1979).

2. A. D. Klose, “The forward and inverse problem in tissue optics based on the radiative transfer equation: a brief review,” J. Quant. Spectrosc. Radiat. Transfer 111(11), 1852–1853 (2010). [CrossRef]

3. G. Eason, A. R. Veitch, R. M. Nisbet, and F. W. Turnbull, “The theory of the back-scattering of light by blood,” J. Phys. D: Appl. Phys. 11(10), 1463–1479 (1978). [CrossRef]

4. S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37(7), 1531–1560 (1992). [CrossRef]

5. W. H. Reed, “New difference schemes for the neutron transport equation,” Nucl. Sci. Eng. 46(2), 309–314 (1971). [CrossRef]

6. P. González-Rodríguez and A. D. Kim, “Comparison of light scattering models for diffuse optical tomography,” Opt. Express 17(11), 8756–8774 (2009). [CrossRef]

7. B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. 10(6), 824–830 (1983). [CrossRef]

8. J. Heiskala, M. Pollari, M. Metsäranta, P. E. Grant, and I. Nissilä, “Probabilistic atlas can improve reconstruction from optical imaging of the neonatal brain,” Opt. Express 17(17), 14977–14992 (2009). [CrossRef]

9. L. Kocsis, P. Herman, and A. Eke, “The modified Beer–Lambert law revisited,” Phys. Med. Biol. 51(5), N91–N98 (2006). [CrossRef]

10. W. B. Baker, A. B. Parthasarathy, D. R. Busch, R. C. Mesquita, H. Joel, and A. G. Yodh, “Modified Beer-Lambert law for blood flow,” Biomed. Opt. Express 5(11), 4053–4075 (2014). [CrossRef]

11. D. A. Boas, T. Gaudette, G. Strangman, X. F. Cheng, J. J. A.. Marota, and J. B. Mandeville, “The accuracy of near infrared spectroscopy and imaging during focal changes in cerebral hemodynamics,” NeuroImage 13(1), 76–90 (2001). [CrossRef]

12. M. Bhatt, K. R. Ayyalasomayajula, and P. K. Yalavarthy, “Generalized Beer–Lambert model for near-infrared light propagation in thick biological tissues,” J. Biomed. Opt. 21(7), 076012 (2016). [CrossRef]

13. D. A. Boas, “Diffuse photon probes of structural and dynamical properties of turbid media: theory and biomedical applications,” Doctoral dissertation, Graduate School of Arts and Sciences, University of Pennsylvania (1996).

14. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Meth. Engng. 25(6), 711–732 (2009). [CrossRef]

15. T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19(4), 879–888 (1992). [CrossRef]

16. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]

17. V. Allen and A. L. McKenzie, “The modified diffusion dipole model,” Phys. Med. Biol. 36(12), 1621–1638 (1991). [CrossRef]

18. D. R. Lynch, Numerical Partial Differential Equations for Environmental Scientists and Engineers: A First Practical Course (Springer Science & Business Media, 2004).

19. W. Q. Lu, J. M. Duan, D. Orive-Miguel, L. Herve, and I. B. Styles, “Graph-and finite element-based total variation models for the inverse problem in diffuse optical tomography,” Biomed. Opt. Express 10(6), 2684–2707 (2019). [CrossRef]

20. A. Buades, B. Coll, and J. M. Morel, “A non-local algorithm for image denoising,” 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05)2, 60–65 (2005).

21. G. Gilboa and S. Osher, “Nonlocal operators with applications to image processing,” Multiscale Model. Simul. 7(3), 1005–1028 (2009). [CrossRef]

22. J. M. Duan, Z. K. Pan, W. Q. Liu, and X. C. Tai, “Color texture image inpainting using the non local CTV model,” JSIP 04(03), 43–51 (2013). [CrossRef]

23. J. M. Duan, Z. K. Pan, B. C. Zhang, W. Q. Liu, and X. C. Tai, “Fast algorithm for color texture image inpainting using the non-local CTV model,” J Glob Optim 62(4), 853–876 (2015). [CrossRef]

24. M. Gunzburger and R. B. Lehoucq, “A nonlocal vector calculus with application to nonlocal boundary value problems,” Multiscale Model. Simul. 8(5), 1581–1598 (2010). [CrossRef]

25. H. Zhang, D. Zeng, H. Zhang, J. Wang, Z. R. Liang, and J. H. Ma, “Applications of nonlocal means algorithm in low-dose X-ray CT image processing and reconstruction: A review,” Med. Phys. 44(3), 1168–1185 (2017). [CrossRef]

26. A. L. Bertozzi and A. Flenner, “Diffuse interface models on graphs for classification of high dimensional data,” Multiscale Model. Simul. 10(3), 1090–1118 (2012). [CrossRef]

27. E. Merkurjev, T. Kostic, and A. L. Bertozzi, “An MBO scheme on graphs for classification and image processing,” SIAM J. Imaging Sci. 6(4), 1903–1930 (2013). [CrossRef]

28. X. Bresson, X. C. Tai, T. F. Chan, and A. Szlam, “Multi-class transductive learning based on l1 relaxations of Cheeger cut and Mumford-Shah-Potts model,” J Math Imaging Vis 49(1), 191–201 (2014). [CrossRef]

29. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995). [CrossRef]

30. A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van Den Bergh, “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37(4), 779–791 (1998). [CrossRef]

31. W. Q. Lu, D. Lighter, and I. B. Styles, “L 1-norm based nonlinear reconstruction improves quantitative accuracy of spectral diffuse optical tomography,” Biomed. Opt. Express 9(4), 1423–1444 (2018). [CrossRef]

32. W. Q. Lu, J. M. Duan, Z. W. Qiu, Z. K. Pan, Q. W. Liu, and L. Bai, “Implementation of high-order variational models made easy for image processing,” Math. Meth. Appl. Sci. 39(14), 4208–4233 (2016). [CrossRef]

33. A. T. Eggebrecht, B. R. White, S. L. Ferradal, C. X. Chen, Y. X. Zhan, A. Z. Snyder, H. Dehghani, and J. P. Culver, “A quantitative spatial comparison of high-density diffuse optical tomography and fmri cortical mapping,” NeuroImage 61(4), 1120–1128 (2012). [CrossRef]

34. D. A. Boas, A. M. Dale, and M. A. Franceschini, “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” NeuroImage 23, S275–S288 (2004). [CrossRef]

35. C. Ichai, H. Quintard, and J. C. Orban, Metabolic Disorders and Critically Ill Patients: From Pathophysiology to Treatment (Springer, 2017).

36. H. Dehghani, B. R. White, B. W. Zeff, A. Tizzard, and J. P. Culver, “Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography,” Appl. Opt. 48(10), D137–D143 (2009). [CrossRef]

	0% noise		1% noise (Mean $\pm$ SD)
	FEM	GNM	FEM	GNM
AC	1.1	1.0	1.1 $\pm$ 0.1	1.1 $\pm$ 0.1
PSNR	54.5	55.8	54.3 $\pm$ 0.5	54.3 $\pm$ 0.4
SSIM	99.6e-2	99.7e-2	99.6e-2 $\pm$ 7.6e-4	99.6e-2 $\pm$ 4.9e-4
RMSE	1.9e-3	1.6e-3	1.9e-3 $\pm$ 1.2e-4	1.9e-3 $\pm$ 8.4e-5

	Mean $\pm$ SD
	FEM	GNM
AC	0.8 $\pm$ 0.0	0.9 $\pm$ 0.1
PSNR	78.9 $\pm$ 0.0	79.0 $\pm$ 0.0
SSIM	99.9e-2 $\pm$ 9.1e-7	99.9e-2 $\pm$ 3.8e-7
RMSE	1.1e-4 $\pm$ 5.1e-7	1.1e-4 $\pm$ 1.8e-6

	0% noise		1% noise (Mean $\pm$ SD)
	FEM	GNM	FEM	GNM
AC	1.1	1.0	1.1 $\pm$ 0.1	1.1 $\pm$ 0.1
PSNR	54.5	55.8	54.3 $\pm$ 0.5	54.3 $\pm$ 0.4
SSIM	99.6e-2	99.7e-2	99.6e-2 $\pm$ 7.6e-4	99.6e-2 $\pm$ 4.9e-4
RMSE	1.9e-3	1.6e-3	1.9e-3 $\pm$ 1.2e-4	1.9e-3 $\pm$ 8.4e-5

	Mean $\pm$ SD
	FEM	GNM
AC	0.8 $\pm$ 0.0	0.9 $\pm$ 0.1
PSNR	78.9 $\pm$ 0.0	79.0 $\pm$ 0.0
SSIM	99.9e-2 $\pm$ 9.1e-7	99.9e-2 $\pm$ 3.8e-7
RMSE	1.1e-4 $\pm$ 5.1e-7	1.1e-4 $\pm$ 1.8e-6

New nonlocal forward model for diffuse optical tomography

Abstract

1. Introduction

2. Methodology

3. Experimental results

3.1 Forward modelling on a 3D rectangular-slab model

3.1.1 3D homogeneous rectangular-slab model

3.1.2 3D heterogeneous rectangular-slab model

3.1.3 Computational time

3.2 Image reconstruction using different forward models

3.2.1 Image reconstruction on a homogeneous circle model

3.2.2 Image reconstruction on a heterogeneous head model

4. Conclusion

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (10)

Tables (2)

Equations (15)

Biomedical Optics Express