Structural information within regularization matrices improves near infrared diffuse optical tomography

Near-Infrared (NIR) tomographic image reconstruction is a non-linear, ill-posed and ill-conditioned problem, and so in this study, different ways of penalizing the objective function with struc ural information were investigated. A simple framework to incorporate struc tural priors is presented, using simple weight matrices that have either Laplacian or Helmholtz-type structures. Using both MRI-derived breast geometry and phantom data, a systematic and quantitative comparison was performed with and without spatial priors. The Helmholtz-type struct ure can be seen as a more generalized approach for incorporating spatial pr ors into the reconstruction scheme. Moreover, parameter reduction (i. e. hard prior information) in the imaging field through the enforcement of spatially explicit regions may lead to erroneous results with imperfe ct spatial priors. © 2007 Optical Society of America OCIS codes: (170.0170) Medical optics and biotechnology; (170.0110) I maging systems; (170.3010) Image reconstruction techniques; (170.3880) Me dical and biological imaging; (170.6960) Tomography; (100.3190) Inverse problems References and links 1. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kil mer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Sig. Proc. Mag. 18, 57–75 (2001). 2. B. W. Pogue, M. Testorf, T. McBride, U. Osterberg, and K. D. Paulsen, “Instrumentation and design of a frequency-domain diffuse optical tomography imager for breas t c ncer detection,” Opt. Express 1, 391–403 (1997). 3. 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 Concen tration, Oxygen Saturation and Scattering Measured in Vivo by Near-Infrared Breast Tomography,” Proc. Nat. Acad. Sci. U.S.A.100, 12349–12354 (2003). 4. Q. Zhu, N. G. Chen, and S. C. Kurtzman, “Imaging tumor angiogen esis by use of combined near-infrared diffusive light and ultrasound,” Opt. Lett. 28, 337–339 (2003). 5. 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 tomograp hic x-ray and optical breast imaging: initial results,” J Biomed. Opt. 10, 024033:1–9 (2005). 6. B. Brooksby, S. Jiang, C. Kogel, M. Doyley, H. Dehghani, J. B. Weaver, S. P. Poplack, B. W Pogue, and K. D Paulsen, “Magnetic resonance guided near infrared tomograp hy of the breast,” Rev. Sci. Instrum. 75, 5262–5270 (2004). 7. A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chav es, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical br east imaging guided by three-dimensional mammography,” Appl. Opt. 42, 5181–5190 (2003). #81441 $15.00 USD Received 26 Mar 2007; revised 11 Jun 2007; accepted 11 Jun 2007; published 13 Jun 2007 (C) 2007 OSA 25 June 2007 / Vol. 15, No. 13 / OPTICS EXPRESS 8043 8. A. Li, G. Boverman, Y. Zhang, D. Brooks, E. L. Miller, M. E. Ki lmer, Q. Zhang, E. M. C. Hillman, and D. A. Boas, “Optimal linear inverse solution with multiple prior s in diffuse optical tomography,” Appl. Opt. 44, 1948–1956 (2005). 9. V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “C oncurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc. Nat. Aca d. Sci. U.S.A.97, 2767-2772 (2000). 10. G. Boverman, E. L. Miller, A. Li, Q. Zhang, T. Chaves, D. H. B rooks, and D. Boas, “Quantitative spectroscopic diffuse optical tomography of the breast guided by imperfect a priori structural information,” Phys. Med. Biol. 50, 3941–3956 (2005). 11. R. L. Barbour, H. L. Graber, J. W. Chang, S. L. S. Barbour, P . C. Koo, R. Aronson, “MRI-guided optical tomography: Prospects and computation for a new imaging method,” IEE E Comp. Sci. Eng. 2, 63–77 (1995). 12. A. H. Hielscher and S. Bartel, “Use of penalty terms in grad ient-based iterative reconstruction schemes for optical tomography,” J. Biomed. Opt. 6, 183-192 (2001). 13. V. Ntziachristos, X. H. Ma, and B. Chance, “Time-correlat d single photon counting imager for simultaneous magnetic resonance and near-infrared mammography,” Rev. Sci. I nstrum.69, 4221–4233 (1998). 14. B. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “ Near infrared (NIR) tomography breast image reconstruction witha priori structural information from MRI: algorithm development for r econstructing heterogeneities,” IEEE J. Sel. Top. Quantum Electron. 9, 199–209 (2003). 15. B. Brooksby, S. Jiang, H. Dehghani, B. W. Pogue, K. D. Paul sen, J. Weaver, C. Kogel and S. P. Poplack, “Combining near infrared tomography and magnetic resonance imagin g to study in vivo breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic r sonance structure,” J. Biomed. Opt. 10, 051504:1–10 (2005). 16. B. Brooksby, B. W. Pogue, S. Jiang, H. Dehghani, S. Sriniv asan, C. Kogel, T. D. Tosteson, J. Weaver, S. P. Poplack and K. D. Paulsen, “Imaging breast adipose and fibrogl andular tissue molecular signatures using Hybrid MRI-Guided Near-Infrared Spectral Tomography,” Proc. Nat. Acad. Sci. U.S.A.103, 8828–8833 (2006). 17. M. F. Ernst and J. A. Roukema, “Diagnosis of non-palpable b reast cancer: a review,” The Breast 11, 13 (2002). 18. H. Dehghani, B. W. Pogue, S. P. Poplack, and K. D. Paulsen, “Multiwavelength three-dimensional near-infrared tomography of the breast: initial simulation, phantom, and cli ni al results,” Appl. Opt. 42, 135–145 (2003). 19. X. Wang, B.W. Pogue, S. Jiang, X. Song, K.D. Paulsen, C. Ko gel, S.P. Poplack, and W.A. Wells, , “Approximation of Mie scattering parameters in near-infrared tomography of n ormal breast tissue in vivo,” J. Biomed. Opt. 10, 051704:1–8 (2005). 20. S. R. Arridge, “Optical tomography in medical imaging,” Inv . Problems15, R41–R93 (1999). 21. S. Srinivasan, B. W. Pogue, H. Dehghani, S. Jiang, X. Song , a d K. D. Paulsen, “Improved quantification of small objects in near-infrared diffuse optical tomography,” J. Biomed. Opt. 9, 1161–1171 (2004). 22. H. Dehghani, B. W. Pogue, J. Shudong, B. Brooksby, and K. D . Paulsen, “Three-dimensional optical tomography: resolution in small-object imaging,” Appl. Opt. 42, 3117–3126 (2003). 23. B. Brooksby, S. Srinivasan, S. Jiang, H. Dehghani, B. W. P ogue, K. D. Paulsen, J. Weaver, C. Kogel, and S. P. Poplack, “Spectral priors improve near-infrared diffuse tomography more than spatial priors,” Opt. Lett. 30, 1968–1970 (2005). 24. H. Jiang, K. D. Paulsen, and U. Osterberg, B. W. Pogue, and M. S. Patterson, “Optical image reconstruction using frequency domain data: simulations and experiments,” J. Opt. Soc. Am. A13, 253–266 (1996). 25. M. Schweiger, S. R. Arridge, M. Hiroaka, and D. T. Delpy, “ The Finite Element Model for the Propagation of Light in Scattering Media: Boundary and Source Conditions, ” Med. Phys.22, 1779–1792 (1995). 26. K. D. Paulsen and H. Jiang, “Spatially varying optical pr operty reconstruction using a finite element diffusion equation approximation,” Med. Phys. 22, 691-701 (1995). 27. A. N. Tikhonov, “Regularization of mathematically incorr ectly posed problems,” Soviet Math. Dokl. 4, 16241627 (1963). 28. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, and K. D. Pauls en, “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med. Phys. 34(6), 2085–2098 (2007). 29. K. Levenberg, “A method for the solution of certain nonlin ear problems in least squares,” Q. Appl. Math. 2, 164-168 (1944). 30. D. W. Marquardt, “An algorithm for least squares estimati on of nonlinear parameters,” J. Soc. Ind. Appl. Math. 11, 431-441 (1963). 31. H. Jiang, K. D. Paulsen, U. Osterberg, B. W. Pogue, and M. S . Patterson, “Simultaneous reconstruction of optical absorption and scattering maps in turbid media from near-infr ared frequency-domain data,” Opt. Lett. 20, 21282130 (1995). 32. B. W. Pogue, K. D. Paulsen, H. Kaufman, and C. Abele, “Calib r t on of near-infrared frequency-domain tissue spectroscopy for absolute absorption coefficient quantita tion in neonatal head-simulating phantoms,” J. Biomed. Opt.5, 185-193 (2000). 33. S. Jiang, B. W. Pogue, T. O. McBride, M. M. Doyley, S. P. Pop lack, and K. D. Paulsen, “Near-infrared breast tomography calibration with optoelastic tissue simulating p hantoms,” J. Electron. Imaging 12, 613–620 (2003). 34. D. R. Lynch,Numerical partial differential equations for environment al scientists and engineers – A first practi#81441 $15.00 USD Received 26 Mar 2007; revised 11 Jun 2007; accepted 11 Jun 2007; published 13 Jun 2007 (C) 2007 OSA 25 June 2007 / Vol. 15, No. 13 / OPTICS EXPRESS 8044 cal course, (Springer, Edition 1, 2005). 35. T. O. Mcbride, B. W. Pogue, S. Jiang, U. L. Osterberg, and K . D. Paulsen, “A parallel-detection frequencydomain near-infrared tomography system for hemoglobin imaging of the breast in vivo,” Rev. Sci. Instrum. 72, 1817–1824 (2001). 36. P. K. Yalavarthy, H. Dehghani, B. W. Pogue, and K. D. Pauls en, “Critical computational aspects of near infrared circular tomographic imaging: Analysis of measurement number, me sh resolution and reconstruction basis,” Opt. Express14, 6113–6127 (2006). 37. B. Brooksby,Combining near infrared tomography and magnetic resonance imaging to improve breast tissue chromophore and scattering assessment , Ph.D. Thesis, Dartmouth College, May 2005. 38. S. Srinivasan, B. W. Pogue, H. Dehghani, S. Jiang, X. Song , a d K. D. Paulsen, “ Effect of image reconstruction bias upon spectroscopy-based quantification of chromophore s in near infrared tomography,” Proc. OSA Biomedical Topical Meetings, OSA Technical Digest, WB3:1-3, Optic al Society of America, Washington, DC (2004).


