Data subset algorithm for computationally efficient reconstruction of 3-D spectral imaging in diffuse optical tomography

Three-dimensional (3-D) models of light propagation in diffuse optical tomography provide an accurate representation of scattering in tissue. Here the use of spectral priors, shown to improve quantification of functional parameters in 2-D, has been extended to 3-D. To make 3-D spectral imaging computationally tractable, a novel technique is presented to deal with the large data set. The basic principle consists of using a dynamic criterion to select optimal data subsets that capture the major changes in the imaging domain. Results from three test cases showed comparable image quality and accuracy with less than 4% difference between the uses of data subset approach versus the entire dataset. Tested on simulated data from two different models, the algorithm was able to discern multiple objects successfully with an average error of 30% in quantifying multiple regions and less than 1% in quantifying the background. ©2006 Optical Society of America OCIS codes: (170.6960) Tomography; (170.3010) Image reconstruction technologies References and Links 1. Q. Zhang, T.J. Brukilacchio, A. Li, J.J. Stott, T. Chaves, E. Hillman, T. Wu, M. Chorlton, E. Rafferty, R.H. Moore, D.B. Kopans, and D.A. Boas, "Coregistered tomographic x-ray and optical breast imaging: initial results,” J Biomed. Opt. 10, 024033-0240339 (2005). 2. Zhu, Q., E.B. Cronin, A.A. Currier, H.S. Vine, M. Huang, N. Chen, and C. Xu, "Benign versus malignant breast masses: optical differentiation with US-guided optical imaging reconstruction,” Radiology, 237(1): p. 57-66 (2005). 3. B. Brooksby, B.W. Pogue, S. Jiang, H. Dehghani, S. Srinivasan, C. Kogel, J. Weaver, S.P. Poplack, and K.D. Paulsen, "Imaging Breast Adipose and Fibroglandular Tissue Molecular Signatures using Hybrid MRI-Guided Near-Infrared Spectral Tomography,” Proceedings of the National Academy of Sciences (in press), (2006). 4. F. F. Jobsis, "Non-invasive, infra-red monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science 198, 1264-1267 (1977). 5. J. S. Wyatt, "Cerebral oxygenation and haemodynamics in the foetus and newborn infant,” Phil. Trnas. R. Soc. Lond. B, 352, 697-700 (1997). 6. G. Strangman, D. A. Boas, and J. P. Sutton, "Non-invasive neuroimaging using near-infrared light,” Biol. Psychiatry 52, 679-693 (2002). 7. J. C. Hebden, A. Gibson, R. M. Yusof, N. Everdell, E. M. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek, J. S. Wyatt, "Three-dimensional optical tomography of the premature infant brain,” Phys. Med. Biol. 47, 4155-66 (2002). #69724 $15.00 USD Received 5 April 2006; revised 19 May 2006; accepted 25 May 2006 (C) 2006 OSA 12 June 2006 / Vol. 14, No. 12 / OPTICS EXPRESS 5394 8. B. J. Tromberg, N. Shah, R. Lanning, A. Cerussi, J. Espinoza, T. Pham, L. Svaasand, and J. Butler, "Noninvasive in vivo characterization of breast tumors using photon migration spectroscopy,” Neoplasia (New York), 2, 26-40 (2000). 9. B. W. Pogue, S. P. Poplack, T.O. McBride, W.A. Wells, O.K. S., U.L. Osterberg, and K.D. Paulsen, "Quantitative Hemoglobin Tomography with Diffuse Near-Infrared Spectroscopy: Pilot Results in the Breast,” Radiology 218, 261-6 (2001). 10. Gratton, E., Fantini, S., Franceschini, M. A., Gratton, G. and Fabiani, M., "Measurements of scattering and absorption changes in muscle and brain,” Phil. Trans. R. Soc. Lond. B 352, 727-735 (1997). 11. P. Van der Zee, M. Cope, S. R. Arridge, M. Essenpreis L. A. Potter, A. D. Edwards, J. S. Wyatt, D. C. McCormick, S. C. Roth, E. O. Reynolds, et al, "Experimentally measured optical pathlengths for the adult head, calf and forearm and the head of the newborn infant as a function of inter-optode spacing,” Adv. Expt. Med. Biol. 316, 143-53 (1992). 12. I. V. Meglinski and S. J. Matcher, "Computer simulation of the skin reflectance spectra,” Computer Methods and Programs in Biomedicine 70, 179-186 (2003). 13. D. B. Jakubowski, A.E. Cerussi, F. Bevilacqua, N. Shah, D. Hsiang, J. Butler, and B.J. Tromberg, "Monitoring neoadjuvant chemotherapy in breast cancer using quantitative diffuse optical spectroscopy: a case study,” J Biomed. Opt. 9, 230-8 (2004). 14. P. Vaupel, K. F., and O. P., "Blood Flow, Oxygen and Nutrient Supply, and Metabolic Microenvironment of Human Tumors: A Review,” Cancer Research 49, 6449-6465 (1989). 15. 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,” PNAS, 100(21): p. 12349-12354 (2003). 16. S. P. Poplack, A.N. Tosteson, M.R. Grove, W.A. Wells, and P.A. Carney, "Mammography in 53,803 women from the New Hampshire mammography network,” Radiology, 217: p. 832-840 (2000). 17. B. W. Pogue and K.D. Paulsen, "High resolution near infrared tomographic imaging simulations of rat cranium using apriori MRI structural information,” Opt. Lett. 23, 1716-8 (1998). 18. B. Brooksby, H. Dehghani, B. W. Pogue, K. D. Paulsen, "Near infrared (NIR) tomography breast image reconstruction with apriori structural information from MRI: algorithm development for reconstructing heterogeneities" IEEE J. Sel. Top. Quantum Electron. 9, 199-209 (2003). 19. M. Schweiger, S. R. Arridge, "Optical tomographic reconstruction in a complex head model using apriori region boundary information,” Phys. Med. Biol. 44, 2703-2721 (1999). 20. H. Dehghani, B.W. Pogue, J. Shudong, B. Brooksby, and K.D. Paulsen, "Three-dimensional opticaltomography: resolution in small-object imaging,” Appl. Opt. 42, 3117-3128 (2003). 21. A. Corlu, T. Durduran, R. Choe, M. Schweiger, E.M. Hillman, S.R. Arridge, and A.G. Yodh, "Uniqueness and wavelength optimization in continuous-wave multispectral diffuse optical tomography,” Opt. Lett. 28, 2339-41 (2003). 22. A. Li, Q. Zhang, J.P. Culver, E.L. Miller, and D.A. Boas, "Reconstructing chromosphere concentration images directly by continuous-wave diffuse optical tomography,” Opt. Lett. 29, 256-8 (2004). 23. 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-69 (2005). 24. J. –L. Boulnois, "Photophysical processes in recent medical laser developments: a review,” Lasers in Medical Science 1, 47-66 (1986). 25. G. M. Hale, and M.R. Querry, "Optical constants of water in the 200-nm to 200-um wavelength region,” Appl. Opt. 12, 555-563 (1973). 26. J. R. Mourant, T. Fuselier, J. Boyer, T.M. Johnson, and I.J. Bigio, "Predictions and measurements of scattering and absorption over broad wavelength ranges in tissue phantoms,” Appl. Opt. 36, 949-957 (1997). 27. H. J. van Staveren, C.J.M. Moes, J. van Marle, S.A. Prahl, and J.C. van Gemert, "Light scattering in Intralipid 10% in the wavelength range of 400-1100nm,” Appl. Opt. 30, 4507-4514 (1991). 28. A. Corlu, R. Choe, T. Durduran, K. Lee, M. Schweiger, S.R. Arridge, E.M. Hillman, and A.G. Yodh, "Diffuse optical tomography with spectral constraints and wavelength optimization,” Appl. Opt. 44, 20822093 (2005). 29. M. Schweiger, and S.R. Arridge, "Comparison of twoand three-dimensional reconstruction methods in optical tomography,” Appl. Opt. 37, 7419-7428 (1998). 30. J. C. Hebden, H. Veenstra, H. Dehghani, E.M. Hillman, M. Schweiger, S.R. Arridge, and D.T. Delpy, "Three-dimensional time-resolved optical tomography of a conical breast phantom,” Appl. Opt. 40, 32783287 (2001). 31. J. P. Culver, R. Choe, M.J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A.G. Yodh, "Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: Evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys. 30, 235-247 (2003). 32. K. D. Paulsen, P.M. Meaney, M.J. Moskowitz, and J.M. Sullivan, "A dual mesh scheme for finite element based reconstruction algorithms,” IEEE Trans Med, Imaging 14, 504-514 (1995). #69724 $15.00 USD Received 5 April 2006; revised 19 May 2006; accepted 25 May 2006 (C) 2006 OSA 12 June 2006 / Vol. 14, No. 12 / OPTICS EXPRESS 5395 33. Intes, X., S. Djeziri, Z. Ichalalene, N. Mincu, Y. Wang, P. St-Jean, F. Lesage, D. Hall, D. Boas, M. Polyzos, P. Fleiszer, and B. Mesurolle, "Time-Domain Optical Mammography SoftScan: Initial Results ,” Acad. Radiology 12, 934-947 (2005). 34. A. Ishimaru, Wave propagation and scattering in random media. Vol. 1. 1978: Academic Press, Inc., New York. 35. M. S. Patterson, B.C. Wilson, and D.R. Wyman, "The propagation of optical radiation in tissue II. Optical properties of tissues and resulting fluence distributions,” Lasers Med. Sci. 6, 379-390 (1990). 36. T. J. Farrell, M. S. Patterson, B. C. Wilson, "A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties,” Med. Phys. 19, 879-888 (1992). 37. K. D. Paulsen, and H. Jiang, "Spatially varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys. 22, 691-701 (1995). 38. S. R. Arridge, and M. Schweiger, "Image reconstruction in optical tomography,” Phil. Trans. R. Soc. Lond. B 352, 717-726 (1997). 39. M. Schweiger, S.R. Arridge, M. Hiraoka, and D.T. Delpy, "The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22, 1779-1792 (1995). 40. 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, 135145 (2003). 41. X. Wang, B.W. Pogue, S. Jiang, X. Song, K.D. Paulsen, C. Kogel, S.P. Poplack, and W.A. Wells, "Approximation of Mie scattering parameters in near-infrared tomography of normal breast tissue in vivo,” J. of Biomed. Opt. 10, 051704-1:051704-8 (2005). 42. J. R. Mourant, A.H. Hielsche


