Comparison of diffusion approximation and higher order diffusion equations for optical tomography of osteoarthritis

Zhen Yuan; Qizhi Zhang; Eric S. Sobel; Huabei Jiang

doi:10.1117/1.3233655

1 September 2009 Comparison of diffusion approximation and higher order diffusion equations for optical tomography of osteoarthritis

Zhen Yuan, Qizhi Zhang, Eric S. Sobel, Huabei Jiang

Author Affiliations +

Journal of Biomedical Optics, Vol. 14, Issue 5, 054013 (September 2009). https://doi.org/10.1117/1.3233655

Abstract

In this study, a simplified spherical harmonics approximated higher order diffusion model is employed for 3-D diffuse optical tomography of osteoarthritis in the finger joints. We find that the use of a higher-order diffusion model in a stand-alone framework provides significant improvement in reconstruction accuracy over the diffusion approximation model. However, we also find that this is not the case in the image-guided setting when spatial prior knowledge from x-rays is incorporated. The results show that the reconstruction error between these two models is about 15 and 4%, respectively, for stand-alone and image-guided frameworks.

1. Introduction

3-D diffuse optical tomography (DOT), as an emerging biomedical imaging modality, has made considerable advances in recent years. ^{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27} So far DOT has been primarily applied to image brain and breast tissues.^{2, 3, 4, 5, 6} Recent phantom and clinical studies show that DOT can also provide quantitative optical images of hand joints for early detection of joint-related diseases.^{7,8,9,10,11,12,13} In DOT, an efficient reconstruction algorithm is critical. Most of the DOT algorithms developed to date make use of so-called model-based image reconstruction schemes. Due to computational complexity, the models used for describing photon migration in biological tissues and image reconstruction have usually been limited to the diffusion approximation (DA) to the Boltzmann radiative transport equation (RTE). However, the DA is not accurate enough to model light transport in some particular tissue regions, including low scattering fluid-filled and high absorption vascular tissues.^{14, 15, 16, 17, 18} In particular, using the DA is very challenging to image small volumes including arthritis in finger joints and body parts in small animals. In such cases, the DA is not able to describe photon migration accurately due to the small optical distance between sources and detectors. To overcome these limitations associated with the DA, several groups have employed the discrete ordinate $(s_{n})$ approximated RTE to recover the optical parameters of small volume tissues.^{10, 11, 12} However, this type of reconstruction algorithm is time consuming. Typically, the $s_{n}$ approximated-RTE-based reconstruction is 40 to 60 times slower than DA based.¹¹

We have developed a 3-D reconstruction method based on simplified spherical harmonics $(p_{3})$ approximated-RTE or higher order diffusion equations.²¹ Our simulations and phantom tests²¹ have demonstrated that the developed higher order diffusion (HD) equations are able to provide solutions that are significantly more accurate than the DA for simulating photon migration and imaging reconstruction in finger joints. We also found that the HD-model-based reconstruction algorithm is only 1 to 2 times slower than the DA based reconstruction.²¹

In this study, we focus on 3-D DOT imaging of finger joints for the diagnosis of osteoarthritis (OA) using the developed HD model. One goal of this study is to investigate if the HD model has the potential to improve the reconstruction accuracy for in-vivo data as compared to that using the DA.⁸ In addition, we evaluate the performance of the HD model using an x-ray/DOT combined multimodality approach.²⁸ The idea of this multimodality approach is to incorporate the high resolution x-ray images of joint tissues into the DOT reconstruction so that both the resolution and accuracy of optical image reconstruction are enhanced. We found from our previous investigations that the improved quantitative capability of imaging absorption and scattering coefficients of the joint tissues allows for more accurate diagnosis of osteoarthritis over the DA-based DOT alone.²⁸ When anatomical a-priori information becomes available, we would like to study if a more accurate modeling of photon migration in tissue can bring in enhanced accuracy for image reconstruction in the finger joints.

2. Reconstruction Methods

The migration of photons within the absorbing and scattering media can be described by the time- and energy-independent form of the Boltzmana transport equation as,^{14, 15, 16, 17}

Eq. 1

