Image reconstruction in diffuse optical tomography based on simplified spherical harmonics approximation

Michael Chu; Hamid Dehghani

doi:10.1364/OE.17.024208

1. Introduction

Near infrared (NIR) optical tomography is a non-invasive imaging modality in which the optical properties within a volume of interest can be reconstructed using measured transmission and reflectance NIR data. These calculated optical maps can then be used to derive functional and structural information about the tissue being imaged [1–3].

In recent years, several NIR imaging systems have been developed with specific applications such as the detection and characterization of breast cancer [1,4–7] and functional imaging of the brain [8–12]. More recently the use of biological markers for function specific optical imaging has led to the development of novel imaging systems and algorithms that aim to recover images of function specific activity within small animals using, for example, fluorescence [13–15] or bioluminescence markers [16–20].

Most NIR imaging systems rely on the use of a model based image reconstruction algorithm in which the light distribution within the imaging domain is approximated to provide a match to measured boundary data [21–23] and therefore providing an estimate to the optical parameters within the medium being imaged. The most commonly applied model for the estimation of NIR light propagation in tissue is achieved through the Diffusion Approximation (DA) to the Radiative Transport Equation (RTE) through its robustness in implementation and computational speed and flexibility.

In order to accurately apply the DA, scattering effects within the medium must be dominant over the absorption, i.e. μ_s’>>μ_a, where μ_s’ is the reduced scattering coefficient and μ_a is the absorption coefficient. Whilst these conditions are generally applicable in most soft tissue, the presence of non-scattering regions, such as the cerebrospinal fluid found within the brain, or small source-detector separation, as encountered in small animal imaging, can mean that the DA is not valid for all cases [24–28].

To accurately model NIR light propagation in small geometries, or models with high absorption (such as for bioluminescence imaging), higher ordered forward models are needed whereby the problem is no longer limited by the ‘diffuse’ approximation. Whilst Monte Carlo simulations which are based on probabilistic models can produce accurate results, they tend to be computationally slow to be of clinical use. Alternatively, numerical approximations to the RTE can be used whereby light propagation is estimated through well established particle transport models based on complex integral-differential equations, which can account for both spatial and directional component of light propagation. The spherical harmonics (P_N) approximation of the RTE [29], for example, expands the angular components of radiance into a series of spherical harmonics. This method has been successfully applied in optical imaging but although more appropriate in conditions where the DA is not valid, it is still less than ideal due to the computational cost incurred [29–32]. The discrete ordinates method (S_N) is also a common implementation of the RTE often encountered in nuclear physics which discretises the angular components into a number of discrete directions. This method has also been used successfully in optical imaging application but as in the latter case results in a much higher computational cost as compared to the DA [27,33,34].

More recently, the application of Simplified Spherical Harmonics (SP_N) approximation to the RTE has been applied and studied [35,36]. The use of SP_N methods have been demonstrated to give accurate solutions for small geometries and in cases where the source/detector separation is small. The SP_N method has the advantage of requiring (N + 1)/2 equations where N is the number of Legendre Polynomials. This is in contrast to (N + 1)² equations for the P_N approximation, and N(N + 2) equations for the S_N method, where N is the number of allowed directions.

Due to the computational complexities outlined above, the vast majority of existing model based reconstruction algorithms make use of the DA. The image reconstructions therefore aim to form images of the absorption coefficient, μ_a and the diffusion coefficient κ = 1/3(μ_a + μ_s’) where μ_s’ = (1-g)μ_s and μ_s and g are the scattering coefficient and anisotropic factor respectively. One drawback of the application of the DA is that through the introduction of the reduced scattering coefficient, information about the scattering phase function is lost [32] and assumed to be constant at g = 0.9 which is accepted as valid for soft tissue. Little work has been done to investigate the effect of varying anisotropy factor within image reconstruction and whether it is at all possible to approximate its value within the inverse problem. The effect of anisotropic structures that affect the propagation of near-infrared photons according to their direction is of current interest and structures have been experimentally shown to cause direction dependence of light propagation in cases of chicken breast tissue and human skin [37, 38]. The theoretical considerations have so far concentrated on developing ways to model the anisotropic propagation of NIR light, and less emphasis has been put on exploring the implications of anisotropic effects in image reconstruction.

In this study, SP_N based models are used to study the effects of anisotropic scattering factor on image reconstruction. More specifically the use of SP_N methods have been utilized to evaluate the sensitivity of measured boundary data from a frequency domain (FD) measurement system (where amplitude and phase of the fluence are used at typically 100 MHz modulation frequency) to not only optical absorption and scatter but also anisotropy factor. Through calculation of error-norm (residual) from a theoretical model, the uniqueness of this problem is demonstrated. Furthermore, a framework for the use of SP_N methods in image reconstruction is presented.

2. Theory

2.1 Forward model

The forward problem is modeled using the SP_N approximation to the RTE. The SP_N equations can be derived by taking the planar spherical harmonics approximation and replacing the 1D derivatives with their 3D counterparts [35,36]. The SP₇ approximation results in a set of four coupled equations:

\begin{array}{l} - \nabla . \frac{1}{3 μ_{a 1}} \nabla ϕ_{1} + μ_{a} ϕ_{1} = \\ Q + \frac{2}{3} μ_{a} ϕ_{2} - \frac{8}{15} μ_{a} ϕ_{3} + \frac{16}{35} μ_{a} ϕ_{4} \end{array}