Introduction
The field of diffuse optical tomography (DOT) has gained considerable momentum over the past two decades due to its ability to provide unique information relating to tissue metabolic status.DOT derives its signal from molecular absorption signatures of dominant tissue chromophores such as oxy-hemoglobin, decoy-hemoglobin and water.Coupled with scattering mechanisms relating to cell organelle, collagen and tissue density, optical tomography offers the possibility for non-invasive functional imaging without extrinsic contrast agents.Future growth in the field is expanding quickly with coupling to MRI, Ultrasound and x-ray tomography all showing applications [1][2][3].The key to all of these functional imaging applications is the need to do multiple wavelengths imaging, such that direct chromophore and scattering reconstruction can be implemented.
The more wavelengths that are used, the better the reconstruction that can be achieved, yet there is an added computational burden in this which cannot be ignored.As more wavelengths are added, eventually the problem becomes too large to solve, and this is especially aggravated in 3-D reconstruction.This paper examines one approach to minimizing this memory requirement and making the problem of multispectral reconstruction tractable, by using datasubsets which have the maximum impact upon image reconstruction.
In DOT, near-infrared light is passed through tissue at multiple wavelengths, and using measurements of reflectance and/or transmittance along the surface/periphery of tissue and an appropriate model for light propagation, tissue optical properties can be recovered.The ability to model light propagation and de-couple the contributions of different chromophores from the absorption spectra [4] provides the key to maximal information from the imaging procedure.The functional parameters typically obtained as a result are images of oxyhemoglobin (HbO 2 ), de-ox hemoglobin(Hb) and water concentrations, and scattering, and extended indices such as total hemoglobin concentration ([Hb T ] = [HbO 2 ] + [Hb]) and oxygen saturation (S t O 2 = [HbO 2 ]/[Hb T ]×100%).Researchers have studied these functional images in brain activation studies [5][6][7], diagnosis of breast cancer [8,9], muscle [10], forearm imaging [11], skin [12] as well as in therapeutic applications [13].In the context of optical imaging for breast cancer, tumors, show a localized increase in total hemoglobin and water arising from neovascularization, separating them from healthy tissue [ 8,9].Malignancies can also exhibit an aggressive metabolic status resulting in hypoxic tissue, which can be seen as a decrease in oxygen saturation in the tumors compared to surrounding healthy tissue and benign lesions in optical images [14].Additional factors such as dense fibrosis and cell proliferation show contrast in scatter [8]; and water may give a measure of the extravascular space [15].The current clinical modality of choice is mammography and suffers from high-false positive rates and community practice has shown that approximately 80% of the women undergoing biopsy as a result of screening mammogram-detected abnormality do not have cancer [16].Overall, NIR may have a niche in the evaluation of cancers, due to its ability to complement information from other modalities such as mammography and MRI, which are structurebased.Recent studies have also indicated that DOT can be used to monitor the response of breast cancer patients to chemotherapy by following the changes observed in the functional parameters [13].
DOT images present high levels of intrinsic contrast [8,17].However due to the dominance of scattering over absorption the resulting images , show poor spatial resolution.The inverse problem leading to reconstructed optical parameters images is inherently illposed.The use of priors reduces the solution space by making it less susceptible to noise in the measured data.Previous studies have shown that priors relating to the tissue structure [17,18] or spatial information associated with the lesion [19,20] can improve image quality and accuracy of the recovered optical coefficients.Prior information considered in this work consists of the known absorbance spectrum of the chromophores as well as an empirical power-law for the scattering based on Mie theory.The use of these priors in the inverse problem leads to significant improvements in image quality and accuracy [21][22][23].In this 'spectral' approach multi-wavelength measurements are used simultaneously bypassing recovery of optical coefficients.The concentrations of oxy-hemoglobin, decoy-hemoglobin, water [24,25] and the scattering parameters [26,27] are recovered directly.This type of parameter reduction and use of constraints yields reduced cross-talk between chromophores having similar spectral behavior (e.g., HbO 2 and water) and shows a robust response to noisy measurements [23,28].
In this paper, the use of these priors has been extended from two-dimensional (2-D) studies, to three-dimensional (3-D) imaging.In a study comparing the use of 2-D and 3-D models [29] using absolute data generated from a 3-D forward model, mimicking measurements by the DOT systems, the images from a 2-D model, while showing higher accuracy in the anomaly region, showed considerable artifacts obscuring detection.The 3-D images were qualitatively superior, and showed accurate localization of inclusions.Other studies using 3-D models [20,30,31] have also shown a similar trend.Quantitatively, these studies in 3-D have shown underestimation in reconstructed parameters and illustrated the need for priors [20].Spectral priors showed a significant improvement in quantitative accuracy in 2-D and a natural extension of the approach is its use to enhance 3-D imaging.However, the use of spectral priors in 3-D is a difficult problem computationally, since the inverse problem involves use of multi-wavelength data simultaneously.In addition, the inversion recovers images of five parameters (chromophore concentrations and scatter) instead of two (optical coefficients at each wavelength) simultaneously.Use of an intermediate mesh for reducing the unknowns (governed by the number of pixels in the 3-d mesh) by interpolating from a fine mesh used for calculation of the forward model to a coarser basis for update of optical coefficients has been used for efficiency [32].However, the spectral method still requires a coarser basis than the conventional method for a comparable memory requirement.Beyond this, the main issue is the large number of measurements available at each wavelength.Fig. 1(b) shows the source-detector configuration for a slab geometry implemented by ART Inc. in a clinical system currently under multi-center clinical trials [33], and used in this study which allows up to 5000 scan-points for intensity and meantime each, at four wavelengths.Since the spectral approach uses multi-wavelength measurements simultaneously, this quickly becomes infeasible computationally in the spectral domain without a more systematic way to solve the inverse problem in stages.
Presented here is a novel approach for the spectral 3-D inverse problem, where measurements are selected preferentially based on their information content relating to the major heterogeneities in the tissue.First the measurement-points are separated into slices (see Fig. 3(c)) and using projection error, a least squares norm that gives the difference between measured data and calculated data from a discretized model on a homogeneous medium, those lines/slice of data corresponding to maximum error are then identified.The projection error is the objective functional minimized in the image reconstruction, and hence this criterion permits using data subsets, rather than all available data, to be used and yet has a maximum effect on the reconstruction.This algorithm has been tested in simulations and is applied here with frequency domain simulations of simple phantom set-ups in the geometry used by ART Inc. for a prototype breast imaging system.The results show that major changes in the tissue can be imaged successfully by using optimal data-subsets, yielding a novel computationally efficient method for 3-d spectral imaging.

