A coupled finite element-boundary element method for modeling Diffusion equation in 3D multi-modality optical imaging

Subhadra Srinivasan; Hamid R. Ghadyani; Brian W. Pogue; Keith D. Paulsen

doi:10.1364/BOE.1.000398

1. Introduction

1.1 Introduction to 3D diffuse optical imaging

Diffuse optical imaging provides functional information related to the physiological status of tissue non-invasively. Absorption, fluorescence and Raman optical imaging have demonstrated ability to provide molecular fingerprints of tissues in healthy and diseased states [1–5]. These optical techniques require a model for image reconstruction from boundary measurements of tissues when used in tomographic applications in-vivo. Image reconstruction involves solving a model for light propagation (called the forward model) iteratively to fit the measured data and recover optical parameters. Traditionally, image reconstruction techniques have used the approximation that light propagation is two-dimensional. However, more recently interest in 3D image reconstruction has grown because it is more accurate than 2D models given that light propagation is inherently three-dimensional [6].

Three-dimensional models have been successfully applied to simple geometries such as cylinders, slabs and spheres where algorithms have been explored for better localization and quantification. For example, Yalavarthy et al [7] used a generalized least squares minimization incorporating data and parameter variances to accelerate 3D image reconstruction for under-determined problems. Using a level-set technique for image reconstruction, Schweiger et al [8] showed that detection and localization of small objects could be improved in 3D. Boverman et al [9] used a parametric approach to reconstruct shape and contrast of piece-wise constant regions in 3D with spherical harmonics for modeling sharp boundaries in tissue and demonstrated quantitative results in a domain with a single inclusion. Zacharopoulos et al [10] used a similar strategy and showed that they could accurately recover location and contrast of an anomaly in experiments on a domain with single inclusion. Srinivasan et al [11] used a dynamic criterion based on the least squares error norm of model-data mismatch to reduce the size of large data sets and speed up 3D image reconstruction. However, applications of 3D image reconstruction to arbitrary shaped geometries such as breast and brain have been more limited, especially as in the setting of multi-modality imaging.

1.2 Multi-modality optical imaging reconstruction techniques

Multi-modality imaging has gained interest as an approach for improving the contrast recovery of diffuse optical imaging and fluorescence [12–15]. Multi-modality imaging uses prior anatomical structure to guide the diffuse optical reconstruction spatially, making it less ill-posed and the images better resolved. In this reconstruction process, the optical imaging domain is typically defined by segmentation and volume meshing of conventional medical images (MRI, X-Ray or CT). Image reconstruction techniques involving multimodal data have generally evolved in two categories of implementation of the spatial data, including: (1) soft prior information and (2) hard prior information. Soft prior info refers to the application of anatomical constraints, which allow for optical property variations to occur within segmented regions. Studies have used algorithms based on total variation minimization [16], sparsity regularization [15], Laplacian and Helmholtz regularizations [14,17,18], data-specific spatially varying regularization [19], with all predominantly in the finite element method (FEM) framework. Hard prior info strictly enforce the tissue boundaries to represent homogeneous or piece-wise constant optical property regions. This has been implemented using FEM [20,21] and the boundary element method (BEM) [22]. Many of these studies have been on simulations with couple of case studies resulting from experimental or clinical data; extensive testing in experimental or clinical data is still to be demonstrated.

1.3 Need for efficient 3D technique for complex 3D domains

In our experience, one of the key challenges in adopting 3D multimodal optical imaging for large clinical studies is in image segmentation and meshing of arbitrary shapes. Figure 1 shows a schematic of a typical workflow before image reconstruction. The process involves segmentation of medical image data, surface rendering (which produces a surface mesh as output) and volumetric meshing. The last step of obtaining a volume mesh for 3D image reconstruction can be time-consuming and difficult to automate in a clinical workflow. Studies in brain and small-animal imaging have used a standard anatomical atlas to by-pass the problem of obtaining subject-specific meshes [23,24]. However, some tissues such as the breast and the prostate show considerably larger heterogeneity between subjects [25] where a subject specific mesh is imperative to the imaging process. Use of a BEM approach as an alternative to FEM for hard priors alleviated the meshing complexity by requiring only surface discretization as compared to volume meshing for modeling light diffusion in 3D [22,26]. BEM showed promise for multimodal image guided diffuse optical spectroscopy of piece-wise constant regions (hard priors) by simplifying the meshing process and implementing the assumption in the forward model itself [22].

Fig. 1 A schematic showing steps from medical image data to obtaining a volumetric mesh for computation with examples from breast data. These steps have to be routinely performed before image reconstruction can be done for 3D multi-modality optical imaging. Methods for image segmentation vary between applications; here thresholding and region-growing techniques were applied for breast tissue. Surface rendering is automatically generated by many open source softwares, but getting a reliable volume mesh can be time-consuming and more difficult to automate.

Download Full Size | PDF

However, using piece-wise constant optical property approximations has limitations: (1) it cannot model tissues which are known to have spatially varying optical property distributions such as large solid tumors [27] (2) results are affected when the prior information on tissue boundaries is imperfect [17,21], and (3) insufficient information exists when the boundary data is simply not available as in the case of false-negative findings in MRI. An efficient method to counter these limitations is needed without the complexity of creating a full 3D volume mesh.

1.4 Coupled finite element – boundary element method (FE-BEM)

Here, we present a hybrid method for modeling the diffusion equation, which combines the strengths of BEM in terms of reduced meshing dimensionality with FEM in terms of modeling optical property heterogeneity. The approach is akin to a tailored method for incorporating soft priors in a modified form in the forward model, itself, i.e. in modeling the light diffusion equation instead of within the image reconstruction formulation. The coupled FE-BEM scheme introduced here assumes homogeneous regions in certain tissue types, which are known to have low variation in functional parameters (e.g. fat) and heterogeneous distributions for other tissues such as tumors, which are known to have large variations in optical properties. The advantage of this technique over FEM is that it does not require volume discretization of the entire 3D domain, but only for tissues with known heterogeneity; surfaces will suffice for the rest of the tissues within the domain of interest. The advantage over BEM is that it can model heterogeneity in certain tissues whereas BEM assumes only piece-wise constant regions. We present an implementation of the coupled FE-BEM system for modeling light diffusion in 3D. Results are reported for light fluence distributions and frequency domain boundary measurements of intensity and phase as well as computational times for realistic tissue geometries and are compared to existing numerical models. The examples presented correspond to breast imaging, although the concept can be readily extended to other sites and applications such as brain and small animal imaging.