\begin{array}{l} - \nabla . \frac{1}{7 μ_{a 3}} \nabla ϕ_{2} + (\frac{4}{9} μ_{a} + \frac{5}{9} μ_{a 2}) ϕ_{2} = \\ - \frac{2}{3} Q + \frac{2}{3} μ_{a} ϕ_{1} + (\frac{16}{45} μ_{a} + \frac{4}{9} μ_{a 2}) ϕ_{3} \\ - (\frac{32}{105} μ_{a} + \frac{8}{21} μ_{a 2}) ϕ_{4} \end{array}

\begin{array}{l} - \nabla . \frac{1}{11 μ_{a 5}} \nabla ϕ_{3} + (\frac{64}{225} μ_{a} + \frac{16}{45} μ_{a 2} + \frac{9}{25} μ_{a 4}) ϕ_{3} = \\ \frac{8}{15} Q - \frac{8}{15} μ_{a} ϕ_{1} + (\frac{16}{45} μ_{a} + \frac{4}{9} μ_{a 2}) ϕ_{2} \\ + (\frac{128}{525} μ_{a} + \frac{32}{105} μ_{a 2} + \frac{54}{175} μ_{a 4}) ϕ_{4} \end{array}

\begin{array}{l} - \nabla . \frac{1}{15 μ_{a 7}} \nabla ϕ_{4} + (\frac{256}{1225} μ_{a} + \frac{64}{245} μ_{a 2} + \frac{324}{1225} μ_{a 4} + \frac{13}{49} μ_{a 6}) ϕ_{4} = \\ - \frac{16}{35} Q + \frac{16}{35} μ_{a} ϕ_{1} - (\frac{32}{105} μ_{a} + \frac{8}{21} μ_{a 2}) ϕ_{2} \\ + (\frac{128}{525} μ_{a} + \frac{32}{105} μ_{a 2} + \frac{54}{175} μ_{a 4}) ϕ_{3} \end{array}

where Q is the isotropic source term,

ϕ_{n}

are composite moments of the fluence and μ_an is the n^th order absorption coefficient given by:

μ_{a n} (x) = μ_{t} (x) - μ_{s} (x) g^{n}

where the total attenuation coefficient μ_t is given by the sum of the absorption and scattering coefficients (μ_t = μ_a + μ_s). In practice, however, it is only possible to experimentally record the total fluence which is calculated from the composite moments as

φ = ϕ_{1} - \frac{2}{3} ϕ_{2} + \frac{8}{15} ϕ_{3} - \frac{16}{35} ϕ_{4}

If we assume that the composite moments are in fact measureable, the SP_N approximation introduces the possibility of reconstructing a wider range of optical parameters. Whereas the DA allows the recovery of just μ_a and μ_s’, the SP_N approximation allows the recovery of μ_a, μ_a1, μ_a2, μ_a3, … μ_an and, from these values μ_s and g can be calculated through Eq. (2). It is important to note that although it is possible to measure the angular dependence of the fluence at the boundary through the use of optical fibers whose numerical aperture can be adjusted, it is not clear how the composite moments as shown in Eq. (1) and 3 may be measured experimentally. Nonetheless, for the purpose of the theoretical study presented here, it will be assumed that these moments can be measured. The SP_N equations have been implemented using the NIRFAST package and previously validated using Monte Carlo models [36].

2.2 Image reconstruction

The goal of the inverse problem is the recovery of the optical properties, μ, by minimizing the difference between the measured data Φ^M and data calculated by the forward model Φ^C using a modified Tinkhonov minimization approach given by [22]

χ^{2} = {}_{μ}^{\min}{\sum_{i = 1}^{N M} {(Φ^{M} - Φ^{C})}^{2} + λ \sum_{j = 1}^{N N} {(μ_{j} - μ_{0})}^{2}}

where NM is the number of measurements, NN is the number of nodes in the FEM, μ₀ is an initial estimate of the optical properties and λ is the Tikhonov regularization parameter. By considering only the first order terms and applying the Levenberg-Marquardt procedure, this leads to an update vector

δ μ = {(J^{T} J + \bar{λ} I)}^{- 1} J^{T} δ Φ

where

\bar{λ} = 2 λ

and J is the Jacobian matrix

(δ Φ / δ μ)

, often referred to as the weight matrix or sensitivity matrix. The Jacobian defines the relationship between changes in calculated boundary data due to small changes in the optical properties. If both phase and amplitude data is available, the structure of the Jacobian becomes