Imaging system geometry
SoftScan ® is a breast imaging clinical device developed by ART Advanced Research Technologies Inc.It is a multi-wavelength transmission system acquiring time-domain signals at 690nm, 730nm, 780nm and 830nm.It uses a raster scan in step-and-shoot mode with a 3mm step size along the x and y directions.The detection geometry consists of one source and five detectors.For a source at (x,y,z)=(0mm, 0mm, 0mm) the relative detector locations are : (-25mm, 5mm, 60mm), (25mm,5mm, 60mm), (0mm, 0mm, 60mm), (-25mm, -15mm, 60mm), (25mm,-15mm, 60mm).This geometry was utilized in the current study.Fig. 1(a) shows (bottom of the slab) the spatial distribution of the detectors relative to a given source location.Fig. 1(b) illustrates the same geometry using a finite-element mesh on a slab with dimensions 96mm × 96mm × 60mm corresponding to a total of 1024 source positions (32 along y by 32 along x).On the same Fig.5120 detector locations are represented at z = 60mm.Fig. 1(c) shows the separation of the scan points into lines along the x-axis with each line containing all scan points along the y-axis.

Forward diffusion modeling
Light photons undergo absorption and scattering processes when passed through tissue and the diffusion equation approximates the bulk light propagation under the assumption that the diffuse fluence behaves as though the scattering is uniformly isotropic with a reduced scattering coefficient, ' S μ , when measured over long distances [34].This condition exists under the assumption that scatter dominates over absorption which is true in the case of (a) (b) (c) several tissue types, including the human breast, in the wavelength region of 650-1350 nm [35].This differential equation is written as [34,36] : where ( , ) r ω Φ is the isotropic fluence at modulation frequency ω and position r, ( ) is the absorption coefficient, c is the speed of light in the medium and 0 ( , ) q r ω is an isotropic source.The diffusion coefficient can be written as , where ' S μ is the reduced scattering coefficient.Equation 1 can be solved using standard numerical techniques, and here a finite element model (FEM) for equation 1, has been applied [37,38].The forward solver obtains the fluence for a given distribution of optical properties by applying suitable boundary conditions, type III (Robintype) in this case given as [39,40]: where α incorporates the reflection at the boundary due to refractive index change and n is the outer normal at the boundary at a point γ .In order to approximate differences in the models for obtaining fluence data (typically experimental) and reconstructing optical images, different finite element meshes were used for the forward and inverse problems, varying in size of elements and nodes for the particular problem.It is known that the diffusion approximation holds only for the region far from the boundary and the source(at least one scattering distance away) [34] .For this reason, the source (Gaussian) in the current application is modeled one scattering distance from the boundary and inside the domain so that the scattering in the tissue can be approximated to be isotropic everywhere.