Ω \cdot \nabla Φ (r, Ω) + σ_{t} (r) Φ (r, Ω) = \int σ_{s} (Ω^{'}, Ω, r) Φ (r, Ω^{'}) d Ω^{'} + s (r, Ω),

in which

r

is the position vector of a photon propagation along the unit direction Ω,

Φ (r, Ω)

is the energy radiance,

s (r, Ω)

is the source term,

σ_{t} (r)

is the total position-dependent macroscopic transport cross section (absorption + scattering), and

σ_{s} (r, Ω^{'} - Ω)

is the differential macroscopic scattering cross section. The radiance using spherical harmonics at position

r

in the direction of unit vector Ω can be expanded as a series of the form,

Eq. 2

Φ (r, Ω) = \sum_{l = 0}^{N} \sum_{m = 0}^{l} (2 l + 1) p_{l}^{m} (\cos θ) [ψ_{l m} (r) \cos (m ϕ) + γ_{l m} (r) \sin (m ϕ)],

where

p_{l}^{m} (\cos θ)

are associated Legendre polynomials of order

l

,

m

, with

ψ_{l m} (r)

and

γ_{l m} (r)

as the coefficients of moments of the series, and θ and ϕ as the axial and azimuthal angles of Ω, respectively.²¹

p_{3}

,

p_{5}

, and

p_{7}

approximations have been performed in various geometries for 2-D problems, and the results validated that the

p_{3}

approximation is sufficiently accurate.²⁰ As such, only the

p_{3}

approximation is used here for an efficient 3-D reconstruction. Based on Eqs. 1, 2, the Eqs. 3, 4, 5, 6, 7, 8 are the derived 3-D

p_{3}

approximated-HD diffusion equations,²¹

Eq. 3

D (\frac{\partial^{2} ψ_{1}}{\partial^{2} x} + \frac{\partial^{2} ψ_{1}}{\partial^{2} y} + \frac{\partial^{2} ψ_{1}}{\partial^{2} z}) - μ_{a} ψ_{1} + D (2 \frac{\partial^{2} ψ_{2}}{\partial^{2} z} - \frac{\partial^{2} ψ_{2}}{\partial^{2} x} - \frac{\partial^{2} ψ_{2}}{\partial^{2} y}) + 6 D (\frac{\partial^{2} ψ_{3}}{\partial^{2} x}) + 6 D (- \frac{\partial^{2} ψ_{3}}{\partial^{2} y}) + 6 D \frac{\partial^{2} ψ_{5}}{\partial z \partial x} + 6 D \frac{\partial^{2} ψ_{6}}{\partial z \partial y} + 6 D \frac{\partial^{2} ψ_{4}}{\partial x \partial y} + 6 D \frac{\partial^{2} ψ_{4}}{\partial y \partial x} = - S (r),

Eq. 4

\frac{25}{7} D (\frac{\partial^{2} ψ_{2}}{\partial^{2} x} + \frac{\partial^{2} ψ_{2}}{\partial^{2} y}) + 2 D (\frac{\partial^{2} ψ_{1}}{\partial^{2} z} - \frac{1}{2} \frac{\partial^{2} ψ_{1}}{\partial^{2} x} - \frac{1}{2} \frac{\partial^{2} ψ_{1}}{\partial^{2} y}) + \frac{60}{7} D (\frac{\partial^{2} ψ_{3}}{\partial^{2} y} - \frac{\partial^{2} ψ_{3}}{\partial^{2} x} - \frac{\partial^{2} ψ_{4}}{\partial x \partial y} - \frac{\partial^{2} ψ_{4}}{\partial y \partial x}) + \frac{55}{7} D \frac{\partial^{2} ψ_{2}}{\partial^{2} z} + \frac{30}{7} D (\frac{\partial^{2} ψ_{5}}{\partial z \partial x} + \frac{\partial^{2} ψ_{6}}{\partial z \partial y}) - 5 μ_{t}^{'} ψ_{2} = 0,

Eq. 5

D (\frac{\partial^{2} ψ_{1}}{\partial^{2} x} - \frac{\partial^{2} ψ_{1}}{\partial^{2} y}) - \frac{10}{7} D (\frac{\partial^{2} ψ_{2}}{\partial^{2} x} - \frac{\partial^{2} ψ_{2}}{\partial^{2} y}) + \frac{30}{7} D (\frac{\partial^{2} ψ_{5}}{\partial x \partial z} - \frac{\partial^{2} ψ_{6}}{\partial y \partial z}) + \frac{90}{7} D (\frac{\partial^{2} ψ_{3}}{\partial^{2} x} + \frac{\partial^{2} ψ_{3}}{\partial^{2} y} + \frac{1}{3} \frac{\partial^{2} ψ_{3}}{\partial^{2} z}) - 10 μ_{t}^{'} ψ_{3} = 0,

Eq. 6

\frac{1}{2} D \frac{\partial^{2} ψ_{1}}{\partial x \partial y} + \frac{1}{2} D \frac{\partial^{2} ψ_{1}}{\partial y \partial x} - \frac{5}{7} D \frac{\partial^{2} ψ_{2}}{\partial x \partial y} - \frac{5}{7} D \frac{\partial^{2} ψ_{2}}{\partial y \partial x} + \frac{45}{7} D (\frac{\partial^{2} ψ_{4}}{\partial^{2} x} + \frac{\partial^{2} ψ_{4}}{\partial^{2} y} + \frac{1}{3} \frac{\partial^{2} ψ_{4}}{\partial^{2} z}) + \frac{15}{7} D (\frac{\partial^{2} ψ_{5}}{\partial y \partial z} + \frac{\partial^{2} ψ_{6}}{\partial x \partial z}) - 5 μ_{t}^{'} ψ_{4} = 0,

Eq. 7

2 D \frac{\partial^{2} ψ_{1}}{\partial x \partial z} + \frac{10}{7} D \frac{\partial^{2} ψ_{2}}{\partial x \partial z} + \frac{45}{7} D (\frac{\partial^{2} ψ_{5}}{\partial^{2} z} + \frac{\partial^{2} ψ_{5}}{\partial^{2} x} + \frac{1}{3} \frac{\partial^{2} ψ_{5}}{\partial^{2} y}) + \frac{60}{7} D (\frac{\partial^{2} ψ_{3}}{\partial x \partial z} + \frac{\partial^{2} ψ_{4}}{\partial y \partial z} + \frac{1}{2} \frac{\partial^{2} ψ_{6}}{\partial z \partial y}) - 5 μ_{t}^{'} ψ_{5} = 0,

Eq. 8

2 D \frac{\partial^{2} ψ_{1}}{\partial y \partial z} + \frac{10}{7} D \frac{\partial^{2} ψ_{2}}{\partial y \partial z} + \frac{45}{7} D (\frac{\partial^{2} ψ_{6}}{\partial^{2} z} + \frac{1}{3} \frac{\partial^{2} ψ_{6}}{\partial^{2} x} + \frac{\partial^{2} ψ_{6}}{\partial^{2} y}) + \frac{60}{7} D (\frac{\partial^{2} ψ_{4}}{\partial x \partial z} - \frac{\partial^{2} ψ_{3}}{\partial y \partial z} + \frac{1}{2} \frac{\partial^{2} ψ_{5}}{\partial x \partial y}) - 5 μ_{t}^{'} ψ_{6} = 0,

in which

μ_{a}

is the absorption coefficient,

μ_{t}^{'} = μ_{a} + μ_{s}

, where

μ_{s}

is the scattering coefficient,

D = 1 ∕ [3 μ_{a} + 3 μ_{s} (1 - g)]

is the diffusion coefficient, and

g

is the scattering anisotropic factor. It should be pointed out that the prior higher order diffusion equations can be further reduced to the diffusion approximation by setting

ψ_{2} = ψ_{3} = ψ_{4} = ψ_{5} = ψ_{6} = 0

,

Eq. 9

\nabla \cdot D (r) \nabla ψ_{1} (r) - μ_{a} (r) ψ_{1} (r) = - S (r) .

Appropriate boundary conditions should be specified for the HD equations. In this study, type 3 boundary conditions (BCs) are used for the first component

ψ_{1}

, while type 1 BCs are assumed for the other five components, i.e.,

Eq. 10

- D \nabla ψ_{1} \cdot n = α ψ_{1}, ψ_{2} = ψ_{3} = ψ_{4} = ψ_{5} = ψ_{6} = 0,

where α is the boundary condition coefficient (α = 0.467 for vacuum case). The choice of the boundary conditions have been discussed in detail in our previous publications.^{21, 25} Generally, type 3 BCs are applied in diffuse optical imaging and contain the linear combination of the photon density and the current along the boundary, while type 1 BCs are the commonly used boundary conditions for the RTE model. The primary reason here for using type 3 BCs for the first component and type 1 BCs for other components is for mathematical simplicity. In fact, type 3 BCs can also be applied to components 2 through 6, but we believe this may not add any benefit, as the first component is the measured quantity for image reconstruction. In finite element discretization, the quantities

ψ (r)

,

μ_{t}^{'}

,

D

, and

μ_{a}

are spatially discretized as,

Eq. 11

ψ_{k} = \sum_{i = 1}^{N} ϕ_{i} {(ψ_{k})}_{i}, μ_{t}^{'} = \sum_{i = 1}^{N} ϕ_{i} {(μ_{t}^{'})}_{i}, D = \sum_{i = 1}^{N} D_{i} ϕ_{i}, μ_{a} = \sum_{i = 1}^{N} {(μ_{a})}_{i} ϕ_{i},

in which

N

is the node number of the 3-D finite element mesh used and

ϕ_{i}

is the basis function. As such, the discretized form of Eqs. 3, 4, 5, 6, 7, 8, 10 is written as

Eq. 12

[A] {ψ} = {b},

in which the elements of the matrix [A] are integrated over the problem domain (V) and boundary domain (Γ), and {b} is the source vector. The inverse solution can be obtained through the following equations:

Eq. 13

[A] {\partial ψ ∕ \partial χ} = {\partial b ∕ \partial χ} - [\partial A ∕ \partial χ] {ψ},

Eq. 14

(I^{T} I + λ I) Δ χ = I^{T} (ψ^{o} - ψ^{c}),

in which χ expresses

D

or

μ_{a}

,

I

is the jacobian matrix formed by

\partial ψ ∕ \partial χ

at the boundary measurement sites; λ, a scalar, and I, the identity matrix, are used to realize the invertible system [Eq. 14];

Δ χ = {(Δ D_{1}, Δ D_{2}, \dots, Δ D_{N}, Δ μ_{a, 1}, Δ μ_{a, 2}, \dots, Δ μ_{a, N})}^{T}

;

ψ^{o} = {(ψ_{1}^{1 o}, ψ_{1}^{2 o}, \dots, ψ_{1}^{M o})}^{T}

, and

ψ^{c} = {(ψ_{1}^{1 c}, ψ_{1}^{2 c}, \dots, ψ_{1}^{M c})}^{T}

, where

ψ_{1}^{o}

and

ψ_{1}^{c}

are observed and computed photon density (the first component) for

i = 1, 2 \dots, M

boundary locations. It should be noted that only the first component, as the measurements, is used for the inverse computation due to the fact that the other five moment variables are equal to zero at the boundary measurement positions.

For x-ray guided 3-D DOT reconstruction, the following updating equation is used,

Eq. 15

Δ χ = {(I^{T} I + I^{T} I + λ I + L^{T} L)}^{- 1} [I^{T} (ψ^{o} - ψ^{c})],

where the x-ray structural a-priori information is incorporated into the iterative process using the spatially variant filter matrix L. In this study, the Laplacian-type filter matrix was used, and its elements

L_{i j}

were constructed according to the visible region or tissue type it was associated with, with x-ray derived priors as follows:^{28, 29, 30}

Eq. 16

L_{i j} = {\begin{matrix} 1 & when i = j \\ - 1 ∕ n n & when i, j \subset one region \\ 0 & when i, j \subset different region \end{matrix},

where

n n

is the finite element node number within a tissue type. Since the finger joint tissues are highly heterogeneous, we have previously shown that a modified Newton method with excellent convergent property is required.¹⁸ Thus the final updating equation is modified as follows:

Eq. 17

χ_{new} = χ_{old} + ζ Δ χ (0 < ζ ⩽ 1),

where ζ is computed from a backtracking line search. Thus the image formation task here is to update optical property distributions via iterative solution of Eqs. 12, 14, or 15 so that a weighted sum of the squared difference between computed and measured photon density can be minimized.

3. Results

In this section, we present reconstruction results using the HD model with phantom and in-vivo data. For comparison, reconstructed optical images using the DA-based algorithm are also provided.

Both phantom and in-vivo data were collected using our homemade 64 × 64-channel photodiodes-based stand-alone DOT system or the hybrid x-ray/DOT imaging system that integrates a modified mini-C-arm x-ray system with the 64 × 64-channel DOT system, which were all detailed elsewhere.^{8, 22, 28} This DOT system consists of several laser modules, a hybrid light delivery subsystem, a fiber optics/tissue interface, a data acquisition module, and light detection modules.²² The cylindrical fiber optics/tissue interface is composed of 64 source and 64 detector fiber bundles that are positioned in four layers along the surface of a plexiglass container and cover a volume of 15 (height) by 30 mm (diameter). In each layer, 16 source and 16 detector fiber bundles are alternatively arranged. The space between the finger and the wall of the plexiglass container is filled with tissue-like phantom materials as coupling media, consisting of distilled water, agar powder, Indian ink, and Intralipid, giving an absorption coefficient of $0.014 {mm}^{- 1}$ and a reduced scattering coefficient of $1.0 {mm}^{- 1}$ .^{8, 22} For in-vivo study, we focused on imaging the distal interphalangeal (DIP) finger joints from healthy volunteers and OA patients. The study was approved by the Institutional Review Board at the University of Florida. In the optical exam, the index finger of the participant was inserted into the 30-mm-diam cylindrical solid phantom medium.⁸

For DOT-alone imaging, we just present the in-vivo results as we have previously reported the phantom study using both the HD and DA models.²¹ Reconstructions were performed using a 3-D mesh with 3009 nodes and 12,800 tetrahedral elements. The mesh used is uniform in the height direction and the average distance of nodes is actually about 1.25 mm. In addition, in the region around the joint cavity/spacing, the mesh was discretized in terms of number of nodes per ${mm}^{3}$ to ensure that the average distance of nodes is around 1 mm in the joint spacing. The relatively coarse mesh density generated inside the bone domain has no significant influence on our reconstruction accuracy, which has been demonstrated and validated by our previous simulation tests and phantom experiments.^{8, 18, 21, 28} As such, the mesh size used is good enough for diffuse optical imaging that only has a resolution limit of 2 to 3 mm. For the imaging systems, 64 sources and 64 detectors were distributed uniformly along the surface of the phantom at four planes ( $z = 2.5$ , $z = 7.5$ , $z = 12.5$ , and $z = 17.5 mm$ ; 16 sources and 16 detectors at each plane). The initial guesses used are $μ_{a} = 0.03$ and $μ_{s}^{'} = 1.0$ for all the tests. For the reconstruction using the HD model, the scattering anisotropic factor $g = 0.7$ . Figures 1a and 1b and Figs. 2a and 2b provide reconstructed two selected coronal ( $y - z$ ) slices for a healthy and OA joint, respectively. The recovered quantitative optical properties are also given in Table 1 . We note that the results provided for both DA and HD models required six iterations for the inverse computation. If we use a LU decomposition-based forward solution solver, the DA and HD models require, respectively, 9.0 and 15.0 min per iteration when a FE mesh of 3009 nodes with 12,800 tetrahedral elements is employed in a high performance PC. If a gradient-based solver is used, the computation time can be reduced to 12 min per iteration for the HD model.

Fig. 1

(a) Reconstructed absorption and (b) scattering images at selected coronal planes for an OA joint using the DA (left column) and HD model (right column).

Fig. 2

(a) Reconstructed absorption and (b) scattering images at selected coronal planes for a healthy joint using the DA (left column) and HD model (right column).

Table 1

Averaged values of absorption and scattering coefficients of bone and joint tissues using the HD and DA models without x-ray guidance.

Test	μa (mm−1) (bone)			μa (mm−1) (joint)			μs′ (mm−1) (bone)			μs′ (mm−1) (joint)
Test	HD	DA	Difference	HD	DA	Difference	HD	DA	Difference	HD	DA	Difference
OA	0.049	0.053	6.3%	0.041	0.0452	10%	2.20	2.30	4.3%	1.85	1.94	4.6%
Healthy	0.048	0.051	6.2%	0.029	0.0337	15%	1.91	1.96	2.6%	1.08	1.16	6.9%

In the phantom experiments with x-ray guidance, a 30-mm-diam cylindrical solid phantom was used as the background medium. Two 20-mm-diam cylindrical solid objects (3-mm off $z$ -axis) mimicking bones were embedded in the background medium. The spacing (“cartilage”) between the two “bones” was 2.5 mm. The optical properties for the background were $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ , while the optical properties of the “cartilage” were assumed to be the same as the background medium. The optical properties for the mimicking bones are $μ_{a} = 0.07 {mm}^{- 1}$ and $μ_{s}^{'} = 4.0 {mm}^{- 1}$ . The recovered images using the DA (left column) and HD equations (right column) at selected longitudinal planes are provided in Figs. 3a and 3b .

Fig. 3

(a) Reconstructed absorption and (b) scattering images at selected planes for the phantom test using the DA (left column) and HD model (right column) with x-ray guidance.

For in-vivo imaging with x-ray guidance, image reconstruction of the DIP finger joints with the optical coupling phantom/media (30 × 20 mm in volume) were performed with a finite element mesh of 2705 nodes and 13,440 tetrahedral elements. This generated mesh is not uniform—the prior x-ray image of bones allows us to use a nodal distance of ∼0.4 mm in the joint spacing, which is accurate enough for joint space imaging. In particular, we are interested in imaging the optical properties of the joint region, and the fine mesh around the cavity is sufficient for our study. Again, the coarse mesh density for the bones has no significant effect on the reconstruction accuracy. The initial optical properties used were optimized based on an x-ray guided forward-fitting algorithm for the approximated joint and bone regions/tissues.²⁴ The optical properties for the background were $μ_{a} = 0.014 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ , which also served as initial guesses for the background medium. Figures 4a and 4b plot the selected absorption and scattering slices of the recovered 3-D images for an OA joint using the DA and HD models, while Figs. 5a and 5b display the reconstructed slices for a healthy joint. The optical absorption and scattering coefficient values are provided in Table 2 for both the phantom and in-vivo data.

Fig. 4

(a) Reconstructed absorption and (b) scattering images at selected coronal planes for an OA joint using the DA (left column) and HD model (right column) with x-ray guidance.

Fig. 5

(a) Reconstructed absorption and (b) scattering images at selected coronal planes for a healthy joint using the DA (left column) and HD model (right column) with x-ray guidance.

Table 2

Averaged values of absorption and scattering coefficients of bone and joint tissues using the HD and DA models with x-ray guidance. Note that difference = (DA-HD) × 100%/DA.

Test	μa (mm−1) (bone)			μa (mm−1) (joint)			μs′ (mm−1) (bone)			μs′ (mm−1) (joint)
Test	HD	DA	Difference	HD	DA	Difference	HD	DA	Difference	HD	DA	Difference
Phantom	0.056	0.054	− 3.8%	0.0098	0.0095	− 3.2%	3.51	3.53	0.6%	1.003	1.02	1.3%
OA	0.046	0.047	1.9%	0.023	0.0230	2.1%	2.25	2.27	0.7%	1.527	1.54	0.7%
Healthy	0.040	0.041	2.3%	0.0098	0.0100	2.0%	2.15	2.16	0.6%	0.952	0.97	1.4%

4. Discussion

From the reconstructed optical images shown in Figs. 1 and 2, we note that the bones are clearly delineated for both the normal and OA joints using DA and HD models. Importantly, compared with the optical property value of the bones, we observe a large drop in the magnitude of absorption and scattering coefficients within the joint space for the healthy case, while only a small drop is seen for the OA joint. These observations also agree well with the x-ray guided optical findings, as displayed in Figs. 4 and 5. In particular, the image resolution from the hybrid imaging method is significantly improved compared to that from DOT alone. We see that x-ray guided optical images exhibit apparent narrowing of the joint space for the OA joint (Fig. 4) relative to those for the healthy joint (Fig. 5).

It is noted from Figs. 1 and 2 that model errors do exist when the DA is used to handle cases involving low scattering, small source/detector distances or large variations of optical properties. As shown in Table 1, the difference in recovered quantitative optical properties between the two models can be as large as 15% (for the absorption coefficient of joint tissue). This demonstrates that the $p 3$ approximated HD model can produce significantly improved reconstruction of both absorption and scattering images for joint imaging. In addition, it should also be pointed out that while a transport or HD model can certainly reduce a large amount of errors in forward computation,²¹ the model errors existing in the forward calculation do not appear to propagate much into the image reconstruction. This is understandable, because the accuracy of an inverse solution depends not only on the accuracy of the forward model, but also on the quality of the experimental data and the use of robust regularization techniques.

More interestingly, model errors appear to lead to even smaller error in reconstruction when the x-ray a-priori structural information is incorporated into DOT reconstruction, as shown in Figs. 3, 4, 5. The difference for the recovered optical properties is less than 4% between the two models, which is evident from Table 2. This is due to the fact that a more accurate modeling of photon migration in tissue can be achieved using the DA when anatomical a-priori information becomes available. The use of prior structural information eliminates the need to look for spatial/anatomy information in the optical inversion, which ensures that optical property profiles are the only parameter(s) that need to be recovered.

In summary, we demonstrate in this study that the improvement resulted from HD-based reconstruction is significant over the DA model for DOT-alone imaging of finger joints. The improvement is moderate using the HD model when x-ray guidance is used because of the high quality of reconstruction already available from the DA model under such a situation.

Acknowledgments

This research was supported in part by the National Institutes of Health (R01 AR048122).

References

1.

A. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today, 48 34 –40 (1995). https://doi.org/10.1063/1.881445 0031-9228 Google Scholar

2.

B. J. Tromberg, N. Shah, R. Lanning, A. Gerussi, J. Espinoza, T. Pham, L. Svaasand, and J. Beuthan, “Noninvasive in vivo characterization of breast cancer using photon migration spectroscopy,” Neoplasia, 2 26 –40 (2000). https://doi.org/10.1038/sj.neo.7900082 1522-8002 Google Scholar

3.

Q. Zhang, T. J. Brukilacchio, A. Li, J. Scott, T. Chaves, E. Hillman, T. Wu, M. Chorlton, E. Rafferty, R. H. Moore, D. B. Kopans, and D. A. Boas, “Coregistered tomography x-ray and optical breast imaging: initial results,” J. Biomed. Opt., 10 024033 (2005). https://doi.org/10.1117/1.1899183 1083-3668 Google Scholar

4.

R. L. Barbour, H. L. Graber, Y. Wang, J. H. Chang, and R. Aronson, “A perturbation approach for optical diffusion tomography using continuous-wave and time-resolved data,” IS11 87 –120 (1993) Google Scholar

5.

J. Hebden S. Arridge, and D. Delpy, “Optical imaging in medicine: I: experimental techniques,” Phys. Med. Biol., 42 825 –840 (1997). https://doi.org/10.1088/0031-9155/42/5/007 0031-9155 Google Scholar

6.

H. Dehghani, D. Delpy, and S. Arridge, “Photon migration in non-scattering tissue and the effects on image reconstruction,” Phys. Med. Biol., 44 2897 –2906 (1999). https://doi.org/10.1088/0031-9155/44/12/303 0031-9155 Google Scholar

7.

A. Pifferi, A. Torricelli, P. Taroni, A. Bassi, E. Chikoidze, E. Giambattistelli, and R. Cubeddu, “Optical biopsy of bone tissue: a step toward the diagnosis of bone pathologies,” J. Biomed. Opt., 9 (3), 474 –480 (2004). https://doi.org/10.1117/1.1691029 1083-3668 Google Scholar

8.

Z. Yuan, Q. Zhang, E. Sobel, and H. Jiang, “3D diffuse optical tomography imaging of osteoarthritis: Initial results in finger joints,” J. Biomed. Opt., 12 034001 (2007). https://doi.org/10.1117/1.2737420 1083-3668 Google Scholar

9.

Y. Xu, N. Iftimia, H. Jiang, L. Key, and M. Bolster, “Three-dimensional diffuse optical tomography of bones and joints,” J. Biomed. Opt., 7 (1), 88 –92 (2002). https://doi.org/10.1117/1.1427336 1083-3668 Google Scholar

10.

A. H. Hielscher, A. D. Klose, A. K. Scheel, B. Moa-Anderson, M. Backhaus, U. Netz, and J. Beuthan, “Sagittal laser optical tomography for imaging of rheumatoid finger joints,” Phys. Med. Biol., 49 1147 –1163 (2004). https://doi.org/10.1088/0031-9155/49/7/005 0031-9155 Google Scholar

11.

K. Ren, G. Bal, and A. H. Hielscher, “Transport and diffuse-based optical tomography in small volumes: a comparative study,” Appl. Opt., 46 6669 –6679 (2007). https://doi.org/10.1364/AO.46.006669 0003-6935 Google Scholar

12.

G. Bal and K. Ren, “Generalized diffuse model in optical tomography with clear layer,” J. Opt. Soc. Am. A, 20 2355 –2364 (2003). https://doi.org/10.1364/JOSAA.20.002355 0740-3232 Google Scholar

13.

V. Prapavat, W. Runge, J. mans, J. Beuthan, and G. Muller, “Development of a finger joint phantom for the optical simulation of early stages of rheumatoid arthritis,” Biomed. Tech., 42 319 –326 (1997). 0013-5585 Google Scholar

14.

A. D. Klose and E. W. Larson, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comp. Physiol., 220 441 –470 (2006). 0373-0859 Google Scholar

15.

E. D. Aydin, C. RE de Oliveira, and A. JH Goddard, “A comparison between transport and diffusion calculations using a finite element-spherical harmonics radiation transport method,” Med. Phys., 29 2013 –2023 (2002). https://doi.org/10.1118/1.1500404 0094-2405 Google Scholar

16.

M. Xu, W. Cai, M. Lax, and R. Alfano, “Photon transport forward model for imaging in turbid media,” Opt. Lett., 26 1066 –1068 (2001). https://doi.org/10.1364/OL.26.001066 0146-9592 Google Scholar

17.

A. Dorn, “Transport-backtransport method for optical tomography,” Inverse Probl., 14 1107 –1130 (1998). https://doi.org/10.1088/0266-5611/14/5/003 0266-5611 Google Scholar

18.

Z. Yuan and H. Jiang, “Image reconstruction schemes that combines modified Newton method and efficient initial guess estimate for optical tomography of finger joints,” Appl. Opt., 46 2757 –2768 (2007). https://doi.org/10.1364/AO.46.002757 0003-6935 Google Scholar

19.

A. D. Kim and J. B. Keller, “Light propagation in biological tissues,” J. Opt. Soc. Am. A, 20 92 –98 (2003). https://doi.org/10.1364/JOSAA.20.000092 0740-3232 Google Scholar

20.

D. A. Boas, “Photon migration within the P3 approximation,” Proc. SPIE, 2389 240 –247 (1995). https://doi.org/10.1117/12.209973 0277-786X Google Scholar

21.

Z. Yuan, X. Hu, and H. Jiang, “A higher-order diffusion model for three-dimensional photon migration and image reconstruction in optical tomography,” Phys. Med. Biol., 54 65 –88 (2009). https://doi.org/10.1088/0031-9155/54/1/005 0031-9155 Google Scholar

22.

Q. Zhang and H. Jiang, “Three-dimensional diffuse optical tomography of simulated hand joints with a 64 × 64-channel photodiodes-based optical system,” J. Opt. A, Pure Appl. Opt., 7 224 –231 (2005). https://doi.org/10.1088/1464-4258/7/5/003 1464-4258 Google Scholar

23.

K. D. Paulsen and H. Jiang, “Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization,” Appl. Opt., 35 3447 –3458 (1996). https://doi.org/10.1364/AO.35.003447 0003-6935 Google Scholar

24.

N. Iftimia and H. Jiang, “Quantitative optical image reconstruction of turbid media by use of direct-current measurements,” Appl. Opt., 39 5256 –5261 (2000). https://doi.org/10.1364/AO.39.005256 0003-6935 Google Scholar

25.

H. Jiang, K. Paulsen, U. Osterberg, B. Pogue, and M. Patterson, “Optical image reconstruction using frequency-domain data: simulations and experiments,” J. Opt. Soc. Am. A, 13 253 –266 (1996). https://doi.org/10.1364/JOSAA.13.000253 0740-3232 Google Scholar

26.

K. D. Paulsen and H. Jiang, “Spatially-varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys., 22 691 –701 (1995). https://doi.org/10.1118/1.597488 0094-2405 Google Scholar

27.

H. Jiang, “Optical image reconstruction based on the third-order diffusion equations,” Opt. Express, 4 241 –246 (1999). https://doi.org/10.1364/OE.4.000241 1094-4087 Google Scholar

28.

Z. Yuan, Q. Zhang, E. Sobel, and H. B. Jiang, “Tomographic x-ray-guided three-dimensional diffuse optical tomography of osteoarthritis in the finger joints,” J. Biomed. Opt., 13 044006 (2008). https://doi.org/10.1117/1.2965547 1083-3668 Google Scholar

29.

B. Brooksby, S. Jiang, H. Dehghani, B. W. Pogue, and K. D. Pausen, “Combining near-infrared tomography and magnetic resonance imaging to study in vivo breast tissues: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. Biomed. Opt., 10 011504 (2005). https://doi.org/10.1117/1.2098627 1083-3668 Google Scholar

30.

B. W. Pogue, T. O. McBride, J. O. Prewitt, U. L. Osterberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt., 38 2950 –2961 (1996). https://doi.org/10.1364/AO.38.002950 0003-6935 Google Scholar

Citation Download Citation

Zhen Yuan, Qizhi Zhang, Eric S. Sobel, and Huabei Jiang "Comparison of diffusion approximation and higher order diffusion equations for optical tomography of osteoarthritis," Journal of Biomedical Optics 14(5), 054013 (1 September 2009). https://doi.org/10.1117/1.3233655

Published: 1 September 2009

Access the abstract

JOURNAL ARTICLE
8 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 14 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Diffusion

Scattering

X-rays

Absorption

Tissues

3D modeling

Bone

1.

Introduction

2.

Reconstruction Methods

Eq. 1

Eq. 2

Eq. 3

Eq. 4

Eq. 5

Eq. 6

Eq. 7

Eq. 8

Eq. 9

Eq. 10

Eq. 11

Eq. 12

Eq. 13

Eq. 14

Eq. 15

Eq. 16

Eq. 17

3.

Results

Fig. 1

Fig. 2

Table 1

Fig. 3

Fig. 4

Fig. 5

Table 2

4.

Discussion

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years