An efficient Jacobian reduction method for diffuse optical image reconstruction

Matthew E. Eames; Brian W. Pogue; Phaneendra K. Yalavarthy; Hamid Dehghani

doi:10.1364/OE.15.015908

1. Introduction

Near-Infrared (NIR) Diffuse Optical Tomography (DOT) is a non-invasive imaging technique whereby NIR light between the wavelengths of 650nm and 900nm is injected into tissue using optical fibers at the surface of the volume of interest and the emergent light is then measured at points along the same surface [1–3]. The measured transmission signals, sometimes known as “boundary data” are then used with a model of light propagation in tissue to derive intrinsic optical properties of the tissue or functional information such as the total hemoglobin, oxygen saturation, water content and scattering properties [4–6].

Diffuse Optical Imaging (DOI) has become an important tool in medical imaging research and has shown potential as a non-invasive technique for quantifying brain function activity [7] as well as breast cancer detection and characterization [2, 8, 9]. Much work has gone into both the generation of accurate forward models to describe the paths taken by the photons in the tissue and the inverse models to accurately recover images of the optical properties of the tissue [1, 3, 10–12].

For the forward problem, depending on the imaging domain, geometry and size, a number of models based on the Radiative Transport Equation (RTE) have been developed [13–15]. However, since the light in soft tissue becomes diffuse within a couple of scattering distances, diffusion based calculations using the Finite Element Method (FEM) have been employed extensively in the case of imaging thick tissues [2, 3, 16]. The use of the Diffusion Approximation for the estimation of light propagation in tissue relies on two main assumptions that (i) the scattering component of light interaction at the applied wavelength is much greater than the absorption and (ii) that the distance between the applied and measured light is much larger than the scattering distance. There have been other proposed hybrid methods that have been used to estimate light propagation in the presence of void regions where scattering is zero [16].

The inverse problem in DOT is concerned with the estimation of the internal intrinsic optical property distribution based on the measured boundary data. Several image reconstruction techniques have been developed. These include the use of gradient based minimization techniques that do not require the explicit calculation of the Jacobian (also known as the sensitivity or weight matrix) [17, 18] but although these are computationally efficient they require the use of an optimization scheme which is not straightforward. Alternatively, the use of full-Newton optimization schemes, whereby a Jacobian (which relates the measured boundary data to the internal optical properties) is calculated and its variant is inverted, have widely been used as these are easy to implement and can easily be generalized for both simple and complex image reconstruction algorithms [12, 19]. Other alternative methods for the reduction of the number of parameters include the use of reconstruction basis [20] as well as the use of adaptive meshing algorithms whereby only the regions of interest are refined to provide numerical accuracy [21]

The problem of image reconstruction in DOT is therefore two fold; there exists a need to ensure that the model used for the prediction of light distribution in tissue is accurate and that the image reconstruction algorithm is both computationally efficient and reliable in estimating the optical properties of the domain in tissue. The use of numerical algorithms such as FEM for the forward model rely on the fact that the volume being imaged and modeled is accurately defined and that the mesh discretization used for the calculation of the FEM solution is adequate and numerically stable. There has been a number of studies that have looked at the effect of mesh resolution in the forward problem in DOT [22, 23] as well as other similar imaging modalities [24, 25]. These studies have essentially demonstrated that in order to accurately calculate the forward problem, the mesh resolution needs to be of a high quality to ensure numerical accuracy. On the other hand, the inverse problem is ill-posed and typically under determined. That is, there typically exist a large number of unknown parameters as compared to the number of available boundary measurements. A number of different strategies have been proposed such as the use of reconstruction basis that addresses both of the problems defined [1, 26].