Image reconstruction: spectral imaging
Image reconstruction is an inverse problem, where optical images are obtained using surface measurements from the tissue surface.Typically, this involves the iterative minimization of an objective function based on the difference between the measured and the model data.The reconstruction procedure used here is based on minimization of the standard sum of squared differences between the measured and calculated optical radiance at specific detector locations.This least squares error norm, called the projection error, is given by: ( ) where M is the total number of measurements at each wavelength, and m j φ and c j φ are the measured and calculated fluence at the boundary for the jth measurement point.Minimization is accomplished with a gradient-based Newton-Raphson method for iteratively updating the optical properties (starting with a homogeneous initial guess) along with Levenberg-Marquardt regularization for stabilizing the inversion.Details of the implementation can be found in earlier work [37,40].Using the recovered optical properties at different wavelengths, the images for chromophore concentrations can be derived using the modified Beer's law, where N c represents the number of chromophores, ε a their extinction coefficients and c a the associated molar concentration.Similarly the spectral behavior of the reduced scattering coefficient is mapped to the empirical power-law [26,27], where the scatter parameters A and b are physical quantities related to the microscopic nature of the turbid medium.Scattering power (b) is governed by the shape or slope of scattering which is predominantly affected by the size [41,42]; the amplitude A relates to the number density of the scatterers.Alternatively spectral priors ( 4) and ( 5) can be incorporated in the image reconstruction procedure [23,43].This allows the recovery of concentrations of the chromophores and scatter parameters directly, using multi-wavelength measurements simultaneously.The use of these priors changes the formation of the inverse problem to minimize a modified functional given by ( ) so that the sum includes measurements over N available wavelengths.The iterative update in the chromophore concentrations and scatter parameters is then governed by the equation: where ( ) , , 1: , , ℑ represent the sensitivity matrices (also known as Jacobian) of the boundary data to each of the chromophore concentrations (c) and scattering parameters (scatter amplitude A and scatter power b, given by Eq. ( 5)).The parameter α in equation 7 represents the regularization parameter and has a starting value in the range 100 to 1 and is reduced successively with iterations by a factor 4 10 found empirically from experiments to work well for the set-up.
The spectral method has been applied in 2-D cases and results have shown substantial improvements in image quality as well as quantification of the functional parameters.It is well-known that while a 3-D model provides a more accurate representation of lightpropagation, quantification of absorption and scatter in the past has shown considerable under-estimation in this model.The application of spectral priors into 3-D should extend the same benefits of reduced cross-talk and improved accuracy observed in 2-D.Computationally, the size of the inverse problem depends on the number of measurements, the size of the 3-D mesh and the number of physical parameters targeted for reconstruction.
For the ART Inc. SoftScan ® imaging geometry shown in Fig. 1(b), a full scan includes 5000 measurements each of amplitude and phase (mean-time was converted to phase) per wavelength.In the spectral method, the Jacobian ( ℑ in Eq. ( 7)) has a size of MN×CP where M is the total number of measurements including amplitude and phase, and N is the number of wavelengths, C is the number of chromophores and scatter parameters (= 5, for HbO 2 , Hb, water, A and b) and P is the size of the reconstruction basis, the unknowns.For a typical 3-d mesh containing 15,000 tetrahedral nodes (= P above) and measurements at four wavelengths (N = 4) each containing 5000 measurements for amplitude and phase separately (M = 10,000) as shown in Fig. 1(b), the size of the Jacobian becomes 40,000 x 75,000, requiring a memory of 48GB.If the pixel basis has 5000 nodes (after interpolation), the size of Jacobian is still 40,000 x 25,000, requiring 16GB for storage.Because of the memory requirements and operations associated with the inversion, the spectral method is not computational feasible, in the current manner.It is evident that a more efficient way to handle the computational requirements of the 3-d problem with spectral priors is desired.