J = [\begin{matrix} \frac{δ \ln I_{1}}{δ μ_{1}} & \frac{δ \ln I_{1}}{δ μ_{2}} & \dots & \frac{δ \ln I_{1}}{δ μ_{N N}} \\ \frac{δ θ_{1}}{δ μ_{1}} & \frac{δ θ_{1}}{δ μ_{2}} & \dots & \frac{δ θ_{1}}{δ μ_{N N}} \\ \frac{δ \ln I_{2}}{δ μ_{1}} & \frac{δ \ln I_{2}}{δ μ_{2}} & \dots & \frac{δ \ln I_{2}}{δ μ_{N N}} \\ \frac{δ θ_{2}}{δ μ_{1}} & \frac{δ θ_{2}}{δ μ_{2}} & \dots & \frac{δ θ_{2}}{δ μ_{N N}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{δ \ln I_{N M}}{δ μ_{1}} & \frac{δ \ln I_{N M}}{δ μ_{2}} & \dots & \frac{δ \ln I_{N M}}{δ μ_{N N}} \\ \frac{δ θ_{N M}}{δ μ_{1}} & \frac{δ θ_{N M}}{δ μ_{2}} & \dots & \frac{δ θ_{N M}}{δ μ_{N N}} \end{matrix}]

where δlnI_i/δμ_j are sub-matrices that define changes in the ith measurement of log amplitude due to a change in optical properties at node j, and δθ_i/δμ_j define changes in the ith measurement of phase due to changes in optical properties at node j.

In this study, the Jacobian as presented in sections 3.1 – 3.3 is constructed using the perturbation method and for section 3.4 these are calculated using the Reciprocity approach which was validated using the perturbation method [39]. As the optical properties of the FEM mesh were defined by absorption coefficient, scattering coefficient and anisotropy factor, the final Jacobian was made up of three separate kernels with the form

J = [\begin{matrix} J_{μ_{a}}; & J_{μ_{s}}; & J_{g} \end{matrix}]

where

J_{μ_{s}}

,

J_{μ_{a}}

and

J_{g}

are the Jacobians due to changes in scattering coefficient, absorption coefficient and anisotropy respectively and have the form

J_{μ_{s}} = [\begin{matrix} \frac{δ \log I}{δ μ_{s}} & ; & \frac{δ θ}{δ μ_{s}} \end{matrix}]

J_{μ_{a}} = [\begin{matrix} \frac{δ \log I}{δ μ_{a}} & ; & \frac{δ θ}{δ μ_{a}} \end{matrix}]

J_{g_{}} = [\begin{matrix} \frac{δ \log I}{δ g} & ; & \frac{δ θ}{δ g} \end{matrix}]

3. Methods and results

3.1 Sensitivity mapping

The sensitivity of the SP₇ models to changes in optical properties were studied using a circular geometry of 43 mm radius with 16 equally spaced and collocated sources and detectors as shown in Fig. 1 . The Jacobian was constructed using the perturbation method before being mapped onto the mesh coordinates to produce a map of sensitivity to changes in optical properties. The reference data was calculated using homogeneous optical properties of μ_a = 0.01 mm⁻¹, μ_s = 10 mm⁻¹, g = 0.9 and refractive index (n) = 1.33 using a modulation frequency of 100MHz. The sources were modeled as point sources located at 1 mm inside the boundary, to correspond with one reduced scattering distance as in DA.

Fig. 1 FEM mesh, of 43mm radius, containing 1785 nodes and 3418 linear triangular elements. Circles represent location of sources, crosses represent location of detectors.

Download Full Size | PDF

The total sensitivity of log amplitude for all sources and detectors pairs to changes in absorption coefficient is mapped in Fig. 2(a) . Areas of high sensitivity exist in the regions surrounding the sources and detectors with a much lower background sensitivity decreasing, as expected, towards the centre of the domain. The negative magnitudes of sensitivity confirm that increases in absorption lead to reductions in measurements of intensity. Figure 2(b) shows the sensitivity of same data-type to scattering coefficient and of similar trend to the previous case. The sensitivity scale for scattering is again negative but has a much lower magnitude than that of absorption. Figure 2(c) shows that the sensitivity of intensity measurements to changes in anisotropy factor and is again highly localized to the sources and detectors. The magnitude of sensitivity of g is several orders lower than that of both absorption and scattering.

Fig. 2 Maps of sensitivity of log Amplitude data to changes in a) absorption, b) scattering and c) anisotropy

Download Full Size | PDF

Figures 2(a)-(b) show that increases in both the absorption and scattering coefficients lead to decreases in intensity measurements. This suggests that it is not possible to distinguish between the two effects. A decrease in intensity, for example, could be caused by either a small increase in absorption or a large increase in scattering.

Figure 3(a) shows the sensitivity of phase measurements to changes in absorption coefficient. It is evident that the sensitivity is greater in the regions surrounding the sources and detectors as compared to the regions in between them. Additionally, the sensitivity although increasing to a maximum at a few scattering depths inside, it then decreases towards the centre of the domain, which is consistent with previous findings [39]. The scattering sensitivity map, Fig. 3(b), shows regions of very high sensitivity near the sources and detectors with much lower sensitivity elsewhere. The magnitude of scattering sensitivity is again much lower than that of absorption but with a positive scale. Sensitivity of phase measurements to changes in anisotropy factor is again limited to the regions surrounding the sources and detectors with a much lower magnitude as compared to absorption and scatter.

Fig. 3 Maps of sensitivity of phase data to changes in a) absorption, b) scattering and c) anisotropy

Download Full Size | PDF

The results show that using log amplitude data alone, it is not possible to separate changes due to perturbations in absorption and scatter properties since both result in the same effect (i.e. increase in either absorption or scatter results in a decrease in log amplitude measurements). With the introduction of phase data, however, it can be seen that changes in amplitude and scatter have opposite effects on such data and so it is possible to distinguish between perturbations in absorption or scatter. However, since scattering and anisotropy factor have opposite sensitivities for both data-types, an increase in one parameter can be compensated by a decrease in the other and vice-versa. Therefore the problem appears non-unique for reconstruction of scatter and anisotropy factor.