Assuming that an optical source detector system consisting of 16 fibers is used to image a three-dimensional (3D) volume, this would give rise to 240 measurements (assuming amplitude only data), whereby each source injects light into the domain sequentially while the other 15 fibers measure the emerging light at the surface in parallel for each source location. In order to model the forward data and assuming a typical breast model, mesh resolutions of 20,000 to 30,000 nodes are required. For image reconstruction, one can either use a regular pixel based reconstruction or a second mesh of the imaging volume, but these contain approximately half of the forward mesh nodes, typically 10,000, unknowns for a single parameter reconstruction when structural prior information is unavailable. This implies that the Jacobian calculated which, in this instance relates intrinsic absorption to a change in the boundary data (log of the amplitude), will have a size of 240×10,000 which will need to be inverted. It becomes apparent, therefore, that the use of a method that would reduce the size of the Jacobian without causing any loss to the accuracy of the imaging algorithm would be beneficial. Moreover, storage of these big matrices reduces the computational efficiency and puts constraint on the resolution of the mesh due to the physical memory available on the computational resources.

In this work we present a method which reduces the size of the Jacobian matrix without compromising the numerical accuracy of both the forward model and the inverse problem. In this technique a fine mesh is used to for the calculation of the forward model. However for the inverse problem, mesh nodes in regions which have a low sensitivity to the boundary data are removed in an efficient manner as to reduce the size of the sensitivity matrix, thus reducing the computation time and memory usage. As an example consider a 2D circular model with a radius of 43 mm and optical properties of 0.01 mm^-1 and 1.0 mm^-1 for absorption and reduced scatter respectively, Fig. 1, for which the Jacobian for a single source and detector has been calculated. In each instance, the contour line which represents the threshold of the magnitude of the Jacobian at different levels is shown. It is evident that given a source detector combination, the total area and therefore the number of nodes within that contribute to the measured data is highly dependent on the amount of sensitivity from the Jacobian which is required.

Previous work has shown that for a given 3D problem where data is measured using a 2D slice rasterizing fashion, only datasets whose projection error (that is the difference between the measured 2D slice data and modeled data) is greater than a tolerance level can be utilized to reconstruct images using a ‘data subset’ algorithm [27]. In this current work however, a general framework is presented that will consider the magnitude of the complete Jacobian, rather than the projection error to provide an update for the non-linear image reconstruction over the whole domain, as compared to specific regions or data subsets.

Fig. 1. The normalized sensitivity for a circular model with a single source and detector. In each case, the contour line shows the threshold of the sensitivity as a percentage of total sensitivity throughout the model.

Download Full Size | PDF

Through the use of simulated data as well as measured phantom data, it is shown that by efficient removal of nodes within the forward mesh which do not contribute more than 1% change in data, no loss in imaging accuracy is seen, whereas the computational speed of the algorithm can be improved by 14 folds. This approach is also shown to be memory efficient, as the size of the Jacobian matrix reduced by 4 fold, leading to 16 fold increase in the available memory. Initial results are, for simplicity, limited to a single data type (log amplitude) and a single parameter reconstruction (absorption). However, to show potential and expandability to more data types and larger parameter (e.g. spectral) reconstruction, the use of frequency domain experimental data from a multi-layered gelatin phantom is also presented whereby both optical absorption and reduced scatter are reconstructed simultaneously using measured log amplitude and phase boundary data at 100 MHz.

2. Theory

For the forward problem we are interested in the fluence rate of photons through the tissue given a set of optical properties. We assume that the propagation of light can be accurately modeled by the diffusion approximation to the RTE in the frequency domain given by

- \nabla \cdot κ (r) \nabla Φ (r, ω) + (μ_{a} (r) + \frac{i ω}{c_{m} (r)}) Φ (r, ω) = q_{0} (r, ω),

where q ₀(r,ω) is an isotropic source, Φ(r,ω) is the photon fluence rate at position r and modulation frequency ω. κ is the diffusion coefficient given by