2. Methods

2.1 Introduction to Diffusion equation

The diffusion equation can be derived from Radiation transport equation under the assumption that light propagation is just linearly anisotropic [28]. This diffusion approximation has been commonly used to model light transport in tissues where scatter dominates over absorption and at distances more than several transport scattering lengths (transport scattering length = $\frac{1}{μ_{s}^{'}}$ , where $μ_{s}^{'}$ is the reduced scattering coefficient) from the source [29]. This model is given in the frequency domain as:

- \nabla . D (r) \nabla Φ (r, ω) + (μ_{a} (r) + \frac{i ω}{c}) Φ (r, ω) = q (r, ω)

where

Φ (r, ω)

is the photon density or fluence at position r in the bounded imaging domain Ω, D is the diffusion coefficient given by:

D (r) = \frac{1}{(3 (μ_{a} (r) + μ_{s}^{'} (r)))}

μ_{a}

is the absorption coefficient, ω is the frequency, and

q (r, ω)

is the isotropic source distribution. The source distribution is modeled as a point source located at a depth of one scattering distance inside the boundary where an optical fiber would be [30]. At the outer boundary of the domain, the relationship between photon fluence and flux is given by a Robin type boundary condition [30]:

Φ (r, ω) + \frac{D}{α} {\frac{\partial Φ}{\partial n} |}_{d Ω} = 0

Where α incorporates refractive index mismatch.

A coupled FE-BEM approach for the diffusion equation in multi-layered media was implemented by assuming homogeneous optical properties in outer layers and heterogeneous optical properties in the innermost tissue layer. Figure 2 shows a schematic of such a layered media illustrated in 2D for simplicity. In this domain, the exterior tissue (labeled I) was homogeneous and bounded by $Γ_{a}$ (containing N_a nodes) and $Γ_{b}$ (containing N_b nodes). $Γ_{b}$ also bounds an interior layer (labeled II) containing N_b nodes on the boundary and N_i nodes on the interior. In the coupled FE-BEM, BEM was used to model the exterior layers and FEM was used for the interior layer. These are discussed below in the context of the coupled system.

Fig. 2 Schematic of a two-layered region in 2D having homogeneous distribution of optical properties in region I and heterogeneous distribution in region II.

Download Full Size | PDF

2.2 Diffusion equation modeled with FEM

The Galerkin formulation was used for FEM where the orthogonality condition $〈 R, W_{i} 〉 = 0$ is satisfied [31]. Here R is the residual of Eq. (1), W_i is the weighting function and symbol $〈〉$ represents integration. Using linear basis functions $ϕ_{j}$ as the weighting function, we obtain the formulation for Eq. (1):

〈 - \nabla D \nabla Φ, ϕ_{j} 〉 + 〈 (μ_{a} + \frac{i ω}{c}) Φ, ϕ_{j} 〉 = 〈 q, ϕ_{j} 〉

The first term in Eq. (4) was integrated using Green’s theorem, to give:

〈 D \nabla Φ, \nabla ϕ_{j} 〉 - \oint D \frac{\partial Φ}{\partial n} ϕ_{j} d s + 〈 (μ_{a} + \frac{i ω}{c}) Φ, ϕ_{j} 〉 = 0

where the integration applies for interior tissues (region II in Fig. 2 bounded by

Γ_{b}

); note that the right hand side source contribution is zero since no source exists within the interior tissue region. Approximating

Φ = \sum_{i = 1}^{N_{v}} Φ_{i} ϕ_{i}

and

D \frac{\partial Φ}{\partial n} = \sum_{i = 1}^{N_{v}} D_{i} \frac{\partial Φ_{i}}{\partial n} ϕ_{i}

, using piece-wise linear basis functions

ϕ_{i}

and nodal values for fluence and flux Eq. (5) becomes:

(〈 D_{i} \nabla ϕ_{i}, \nabla ϕ_{j} 〉 + 〈 {(μ_{a} + \frac{i ω}{c})}_{i} ϕ_{i}, ϕ_{j} 〉) \sum_{i = 1}^{N_{v}} Φ_{i} = (\oint ϕ_{i} ϕ_{j} d s) \sum_{i = 1}^{N_{v}} D_{i} \frac{\partial Φ_{i}}{\partial n}

where N_v is the total number of volume nodes (N_v = N_b + N_i). Equation (6) can be written in matrix form as:

[A] (Φ) = [B] (D \frac{\partial Φ}{\partial n})

where

\begin{array}{l} A_{k l} = 〈 D \nabla ϕ_{k} \nabla ϕ_{l} 〉 + 〈 (μ_{a} + \frac{i ω}{c}) ϕ_{k} ϕ_{l} 〉 \\ B_{k l} = \oint ϕ_{k} ϕ_{l} \end{array}

Separating boundary (b) and interior (i) nodes of the inner region II, Eq. (7) expands as:

\begin{array}{l} (\begin{matrix} A_{b b} & A_{b i} \\ A_{i b} & A_{i i} \end{matrix}) {\begin{matrix} Φ_{b} \\ Φ_{i} \end{matrix}} = (\begin{matrix} B_{b b} & 0 \\ 0 & 0 \end{matrix}) {\begin{matrix} D {\frac{\partial Φ}{\partial n} |}_{Γ}_{b} \\ 0 \end{matrix}} \\ \Rightarrow {\begin{matrix} Φ_{b} \\ Φ_{i} \end{matrix}} = (\begin{matrix} A I_{b b} & A I_{b i} \\ A I_{i b} & A I_{i i} \end{matrix}) (\begin{matrix} B_{b b} & 0 \\ 0 & 0 \end{matrix}) {\begin{matrix} D {\frac{\partial Φ}{\partial n} |}_{b} \\ 0 \end{matrix}} \end{array}

where AI = A⁻¹. Φ_b can be obtained from

Φ_{b} = [A I_{b}_{b}] B_{b b} (D {\frac{\partial Φ}{\partial n} |}_{b})

this relationship between fluence

Φ_{b}

and flux

D {\frac{\partial Φ}{\partial n} |}_{b}

is applied within the BEM integral equation as described in the next section.

2.3 Diffusion equation modeled with BEM

Under the assumption that the tissue contains boundaries known a priori which separate into piece-wise constant homogeneous regions, the diffusion equation can be written in the form of a modified Helmholtz equation given in each region by [22,26]:

\nabla . D_{l} \nabla Φ - k_{l}^{2} Φ = - q_{0} (r, ω)

where

(μ_{a} (r) + \frac{i ω}{c}) = k_{l}^{2}

Here subscript l refers to the region label and applies to homogeneous region I in Fig. 2 bounded by $Γ_{a}$ and $Γ_{b}$ . The fundamental solution given by the Green’s function for Eq. (11) satisfies:

D_{l} \nabla^{2} G (r, r_{i}) - k_{l}^{2} G (r, r_{i}) = - δ (r - r_{i})

where

G_{i} (r, r i) = \frac{\exp (\frac{- k_{l} | r - r_{i} |}{\sqrt{D_{l}}})}{4 π D_{l} | r - r_{i} |}, 3 - D

The boundary integral form of Eq. (11) was derived using weighted residuals, Green’s third identity and the fundamental solutions [32] and appears as:

c_{i} Φ_{i} + \oint D_{l} \frac{\partial G_{i}}{\partial n} Φ - \oint D_{l} \frac{\partial Φ}{\partial n} G_{i} = 〈 q_{0}, G_{i} 〉

for the Green’s function which is singular in node i where , and Ω is the solid angle enclosed by the boundary at node i.

The photon fluence and flux are discretized using linear basis functions $ψ_{i}$ defined on the triangles of the surfaces, as and , where N_s is the number of boundary nodes on the surface (N_s = N_a + N_b). In discretized form, Eq. (15) becomes:

c_{i} Φ_{i} + \oint D_{l} \frac{\partial G_{i}}{\partial n} Φ d s - \oint D_{l} \frac{\partial Φ}{\partial n} G_{i} d s = 〈 q_{0}, G_{i} 〉

which can be written as matrix equation

[\tilde{A}] {Φ_{i}} - [\tilde{B}] {D_{l} \frac{\partial Φ}{\partial n}} = {{\tilde{Q}}_{i}}

where

\begin{array}{l} {\tilde{A}}_{i, j} = c_{i} δ_{i j} + \oint D_{l} \frac{\partial G_{i}}{\partial n} ψ_{j} d s \\ {\tilde{B}}_{i, j} = \oint G_{i} ψ_{j} d s \\ {\tilde{Q}}_{i} = 〈 q_{0}, G_{i} 〉 \end{array}

The Robin boundary condition specified in Eq. (3) is applied for the outer boundary. For multi-region problems, continuity conditions are enforced across the interior boundaries. For a two-region problem, the matrix form was derived by separating nodes on boundaries and as (see Appendix for details).

(\begin{array}{l} \begin{matrix} {\tilde{A}}_{a a} + α {\tilde{B}}_{a a} & {\tilde{A}}_{a b} & - {\tilde{B}}_{a b} \end{matrix} \\ \begin{matrix} {\tilde{A}}_{b a} + α {\tilde{B}}_{b a} & {\tilde{A}}_{b b} & - {\tilde{B}}_{b b} \end{matrix} \end{array}) {\begin{cases} Φ_{a} \\ Φ_{b} \\ {D \frac{\partial Φ}{\partial n} |}_{b} \end{cases}} = {\begin{cases} {\tilde{Q}}_{a} \\ {\tilde{Q}}_{b} \end{cases}}

Note from Eq.s 10 and 19 that both BEM and FEM formulations containing fluence on boundary nodes of interior tissue which couples the FEM and BEM system of equations.

2.4 Coupled FE-BEM for Diffusion equation

To derive the coupled FE-BEM formulation, we note that the fluence has to be the same whether derived from BEM or FEM for interior boundaries and the flux has to be continuous. This can be stated mathematically as:

\begin{array}{l} Φ_{b} (FEM) = Φ_{b} (BEM) \\ {D \frac{\partial Φ}{\partial n} |}_{b} (FEM) = - {D \frac{\partial Φ}{\partial n} |}_{b} (BEM) \end{array}

The negative sign for the flux is because the BEM formulation derived flux going outwards from region I into II, and FEM formulation has flux going into region I from II.

Using these relations and substituting for $Φ_{b}$ from Eq. (10) into Eq. (19) produces

\begin{array}{l} (\begin{array}{l} \begin{matrix} {\tilde{A}}_{a a} + α {\tilde{B}}_{a a} & {\tilde{A}}_{a b} & - {\tilde{B}}_{a b} \end{matrix} \\ \begin{matrix} {\tilde{A}}_{b a} + α {\tilde{B}}_{b a} & {\tilde{A}}_{b b} & - {\tilde{B}}_{b b} \end{matrix} \end{array}) {\begin{cases} Φ_{a} \\ - A I_{b b} B_{b b} {D \frac{\partial Φ}{\partial n} |}_{b} \\ {D \frac{\partial Φ}{\partial n} |}_{b} \end{cases}} = {\begin{cases} {\tilde{Q}}_{a} \\ {\tilde{Q}}_{b} \end{cases}} \\ \Rightarrow (\begin{matrix} {\tilde{A}}_{a a} + α {\tilde{B}}_{a a} & - {\tilde{A}}_{a b} A I_{b b} B_{b b} - {\tilde{B}}_{a b} \\ {\tilde{A}}_{b a} + α {\tilde{B}}_{b a} & - {\tilde{A}}_{b b} A I_{b b} B_{b b} - {\tilde{B}}_{b b} \end{matrix}) {\begin{cases} Φ_{a} \\ {D \frac{\partial Φ}{\partial n} |}_{b} \end{cases}} = {\begin{cases} {\tilde{Q}}_{a} \\ {\tilde{Q}}_{b} \end{cases}} \end{array}

This system was solved for fluence on the outer boundary and flux on inner boundary. The flux was used from this solution to solve the FEM equation [Eq. (9)] for interior field. Also note that matrix A has already been inverted when solving Eq. (10), so this step is straightforward. The size of the matrix to be inverted in Eq. (21) is N_s x N_s. Equation (21) represents a two-region problem but the approach is easily extended to multiple regions as shown in the Appendix. The coupled FE-BEM equations were implemented in Matlab and C and used to generate fluence distributions in the domain.

2.5 Simulation setup

Realistic breast-shaped imaging domains were generated using a clinical MRI data set collected from a female volunteer diagnosed with infiltrating ductal carcinoma as part of an ongoing clinical trial with MRI/optical imaging. A 3T Phillips scanner was used to collect the MRI and contrast-enhanced MR data sets. Using the MR volume, image segmentation of adipose, fibroglandular and tumor tissues was performed with the use of software package Mimics^TM [33]. In addition, spherical inclusions were also simulated within the outer breast region. Using these geometries, six test cases of multiple regions were created for the simulations as shown in Fig. 3 . The volume meshes for interior tissues of interest were generated with the same software. Combining these surfaces and volumes provided meshes for the coupled FE-BEM. The corresponding mesh sizes are given in Table 1 for each of the test cases.

Fig. 3 Surface renderings of the six test cases used in this study are shown, with two-three regions created. Clockwise from top left, the six test cases show cases (1) the outer breast contour and tumor created from clinical MRI (2) outer breast and simulated spherical inclusion (3) Outer breast, sphere and tumor (4) Outer breast, larger sphere and tumor (5) Outer breast and two spherical inclusions and (6) Outer breast, fibroglandular and tumor tissues.

Download Full Size | PDF

Table 1. Mesh sizes for the different test cases used in the simulations. The first two columns of mesh sizes correspond to the coupled FE-BEM and the last two columns correspond to mesh sizes for BEM and FEM

View Table

To compare the results from the coupled FE-BEM, forward data was also generated using BEM and FEM techniques both of which have been validated previously [34,35]. For the BEM, only surfaces were required, and multiple homogeneous regions were simulated. For the FEM, a full 3D volume mesh was required with the interior boundaries preserved for consistency. The volume meshes for each of the test cases were created with a 3D pixel-based mesh generator [36], which used the average edge size from the surfaces for generating the volume mesh. A schematic of such a mesh is shown in Fig. 1 (last step). Mesh sizes used in BEM and FEM only reconstructions are also given in Table 1. The meshes for testing all three models were of comparable mesh resolutions and with interior boundaries preserved. The computer time for volume mesh generation varied from 260 seconds to 323 seconds. The source-detector geometry for the imaging domains contained sixteen sources with fifteen detectors per source in a circular ring around the periphery of the breast, giving a total of 240 measurements [4]. The fiber indentations for the sixteen locations can be seen in the surface rendering (Fig. 3).

3. Results

3.1 Photon fluence distribution from coupled FE-BEM

The coupled FE-BEM was applied to generate the photon fluence in the six test cases shown in Fig. 3. In the simulation, both the exterior and interior tissues had homogeneous distributions of optical properties where $μ_{a}$ = 0.006 mm⁻¹ and $μ_{s}^{'}$ = 1.0 mm⁻¹ for outer region(s) and $μ_{a}$ = 0.02 and $μ_{s}^{'}$ = 2.0 mm⁻¹ in the interior tissue. The logarithm of fluence distribution at the boundaries of the tissues for a single source is shown in Fig. 4 for test cases 1 and 6 where the diffusive pattern typically expected from the diffusion equation is seen.

Fig. 4 Logarithm of photon fluence obtained using coupled FE-BEM for a single source in test cases 1, and 6. Left: Results from test case 1 showing outer boundary; Middle: inner tumor boundary by making outer surface transparent; Right: Results from test case 6 for inner tissues.

Download Full Size | PDF

3.2 Comparison of boundary data using BEM, FEM and Coupled FE-BEM

To compare the results from the coupled FE-BEM with existing models, the boundary data at detector locations were computed. The logarithm of intensity and phase is shown in Fig. 5 at the boundary detector locations for 240 measurements (16 sources x 15 detectors/source) generated using the three models (BEM, FEM and coupled FE-BEM) for test case 1. The measurements show good agreement with RMS differences in logarithm of intensity between BEM and the coupled model of less than 0.1 and in phase of less than 1 degree. The RMS differences between FEM and the coupled model was less than 0.5 for logarithm of intensity and 2.5 degrees for phase. These differences are likely due to the differences in the mesh types and discretization.

Fig. 5 Comparison of (a) logarithm of intensity and (b) phase at the detector locations on the boundary ( = 240 measurement points) obtained from BEM, FEM and coupled FE-BEM for test case 1.

Download Full Size | PDF

3.3 Modeling heterogeneity

One of the drawbacks of BEM is that it cannot model heterogeneity of tissue due to the inherent assumption in the model: the Diffusion equation only reduces to modified Helmholtz in BEM formulation for piece-wise constant or homogeneous regions. For modeling heterogeneity, the coupled model offers an alternative solution. To illustrate the change in fluence with increasing heterogeneity, a cross-section along the center of the inner sphere in test case 2 is shown in Fig. 6 for a single source. The left column indicates the property distribution and right column shows the corresponding logarithm of fluence distribution for a (1) homogeneous domain (sphere to background contrast of 1:1), (2) heterogeneous domain (2:1 sphere to background contrast) and (3) heterogeneous domain with spatially varying contrast in the sphere (2:1 varying with background). As the heterogeneity in the absorption increases, a decrease in fluence is observed in parts of the sphere, as expected. A decrease in intensity also occurred at the boundary as a result of the heterogeneity.

Fig. 6 2-D cross-sections along the center of the interior spherical inclusion in test case 2 for $μ_{a}$ (left column) and logarithm of fluence (right column). The background was always homogeneous. Top row shows cross-section of sphere for a homogeneous domain (1:1 contrast between sphere and background), Middle row shows 2:1 contrast between sphere and background and bottom row shows a spatially varying distribution in the sphere (2:1 varying). As expected the fluence decreases with increasing heterogeneity.

Download Full Size | PDF

3.4 Analysis of computational times between Coupled FE-BEM and FEM

The computational time required by coupled FE-BEM was a function of the surface mesh size and was found to scale as N_s^3.2, where N_s is the number of nodes in the surface mesh. This outcome was expected given that the matrix assembly and solving the BEM component of the coupled model consumed the most time and the BEM was found to scale with surface node size as N_s^3.5. The scaling was obtained for the two region and three region problems in complex domains presented here, but was smaller (N_s^2.7 quoted previously for BEM [22]) in simple two region domains. The FEM component of the coupled model consumed less than 0.5% of the total time.

Since the computational time of coupled FE-BEM scales with surface mesh size, it is reasonable to assume that the speed-up of the coupled model when compared to stand-alone FEM will be a function of the ratio of the number of surface to full 3D volume nodes (N_s/N). Figure 7 (top row) shows a plot of the ratio of computational times of coupled FE-BEM to FEM time, as a function of N_s/N, the values for N_s and N can be found in Table 1 (first column and last columns respectively). The plot shows that for ratio of N_s/N < 20%, coupled FE-BEM was faster (ratio of times < 1) whereas for N_s/N > 20%, stand-alone FEM was faster (ratio of times > 1). This data did not include the computational time for creating a large 3D volume mesh for FEM. It is important to note that when the meshing time for FEM was included, coupled FE-BEM was always faster than FEM (ratio < 1) for the cases presented here (ratio of times ranged from 0.14 to 0.92).

Fig. 7 Ratio of computational time of coupled FE-BEM to stand-alone FEM for the six test cases, plotted as a function of % surface to volume nodes (top) from the respective meshes (N_s/N) where Ns is the number of boundary nodes in the coupled mesh and N is the number of nodes in the FEM mesh and % surface area to volume ratio (bottom) of the total tissue domain.

Download Full Size | PDF

Since the metric (N_s/N) requires a volume mesh to be created, we also chose the physical surface area to volume ratio (SA/V) as another metric for comparing computational times, and can be obtained from image segmentation. Figure 7 (bottom row) shows that the coupled model was faster than FEM (ratio of times < 1) when SA/V < 10%. These plots illustrate that we can use quantitative metrics to determine the most efficient 3D forward model for the imaging domain under consideration.

3.5 Analysis of computational times between Coupled FE-BEM and BEM

A similar comparison was performed for the ratio of computational times of coupled FE-BEM and BEM. Since the number of surface nodes was the same for the coupled FE-BEM and BEM models (See Table 1), the time differences depend on the total number of volume nodes used in the interior tissue region (N_v = N_b + N_i) as compared to the surface nodes (N_b) on the boundary in the same region (see region II in schematic of Fig. 2). For small N_b/N_v, the volume nodes dominate such that coupled FE-BEM was longer to compute than BEM. For larger N_b/N_v, surface nodes dominate and hence coupled FE-BEM was faster than BEM. Overall, the differences in the two models were less than 10% for the test cases presented here (see Fig. 8 , top row). A ratio of 50% N_b/N_v appeared to be the delineating value. Similarly, a ratio of 20% appeared to separate the two models in terms of ratio of interior tissue surface area (ISA) to interior tissue volume (IV), see Fig. 8 (bottom).

Fig. 8 Ratio of computational time of coupled FE-BEM model to BEM for the six test cases, plotted as a function of % surface to volume nodes (top) of the interior tissue (N_b/N_v) where N_b is the number of nodes on boundary of interior tissue and N_v is the number of volume nodes of interior tissue, and % surface area to volume ratio (ISA/IV) (bottom) of the interior tissue domain.

Download Full Size | PDF

4. Discussion and Conclusions

Coupled FE-BEM methods have been used extensively in other fields such as electrostatics [37], electromagnetics [38] and in biomedical applications to model cardiac tissue [39]; among others, Here we present application of this technique to diffuse optical tomography. The coupled FE-BEM method provides an elegant solution to the practical problem in multi-modality optical imaging of how to model heterogeneity in tissues whose boundaries are known, without complex volumetric meshing of the full 3D domain. In this method, the volume meshing has not been eliminated, but rather the size of the domains were reduced for which it is needed. Therefore, this has an impact on both the meshing time as well as the computational time for the forward solver.

Different implementation options exist [40], and we chose one does not change the bandwidth of the matrices involved. Specifically, the sparsity of the FEM matrix, which is a highly desirable aspect of finite elements, was not altered. No increase in the size of dense BEM matrix to be solved occurred as well. The computational time of the coupled method was governed primarily by the BEM matrix size (> 99% of total time) for the domains described here. This will likely change for larger volumetric FEM computations within the domain, or larger areas of heterogeneity, but is not anticipated in the current application. Comparison to existing and validated numerical models based on FEM alone and BEM alone showed good agreement with RMS differences of less than 0.5 in logarithm of intensity and less than 2.5 degrees in phase.

The coupled FE-BEM method incorporates the idea of soft priors directly into the forward model itself, which is different from traditional techniques where regularization is used in the image reconstruction or inverse problem. The choice of numerical technique for the forward model will depend on the problem, the imaging domain and its approximations with respect to homogeneity/heterogeneity. These a priori assumptions when used intelligently can greatly influence the choice of the model to be used. We have shown that the coupled FE-BEM is faster than FEM when the surface to volume node ratio was less than 20% and when the total surface area to volume was less than 10%. However, when meshing time was included, the coupled FE-BEM was always faster and the ratio of computational times (Coupled / FEM) ranged from 0.14 to 0.92. Coupled FE-BEM was comparable to BEM ( ± 10%) for the range of mesh sizes and tissue types examined here. We have presented results from realistic breast-shaped models in these simulations. While the results presented here are from breast geometries, the model can be applied to other tissue regions as well.

In conclusion, a coupled FE-BEM method was implemented for modeling light diffusion in 3D for multi-modality optical imaging systems and the results show good agreement with existing numerical models but utilize a fraction of the volume mesh size required by corresponding FEM techniques.

5. Appendix

Equation (17) describes the matrix form of the BEM for a single region. For an external region consisting of boundaries a and b, in region I, the matrix formulation extension of Eq. (17) is

(\begin{matrix} {\tilde{A}}_{a a} & {\tilde{A}}_{a b} \\ {\tilde{A}}_{b a} & {\tilde{A}}_{b b} \end{matrix}) {\begin{cases} Φ_{a} \\ Φ_{b} \end{cases}} - (\begin{matrix} {\tilde{B}}_{a a} & {\tilde{B}}_{a b} \\ {\tilde{B}}_{b a} & {\tilde{B}}_{b b} \end{matrix}) {\begin{cases} D_{I} {\frac{\partial Φ}{\partial n} |}_{a} \\ D_{I} {\frac{\partial Φ}{\partial n} |}_{a}_{b} \end{cases}} = {\begin{cases} {\tilde{Q}}_{a} \\ {\tilde{Q}}_{b} \end{cases}}

(\begin{array}{l} \begin{matrix} {\tilde{A}}_{a a} & {\tilde{A}}_{a b} & - {\tilde{B}}_{a a} & - {\tilde{B}}_{a b} \end{matrix} \\ \begin{matrix} {\tilde{A}}_{b a} & {\tilde{A}}_{b b} & - {\tilde{B}}_{b a} & - {\tilde{B}}_{b b} \end{matrix} \end{array}) {\begin{matrix} Φ_{a} \\ Φ_{b} \\ D_{I} {\frac{\partial Φ}{\partial n} |}_{a} \\ D_{I} {\frac{\partial Φ}{\partial n} |}_{b} \end{matrix}} = {\begin{cases} {\tilde{Q}}_{a} \\ {\tilde{Q}}_{b} \end{cases}}

Substituting the boundary condition in Eq. (3) for the outer boundary, Eq. (23) becomes:

\begin{array}{l} (\begin{array}{l} \begin{matrix} {\tilde{A}}_{a a} & {\tilde{A}}_{a b} & - {\tilde{B}}_{a a} & - {\tilde{B}}_{a b} \end{matrix} \\ \begin{matrix} {\tilde{A}}_{b a} & {\tilde{A}}_{b b} & - {\tilde{B}}_{b a} & - {\tilde{B}}_{b b} \end{matrix} \end{array}) {\begin{matrix} Φ_{a} \\ Φ_{b} \\ - α Φ_{a} \\ D_{I} {\frac{\partial Φ}{\partial n} |}_{b} \end{matrix}} = {\begin{cases} {\tilde{Q}}_{a} \\ {\tilde{Q}}_{b} \end{cases}} \\ \Rightarrow (\begin{array}{l} \begin{matrix} {\tilde{A}}_{a a} + α {\tilde{B}}_{a a} & {\tilde{A}}_{a b} & - {\tilde{B}}_{a b} \end{matrix} \\ \begin{matrix} {\tilde{A}}_{b a} + α {\tilde{B}}_{b a} & {\tilde{A}}_{b b} & - {\tilde{B}}_{b b} \end{matrix} \end{array}) {\begin{matrix} Φ_{a} \\ Φ_{b} \\ D_{I} {\frac{\partial Φ}{\partial n} |}_{b} \end{matrix}} = {\begin{cases} {\tilde{Q}}_{a} \\ {\tilde{Q}}_{b} \end{cases}} \end{array}

which yields Eq. (19). For successive layers bounded by a, b and c, the matrix for BEM is

(\begin{matrix} {\tilde{A}}_{a a I} + α {\tilde{B}}_{a a I} & {\tilde{A}}_{a b I} & - {\tilde{B}}_{a b I} \\ \begin{array}{l} {\tilde{A}}_{b a} + α {\tilde{B}}_{b a I} \end{array} & \begin{array}{l} {\tilde{A}}_{b b I} \\ {\tilde{A}}_{b b I I} \\ {\tilde{A}}_{c b I I} \end{array} & \begin{array}{l} - {\tilde{B}}_{b b I} \\ {\tilde{B}}_{b b I I} \\ {\tilde{B}}_{c b I I} \end{array} & \begin{array}{l} {\tilde{A}}_{b c I I} \\ {\tilde{A}}_{c c I I} \end{array} & \begin{array}{l} - {\tilde{B}}_{b c I I} \\ - {\tilde{B}}_{c c I I} \end{array} \end{matrix}) {\begin{matrix} Φ_{a I} \\ Φ_{b I} \\ \begin{array}{l} D_{I} {\frac{\partial Φ}{\partial n} |}_{b I} \\ \begin{matrix} Φ_{c I I} \\ D_{I I} {\frac{\partial Φ}{\partial n} |}_{c I I} \end{matrix} \end{array} \end{matrix}} = {\begin{cases} {\tilde{Q}}_{a I} \\ {\tilde{Q}}_{b I} \\ 0 \\ 0 \end{cases}}

and the FEM relationship is given for an interior region as which is used along with continuity conditions to derive the coupled FE-BEM given by:

\begin{array}{l} (\begin{matrix} {\tilde{A}}_{a a I} + α {\tilde{B}}_{a a I} & {\tilde{A}}_{a b I} & - {\tilde{B}}_{a b I} \\ \begin{array}{l} {\tilde{A}}_{b a} + α {\tilde{B}}_{b a I} \end{array} & \begin{array}{l} {\tilde{A}}_{b b I} \\ {\tilde{A}}_{b b I I} \\ {\tilde{A}}_{c b I I} \end{array} & \begin{array}{l} - {\tilde{B}}_{b b I} \\ {\tilde{B}}_{b b I I} \\ {\tilde{B}}_{c b I I} \end{array} & \begin{array}{l} {\tilde{A}}_{b c I I} \\ {\tilde{A}}_{c c I I} \end{array} & \begin{array}{l} - {\tilde{B}}_{b c I I} \\ - {\tilde{B}}_{c c I I} \end{array} \end{matrix}) {\begin{matrix} Φ_{a I} \\ Φ_{b I} \\ \begin{array}{l} D_{I} {\frac{\partial Φ}{\partial n} |}_{b I} \\ \begin{matrix} - A I_{c c} B_{c c} D_{I I} {\frac{\partial Φ}{\partial n} |}_{c I I} \\ D_{I I} {\frac{\partial Φ}{\partial n} |}_{c I I} \end{matrix} \end{array} \end{matrix}} = {\begin{cases} {\tilde{Q}}_{a I} \\ {\tilde{Q}}_{b I} \\ 0 \\ 0 \end{cases}} \\ \Rightarrow (\begin{matrix} {\tilde{A}}_{a a I} + α {\tilde{B}}_{a a I} & {\tilde{A}}_{a b I} & - {\tilde{B}}_{a b I} \\ {\tilde{A}}_{b a} + α {\tilde{B}}_{b a I} & {\tilde{A}}_{b b I} & - {\tilde{B}}_{b b I} \\ {\tilde{A}}_{b b I I} & {\tilde{B}}_{b b I I} & - {\tilde{A}}_{b c I I} A I_{c c} B_{c c} - {\tilde{B}}_{b c I I} \\ {\tilde{A}}_{c b I I} & {\tilde{B}}_{c b I I} & - {\tilde{A}}_{c c I I} A I_{c c} B_{c c} - {\tilde{B}}_{c c I I} \end{matrix}) {\begin{matrix} Φ_{a I} \\ Φ_{b I} \\ D_{I} {\frac{\partial Φ}{\partial n} |}_{b I} \\ D_{I I} {\frac{\partial Φ}{\partial n} |}_{c I I} \end{matrix}} = {\begin{cases} {\tilde{Q}}_{a I} \\ {\tilde{Q}}_{b I} \\ 0 \\ 0 \end{cases}} \end{array}

Acknowledgments

This work was funded by NIH-NIBIB grant R01EB007966 and patient data collected through NCI grant P01CA080139.

References and links

1. R. Weissleder and M. J. Pittet, “Imaging in the era of molecular oncology,” Nature 452(7187), 580–589 (2008). [CrossRef] [PubMed]

2. A. Cerussi, D. Hsiang, N. Shah, R. Mehta, A. Durkin, J. Butler, and B. J. Tromberg, “Predicting response to breast cancer neoadjuvant chemotherapy using diffuse optical spectroscopy,” Proc. Natl. Acad. Sci. U.S.A. 104(10), 4014–4019 (2007). [CrossRef] [PubMed]

3. M. V. Schulmerich, J. H. Cole, K. A. Dooley, M. D. Morris, J. M. Kreider, S. A. Goldstein, S. Srinivasan, and B. W. Pogue, “Noninvasive Raman tomographic imaging of canine bone tissue,” J. Biomed. Opt. 13(2), 020506 (2008). [CrossRef]

4. C. M. Carpenter, B. W. Pogue, S. Jiang, H. Dehghani, X. Wang, K. D. Paulsen, W. A. Wells, J. Forero, C. Kogel, J. B. Weaver, S. P. Poplack, and P. A. Kaufman, “Image-guided optical spectroscopy provides molecular-specific information in vivo: MRI-guided spectroscopy of breast cancer hemoglobin, water, and scatterer size,” Opt. Lett. 32(8), 933–935 (2007). [CrossRef] [PubMed]

5. M. J. Niedre, R. H. de Kleine, E. Aikawa, D. G. Kirsch, R. Weissleder, and V. Ntziachristos, “Early photon tomography allows fluorescence detection of lung carcinomas and disease progression in mice in vivo,” Proc. Natl. Acad. Sci. U.S.A. 105(49), 19126–19131 (2008). [CrossRef] [PubMed]

6. M. Schweiger and S. R. Arridge, “Comparison of two- and three-dimensional reconstruction methods in optical tomography,” Appl. Opt. 37(31), 7419–7428 (1998). [CrossRef] [PubMed]

7. P. K. Yalavarthy, D. R. Lynch, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Implementation of a computationally efficient least-squares algorithm for highly under-determined three-dimensional diffuse optical tomography problems,” Med. Phys. 35(5), 1682–1697 (2008). [CrossRef] [PubMed]

8. M. Schweiger, O. Dorn, A. Zacharopoulos, I. Nissila, and S. R. Arridge, “3D level set reconstruction of model and experimental data in Diffuse Optical Tomography,” Opt. Express 18(1), 150–164 (2010). [CrossRef] [PubMed]

9. G. Boverman, E. L. Miller, D. H. Brooks, D. Isaacson, Q. Fang, and D. A. Boas, “Estimation and statistical bounds for three-dimensional polar shapes in diffuse optical tomography,” IEEE Trans. Med. Imaging 27(6), 752–765 (2008). [CrossRef] [PubMed]

10. A. D. Zacharopoulos, M. Schweiger, V. Kolehmainen, and S. R. Arridge, “3D shape based reconstruction of experimental data in Diffuse Optical Tomography,” Opt. Express 17(21), 18940–18956 (2009). [CrossRef] [PubMed]

11. S. Srinivasan, B. W. Pogue, H. Dehghani, F. Leblond, and X. Intes, “Data subset algorithm for computationally efficient reconstruction of 3-D spectral imaging in diffuse optical tomography,” Opt. Express 14(12), 5394–5410 (2006). [CrossRef] [PubMed]

12. Q. Zhang, T. J. Brukilacchio, A. Li, J. J. Stott, T. Chaves, E. Hillman, T. Wu, M. Chorlton, E. Rafferty, R. H. Moore, D. B. Kopans, and D. A. Boas, “Coregistered tomographic x-ray and optical breast imaging: initial results,” J. Biomed. Opt. 10(2), 024033–0240339 (2005). [CrossRef] [PubMed]

13. S. C. Davis, K. S. Samkoe, J. A. O’Hara, S. L. Gibbs-Strauss, H. L. Payne, P. J. Hoopes, K. D. Paulsen, and B. W. Pogue, “MRI-coupled fluorescence tomography quantifies EGFR activity in brain tumors,” Acad. Radiol. 17(3), 271–276 (2010). [CrossRef] [PubMed]

14. R. B. Schulz, A. Ale, A. Sarantopoulos, M. Freyer, E. Soehngen, M. Zientkowska, and V. Ntziachristos, “Hybrid system for simultaneous fluorescence and x-ray computed tomography,” IEEE Trans. Med. Imaging 29(2), 465–473 (2010). [CrossRef] [PubMed]

15. C. Li, G. Wang, J. Qi, and S. R. Cherry, “Three-dimensional fluorescence optical tomography in small-animal imaging using simultaneous positron-emission-tomography priors,” Opt. Lett. 34(19), 2933–2935 (2009). [CrossRef] [PubMed]

16. K. D. Paulsen and H. Jiang, “Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization,” Appl. Opt. 35(19), 3447–3458 (1996). [CrossRef]

17. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, C. M. Carpenter, S. Jiang, and K. D. Paulsen, “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express 15(13), 8043–8058 (2007). [CrossRef] [PubMed]

18. Z. Yuan, Q. Zhang, E. S. Sobel, and H. Jiang, “Tomographic x-ray-guided three-dimensional diffuse optical tomography of osteoarthritis in the finger joints,” J. Biomed. Opt. 13(4), 044006 (2008). [CrossRef] [PubMed]

19. D. Hyde, E. L. Miller, D. H. Brooks, and V. Ntziachristos, “Data specific spatially varying regularization for multimodal fluorescence molecular tomography,” IEEE Trans. Med. Imaging 29(2), 365–374 (2010). [CrossRef] [PubMed]

20. C. M. Carpenter, S. Srinivasan, B. W. Pogue, and K. D. Paulsen, “Methodology development for three-dimensional MR-guided near infrared spectroscopy of breast tumors,” Opt. Express 16(22), 17903–17914 (2008). [CrossRef] [PubMed]

21. G. Boverman, E. L. Miller, A. Li, Q. Zhang, T. Chaves, D. H. Brooks, and D. A. Boas, “Quantitative spectroscopic diffuse optical tomography of the breast guided by imperfect a priori structural information,” Phys. Med. Biol. 50(17), 3941–3956 (2005). [CrossRef] [PubMed]

22. S. Srinivasan, B. W. Pogue, C. Carpenter, P. K. Yalavarthy, and K. D. Paulsen, “A boundary element approach for image-guided near-infrared absorption and scatter estimation,” Med. Phys. 34(11), 4545–4557 (2007). [CrossRef] [PubMed]

23. X. Zhang, C. T. Badea, and G. A. Johnson, “Three-dimensional reconstruction in free-space whole-body fluorescence tomography of mice using optically reconstructed surface and atlas anatomy,” J. Biomed. Opt. 14(6), 064010 (2009). [CrossRef] [PubMed]

24. A. Custo, D. A. Boas, D. Tsuzuki, I. Dan, R. Mesquita, B. Fischl, W. E. Grimson, and W. Wells 3rd, “Anatomical atlas-guided diffuse optical tomography of brain activation,” Neuroimage 49(1), 561–567 (2010). [CrossRef] [PubMed]

25. B. W. Pogue, S. Jiang, H. Dehghani, C. Kogel, S. Soho, S. Srinivasan, X. Song, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “Characterization of hemoglobin, water, and NIR scattering in breast tissue: analysis of intersubject variability and menstrual cycle changes,” J. Biomed. Opt. 9(3), 541–552 (2004). [CrossRef] [PubMed]

26. J. Sikora, A. Zacharopoulos, A. Douiri, M. Schweiger, L. Horesh, S. R. Arridge, and J. Ripoll, “Diffuse photon propagation in multilayered geometries,” Phys Med Biol (2006).

27. P. Vaupel, F. Kallinowski, and P. Okunieff, “Blood flow, oxygen and nutrient supply, and metabolic microenvironment of human tumors: a review,” Cancer Res. 49(23), 6449–6465 (1989). [PubMed]

28. A. Ishimaru, Wave propagation and scattering in random media (Academic Press, Inc., New York, 1978), Vol. 1.

29. M. S. Patterson, B. C. Wilson, and D. R. Wyman, “The propagation of optical radiation in tissue I. models of radiation transport and their application,” Lasers Med. Sci. 6(2), 155–168 (1991). [CrossRef]

30. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20(2), 299–309 (1993). [CrossRef] [PubMed]

31. D. R. Lynch, Numerical partial differential equations for environmental scientists and engineers (Springer, 2005).

32. C. A. Brebbia and J. Dominguez, Boundary elements: an introductory course.

33. http://materialise.com/materialise/view/en/92078-Mimics.html, retrieved.

34. S. Srinivasan, C. Carpenter, B. W. Pogue, and K. D. Paulsen, “Image-guided near infrared spectroscopy using boundary element method: phantom validation,” in SPIE 2009 BiOS Biomedical Optics Symposium: Multimodal Biomedical Imaging IV, 2009)

35. H. Dehghani, S. Srinivasan, B. W. Pogue, and A. Gibson, “Numerical modelling and image reconstruction in diffuse optical tomography,” Philos. Transact. A Math. Phys. Eng. Sci. 367(1900), 3073–3093 (2009). [CrossRef] [PubMed]

36. H. Ghadyani, J. M. Sullivan Jr, and Z. Wu, “Boundary recovery for delaunay tetrahedral meshes using local topological transformations,” Finite Elem. Anal. Des. 46(1-2), 74–83 (2010). [CrossRef] [PubMed]

37. C. P. Bradley, G. M. Harris, and A. J. Pullan, “The computational performance of a high-order coupled FEM/BEM procedure in electropotential problems,” IEEE Trans. Biomed. Eng. 48(11), 1238–1250 (2001). [CrossRef] [PubMed]

38. H. Ammari, ed., Modeling and computations in electromagnetics (Springer, 2007).

39. G. Fischer, B. Tilg, R. Modre, G. J. Huiskamp, J. Fetzer, W. Rucker, and P. Wach, “A bidomain model based BEM-FEM coupling formulation for anisotropic cardiac tissue,” Ann. Biomed. Eng. 28(10), 1229–1243 (2000). [CrossRef] [PubMed]

40. K. D. Paulsen and W. Liu, “Memory and operations count scaling of coupled finite element and boundary element systems of equations,” Int. J. Numer. Methods Eng. 33(6), 1289–1303 (1992). [CrossRef]

Test Case #	#Surface Nodes	# Volume Nodes	# Nodes BEM Mesh	# Nodes FEM Mesh
1	6471	798	6471	70423
2	6869	2171	6869	61468
3	7346	798	7346	65949
4	8297	798	8297	50203
5	8695	2171	8695	51041
6	11415	798	11415	50243

A coupled finite element-boundary element method for modeling Diffusion equation in 3D multi-modality optical imaging

Abstract

1. Introduction

1.1 Introduction to 3D diffuse optical imaging

1.2 Multi-modality optical imaging reconstruction techniques

1.3 Need for efficient 3D technique for complex 3D domains

1.4 Coupled finite element – boundary element method (FE-BEM)

2. Methods

2.1 Introduction to Diffusion equation

2.2 Diffusion equation modeled with FEM

2.3 Diffusion equation modeled with BEM

2.4 Coupled FE-BEM for Diffusion equation

2.5 Simulation setup

3. Results

3.1 Photon fluence distribution from coupled FE-BEM

3.2 Comparison of boundary data using BEM, FEM and Coupled FE-BEM

3.3 Modeling heterogeneity

3.4 Analysis of computational times between Coupled FE-BEM and FEM

3.5 Analysis of computational times between Coupled FE-BEM and BEM

4. Discussion and Conclusions

5. Appendix

Acknowledgments

References and links

Cited By

Figures (8)

Tables (1)

Equations (26)

Biomedical Optics Express