3.2 Uniqueness of boundary data

Next the uniqueness of the imaging problem is analyzed. The problem of non-uniqueness arises when more than one set of optical property distributions lead to an identical set of boundary data, in which case, the simultaneous recovery of multiple optical properties, such as the absorption and scattering coefficients, is not possible [40].

The previously used circular mesh, Fig. 1, was used to test the uniqueness of the SP₇ equations. A circular anomaly, with radius 10 mm, was inserted into the mesh as shown in Fig. 4 . The optical properties within the anomaly were perturbed with either the absorption coefficient ranging between 0.0125mm⁻¹ to 0.0375mm⁻¹, scattering coefficient ranging between 10mm⁻¹ to 30mm⁻¹ and anisotropy factor ranged between 0 and 1 and corresponding boundary data for all possible combinations were calculated. These wide range of values are chosen in order to allow a comprehensive study of their effect on problem uniqueness. Reference data was also calculated using the optical properties of μ_a = 0.025 mm⁻¹, μ_s = 20 mm⁻¹ and g = 0.9. The absolute error (L1 norm) for each data-type between the reference and each set of perturbed data was the defined as

δ = \sum | (φ_{r e f} - φ_{μ}) |

where

φ_{r e f}

is the reference homogeneous unperturbed boundary data (either phase or log amplitude) and

φ_{μ}

is the perturbed data. Figure 5(a) shows the map of error (Eq. (9) in amplitude data for each combination of μ_a and μ_s. A large band of optical property combinations exists for which the error between the reference data and perturbed data falls to zero. These equivalent combinations lead to identical amplitude data at the boundary and, as such, amplitude data alone is insufficient to recover both absorption and scattering properties simultaneously. By introducing phase measurements, however, the number of non-unique solutions can be minimized. Figure 5(b) shows the error map of phase measurements for the same range of optical properties. In this case, the optical property combinations leading to identical boundary data lie within a much more localized region. By combining the two data types, however, Fig. 5(c), the solution can converge on a smaller range of optical properties that lead to the perturbed boundary data.

Fig. 4 Circular mesh with a single inclusion of 10mm radius used to test uniqueness of boundary measurements.

Download Full Size | PDF

Fig. 5 Error maps of SP₇ data with varying optical properties with arbitrary units of error. (a) map of log(Amplitude) with varying µ_a (y axis) and µ_s (x axis), (b) same as (a) but for phase and (c) error map of combined log(Amplitude) and Phase data.

Download Full Size | PDF

The study was then repeated with varying absorption coefficients and anisotropy factors. Figure 6(a) shows that a large band of non-unique optical property combinations exists for amplitude measurements. The errors in phase measurements, Fig. 6(b), show an even larger range of equivalent combinations. When combining the two data types, however, it can be seen that the number of $μ_{a}$ and g combinations leading to the reference data actually falls within a smaller narrow region, Fig. 6(c).

Fig. 6 Same as Fig. 5, but for (a) map of log(Amplitude) data with varying g and µ_a, (b) same as (a) but for phase, (c) sum of (a) and (b)

Download Full Size | PDF

The map of error in amplitude measurements with varying μ_s and g, are also shown in Fig. 7 . Figure 7(a), shows a large band of equivalent combinations of the two parameters for log amplitude data. The error in phase data actually shows two large regions of non-unique optical property combinations. The band of equivalent combinations in the amplitude error map is parallel to one of the bands in the phase error map. As such, even when combining the two data types, a large band of equivalent combinations still exists, Fig. 7(c), and therefore, even with both phase and amplitude data, it is not possible to distinguish between changes due to $μ_{s}$ and g.

Fig. 7 Same as Fig. 6, but for (a) map of log(Amplitude) data with varying g and µ_s, (b) same as (a) but for phase, (c) sum of (a) and (b)

Download Full Size | PDF

By direct visual comparison between Figs. 6(c) and 7(c), it is evident that there exists a large range of absorption, scatter and anisotropy factor that would give rise to the same boundary measurements, even when considering both data-types and hence indicating the non-uniqueness of the problem to reconstruct all three parameters.

3.3 Uniqueness with the introduction of higher order terms

The use of higher ordered SP_N approximations introduces a number of new variables, μ_an, which are defined as in Eq. (2). These higher ordered absorption coefficients are dependent on both μ_s and g and as such, may make it possible to separate the effects due to changes in the two parameters, if it is assumed that the composite moments of the boundary data can be measured.

To test the uniqueness of the higher absorption orders, the values of μ_an for n = 1:4 (as available for SP₅) were calculated using optical properties of μ_a = 0.01 mm⁻¹, μ_s = 20 mm⁻¹ and g = 0.5 as a reference. The values of μ_an were then re-calculated with values of the scattering coefficient ranging from 1mm⁻¹ to 40mm⁻¹ and the anisotropy factor ranging from 0.1 to 0.9.

Figure 8(a) shows the difference between the reference value of μ_a1 and those calculated for each combination of scattering coefficient and anisotropy factor. There is a large band of optical property combinations for which the residual falls to zero and therefore result in the same value of μ_a1. Figure 8(b) shows the same differences between values of μ_a2. There is still a large number of scattering coefficient / anisotropy factor combinations that lead to the same value in this case, although slightly smaller than for μ_a1. Figures 8(c)-(d) show similar plots for μ_a3 and μ_a4 respectively. It can be seen that as higher order moments of absorption are considered, the number of scattering coefficient / anisotropy factor combinations leading to the same value reduces, but does not converge to a particular combination. As such, the scattering coefficient and anisotropy factors cannot be simultaneously identified.

Fig. 8 Maps of residuals between (a) μ_a1, (b) μ_a2, (c) μ_a3 and (d) μ_a4 for a range of scattering coefficients and anisotropy factors as compared to a set of reference values. The range of anisotropy factors is listed on the x-axes whilst the y-axes represent the range of scattering coefficients.

Download Full Size | PDF

3.4 Multi-parameter image reconstruction using the SP_N approximation

SP_N based image reconstruction algorithms were developed for N = 1, 3, 5 and 7 using the modified Tikhonov minimization method. Test boundary data was generated using the same circular geometry, Fig. 1, with μ_a = 0.01 mm⁻¹, μ_s = 10 mm⁻¹, g = 0.9 and n = 1.33. A highly absorbing anomaly, with μ_a = 0.02 mm⁻¹ and a highly scattering anomaly with μ_s = 20 mm⁻¹, were inserted into the mesh as shown in Fig. 9 . Both anomalies have the same anisotropic factor and refractive index as the background.

Fig. 9 FEM model with one highly absorbing target and one highly scattering target.

Download Full Size | PDF

Each of the SP_N reconstruction algorithms (for N = 1, 3, 5 and 7) was capable of reconstructing different optical properties, depending on the order N. In order to test the capability of each of the SP_N models, the optical properties listed in Table 1 were reconstructed, and from these, values of μ_a and μ_s were extracted. For the inverse model, the same mesh as the forward model was used without added noise to ensure that any errors in the reconstructed optical maps were only due to the complexity of the higher ordered approximations. Each of the algorithms were allowed to continue until there was less than a 0.1% change in successive iterations. The initial regularization parameter λ (Eq. (5) for all algorithms was set to 10 times the maximum diagonal of the matrix J^TJ and was reduced by a factor of 10^0.25 at each iteration [22].

Table 1. Reconstructed values of SP_N image reconstruction algorithms where the diffusion terms κ_n = 1/(Aμ_an) where A is a constant.

View Table

Figure 10(a) shows the reconstructed images generated using an SP₁ based reconstruction algorithm and SP₁ forward data. The location of the two anomalies has been accurately reconstructed. The background optical properties have been accurately recovered although the values of the two anomalies have been over-estimated. The reconstructed image shows a maximum absorption coefficient of 0.0249 mm⁻¹ and a maximum scattering coefficient of 23.9 mm⁻¹. The SP₁ reconstruction required 32 iterations at 3.1 second per iteration to recover the two relevant unknowns.

Fig. 10 Reconstructed optical maps using a) SP₁, b) SP₃, c) SP₅, d) SP₇ based reconstruction algorithms