Criterion for selecting optimal data subsets
While a large number of data measurements exist; for efficient computation, a smaller subset of data can be chosen based on some criteria relating to information content of the data.In the scheme presented here, the measurements are separated into lines, where each line or slice refers to the data collected by raster-scanning the source along a particular y-direction for a fixed x.Fig. 1(c) shows this separation; multiple lines of data are obtained for varying xposition.Depending on the position of the line relative to changes in the imaging field, its effect on the reconstruction will differ.The image reconstruction is dependent on the minimization of least squares functional of the difference between the measured data and the model data calculated initially for a homogeneous domain with approximately background properties.Hence it is reasonable to assume that if the projection error were calculated individually for each of the slices, the slice with the maximum will correspond to spatial areas of most heterogeneity and will have most effect on the images recovered.In this approach, which will be referred to as the 'data subset' algorithm, the projection error is used as a criterion to select fewer slices or subsets of data which capture the most significant changes in the imaging domain.
In the data subset method as applied to spectral imaging in 3-D, starting with an initial estimate for the imaging parameters, the forward model is solved and the projection error giving the least squares norm of the difference between the measured data and the model data, is calculated for all the slices of data individually.After sorting them in descending order, three slices of data were then chosen based on maximum projection errors in descending order, and used for image reconstruction.The number of slices was chosen arbitrarily and can be increased, as shown in Section 3.3.In a particular iteration, each of the slices is used for calculating the objective function and the unknowns (nodes of the 3-D mesh) for all five NIR parameters are updated everywhere in the 3-D volume in the direction of minimization.In this manner, each of the lines is used successively to update the parameters.The regularization is varied with iteration, and stays constant for all the slices within the iteration.However, the regularization is scaled by the maximum along the diagonal of the Hessian (= T ℑ ℑ in Eq. ( 7)) which automatically controls the step-size of the update, for each of the slices.Simulations have shown that using 3-12 lines of data (from a total of 32 lines) show reasonable recovery of images.A log file containing the slice numbers used, as well as the regularization factor and the projection error and error change with iteration is stored for each reconstruction.All simulations presented here were carried out by generating measurements of amplitude and phase using the forward solver on a dense finite element mesh with a different sized mesh for image reconstruction.A starting regularization value of 1 was found to work optimally for the simulations.The background values were used as the initial guess and increasing the sampling of the domain by choosing more slices based on the projection error criterion, accommodates errors in initial estimates of the parameters (values up to 10% off could yield similar images and quantification).The stopping criterion for image reconstruction was chosen empirically as the change in projection error for a particular slice to be less than 1% of the previous iteration value.Once this criterion is reached for a particular slice, it is no longer considered in the minimization procedure, for successive iterations.

Results
The results presented here are divided into three sections.The first section deals with image reconstruction using the data subset approach on forward data from an analytic model.Results from comparison of this approach, with use of all available data are presented in Section 2 using three test phantoms.The behavior of the data subset method on a complex phantom geometry is analyzed in Section 3.