κ = 1 ⁄ 3 (μ_{a} + {μ'}_{s}),

where µ_a and µ_s′ are absorption and reduced scattering (or transport scattering) coefficients respectively and c_m(r) is the speed of light in the medium. The air-tissue boundary is derived from a type III condition given by

Φ (ξ) + 2 A \hat{n} \cdot κ (ξ) \nabla Φ (ξ) = 0,

where ξ is a point on the external boundary and A depends upon the relative refractive index mismatch between the tissue domain and air derived from Fresnel’s law.

The inverse problem is solved with the aim of recovering the optical properties of the tissue µ=[µ_a, µ_s′] for each node with the FEM mesh. The inversion is achieved by finding the minimum between the measured data, Φ^M, and the calculated data, Φ^C, using a modified- Levenberg-Marquardt minimization approach given by

Ω = {}_{μ}^{min}{{\sum_{i = 1}^{NM} {(Φ_{i}^{M} - Φ_{i}^{C})}^{2}},}

where NM is the number of measurements and µ represents the optical properties that need to be reconstructed. For an ill-posed problem like DOT, minimizing Eq. (4) with respect to the optical properties µ, and considering only the linear terms leads to the update equation

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

where δµ is now the update vector for the optical properties. λ is implemented similar to Levenberg-Marquardt approach [28] where it starts at a given value scaled by maximum of diagonal of J^TJ and is systematically reduced at each iteration and J is the Jacobian matrix.

3. Jacobian reduction

The Jacobian matrix J is of the size number of measurements, NM, by the number of FEM nodes, NN, and is calculated by using the Adjoint method. For simplicity, limiting the Jacobian to the measured amplitude data and optical absorption, it links a change in log amplitude, at the boundary, with a change in absorption coefficient, µ_a. The Jacobian thus has the form

J = [\begin{matrix} \frac{δ \ln I_{1}}{δ μ_{a 1}} & \frac{δ \ln I_{1}}{δ μ_{a 2}} & \dots & \frac{δ \ln I_{1}}{δ μ_{aNN}} \\ \frac{δ \ln I_{2}}{δ μ_{a 1}} & \frac{δ \ln I_{2}}{δ μ_{a 2}} & \dots & \frac{δ \ln I_{2}}{δ μ_{aNN}} \\ \frac{δ \ln I_{3}}{δ μ_{a 1}} & \frac{δ \ln I_{3}}{δ μ_{a 2}} & \dots & \frac{δ \ln I_{3}}{δ μ_{aNN}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{δ \ln I_{NM}}{δ μ_{a 1}} & \frac{δ \ln I_{NM}}{δ μ_{a 2}} & \dots & \frac{δ \ln I_{NM}}{δ μ_{aNN}} \end{matrix}] .

For the inversion procedure the inverse of the Hessian matrix J^TJ is needed which has an order of NN×NN.

Due to the nature of tissue (highly scattering) it is unlikely that a photon entering the tissue at a given point will travel far from the source and the greatest sensitivity is found close to the source-detector plane of measurement, see for example Fig. 1. Therefore the sensitivity at distances far from the source or the detector is likely to be small. One option to remove the sensitivity of these (almost) non-contributing regions is to set the associated nodal values to zero. But these nodes will then remain within the matrix and contribute to the calculation of the inverse of the Hessian. Note that inversion of a N×N order matrix, requires an order of N³ operations and N² memory [29]. Thus the usage of full Hessian matrix with less sensitive parts being zero does not reduce the computational burden, unless a sparse matrix solver is used. But usage of sparse matrix solvers requires calculation of half band width of the matrix, which is difficult to measure in these cases where the structure of Hessian largely depends on node numbering and the sequence how the measurements are taken. Alternatively we can look at the total sensitivity within the volume and find regions where the value is very small, typically below a certain threshold. These regions can then be removed from the matrix thus reducing the size of the Hessian. In effect this method reduces the volume of interest for the case of the update equation without decreasing the mesh resolution, but the entire volume is used to accurately calculate the forward model.

In order to reduce the size of the Jacobian, the total sensitivity throughout the imaging domain is calculated and a new Jacobian, J̃_ij, is formed:

{\tilde{J}}_{ij} = {\begin{matrix} J_{ij} & if & \sum_{i = 1}^{NM} J_{ij} > = threshold \\ 0 & if & \sum_{i = 1}^{NM} J_{ij} < threshold \end{matrix},

where j corresponds to the node number within the domain. Given a new Jacobian, J̃_ij, where the total sensitivity for a given node is equal to zero, the entire column corresponding to that node is removed to produce a much smaller Jacobian matrix drastically reducing the number of operations and memory required for the calculation of the update of optical properties (Eq. 5).

4. Method and results

4.1 Simulated breast model

A uniform mesh of a female breast was used which consisted of 18374 nodes corresponding to 75215 linear tetrahedral elements, Fig. 2. To generate forward data (log amplitude) at a single wavelength of 785 nm, the forward model was assumed to contain 3 different regions of Adipose, glandular and tumor whose optical properties are given in Table 1, and cross section of absorption is shown in Fig. 3. Sixteen optical source detector fibers were used as the data collection strategy, which are placed equidistant at the mid-plane of the model. To reconstruct absorption only images using log amplitude data, the reduced Jacobian as defined earlier (Eq. 7), was utilized with varying levels of threshold of 0%, 1% and 10% for the reduction scheme. It should be noted that assuming a 0% threshold for the Jacobian reduction corresponds to using the full Jacobian and will serve as a reference to image detriment. The initial guess for optical absorption for the iterative image reconstruction was assumed to be that of a homogeneous breast model containing Glandular tissue. Images were reconstructed until the projection error, that is the difference between the estimated and actual data, did not improve more than 2%.

Table 1. The chromophore concentration of different regions within the 3D breast model, and the corresponding optical properties at 785 nm are listed.

View Table | View all tables in this article

Fig. 2. The female breast mesh used for generation of simulated data. The indentation represents the modelling of the optical fiber source and detectors.

Download Full Size | PDF

Fig. 3. A sagittal cross-section of the breast model shown in Fig. 2 illustrating the truest fields of optical absorption used at 785 nm.

Download Full Size | PDF

In order to demonstrate the level of total sensitivity (for log amplitude and optical absorption) throughout the whole model, a sagittal cross-section of the total sensitivity, at the first iteration of image reconstruction, is shown in Fig. 4. In each case a contour line is added to demonstrate the limits of sensitivity at different threshold levels. As evident, at high levels of threshold, the total sensitivity is constrained to regions near the sources and detectors and as the threshold in decreased more regions within the breast model contribute to the total sensitivity.

Fig. 4. The sagittal plot of the normalized total sensitivity for a breast model. In each case, the contour line shows the threshold of the sensitivity as a percentage of total sensitivity throughout the model.

Download Full Size | PDF

Reconstructed images using varying levels of threshold of 0%, 1% and 10% for the reduction scheme are shown in Fig. 5. In this proof of concept only optical absorption is reconstructed using the log amplitude of the boundary data. It is seen that as the threshold level increases the qualitative and quantitative accuracy of the reconstructed images decreases with the reconstructed changes seen more locally around the plane of imaging and near the sources and detectors. However, the reconstructed images at 0% threshold (that is a Jacobian that has not been reduced) are the same as those seen when a threshold of 1% is used. Table 2 shows the computational aspects of each reconstruction scheme including the number of nodes for the forward and inverse problem as well as computation time per iteration of image reconstruction. Given that the sensitivity decreases rapidly away from the plane of imaging the Jacobian reduction method removes areas of low sensitivity from the calculation and does not therefore update voxels within that region. As shown with a tolerance of 0%, even by including these regions of low sensitivity, due to the lack of photon sampling, any contribution from these regions is negligible and therefore the Jacobian reduction is a suitable method for removing their contribution within the algorithm.

Table 2. The details of each method used for the calculation of the forward model and the reconstruction of the breast mesh model.

View Table | View all tables in this article

Fig. 5. The sagittal and coronal cross-section of the reconstructed images of the breast model using the reduced and non-reduced (0% threshold) Jacobian.

Download Full Size | PDF

4.2 Experimental phantom

A cylindrical (height=40mm, radius=43mm) gelatin based multi-layered phantom was fabricated with different optical properties using India ink for absorption and TiO₂ for scattering. Layers are constructed by successively hardening gel solutions with different concentrations of ink and TiO₂. The outer layer of thickness 5mm, similar to a typical fatty breast layer, has the optical properties 0.0065 mm^-1 and 0.65 mm^-1 (at 785 nm) for absorption and reduced scatter respectively. The fibroglandular layer has optical properties of µ_a=0.01 mm^-1 and µ_s’=1.0 mm^-1 and has a radius of 38mm. An anomaly, which represents the tumor, of radius 8 mm of cylindrical shape stretching the entire z direction was placed approximately 20 mm from the centre with optical properties µ_a=0.02 mm^-1 and µ_s’=1.2 mm^-1. NIR data was collected using a frequency domain system with a modulation frequency of 100MHz using 16 optical fibers placed circularly along the mid-height of the phantom.

Fig. 6. The coronal cross-section of the reconstructed images for µ_a and µ_s’ of the cylindrical phantom using the reduced and non-reduced (0% threshold) Jacobian.

Download Full Size | PDF

Images were reconstructed using both log amplitude of the measured boundary data and phase using a cylindrical mesh containing 8990 nodes corresponding to 44803 linear tetrahedral elements. For image reconstruction the reduced Jacobian, as defined earlier, was extended to incorporate both the change in measured data (phase and amplitude) due to a change in µ_a and κ. Since both the absorption and scattering parameters determine the sensitivity, each were utilized accordingly using varying levels of threshold of 0%, 1%, 5% and 10% for the reduction scheme. The initial guess for optical absorption and reduced scattering for the iterative image reconstruction was assumed to be that of a homogeneous phantom, together with a uniform pixel basis of 20×20×20 for image reconstruction.

Reconstructed images for µ_a and µ_s’ using varying levels of threshold of 0%, 1%, 5% and 10% for the reduction scheme are shown in Fig. 6. For each method we find a large contrast for the reconstructed images for µ_a and as we would expect the contrast for the reduced scattering is small. There is no difference in the quantitative and qualitative accuracy between the reduction method with a tolerance of 1% and the usual reconstruction method (0%). As the tolerance increases there is a large saving in memory usage for the inversion of the Jacobian, table 3. With a tolerance of just 1% the number of nodes within the Jacobian is reduced by 8% saving 0.15Gb of memory for storage alone (18% reduction). At a tolerance 10%, the number of nodes decreases further to just 51% of the total, saving 0.73Gb of memory (74% reduction). The reconstructed images are extremely sensitive to the initial guess of the optical parameters and in this case the method has been unable to reconstruct the true µ_s’. However, the difference between the standard method (0%) and a tolerance of either 1% or 5% is negligible and demonstrating that the reduction method gives the same result at improved computational speed.

Table 3. The details of each method used for the calculation of the forward model and the reconstruction of the phantom.

View Table | View all tables in this article

5. Discussions and conclusions

A method for the reduction of the Jacobian matrix, which improves the computational speed and efficacy of the image reconstruction algorithm, is outlined and presented. This method utilizes a fine mesh to calculate the forward problem accurately, but considers the magnitude of the total sensitivity to discount any regions which have sensitivity below a given tolerance level and are therefore unlikely to contribute to the inverse problem. Simulated data using a realistic breast mesh were used to demonstrate the proof of concept using a single data type (log amplitude) and reconstruction the optical absorption. It is shown that constraining the sensitivity to different thresholds will only consider those regions whose total sensitivity is above the set threshold values, Fig. 4. For example, given an optical fiber setup which is radially placed around the breast, it is shown that the sensitivity within the plane of imaging is dominant, whereas the sensitivity near the chest wall or the nipple is not. It is shown that assuming a threshold of 1%, produced reconstructed images that are identical, both qualitatively and quantitatively, as compared to using the entire Jacobian, Fig. 5. Additionally, the computation time, per iteration, for the reduced method at a threshold of 1% is dramatically faster (46 seconds) as compared to conventional method (663 seconds). However, setting the threshold at a higher value of 10% results in artefacts within the reconstructed image, and does not improve reconstruction time further as compared to 1% threshold.

In order to validate the methods, measured phantom data were also used to demonstrate the applicability of more data-types (amplitude and phase) as well as larger unknown parameter set (absorption and reduced scatter). It was also demonstrated that assuming a varying threshold of 1%–10% results in reconstructed images which are similar to those using the entire Jacobian. The multi-layered phantom data is showing more robustness to higher threshold levels, since the target is cylindrical in shape and extended in Z-direction through out the domain. Furthermore, the breast model used was a more complicated case, consisting of layers as well irregular shaped target. Therefore although the experimental data has shown to be more robust to the level of assumed threshold, it should be noted that for more complex, non-homogeneous and layered problems, it is likely that a lower threshold of 1% would be more applicable.

The reduction of Jacobian size has improved the computational efficacy by reducing the overhead memory requirement and therefore computational speed. At the same time the quality of the reconstructed optical images have been shown to be uncompromised depending on the amount of reduction used. It is also important to remember that non-linear inverse problems which use Newton-type type algorithms typically require few iterations (typically between 8–15 for DOT problems) for convergence. From the above studies, it is demonstrated that usage of 1% threshold in reducing the Jacobian does not compromise the obtained image quality and can improve the computational speed by up to 14 times per iteration.

The methods demonstrated here are easily expandable to spectral imaging methods as well as methods with larger data sets. It is demonstrated that using the Jacobian reduction method outline here, not only an increase in computational speed is achieved, but the overall accuracy of the reconstructed images remain uncompromised.

Acknowledgments

This work has been sponsored by the National Cancer Institute through grants RO1CA78734, PO1CA80139, as well as the Engineering and Physical Sciences Research Council, UK. P.K.Y. acknowledges the DOD Breast Cancer predoctoral fellowship (BC050309). The authors would like to thank Dr. Shudong Jiang and Colin Carpenter for their help in conducting the phantom experiment.

References and links

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

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

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

4. S. Srinivasan, B. W. Pogue, S. Jiang, H. Dehghani, and K. D. Paulsen, “Spectrally constrained Chromophore and Scattering NIR Tomography provides quantitative and robust reconstruction,” Appl. Opt. 44, 1858–1869 (2005). [CrossRef] [PubMed]

5. S. Srinivasan, B. W. Pogue, B. Brooksby, S. Jiang, H. Dehghani, C. Kogel, S. P. Poplack, and K. D. Paulsen, “Near-infrared characterization of breast tumors in vivo using spectrally-constrained reconstruction,” Technol. Cancer Res. Treat. 5, 513–526 (2005).

6. A. Corlu, R. Choe, T. Durduran, K. Lee, M. Schweiger, S. R. Arridge, E. M. C. Hillman, and A. G. Yodh, “Diffuse optical tomography with spectral constraints and wavelength optimization,” Appl. Opt. 44, 2082–2093 (2005). [CrossRef] [PubMed]

7. 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, 12169–12174 (2007). [CrossRef] [PubMed]

8. S. Srinivasan, B. W. Pogue, S. Jiang, H. Dehghani, C. Kogel, S. Soho, J. J. Gibson, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “Interpreting hemoglobin and water concentration, oxygen saturation and scattering measured in Vivo by near-infrared breast tomography,” Proc. Natl. Acad. Sci. U. S.A. 100, 12349–12354 (2003). [CrossRef] [PubMed]

9. 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, 1128–1139 (2005). [CrossRef] [PubMed]

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

11. A. D. Klose and A. H. Hielscher, “Quasi Newton methods in optical tomographic image reconstruction,” Inverse Probl. 19, 387–409 (2003). [CrossRef]

12. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, and K. D. Paulsen, “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med. Phys. 34, 2085–2098 (2007). [CrossRef] [PubMed]

13. A. D. Klose and A. H. Hielscher, “Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer,” Med. Phys. 26, 1698–1707 (1999). [CrossRef] [PubMed]

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

15. O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl. 14, 1107–1130 (1998). [CrossRef]

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

17. S. R. Arridge and M. Schwieger, “Gradient-based optimisation scheme for optical tomography,” Opt. Express. 2, 212–226 (1998). [CrossRef]

18. 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, 262–271 (1999). [CrossRef] [PubMed]

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

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

21. A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12, 5402–5417 (2004). [CrossRef] [PubMed]

22. P. K. Yalavarthy, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Critical computational aspects of near infrared circular tomographic imaging: Analysis of measurement number, mesh resolution and reconstruction basis,” Opt. Express 14, 6113–6127 (2006). [CrossRef] [PubMed]

23. M. Guven, B. Yazici, K. Kwon, E. Giladi, and X. Intes, “Effect of discretization error and adaptive mesh generation in diffuse optical absorption imaging,” Inverse Probl. 23, 1135–1160 (2007). [CrossRef]

24. M. Molinari, B. H. Blott, S. J. Cox, and G. J. Daniell, “Optimal imaging with adaptive mesh refinement in electrical impedance tomography,” Physiol. Meas. 23, 121–128 (2002). [CrossRef] [PubMed]

25. H. Dehghani and M Soleimani, “Numerical modelling errors in electrical impedance tomography,” Phys Meas 28, S45–S55 (2007). [CrossRef]

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). [CrossRef] [PubMed]

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

28. D. W. Marquardt, “An algorithm for least squares estimation of nonlinear parameters,” J. Soc. Ind. Appl. Math 11, 431–441 (1963). [CrossRef]

29. J. R. Westlake, A handbook of numerical matrix inversion and solution of linear equations (John Wiley & Sons Inc, New York, 1968).

	Adipose	Glandular	Tumour
HbO (mM)	0.012	0.015	0.016
Hb (mM)	0.005	0.006	0.024
Water (%)	47	70	80
Scatter amplitude	1.34	0.94	0.5
Scatter power	0.56	0.79	1
µ_a at 785 nm (mm^-1)	0.0046	0.0059	0.0104
µ_s′ at 785 nm (mm^-1)	1.53	1.1381	0.6369

Threshold level	Number of nodes in forward model	Number of nodes in reconstruction	Memory usage of J ^T J	Computational time per iteration (seconds)
0%	18374	18374	2.70Gb	663
1%	18374	4414	0.16Gb	46
10%	18374	1800	0.03Gb	36

Threshol d level	Number of nodes in forward model	Number of nodes in reconstruction for µ_a	Number of nodes in reconstructi on for κ	Memory usage of J ^T J	Computation al time per iteration (seconds)
0%	17980	5520	5520	0.98 Gb	318
1%	17980	4803	5409	0.83 Gb	299
5%	17980	3552	4703	0.55 Gb	225
10%	17980	2762	2828	0.25 Gb	165

	Adipose	Glandular	Tumour
HbO (mM)	0.012	0.015	0.016
Hb (mM)	0.005	0.006	0.024
Water (%)	47	70	80
Scatter amplitude	1.34	0.94	0.5
Scatter power	0.56	0.79	1
µ_a at 785 nm (mm^-1)	0.0046	0.0059	0.0104
µ_s′ at 785 nm (mm^-1)	1.53	1.1381	0.6369

Threshold level	Number of nodes in forward model	Number of nodes in reconstruction	Memory usage of J ^T J	Computational time per iteration (seconds)
0%	18374	18374	2.70Gb	663
1%	18374	4414	0.16Gb	46
10%	18374	1800	0.03Gb	36

An efficient Jacobian reduction method for diffuse optical image reconstruction

Abstract

1. Introduction

2. Theory

3. Jacobian reduction

4. Method and results

4.1 Simulated breast model

4.2 Experimental phantom

5. Discussions and conclusions

Acknowledgments

References and links

Cited By

Figures (6)

Tables (3)

Equations (7)

Optics Express