Download Full Size | PDF

Reconstruction using the SP₅ model (using SP₅ data) is shown in Fig. 10(c). In this case, the optical properties have again been overestimated with a maximum absorption coefficient of 0.0224mm⁻¹ and a maximum scattering coefficient of 21.4mm⁻¹. The cross talk between the absorption and scattering images is still present. The SP₅ reconstruction required 34 iterations at 54 second per iteration to recover the six relevant unknowns.

The optical maps generated by the SP₇ model (using SP₇ data), Fig. 10(d) contains artifact throughout the domain. The location of the highly absorbing target has been recovered with poor size and contrast accuracy. The reconstruction failed to recover the highly scattering target. It is likely that the failure of the SP₇ reconstruction is due to the highly under-determined nature of the problem. The SP₇ model requires 8 unique variables to be calculated at each node, i.e. 14280 unknowns, based on just 1920 boundary measurements. The SP₇ reconstruction required 20 iterations at 52 second per iteration to recover the eight relevant unknowns.

3.5 Multi-parameter image reconstruction and hard priori information

The under-determined nature of the image reconstruction problems can be minimized with the use of prior information. MRI data can be used, for example, to determine regions of similar tissue types and therefore, regions of homogeneous optical properties [41]. This then reduces the problem of calculating optical properties at every node to just a small number of known regions.

The previous study was repeated using region-based reconstruction for the SP_N methods with N = 1, 3 and 5. The SP₇ model has been omitted from further studies as the increase in accuracy over the SP₅ model has been shown to be very small and neglectable [35].

It can be seen (Fig. 11 ) that the use of prior information, whereby instead of reconstruction the unknown parameters at each spatial variable, homogeneous values are estimated for each unknown region (in this case 3) has enabled the SP₅ reconstruction (when using SP₅ forward data) to accurately recover both targets. This is due to the vast reduction in the number of unknowns from 6 times the number of nodes (10710) to 6 times the number of regions (18). The same findings (not shown) were found for SP₁ and SP₃ cases.

Fig. 11 Recovered optical map generated using SP₅ reconstruction with prior information.

Download Full Size | PDF

3.6 Diffusion based image reconstruction using the SP_N approximation

