1. Introduction

Near Infrared (NIR) optical tomography is a method that uses light in the 650-900nm wave-length range to recover images of the internal spatial distribution of tissue optical properties, absorption (or chromophore concentrations) and scattering parameters[1–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–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–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–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.

2. Methods: mathematical formulation

2.1. 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,ω = 100MHz). 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.

2.2. 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. ∥.∥² 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 optical properties and the initial values, and the data-model misfit. More specifically, it is the ratio of the variances of the data noise and the parameter (λ = σ² _y/σ² _(D,μa))[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 ₀,μ_a0)] leading to [31]

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 LMapproach, λ^* 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 ₀,μ_a0), are obtained from a pre-reconstruction step where the data is calibrated by a homogeneous fitting procedure[32, 33].

2.3. 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.

2.3.1. 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]

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.

2.3.2. 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)²)’ gives

Writing this as a generalized L-matrix form similar to Eq. 12

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.

2.4. 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)×(2^*NR) instead of (2^*NM)×(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)×(2^*NR), having a block form $R=\left[\begin{array}{c}\sum \\ \sum \end{array}\right]$. 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)×1, 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].

3. Methods: simulations and experiments

3.1. 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 T1-weighted 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(a) first column). Sixteen light collection/delivery fibers were arranged 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.

 

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.

Download Full Size | PDF

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)

View Table | View all tables in this article

3.2. 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 FEMmesh 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.

3.3. 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₂ (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₂. 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].

 

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 %).

Download Full Size | PDF

 

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

Download Full Size | PDF

4. 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 (ma 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. 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 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).

5. 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₂ 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].

 

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.

Download Full Size | PDF

 

Fig. 5. (a) Reconstructed μ_a and μ´_s images from the experimental phantom data using Helmholtz-type regularization matrix for different values of κ, which are given at the top of each column. (b) Cross-sectional plots along the dotted line of the actual images in Fig. 4(a) (first column) are shown with the data from reconstructed μ_a and μ´_s images in (a). The best prior estimate (BPE) case (κ = 1/16 mm^-1) is also presented for comparison.

Download Full Size | PDF

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).

View Table | View all tables in this article

Theoretically Helmholtz and Laplacian structures are identical when κ = 0 (equivalently λ 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. Srinivasan 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.

6. 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 represent multi-layered structure which is commonly found in breast tissue with 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_tO₂%) and water fraction (H₂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.