Reconstruction using data subset algorithm
In order to test the data subset method for spectral image reconstruction, time-domain simulated data were used.Amplitude and mean-time measurements were obtained from an analytical model [33,[44][45][46] on a test phantom with dimensions of 96x96x60mm with a single inclusion of size 15×15×20mm centered at (30,30,30)mm.The simulated inclusion has a 4:1 contrast in HbO 2 with background concentrations of HbO 2 = 14 µM, Hb = 6 µM, Water = 30%, A = 0.8638 and b = 0.6.The amplitude data was assumed to be at 100MHz, since this was found to be nearly equivalent to the steady state intensity.The mean-time was converted to phase to be compatible with the finite element model used, developed in the frequency domain.Phase was calculated using the relationship: where f is the frequency (100MHz) and τ is the mean-time.The calibrated log amplitude and phase at 785nm are shown in Figs.2(a) and (b) for the test phantom with inclusion, compared with data generated using the FEM model on a homogeneous domain with identical background concentrations.
The difference between the data-sets or the residual data is illustrated in 2(c) for log amplitude and 2(d) for phase.Figure 2 shows that the general shape and magnitude of the measurements between the two models are comparable.The residual for log Amplitude shows a Gaussian-shape, with center corresponding to the location of the HbO 2 inclusion.The residual of the phase does not show this trend, possibly because phase predominantly contains information relating to scattering, and the inclusion does not contain a contrast in scatter.
For the purpose of image reconstruction, measurement points below a certain threshold (log of amplitude < -16) were considered to be below the noise level and removed from the data before reconstruction.Based on the projection error calculated using measurements shown in Fig. 2, after separation of the measurement points into slices, the projection errors are plotted for all slices in Fig. 3(a).The criterion gave a maximum at slice 10, corresponding to position of x = 30mm, since each slice is with a resolution of 3mm along x-axis.This corresponds well with the location of the inclusion.Three slices of data with corresponding magnitudes of the projection errors in decreasing order were chosen and image reconstruction was performed using the data subset method for spectral imaging.A finite element rectangular mesh of corresponding domain size with 56,753 tetrahedral elements and 11,151 nodes was used for reconstruction with second pixel basis of size 15×15×8 for interpolation of the Jacobian before inversion.The initial estimate for the parameters was assumed to have an average of 6% error relative to the true background concentrations.A starting value of 1 was used for regularization, found to be optimal for the slice approach.The reconstruction was terminated when the projection error change between successive iterations for each slice was less than 1% and the algorithm converged in 17 iterations.Images for all five NIR parameters were recovered and Fig. 3(c) shows the recovered image for HbO 2 compared with the true image in Fig. 3(b).The inclusion is clearly visible in the HbO 2 image in the expected position.The diffuse nature of the recovered image is likely due to dominance of scattering.While object is discernible qualitatively, the quantitative accuracy is limited to 30% of expected contrast.All other NIR parameter images were found to be homogeneous as expected with average values within 1.25% of the true values and standard deviation within 2.5% of the mean.

Comparison of data subset approach with regular method
In order to compare the results from using the data subset approach with the use of entire available data (called 'regular' method here), three test simulations were carried out on a phantom with varying inclusion size.The phantom is of size 60×60×60mm, a smaller size than conventionally used, to make the regular method computationally feasible.A spherical inclusion was centered at (20, 20, 20mm) with diameter varying as 10, 15, 20mm for the three test cases.The background had chromophore concentrations of 14µM HbO 2 , 8 µM Hb, 30% water, A = 1.0 and b = 1.0, based on typical concentrations found in normal breast tissue [15].In order to mimic a typical lesion, the inclusion in all three cases, had nearly 3:1 contrast in HbO 2 (containing 40µM), 2:1 contrast in Hb and scatter amplitude and 60% water.There was 20% contrast in scatter power.The amplitude and phase measurements were generated with the source raster-scanned in 3mm resolution along x and y axis giving 21 source positions along each direction.The forward solver was used with a finite element mesh containing 70,487 tetrahedral elements and 13,946 nodes.A finite element mesh containing 43,182 tetrahedral elements and 8,812 nodes was used for image reconstruction with a pixel basis of 15x15x8 for interpolation.The initial homogeneous estimates for all parameters used in the reconstruction, were set at the background values.
The regular method and the data subset approach were carried out on the generated data.For the latter, the projection error for different lines of data was computed, where each line was at a different value along x-axis.Six slices in the descending order of projection error were used for the slice method for image reconstruction.Both methods used the same meshes and reconstruction basis and value for starting regularization (=1).The stopping criterion was when the change in projection error was less than 1% between successive iterations.The regular method converged in 9 to 11 iterations for the three test cases and the data subset approach involved 11 to13 iterations.In order to illustrate a visual comparison of the images recovered by both method, the HbO2 and Hb images recovered for the second test case (diameter of inclusion = 15mm) are displayed in Fig. 4. Fig. 4 (a) shows a cross-section of the plane containing the center of the inclusion (labeled as region 1).Fig. 4(b) shows the same cross-section of the reconstructed images for HbO 2 using the data subset approach (on the left) and regular method (on the right) respectively.Fig. 4(c) shows the corresponding images for Hb.The image for HbO 2 is comparable qualitatively in both cases.In the case of Hb, a pseudo artifact is visible in the top right corner of the image from the regular method using all data, (not present in the true crosssection) and is reduced in magnitude using the data subset technique.This is possibly because of the localized use of measurements corresponding to the position of the inclusion, which reduces background heterogeneity.The images for water scatter amplitude and power (not shown here) were comparable with both methods.In order to compare the quantitative values recovered, the maximum in the inclusion is plotted for both methods for each of the three test cases, in Figs.5(a), (b) and (c).The chromophore concentrations and scatter parameters show comparable values, indicating that use of reduced measurements in the form of subsets, specifically corresponding to areas of maximum heterogeneity successfully capture the imaging domain in the reconstruction procedure.As expected, the quantification improves as the size of the inclusion increases.For all three cases, the regular approach appears slightly better in recovering scatter parameters, but for the 20mm case, the use of a subset of the data shows better quantification of HbO 2 and water.On an average, the percent difference between the two methods was less than 4% for all three test cases.In terms of reconstruction time, the data subset method converged faster taking less than 45% (on an average over the three test phantoms)of the time required for regular approach.The domain used here is smaller than typical domains used in this geometry, which was necessary for using the regular method which was infeasible in larger domains of sizes typically observed in the clinical setting.In future simulations involving larger imaging domains, only the results from the data subset method are shown.Fig. 5. Quantitative comparison of the data subset method with use of entire available dataset (regular method) for three test cases having dimensions 60×60×60 mm with single inclusion whose size varies as in (a) 10mm (test 1) (b) 15mm (test 2) and (c) 20mm (test 3).The maximum in the region of interest (ROI) is plotted for each of the five NIR parameters.The difference between the two methods was less than 4% overall for all three cases.

Complex test phantom results
Previously, the image reconstruction was applied to test phantoms with a single anomaly.To characterize the behavior of the data subset approach further, a more complicated test phantom having three anomalies, was used to generate amplitude and phase data.The phantom had dimensions of 96×96×60mm, with three spherical inclusions of diameter 20mm.The first inclusion was located at the center of the phantom at (48, 48, 30 mm) and had a contrast in all five NIR parameters with respect to the background.The second inclusion was off-center and closer to the source side, centered at (24, 24, 15mm) and had contrast only in Hb and water.The third inclusion was also off-center and closer to the source side along zdirection centered at (72, 24, 15 mm), and had contrast only in the scatter parameters A and b.This set-up of heterogeneities allows the examination of cross-talk between HbO 2 and water, as well as Hb and scatter, the parameters most susceptible because of their similar spectral shapes.The actual concentrations and contrast in parameters are shown in Table 1.
Amplitude and phase measurements were generated on a fine mesh having 117,000 tetrahedral elements with 22,715 nodes and a second mesh having 56,753 elements and 11,151 nodes was used for image reconstruction with a pixel basis of 15x15x8 for interpolation.The source positions were divided into 32 slices, for varying x-position in 3 mm resolution.Each slice had the source raster-scanned along 32 positions for y-direction, with five detector positions for each source, giving 160 measurement points for a single slice.The true background properties were given as the initial estimate.The projection error was calculated for all slices and sorted in descending order of magnitude.The images reconstructed using 12 slices are shown in Fig. 6.The algorithm converged in 6 iterations.The cross-sections from the true phantom along two planes, depicts the inclusions labeled as regions 1, 2 and 3 in Fig. 6(a) with associated concentrations as in Table 1.The reconstructed image for HbO 2 in Fig. 6(b) shows the central inclusion clearly in the central cross-section as well as in the second cross-section with no visible cross-talk from the other NIR parameters having contrasts in regions 2 and 3. Two inclusions are expected in Hb (Regions 1 and 2) and the cross-sections of the reconstructed Hb image in Fig. 6(c) show both of these inclusions in the expected positions.The central inclusion is well-localized and is not visible on the top of the imaging domain (not shown here).The peak in region 2 is higher in second cross-section as expected, since it is centered lower in z-direction, compared to region 1.The image for water was found to mostly homogeneous (not shown here) with region 1 faintly visible.Regions 1 and 3 have a contrast in scatter amplitude and in the reconstructed scatter amplitude image in Figs.6(d), both inclusions are clearly visible in both cross-sections.In addition inclusion labeled region 2 is also visible due to cross-talk from Hb.In order to quantify the chromophore concentrations and background values for the different regions, the maximum for each parameter in the three regions as well as the average in the background was calculated and the difference between the two is plotted for each region in Fig. 7(a).The HbO 2 increase in region 1 is evident quantitatively, when compared to regions 2 and 3. Contrasts in Hb appear as comparable increases in regions 1 and 2 compared to region 3, which did not posses a contrast in this parameter.Water shows a small increase in region 1, while the quantification in region 2 is comparable to region 3. Scatter amplitude shows an increase in all three regions over the background.Regions 1 and 3 show expected trend, the increase in region 2 is possibly due to cross-talk from de-ox hemoglobin.Scatter power shows a small increase in all three regions over the background, showing similar behavior as scatter amplitude.The mean error in quantifying the parameters with contrast was 42.5% and mean error over all three regions and five parameters was 30%, with the minimum being in the quantification of region 3 (< 15%) containing scatter contrast.Further, for this complex phantom set-up, increasing the number of slices to be used for the image reconstruction was examined by reconstructing using 6, 9 and 12 slices of measurements.Increasing the number of slices will provide more uniform sampling, but increase the computational burden.12 slices was assumed to provide sufficient sampling of the heterogeneities in the imaging domain.Using the projection error criterion, the slices were sorted in descending order and different number of slices (6, 9 and 12) was chosen for independent spectral reconstructions with all other reconstruction parameters kept constant.To illustrate the trend in quantification, the difference between the maximum of each parameter in region 1 (possessing a contrast in all five NIR parameters) and the corresponding average in the background are plotted in Fig. 7(b).The quantification of HbO 2 and water clearly improve with increasing number of slices.The recovered values for Hb are comparable for 9 and 12 slices.For scatter parameters, the use of 9 slices shows the best quantification.The mean error over all three regions and five NIR parameters was less than 33%.The quantification of the background was accurate with less 1% mean error overall.The difference between the maximum in region of interest (ROI) and the average in the background is plotted for each region and each parameter, from reconstructed images shown in Fig. 6 using 12 lines of data.HbO 2 shows an increase in region 1 higher than regions 2 and 3; Hb shows comparable increases in regions 1 and 2 (as expected).Scatter amplitude and power show contrasts in all three regions, with maximum in region 3. (b) Difference between maximum in region 1 and average in background is plotted for all five NIR parameters from images reconstructed using 6, 9 and 12 lines of data separately.Quantification of HbO 2 and water improves with increasing data subsets.

Discussion
Three dimensional imaging in the context of reconstructing optical properties has clearly shown benefits over 2-D imaging [29], but there is a strong need for priors to improve quantitative accuracy [20].The use of spatial priors improved quantitative accuracy substantially in 3-D by reducing the size of the problem [20].However, the sensitivity of the inverse problem to errors in the spatial priors are yet to be investigated in detail and 2-D studies indicate that erroneous perturbations in the spatial priors introduce severe bias in the reconstructed images [47].Spectral priors offer physiological constraints based on expected behavior of absorbers and scatterers with wavelength and with reasonable assumptions regarding the composition of tissue, they extend the benefits of a parametric-type reconstruction [48] to DOT.Researchers have shown using different types of DOT measurements (CW and frequency domain) that spectral reconstructions were robust to higher levels of noise in measurements [23], possessed lesser cross-talk between parameters [22,28], and were quantitatively superior to spatial priors [49].Corlu et al. [28] showed a 3-D application of spectral imaging in clinical data, however their tests regarding the behavior of spectral priors were limited to 2-D situations.Presented here is one of the preliminary studies in the application of spectral priors in 3-D, with a novel technique using an optimal data subset rather than an entire data-set for computational feasibility.
As 3-D imaging advances, research into alternate techniques for solving the inverse problem, preserving qualitative advantages and improving accuracy while placing a premium on efficiency and memory requirements come into place.Arridge et al. [50] presented an alternate approach of using conjugate gradient method instead of Newton's method for inversion to reduce computational costs.We are exploring methods that can be incorporated into our existing solver based on Newton's method, which has been validated in studies previously [37], hence the need for data subset approach.Schweiger et al. [51] presented an alternate formulation for the Hessian in the Gauss-Newton setting, in order to exclude the explicit computation and storage of the Hessian.Demonstrating in the framework of optical properties recovery without spectral priors, they also acknowledged that the formulation was more efficient when M<P/2 where M is the total number of measurements and P is the total number of unknowns in the reconstruction.The current application is an unlikely scenario for application of this algorithm due to the large dataset available.Dierkes et al. [52] implemented a fast Fourier space method for paraxial data to reduce the computational burden from a large data set.However their method was in the framework of linear reconstruction using Rytov approximation.Eppstein et al. [53] used a Bayesian approach to 3-D image reconstruction employing a domain decomposition scheme where the entire 3-D volume was assumed to comprise of non-overlapping smaller sub-domains.The complete domain was used for the forward model, but inversion was carried out in the sub-domains successively thus recovering the entire volume.A sub-zone based approach has also been implemented in Magnetic Resonance Elastography (MRE) studies [54,55] where model-based optimization was performed on small overlapping sub-zones of the total tissue region of interest that are processed in an hierarchical order determined by progressive error minimization.The objective function was recast as the sum of minimizations and results showed recovery of high quality distribution images and robust behavior to noise supported by a significant reduction in processing time.Doyley et al. [56] have presented clinical studies using the algorithm and demonstrate that computational time is linearly related to the number of subzones used during image reconstruction, and accuracy is comparable to full-volume studies if the domain is large enough to negate effects of boundary conditions.
In the algorithm presented here, a similar and yet distinct approach of separating measurements and replacing the objective function to contain a sum of minimizations of optimal lines of data has been offered.At each step of minimization, the global image is updated, instead of sub-zones.The primary computational issue handled here is the large measurement data set; and the projection error criterion offers a potential scheme for selecting optimal data subsets.The algorithm can also handle the minimization of all available lines of data but the use of the particular slices makes the reconstruction efficient by identifying key changes in the imaging domain and reducing redundant information.This translates to a balance of capturing the dominant changes in tissue such as a high contrast tumor, at the cost of smaller less absorbing in-homogeneities, depending on the number of slices chosen.The algorithm can be implemented for any raster-scanning geometry.It can also be used for systems involving multi-plane measurements, so that data from planes closest to the tumor can be used for image reconstruction.
Section 1 dealt with calibration and image reconstruction with the data subset approach, using time domain measurements of mean-time and intensity from an analytical model.After suitable calibration, the projection error criterion when applied on the data, indicated that the optimal subsets to be used were along the lines closest to the inclusion.The reconstructed images using three slices of data and the 3-D spectral method recovered the HbO 2 anomaly in the location expected (Fig. 3(c)), with the other images remaining homogeneous.Quantitatively, only 30% of the true contrast could be recovered in the inclusion, though the background values were accurate with less than 1.25% error.However, this under-estimation may be due to the raster scanning imaging geometry which limits the angular sampling [57], rather than the 3-D model.A comparison of the data subset method with use of the entire available dataset (regular method) was analyzed using three test cases in Section 2. The results demonstrated that the images recovered with both methods were comparable both qualitatively and quantitatively.The data subset approach showed suppression of pseudo artifacts in the background of Hb image (see Fig. 4(d) and (e)), possibly due to the selection of the optimal measurements.As the size of the inclusion increased, the quantification improved and for the test case having the 20mm inclusion, the data subset approach showed slightly higher accuracy in HbO 2 and water compared to the regular method.The difference between the two methods was quantified using the maximum in the region of interest, to be less than 4% overall in all three test cases.When applied on a complex test phantom, the use of selected optimal data subsets could recover the multiple anomalies with differing contrasts.HbO 2 contrast was expected only in the central anomaly and the image in Fig. 5 reconstructed this with a localized increase; the absence of cross-talk from the other parameters is encouraging.The reconstructed Hb image showed two regions of contrast in agreement with the expected locations.The scatter images showed accurate recovery in anomaly regions 1 and 3, also showing some cross-talk from Hb in region 2. The cross-talk occurs because the spectrum of Hb is similar in shape to Mie scattering spectra and this can be minimized by the availability of shorter wavelengths.A comparison of the accuracy in the reconstructed images, with different lines of data (6.9 and 12 lines) showed that quantification of HbO 2 and water improved with increasing number of slices.

Conclusions
Overall, use of the projection error criterion in selecting optimal data subsets has made 3-D spectral imaging computationally tractable allowing direct reconstruction of chromophore concentrations and scattering without significant cost on image quality and accuracy in characterization of the main heterogeneities in the imaging domain.This will be employed in future studies to analyze effects of increasing background heterogeneity and applied on clinical data.These spectral priors in 3-D can be coupled with spatial priors relating to structure of the tissue to obtain the best possible images from DOT.

Fig. 2 .Fig. 3 .
Fig. 2. (a) Comparison of log of calibrated amplitude data from analytical model versus forward data from finite element model.The analytical model data was generated on a test phantom with a single anomaly containing 4:1 contrast in HbO2 and the FEM data was generated on a homogeneous medium with identical background properties as the former.(b) Same as (a) for phase data, (c) the difference in the log amplitude data between the two models shown in (a), also the residual to be minimized in the image reconstruction (d) the residual for phase.

Fig. 4 .
Fig. 4. (a) Cross-section showing the location of a 15mm spherical inclusion in a test phantom(size: 60mm cube), with region label 1 compared to 0 in the background.The crosssection is along the center of inclusion.Frequency domain measurements were generated and 3-D spectral images were recovered.(b) reconstructed HbO 2 image (in mM) using 6 lines of optimal data subsets (left) and using entire available dataset (regular method) on the right (c) same as (b) for Hb.The images are comparable using both methods, though the data-subset method shows reduced artifacts in the background.

Fig. 6 .
Fig. 6.(a) two cross-sections of a test phantom showing three inclusions labeled with different region numbers (I cross-section is along center of region 1 and II cross-section is along plane containing centers of regions 2 and 3).Region 1 has contrast in all five NIR parameters.Region 2 has contrast only in Hb and water and Region 3 has contrast only in the scatter parameters (see Table 1 for actual concentrations).(b) Cross-sections from reconstructed HbO 2 image (in mM): the single central anomaly is visible (c) cross-sections from the Hb image (in mM) showing contrasts in the two regions expected and (d) cross-sections from reconstructed scatter amplitude Both scatter amplitude and power (not shown here) show contrast in regions 1 and 3 (expected) as well as in region 2, due to cross-talk from Hb.

Fig. 7 .
Fig. 7. (a)  The difference between the maximum in region of interest (ROI) and the average in the background is plotted for each region and each parameter, from reconstructed images shown in Fig.6using 12 lines of data.HbO 2 shows an increase in region 1 higher than regions 2 and 3; Hb shows comparable increases in regions 1 and 2 (as expected).Scatter amplitude and power show contrasts in all three regions, with maximum in region 3. (b) Difference between maximum in region 1 and average in background is plotted for all five NIR parameters from images reconstructed using 6, 9 and 12 lines of data separately.Quantification of HbO 2 and water improves with increasing data subsets.

Table 1 .
True concentrations for HbO2, Hb, water, scatter amplitude and power in the background and contrasts with respect to background in each region, for the test phantom shown in Fig.6.