As discussed earlier, the SP_N approximations, Eqs. (1)(a)-(d), are based on composite moments of fluence with the total flux calculated using Eq. (2). The SP_N reconstruction algorithms introduced so far were therefore based on various moments of the phase and amplitude boundary data which attempted to recover all moments of absorption and scattering coefficients. In practice, however, there are no experimental methods to easily measure the angular components of amplitude and phase and so image reconstructions must be based on total values of fluence. A new reconstruction algorithm was therefore developed that used the SP₅ model but only required measurements of phase and amplitude as available from the total fluence and not composite moments.

The new reconstruction algorithm uses the SP₅ forward model to calculate the composite moments of fluence. The total fluence is then calculated using Eq. (3). The Jacobian is then constructed using Eqs. (8a) and (8b) where log of amplitude and absolute phase are extracted from the calculated total fluence. By eliminating the composite moments of fluence, however, it is only possible to reconstruct for μ_a. and μ_s (assuming g = 0.9). This new reconstruction algorithm is tested using SP₅ forward data and the resulting optical map is shown in Fig. 12(a) . The reconstruction has performed well recovering the optical parameters within 10% of the target values. For comparison, the same forward data is also used to reconstruct optical parameters using the DA based algorithm, Fig. 12(b). In this case, the target values have been over-estimated and boundary artifact has been introduced, indicating the maximum errors seen from high order model mismatch, is as expected near the source / detector positions, which lead to image inaccuracy and artifacts.

Fig. 12 Diffusion based images from SP₅ data using (a) SP₅ based Jacobian and (b) SP₁ based Jacobian.

Download Full Size | PDF

4. Discussions and conclusions

Forward models for image reconstruction based on the simplified spherical harmonic approximation have been presented. Sensitivity maps for amplitude and phase data are calculated, Figs. 2 and 3, using SP₇ and it is shown that the sensitivity of measured boundary data for changes in the anisotropy factor are many orders of magnitude lower than those of absorption and scattering coefficients. It is also shown that the sensitivity to both scattering coefficient anisotropy factor is highly localized to the regions adjacent to the source and detector (Figs. 2 and 3). This suggests that beyond a few scattering lengths, the light distribution is fully diffuse and as such, distant changes in scattering and anisotropy factors have little effect on measured data. Although it is possible to normalize the sensitivity maps prior to use in image reconstruction [42], it is also evident that an increase in scattering parameter will have an effect on boundary data which is opposite of that due to anisotropy factor. This in turn indicates that an increase in one parameter can be compensated by a decrease on the other and vice-versa. This effect is unlike that seen for absorption and scattering, whereby the two parameters can be separated using log amplitude and phase data, but it is unclear whether all 3 parameters can be resolved using these 2 data-types.

Through the use of error maps, whereby the difference in some reference data with respect to data in the presence of an anomaly is calculated, it is shown that using frequency domain data, it is possible to distinguish, with some accuracy, the absorption and scattering properties of the medium, Fig. 5, which has also been presented previously [21]. The same principles have been applied to higher order models and additionally the effect of the anisotropy factor has been studied as well as the optical absorption and scatter. It is demonstrated, Fig. 6 and 7, that using frequency domain data, a large and non-unique range of all three optical properties can lead to the same measured data, indicating that it is not possible to reconstruct all three parameters using only 2 data-types. This finding is in line with earlier results, Fig. 3, and demonstrate that it is difficult to separate for both scattering and anisotropic factors simultaneously using frequency domain data. Further work is needed to investigate whether other data-types, for example using time resolved data, and / or spectral data, would allow separation of these variables.

Residual maps for different orders of absorption coefficient (μ_an) as available using SP_N models have been calculated to investigate the possibility of resolving scattering and anisotropy factors, if μ_an coefficients can be measured and calculated. It is shown, Fig. 8, that even if all higher moments of the absorption coefficients can be calculated, there still exists a large range of scattering and anisotropy factors that would lead to the same μ_an values, thus indicating that these cannot be easily resolved.

Reconstruction algorithms based on the SP_N approximation have been developed and tested. For image reconstruction, the same FEM model as for the forward data was used and no noise was added. The inverse crime has been strictly committed, since we are concerned with the accuracy of each model in determining the unknown parameters associated with each model. Additionally, noise free data have only used to only highlight the error seen due to model mismatch, rather than the effects of noise within the data, which can show similar trends. It was shown, Fig. 10(a)-(c), that for N = 1, 3 and 5 the reconstructions performed well, with the reconstructed values being within 24% of expected values with the worse results being obtained from SP₁ and absorption coefficient. This accuracy could be further improved by the optimization of the regularization parameter and stopping criteria. The N = 7 model, however, failed to accurately recover the target values, Fig. 10(d). The absorption coefficient, although located, is extremely over estimated, with some cross-talk from the scattering value. The reconstructed scatter image contains high valued artifacts and can be considered inaccurate. By limiting the presented results to these larger domain, it is also worth noting that the majority of image artifacts are seen at the boundary, near the source locations, where SP₁ solution is known the be less accurate. It is therefore expected that the errors seen due to the lack of higher order approximations to be substantially more significant in small geometry imaging experiments. The poor performance of the SP₇ model was most likely due to the under-determined nature of the problem. As stated earlier, the SP₇ model contains 8 unique unknown variables which need to be calculated at each spatial location (FEM node), i.e. 14280 unknowns, based on just 1920 boundary measurements (240 log amplitude and 240 Phase, based on 16 co-located source and detectors and 4 composite moments of fluence). In order to eliminate the large degree of freedoms, the accuracy of the reconstruction algorithms was further improved by creating a better determined problem. This was achieved for the SP₁, SP₃ and SP₅ approximations with the use of prior information. The prior information was used to identify a number of homogeneous regions, reducing the number of properties to be recovered. This simplification of the problem led to improvements in all three of the reconstructions (for N = 1, 3 and 5). Both the absorbing and scattering targets were accurately recovered by all models to within 2%, Fig. 11.

