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 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 approximated-RTE-based reconstruction is 40 to 60 times slower than DA based.11
We have developed a 3-D reconstruction method based on simplified spherical harmonics approximated-RTE or higher order diffusion equations.21 Our simulations and phantom tests21 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.21
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.8 In addition, we evaluate the performance of the HD model using an x-ray/DOT combined multimodality approach.28 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.28 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.
Reconstruction Methodsis the position vector of a photon propagation along the unit direction Ω, is the energy radiance, is the source term, is the total position-dependent macroscopic transport cross section (absorption + scattering), and is the differential macroscopic scattering cross section. The radiance using spherical harmonics at position in the direction of unit vector Ω can be expanded as a series of the form, are associated Legendre polynomials of order , , with and as the coefficients of moments of the series, and θ and ϕ as the axial and azimuthal angles of Ω, respectively.21 , , and approximations have been performed in various geometries for 2-D problems, and the results validated that the approximation is sufficiently accurate.20 As such, only the 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 approximated-HD diffusion equations,21 is the absorption coefficient, , where is the scattering coefficient, is the diffusion coefficient, and 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 , , while type 1 BCs are assumed for the other five components, i.e.,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 , , , and are spatially discretized as, is the node number of the 3-D finite element mesh used and is the basis function. As such, the discretized form of Eqs. 3, 4, 5, 6, 7, 8, 10 is written as or , is the jacobian matrix formed by at the boundary measurement sites; λ, a scalar, and I, the identity matrix, are used to realize the invertible system [Eq. 14]; ; , and , where and are observed and computed photon density (the first component) for 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,were constructed according to the visible region or tissue type it was associated with, with x-ray derived priors as follows:28, 29, 30 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.18 Thus the final updating equation is modified as follows:12, 14, or 15 so that a weighted sum of the squared difference between computed and measured photon density can be minimized.
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.22 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 and a reduced scattering coefficient of .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.8
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.21 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 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 ( , , , and ; 16 sources and 16 detectors at each plane). The initial guesses used are and for all the tests. For the reconstruction using the HD model, the scattering anisotropic factor . Figures 1a and 1b and Figs. 2a and 2b provide reconstructed two selected coronal ( ) 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.
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)|
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 -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 and , 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 and . The recovered images using the DA (left column) and HD equations (right column) at selected longitudinal planes are provided in Figs. 3a and 3b .
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.24 The optical properties for the background were and , 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.
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)|
|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%|
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 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,21 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.
This research was supported in part by the National Institutes of Health (R01 AR048122).