Introduction
Near Infrared (NIR) optical tomography is a method that uses light in the 650-900nm wavelength range to recover images of the internal spatial distribution of tissue optical properties, absorption (or chromophore concentrations) and scattering parameters [1][2][3].When imaged at multiple wavelengths, these quantities can be used to estimate tissue hemoglobin and water concentration [3] and have the advantage of being acquired non-invasively and without ionizing radiation.The imaging procedure can be rapidly or repeatedly applied to investigate physiological state, and systems can be integrated into conventional imaging platforms such as X-ray mammography, Ultrasound and MRI [4][5][6][7][8][9].These hybrid systems have been shown to achieve superior performance in terms of resolution and quantitative accuracy which should provide more accurate physiological data from the tissue under investigation [4][5][6][7][8][9][10][11][12][13][14][15][16].However, a fundamental question is how to utilize the spatial information from the clinical system to maximize the accuracy of NIR tomography.In this study, the ability to improve the quantitative accuracy of regions imaged with NIR tomography was investigated, in the setting where prior spatial information is readily available.This work explores image reconstruction strategies that take advantage of multimodality image data, in particular, the combination of MRI with NIR optical tomography for breast cancer imaging.MRI provides structural information at high spatial resolution (near 1 mm), whereas NIR imaging has relatively poor spatial resolution (near 4-7 mm).Yet MR imaging would benefit from the molecular-specific signatures available through NIR [3,[16][17][18], specifically tissue hemoglobin content, oxygenation level, and water, as well as scattering particle size and number density [19].The inverse problem (image reconstruction procedure) in NIR imaging is known to be a non-linear, ill-posed and ill-conditioned problem [20].Use of structural information in NIR reconstruction schemes has been explored by several research groups.For example, Li et.al [7] have used the data derived from X-ray mammography for choosing different regularization parameters for the region of interest (ROI) and surrounding tissue, and have shown that the contrast and resolution of the reconstructed images can be improved.Srinivasan et al. [21] have developed a three-step reconstruction process for improving the quantification accuracy of small-objects in NIR tomography, where they use the conventional NIR reconstructed images (first step) as a structural prior for the last two steps.Earlier papers have shown that optical contrast can be correlated to MR contrast [6,9,13] and that structural MRI images can be used to reduce the number of unknown parameters to be estimated [14].The difficulty with parameter reduction approaches (referred to as hard-priors [22]) is the potential of introducing error by imposing incorrect model assumptions and introducing variation due to uncertainty in the prior information (even when the underlying model is appropriate).For example, the features which lead to contrast in one imaging system may not be spatially coregistered with those that produce contrast in another imaging system.Further, segmentation of congruent features always includes classification errors resulting from digitization.Recently, Boverman et.al [10] showed that even imperfect priors which encode breast background structure improves anomaly localization, but at the expense of biasing the spectroscopic dimension of the image reconstruction.Another type of approach for constraining the problem, often described as soft-priors, penalizes the variation within regions which are assumed to have the same properties by controlling regularization.Brooksby et.al [15,16,23] have developed a Laplacian type of regularization that allows intra-region variability.This is a method which works well even if the confidence in the prior structural information is low.
This paper develops a more generalized framework, known as soft-priors, for incorporating the structural priors into the NIR image reconstruction process, and explores a covariance-based constraint scheme adopted from finite differencing of the Helmholtz equation in addition to the soft and hard prior approaches noted above.The soft-prior approach allows optical property variation with a given region, reducing biases caused by the use of imperfect prior information.The results indicate that imperfect structural information can generate errors in the hard-priors case, whereas the soft-priors are able to quantify regions more appropriately.Simulation and experimental studies are performed to demonstrate the superior reconstruction image quality.These types of procedures are needed to improve NIR imaging both in terms of high spatial resolution available from MRI and high contrast inherent in the NIR signal.

Diffusion-based Light Transport Model
Light transport in breast tissue can be described accurately by the Diffusion Equation (DE), which is an approximation to the Radiative Transport Equation (RTE) [24].In the frequencydomain, the DE is given by where Φ(r, ω) is the photon density at position r and light modulation frequency is given by ω (in this work, ω = 100 MHz).The isotropic source term is represented by Q o (r, ω) and speed of light in tissue by c. µ a (r) is the optical absorption coefficient and D(r) is the optical diffusion coefficient, which is defined as where µ ′ s (r) is the reduced scattering coefficient, which is defined as µ ′ s = µ s (1 − g).µ s is the scattering coefficient and g is the anisotropy factor.A Robin-type (Type-III) boundary condition is applied to exactly model the refractive-index mismatch at the boundary [25].The boundary data for a frequency domain system are the amplitude and phase of the measured signal, which is used with a Finite Element Method (FEM) based reconstruction procedure to obtain the internal spatial distributions of µ a and µ ′ s .

Standard image reconstruction
The objective function (Ω) for this procedure can be written as Where, F is the forward operator that generates the model response and y is the experimental measured data. . 2 represents L2-Norm of the vector.This is also known as the Tikhonov approach [27], where λ is the regularization parameter that balances the current estimate of ) [28].Minimizing Eq. (3) (by setting first derivatives with respect to µ a and D to zero) leads to an update equation Note that deriving this update equation is not so straight forward, the full derivation is given in Ref. [28].In Eq. 4, J is the Jacobian matrix and I is the identity matrix.Note that J T J is ill-conditioned; λ I stabilizes the matrix.However, a slight deviation from this update equation is generally employed, which is also known as the Levenberg-Marquardt (LM) regularization method [29,30], assuming [(δ D, δ µ a ) = (D, µ a ) − (D 0 , µ a0 )] leading to [31]   (J Most of the literature reports λ * ≡ 2λ [20,31], which is true only for the Levenberg-Marquardt minimization which does not involve the parameter field in the objective function (Eq. 3) [29,30].A detailed discussion about different least-squares minimization methods is presented in Ref. [28].In this LM approach, λ * typically starts being the ratio of the variances and is reduced at each of the iterations by a small factor (in here, it is √ 10 and also multiplied by the maximum of the diagonal values of J T J [21]).The iterative procedure is repeated until experimental data matches with modeled data within a preset value ε (≈ data noise level).In general, the initial values, (D 0 , µ a0 ), are obtained from a pre-reconstruction step where the data is calibrated by a homogeneous fitting procedure [32,33].

Inclusion of a priori information: Soft-Priors
The objective function with inclusion of prior information is given as [15] Here also λ is the ratio of the variance of the data noise to parameter field and L is a penalty matrix (dimensionless in all the cases considered in this work) which can be derived from MRI structural priors as indicated below.The update equation resulting from this procedure is: In this work, each location in the computational discretized model is labeled according to tissue type (fatty, fibroglandular or tumor) determined from MRI T1-weighted images [15,16,23].It was also assumed that there is no covariance between the different regions of the imaging domain.Since the domain model does not itself change throughout the iterative reconstruction algorithm, the L-matrix is calculated before the reconstruction procedure and it is used through out the process to penalize the solution.In implementation of Eq. 7 for simultaneous reconstruction of D and µ a , using the Jacobian form described in Ref. [18] leads to the block matrices where H represents Hessian matrix (J T J).Since D and µ a are considered to be independent of each other in the estimation procedure, the cross-terms in the combined regularization matrix (λ L T L) become zero.Two forms for the L-matrix are discussed in the subsections below, including the Laplacian and Helmholtz structures.

a Laplacian-type structured regularization matrix
The Laplace equation in parameter u(r) can be written as A finite difference approximation to this equation for 'N' number of equi-space (step size = h) nodes can be written as [34] Dividing the whole equation by '-N' leads to The L matrix is a matrix that relates each nodal property of the numerical model to all other nodes.Therefore given a node i within the mesh, its relationship to another node j having Laplacian structure (Eq.11) within the same mesh can be given as [15,16] L i j =    0 if i and j are not in the same region −1/N if i and j are in the same region where N is the number of finite element mesh nodes comprising a given region.In this case, L T L approximates a second-order Laplacian smoothing operator within each region, and functionally works to average the update within a region, while allowing discontinuity between different regions.

a Helmholtz-type structured regularization matrix
The Helmholtz equation in parameter u(r) for a damped wave can be written as where κ is the wave number, specifically κ = 1/l, where l covariance length [34].l also represents the decay length scale over which the parameter u(r) has correlation.Making κ = 0, will give the Laplace equation (Eq.9).A finite difference approximation to this equation for 'N' number of equi-space (step size = h) nodes can be written as [34] Dividing the whole equation by '−(N + (κh) 2 )' gives = 0 (15) Writing this as a generalized L-matrix form similar to Eq. 12 if i and j are not in the same region − 1 N+(κh) 2   if i and j are in the same region For the FEM nodes case, h is chosen to be the distance between the nodes.Moreover, κ = 1/l is generally chosen to be the inverse of the size of the feature (tumor in this case) in the imaging domain.This can be seen as using the best priors for the estimation of optical properties [28], termed as best priors estimate (BPE).As the prior structural information is available through MRI, l is chosen to be the diameter of the target (tumor) in the BPE case.In this case, L T L (L given by Eq. 16) approximates a second-order Helmholtz smoothing operator.To determine the effect of κ on the parameter reconstruction, different values for κ are chosen.It is shown that for small values of κ, which corresponds to a large correlation length (l), both Laplacian and Helmholtz structures recover the same optical property distribution.

Inclusion of a priori information: hard-priors
Reduction of parameter space to the number of regions segmented from a high resolution imaging modality is known as hard-priors [28].In this section, a detailed formulation for this approach is given.The estimation optical properties (D and µ a ) in this procedure is performed through LM minimization, which was discussed in detail in Ref. [28].The update equation in this case is given by Eq. 5, with Jacobian having dimension of (2*NM)x(2*NR) instead of (2*NM)x(2*NN).In here, NM, NN and NR represent the number of measurements, number of FEM nodes, and number of regions, respectively.The multiplication factor 2 in the dimensions results from the treatment of amplitude and phase separately for the frequency-domain signal, and the estimated parameters being D and µ a .This transformation requires multiplication of original Jacobian, J, by region mapper matrix (R) [22].
where the R has dimension of (2*NN)x(2*NR), having a block form R = Σ Σ .J can be also seen as a block matrix with J = [J D J µ a ].If RL j represented a segmented region (j) in the FEM mesh for j = 1, 2, ...NR, then Σ has the from Overall, the new Jacobian J matrix elements were produced by adding the sensitivity of nodes belonging to the same region [22].Using J instead of J, in Eq. 5 gives the update of optical properties of the segmented regions δ D, δ µ a having dimension of (2*NR)x1.To map back this solution vector to the original dimensions (2*NN)x1, requires It is also important to relaize that even though, NR≪NM, solving this might still be illposed [20], leading to LM update equation [28].

Breast geometry -Effect of imperfect a priori information
The techniques described in Sec. 2 were used to reconstruct images from synthetic data with 1% random noise added.Numerical experiments using synthetic data generated on MRI T1weighted breast images with incorrect priors to show the effectiveness of soft-priors.Figure 1(a) (first column) shows the original distribution of three tissue layers, namely fatty (µ a = 0.006 mm −1 and µ ′ s = 0.6 mm −1 ), fibro glandular (µ a = 0.012 mm −1 and µ ′ s = 1.2 mm −1 ), and tumor (µ a = 0.018 mm −1 and µ ′ s = 1.8 mm −1 ) for the breast geometry (labeled as 0, 1 and 2 respectively, in Fig. 1 equally spaced on a circle (indentions in Fig. 1(a)).In succession, one fiber was used as the source while all other fibers served as detectors which provided a total of 240 measurements.In these studies, the source was modeled as a Gaussian profile with a Full Width Half Maximum (FWHM) of 3 mm to most accurately represent the light applied in the clinical system[35], and was placed at a depth of one transport scattering distance from the tissue boundary.A mesh of 1785 nodes (corresponding to 3418 linear triangular elements) with uniform nodal density was used for the diffusion model and reconstruction calculations [36].A total of 7% of the glandular layer (label-1) FEM nodes were labeled (relative to the original glandular layer nodes) as fat (label-0) to introduce imperfect structural priors.
The same initial estimates (optical properties of region '0') were used as homogeneous starting conditions.The iterative procedure was stopped once the data-model misfit (residual) did not improve by more than 2% when compared with the previous iteration.The starting value for λ is chosen to be 25000 and 75 for µ a and D respectively, derived from the noise characteristics [28], for Eq. 5.For the case of Eq. 7, a starting value of 10 was chosen and was decreased as the iterations progressed, this procedure is outlined in Ref. [21].
Table 1.Mean and standard deviation of the reconstructed (a) µ a and (b) µ ′ s values in different regions (labeled in first column of Fig. 1(a)) recovered with simulated data having 1% random noise and imperfect structural information defining region '1' (7% reduction compared to the original segmentation).The corresponding reconstructed images are shown in Fig. 1(a

Breast geometry -effect of data noise level
The techniques described in Sec. 2 were used to reconstruct images from synthetic data with 1, 3, 5 and 10% Gaussian distributed noise to see the effect of data noise level on the reconstruction techniques.The breast geometry (and the optical properties) were equivalent to the previous section (Fig. 1(a)), however, perfect spatial priors were used.The same FEM mesh as described above was employed in the forward and reconstruction problems.Optical properties of region '0' were used as initial guess for the iterative procedure.The regularization parameter (λ ) and stopping criterion was chosen according to each data noise level.

Phantom studies
A multi-layered gelatin phantom (86 mm diameter, 25 mm height) was fabricated with different optical properties using heated mixtures of water (80%), gelatin (20%) (G2625, Sigma Inc.), India ink for absorption, and TiO 2 (titanium oxide powder, Sigma Inc.) for scattering [15] (see Fig. 3).Different layers of gelatin were constructed by successively hardening gel solutions containing different amounts of ink and TiO 2 .A cylindrical hole (diameter of 16 mm and height of 24 mm) was filled with liquid to mimic the tumor.The first column in Fig. 4 shows the axial cross-section of three-layers of the phantom (Fig. 3) where the region labeled '0' has the optical properties (µ a = 0.0065 mm −1 and µ ′ s = 0. 65 mm −1 ), similar to the typical fatty layer in the breast [6] and a thickness of 10 mm.The fibroglandular layer (diameter 76mm) has optical properties (region labeled '1') of µ a = 0.01 mm −1 and µ ′ s = 1.0 mm −1 .The tumor (represented by the region labeled '2') has a diameter of 16 mm with optical properties of µ a = 0.02 mm −1 and µ ′ s = 1.2 mm −1 .The optical properties were validated by measuring large cylindrical samples of each layer.Appropriate mixtures of Intralipid and India ink were used to achieve the desired optical parameters of the tumor.Data was acquired using a clinical NIR system [35] where the fibers were marked and photographed to extract region information (analogous to MRI images).This regional information was used to label the corresponding regions in the FEM mesh [15].A mesh of 1785 nodes (corresponding to 3418 linear triangular elements) was used for the diffusion model calculations and a mesh having 1360 nodes was used in the reconstruction [36].The meshes considered in this work have uniform nodal density across the total computation volume.NIR data was calibrated using a reference homogenous phantom to obtain initial optical properties estimates and minimize the variation between the 16 optical channels according to standard practice in human imaging studies [32,33].

Results
Reconstructed µ a and µ ′ s images obtained from the noisy simulated data with imperfect (7% error) glandular layer priors using the methods described in Sec. 2 are shown in Fig. 1.Using hard priors, the total number of unknowns is reduced to 6 parameters (µ a and µ ′ s for each of the 3 regions) and these images are presented in the last column of Fig. 1(a).The images from two different approaches of soft-priors are shown in the middle 2 columns; the first column shows the expected results.For the Helmholtz case, κ = 1/8 mm −1 (BPE) was used, where 8 mm is the diameter of the tumor.Cross-sectional plots of the reconstructed µ a and µ ′ s distributions along the dotted line in Fig. 1(a) (first column) are provided in Fig. 1(b).Table 1(a) and (b) show the mean and standard deviation of the optical property estimations in each region of the reconstructed images.Note that the NIR reconstruction procedure without prior information (not shown) did not generate meaningful images in this complex case.
Figure 2(a) shows the reconstructed images from the data with 5% noise using all four techniques described in Sec. 2. The first column of Fig. 1(a) shows the actual distribution of optical properties.gives lesser standard deviation for the absorption coefficient (µ a ) reconstruction compared to the Helmholtz structure.On the other hand, Helmholtz structure performs better than Laplacian in the case of scattering coefficient (µ ′ s ).A photograph of the phantom used to collect data at 785nm with 16 fibers in a single plane (giving 240 measurements) is shown in Fig. 3. Images obtained from the procedures described in Sec. 2 are presented in Fig. 4(a) along with cross-sectional plots of the reconstructed µ a and µ ′ s distributions in Fig. 4(b).Table 2(a) and (b) show the mean and standard deviation of the optical property estimations in each region of the reconstructed images.Here also, for Helmholtz case, the BPE (κ = 1/16 mm −1 ) was used.Figure 5(a) gives the results for different values of κ (given on the top of the each column, true distribution is shown as first column in Fig. 4(a)) in the Helmholtz case.Corresponding cross-sections are plotted in Fig. 5(b).The mean and standard deviation values for each of this case are also given in Table 2(a) and (b).

Discussion
The reconstructed results (Fig. 1, 2, 4 and 5) show that the structural priors improve the reconstructed image quality dramatically.The penalized problem formulation (different type of regularizations) generates smoother images resulting in smaller standard deviations from the mean values (see Table -2) as compared to the generalized problem that does not incorporate prior information.
The hard-prior case produces significant optical property value error when the structural a priori information is imperfect in the breast geometry (Fig. 1 and Table-1).In this case, a 7% variation in the definition of the glandular layer caused false estimates of the optical properties.On the other hand, soft-priors (Laplacian and Helmholtz) yield good estimates for each layer.Hard-priors over-estimate the tumor absorption coefficient by 360% and under-estimate its scattering coefficient by 88% (Fig. 1).Soft-priors are within 6% of the expected values even with the error in the structural prior.Note that in this particular case, hard-priors failed when 7% of glandular layer made as a fatty layer, yet below this error value it gave reasonable estimates of optical properties.It is also important to realize that this is only a test case of the spatial-prior error progation.There is not enough in-vivo data for a systematic investiagion of spatial-error propagation in the reconstruction procedures, especially in the hard-priors case.
With perfect structural priors, increasing the data noise level also increases the quantification error in the reconstructed images (Fig. 2(a) and(b)).Hard-Priors clearly fail in quantifying the tumor optical properties as the data noise level increases.Soft-Priors does better in the quantification than Hard-Priors.It is also evident that incorporation of structural information is key for accurate quantification of the optical properties.The experimental results (Fig. 4) from the three-layer gelatin phantom (Fig. 3) also show that incorporation of perfect structural information improves the quantification and quality of the reconstructed images.Specifically, the mean and standard deviation of the reconstructed optical property values in each region are both more accurate and more precise where the priors are included.The mean values (Table -2) show that the absorption coefficient for the Region-0 (fat) is under-estimated by 77% in the case of the Helmholtz regularization (BPE).This error can be explained by the fact that the Euclidean distance between the nodes was used rather than the distance between the nodes along the boundary of a particular region.Both Laplacian and Hard-Priors over estimate the scattering coefficient of the tumor region by a factor of 2. It is known that the photon path length is affected by the scattering coefficient.By constraining the problem based on the distance, one can expect to estimate the scattering better.In this experimental case, table-2 indicates that the Helmholtz technique produces more quantitative accuracy for the scattering coefficient estimation of the tumor and the Laplacian technique is best for the absorption coefficient estimation (as well as Fig. 2 and 4).Usage of the Helmholtz-type form for the diffusion part and the Laplacian- type form for the absoprtion part (in Eq. 8), may lead to a better estimation of both optical properties, which will be considered as a future work.Under estimation of optical properties in the regions '0' and '1' can be attributed to the random fluctuations in heating/cooling of the gelatin solution especially in making layered phantom, which eventually changes optical properties.Moreover, ink and TiO 2 particles settling to the bottom, making the mixture less homogenous also causes the change optical properties.These type of errors are more common for layered gelatin phantoms[37].

Methods
Region-0 Region- Theoretically Helmholtz and Laplacian structures are identical when κ = 0 (equivalently l is large).Figures 5(a) and (b) (as well as Table -2) show that when κ = 1/86 mm −1 (l = 86 mm is the diameter of the domain), Laplacian and Helmholtz structures give reasonably close reconstruction values of optical properties, which indicates the expected trend presented in this paper for the two methods.It also indicates that the BPE case (κ = 1/16 mm −1 ) gives the best results, as the priors applied are close to the true structure of the feature (tumor).These results also indicate that unreasonable constraints (like κ = 1/5 mm −1 , Fig. 5(a) first column) makes the estimation problem amplify the noise resulting in physiologically invalid (scattering coefficient is greater for fatty layer compared to fibroglandular layer) estimates of optical properties.In the tumor region, as κ decreases, the estimated values of optical parameters become closer to expected values (Table-2 and Fig. 5, κ = 1/43 and κ = 1/86 cases).This is due to the correlation length becoming larger, making the covariance in the neighboring nodes larger.
In this study, it is shown that imperfect priors (commonly caused by improper image segmentation and image artifacts in MRI or X-ray) can lead to error-prone results in the hard-prior case whereas soft-priors are more immune to uncertainty in the prior data.It is also shown that the techniques used to incorporate the soft structural prior information influences the image outcome, which may lead to improvements in image accuracy if properly implemented.Srini-vasan et.al [38] have found that 5% error in optical properties introduces approximately 45% error in the chromophore concentration when a limited number of wavelengths are used for imaging.The correct "soft" inclusion of a priori information therefore can be expected to lead to a more accurate quantification of chromophore concentrations as well.It should be noted that over-weighting of the penalty term in the problem formulation may make the solution ignore the data-model misfit and emphasize smooth feature extraction.The techniques developed in this work were applied for two-dimensional test objects, and can be easily extended to threedimensional case.A more extensive study of this is left for future investigations.

Conclusions
This work has investigated several ways of incorporating structural information into an iterative image reconstruction.The results have been supported by gelatin phantom experiments that multi-layered structure which is commonly found in breast tissue adipose (fatty) tissue on the exterior and fibroglandular tissue nearer to the interior.Soft-priors allow the tissue optical properties to vary within predefined regions, unlike hard-priors which constrain these zones to be homogeneous.Hard-priors were found to perform poorly when the prior information contained area errors as small as 7% which can easily be produced by most segmentation algorithms.True boundary extraction from MRI images introduces unavoidable segmentation and discretization errors that are better tolerated when the structural information is encoded through the soft-prior approach involving a penalty matrix.
The results reported here indicate that the optical properties of different tissue types can be quantified more accurately when their estimation is properly guided by "soft" structural a priori information.The problem formulation and results presented in this work indicate that data from other imaging modalities such as ultrasound or x-ray tomography, could also be used as the source of the structural prior.In the experimental cases investigated, the Helmholtz structure using best priors gave a better estimation of scattering coefficient of the target (tumor).However, the Laplacian type of regularization leads to more superior absorption coefficient estimate.So it is reasonable to conclude that Laplacian structure gives the best estimates of total hemoglobin concentration (Hb T ), hemoglobin oxygen saturation (S t O 2 %) and water fraction (H 2 O) (which are the main absorbers).Helmholtz structure gives the best estimates of the scattering power and scattering amplitude (scattering parameters).The framework presented here can also be extended to other regularization terms such as total variation minimization or spectral prior constraints, which may be studied in future work.

Fig. 1 .
Fig. 1.(a) Simulated µ a and µ ′ s distributions from a breast (obtained from a volunteer) are shown in the first column.Optical properties for the region labeled '0' (fat) are: µ a = 0.006 mm −1 and µ s = 0.6 mm −1 .Region '1' (fibroglandular) values are: µ a = 0.012 mm −1 and µ ′ s = 1.2 mm −1 .Region '2' (tumor) values are: µ a = 0.018 mm −1 and µ ′ s = 1.8 mm −1 .Reconstructed µ a and µ ′ s images from different techniques with simulated data having 1% random noise and imperfect structural information in defining region '1' (7% reduction compared to the original segmentation) are shown in the rest of the columns.The middle two columns use soft prior structural information while the last column shows the result with hard prior information.In the Helmholtz case, κ = 1/8 mm −1 (BPE) was used.(b) Cross-sectional plots, along the dotted line in the actual image (see first column of (a)), of true and reconstructed µ a and µ ′ s distributions.#81441 -$15.00USD Received 26 Mar 2007; revised 11 Jun 2007; accepted 11 Jun 2007; published 13 Jun 2007

#Fig. 2 .
Fig. 2. (a) Reconstructed µ a and µ ′ s images from different techniques with simulated data having 5% random noise and perfect structural priors (actual images are shown in the first column of Fig. 1(a)).The first column shows the reconstruction results without the use of prior information.The middle two columns use soft prior structural information while the last row shows the result with hard prior information.In the Helmholtz case, κ = 1/8 mm −1 (BPE) was used.(b) The mean values and standard deviations (plotted as error bars) in µ a and µ ′ s for different regions of breast geometry (labeled in actual image) with increasing noise level (1% to 10 %).

Fig. 3 .
Fig. 3. Photograph for gelatin phantom (representing the idealized two-dimensional crosssectional geometry shown as first column in Fig. 4(a)) used in the experimental studies.

Figure 2 (
b) shows the mean and standard deviation values (as error bars) of reconstructed images using different techniques for different regions of the breast geometry with increasing data noise level.In the Helmholtz case, κ = 1/8 mm −1 (BPE) was used.The actual values are also plotted for the comparison.It can be clearly seen Laplacian regularization #81441 -$15.00USD Received 26 Mar 2007; revised 11 Jun 2007; accepted 11 Jun 2007; published 13 Jun 2007

Fig. 4 .Fig. 5 .
Fig. 4. (a) Actual µ a and µ ′ s distributions (axial cross-section) of phantom (Fig. 3) case are shown in the first column.Optical properties for the region labeled '0' are: µ a = 0.0065 mm −1 and µ ′ s = 0. 65 mm −1 .Region '1' values are: µ a = 0.01 mm −1 and µ ′ s = 1.0 mm −1 .Region '2' (tumor) values are: µ a = 0.02 mm −1 and µ ′ s = 1.2 mm −1 .Reconstructed µ a and µ ′ s distribution from different techniques (discussed in Sec. 2) from the experimental phantom data.Second column of images does not use prior information.The middle rows use soft prior structural information and the last row of images were recovered with hard priors.In the Helmholtz case, κ = 1/16 mm −1 (BPE) was used.(b) Cross-sectional plots along the dotted line in the actual image (see first column of (a)) of the true and reconstructed µ a and µ ′ s distributions.

Table 2 .
Mean and standard deviation of the reconstructed (a) µ a and (b) µ ′ s values in different regions (labeled in first column of Fig. 4(a)) recovered from the experimental phantom data.The corresponding reconstructed images are shown in Fig. 4(a) and 5(a). )