The SP_N equations are based on composite moments of fluence, Eq. (2). In reality, however, experimental systems can only measure the total fluence at the boundary of the domain and so the full SP_N reconstruction algorithms are of limited use in image reconstruction where all composite moments are needed. An alternative method has been proposed in which the forward models are based on absolute fluence but the Jacobian matrix was calculated using the SP₅ model to allow the calculation of absorption (μ_a) and scattering (μ_s) coefficient only. The diffusion parameter based images reconstructed from simulated SP₅ data whereby the Jacobian is based on either SP₁ or SP₅, Fig. 12, indicate that although both models can be used, the higher order model is more accurate both in terms of quantitative and qualitative analysis. The target values of the test problem calculated using SP₅ are recovered with more accuracy as compared to the SP₁ based model with much less artifact. Additionally, eliminating the use of complex moments also results in decreased computation time and memory requirements.

The results presented in this work indicate that for image reconstruction whereby the DA is less valid, in for example, small animal imaging and / or where the absorption coefficient is more dominant, the higher order models based on simplified spherical harmonics can be used to generate the sensitivity matrix for diffusion based image reconstruction, without the additional computation complexity in terms of the number of unknown parameters. The incorporation of these more accurate models can however allow for a better accuracy in terms of light propagation models.

Acknowledgements

This work has been sponsored by the Engineering and Physical Sciences Research Council, UK.

References and links

1. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50(4), R1–R43 (2005). [CrossRef] [PubMed]

2. D. R. Leff, O. J. Warren, L. C. Enfield, A. Gibson, T. Athanasiou, D. K. Patten, J. Hebden, G. Z. Yang, and A. Darzi, “Diffuse optical imaging of the healthy and diseased breast: a systematic review,” Breast Cancer Res. Treat. 108(1), 9–22 (2008). [CrossRef] [PubMed]

3. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18(6), 57–75 (2001). [CrossRef]

4. R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32(4), 1128–1139 (2005). [CrossRef] [PubMed]

5. H. Dehghani, B. W. Pogue, S. P. Poplack, and K. D. Paulsen, “Multiwavelength three-dimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results,” Appl. Opt. 42(1), 135–145 (2003). [CrossRef] [PubMed]

6. H. Jiang, Y. Xu, N. Iftimia, J. Eggert, K. Klove, L. Baron, and L. Fajardo, “Three-dimensional optical tomographic imaging of breast in a human subject,” IEEE Trans. Med. Img. 20(12), 1334–1340 (2001). [CrossRef]

7. L. C. Enfield, A. P. Gibson, N. L. Everdell, D. T. Delpy, M. Schweiger, S. R. Arridge, C. Richardson, M. Keshtgar, M. Douek, and J. C. Hebden, “Three-dimensional time-resolved optical mammography of the uncompressed breast,” Appl. Opt. 46(17), 3628–3638 (2007). [CrossRef] [PubMed]

8. B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar, and J. P. Culver, “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Sci. U.S.A. 104(29), 12169–12174 (2007). [CrossRef] [PubMed]

9. D. A. Boas, K. Chen, D. Grebert, and M. A. Franceschini, “Improving the diffuse optical imaging spatial resolution of the cerebral hemodynamic response to brain activation in humans,” Opt. Lett. 29(13), 1506–1508 (2004). [CrossRef] [PubMed]

10. J. P. Culver, A. M. Siegel, J. J. Stott, and D. A. Boas, “Volumetric diffuse optical tomography of brain activity,” Opt. Lett. 28(21), 2061–2063 (2003). [CrossRef] [PubMed]

11. J. C. Hebden, A. Gibson, T. Austin, R. M. Yusof, N. Everdell, D. T. Delpy, S. R. Arridge, J. H. Meek, and J. S. Wyatt, “Imaging changes in blood volume and oxygenation in the newborn infant brain using three-dimensional optical tomography,” Phys. Med. Biol. 49(7), 1117–1130 (2004). [CrossRef] [PubMed]

12. A. Y. Bluestone, G. Abdoulaev, C. Schmitz, R. L. Barbour, and A. H. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express 9(6), 272–286 (2001). [CrossRef] [PubMed]

13. A. D. Klose and A. H. Hielscher, “Fluorescence tomography with simulated data based on the equation of radiative transfer,” Opt. Lett. 28(12), 1019–1021 (2003). [CrossRef] [PubMed]

14. V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8(1), 1–33 (2006). [CrossRef] [PubMed]

15. E. M. Sevick-Muraca, G. Lopez, J. S. Reynolds, T. L. Troy, and C. L. Hutchinson, “Fluorescence and absorption contrast mechanisms for biomedical optical imaging using frequency-domain techniques,” Photochem. Photobiol. 66(1), 55–64 (1997). [CrossRef] [PubMed]

16. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(4225), 4241 (2005). [CrossRef]

17. C. H. Contag and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002). [CrossRef] [PubMed]

18. C. Kuo, O. Coquoz, T. L. Troy, H. Xu, and B. W. Rice, “Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging,” JBO 12, 24007 (2007).

19. G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. McLennan, and M. Henry, “In vivo mouse studies with bioluminescence tomography,” Opt. Express 14(17), 7801–7809 (2006). [CrossRef] [PubMed]

20. H. Dehghani, S. C. Davis, and B. W. Pogue, “Spectrally resolved bioluminescence tomography using the reciprocity approach,” Med. Phys. 35(11), 4863–4871 (2008). [CrossRef] [PubMed]

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

22. 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: Algorithms for Numerical Model and Image Reconstruction Algorithms,” Communications in Numerical Methods in Engineering (2008).

23. A. H. Hielscher, A. D. Klose, and K. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imaging 18(3), 262–271 (1999). [CrossRef] [PubMed]

24. J. Riley, H. Dehghani, M. Schweiger, S. Arridge, J. Ripoll, and M. Nieto-Vesperinas, “3D optical tomography in the presence of void regions,” Opt. Express 7(13), 462–467 (2000). [CrossRef] [PubMed]

25. S. R. Arridge, H. Dehghani, M. Schweiger, and E. Okada, “The finite element model for the propagation of light in scattering media: a direct method for domains with nonscattering regions,” Med. Phys. 27(1), 252–264 (2000). [CrossRef] [PubMed]

26. H. Dehghani, S. R. Arridge, M. Schweiger, and D. T. Delpy, “Optical tomography in the presence of void regions,” J. Opt. Soc. Am. 17(9), 1659–1670 (2000). [CrossRef]

27. A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. 43(5), 1285–1302 (1998). [CrossRef] [PubMed]

28. J. Ripoll, M. Nieto-Vesperinas, S. R. Arridge, and H. Dehghani, “Boundary conditions for light propagation in diffusive media with nonscattering regions,” J. Opt. Soc. Am. A 17(9), 1671–1681 (2000). [CrossRef]

29. K. M. Case, and P. F. Zweifel, Linear Transport Theory (Addidon-Wesley, Reading, MA, 1967).

30. H. B. Jiang, “Optical image reconstruction based on the third-order diffusion equations,” Opt. Express 4(8), 241–246 (1999). [CrossRef] [PubMed]

31. S. Wright, M. Schweiger, and S. R. Arridge, “Reconstruction in optical tomography using the P-N approximations,” Meas. Sci. Technol. 18(1), 79–86 (2007). [CrossRef]

32. E. D. Aydin, C. R. E. de Oliveira, and A. J. H. Goddard, “A comparison between transport and diffusion calculations using a finite element-spherical harmonics radiation transport method,” Med. Phys. 29(9), 2013–2023 (2002). [CrossRef] [PubMed]

33. A. D. Klose, U. Netz, J. Beuthan, and A. H. Hielscher, “Optical tomography using the time-independent equation of radiative transfer — Part 1: forward model,” J. Quant. Spectrosc. Radiat. Transf. 72(5), 691–713 (2002). [CrossRef]

34. S. Cahandrasekhar, Radiative Transfer (Clarendon Press, London, 1950).

35. A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys. 220(1), 441–470 (2006). [CrossRef]

36. M. Chu, K. Vishwanath, A. D. Klose, and H. Dehghani, “Light transport in biological tissue using three-dimensional frequency-domain simplified spherical harmonics equations,” Phys. Med. Biol. 54(8), 2493–2509 (2009). [CrossRef] [PubMed]

37. G. Marquez, L. V. Wang, S. P. Lin, J. A. Schwartz, and S. L. Thomsen, “Anisotropy in the absorption and scattering spectra of chicken breast tissue,” Appl. Opt. 37(4), 798–804 (1998). [CrossRef] [PubMed]

38. S. Nickell, M. Hermann, M. Essenpreis, T. J. Farrell, U. Krämer, and M. S. Patterson, “Anisotropy of light propagation in human skin,” Phys. Med. Biol. 45(10), 2873–2886 (2000). [CrossRef] [PubMed]

39. S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part2: Finite-element-method calculations,” Appl. Opt. 34(34), 8026–8037 (1995). [CrossRef] [PubMed]

40. S. R. Arridge and W. R. B. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. 23(11), 882–884 (1998). [CrossRef] [PubMed]

41. H. Dehghani, B. W. Pogue, J. Shudong, B. Brooksby, and K. D. Paulsen, “Three-dimensional optical tomography: resolution in small-object imaging,” Appl. Opt. 42(16), 3117–3128 (2003). [CrossRef] [PubMed]

42. M. E. Eames and H. Dehghani, “Wavelength dependence of sensitivity in spectral diffuse optical imaging: effect of normalization on image reconstruction,” Opt. Express 16(22), 17780–17791 (2008). [CrossRef] [PubMed]

Image reconstruction in diffuse optical tomography based on simplified spherical harmonics approximation

Abstract

1. Introduction