Near-infrared spectral tomography (NIRST) has demonstrated its potential for breast lesion characterization and treatment monitoring.^1,2 However, most studies indicate that NIRST alone may not be capable of providing the quantitative accuracy required for clinical applications because of spatial resolution limitations.³ High-spatial resolution guidance from structural images, such as x-ray or MRI, has improved the quantitative accuracy of NIRST.^1,2,4,5 MRI-guided NIRST is especially suitable for breast cancer imaging because the approach leverages the high sensitivity of MRI, but improves specificity through recovery of NIRST parameters.²

To incorporate prior information from MRI into NIRST, previous approaches have segmented the MRI images into tumor and nontumor regions, each of which is assumed to be optically homogeneous. This “hard-prior” strategy reduces the size of the inversion problem by lumping regions together into a few super pixels or zones.^4,5 An alternative, so called “soft-priors,” approach incorporates Laplacian- or Helmholtz-type structures into the regularization to encode spatial information that retains the inversion parameters at each pixel over the entire domain.⁶ Both techniques require manual segmentation to identify regions within the MRI image volume. This segmentation step adds further complexity to the processing and reduces objectivity when combining the image information. Additionally, segmentation can be time-consuming and requires sufficient experience to avoid introducing bias or error. To overcome the need for segmentation, a direct regularization method of encoding spatial priors has been developed.⁷

To date, spatial information incorporated into the NIRST image reconstruction has been from a single MRI sequence. However, T1-weighted (T1W), T2-weighted (T2W), T1-weighted dynamic contrast enhanced (DCE), and diffusion weighted (DW) images are often acquired in standard clinical breast MRI. Predominantly, only DCE images have been used to improve the accuracy of NIRST; however, the type of MRI sequences and the way in which it is used influences the resulting values.

DCE-guided NIRST image reconstruction is based on vascular-delivered contrast, but few studies have correlated DCE images with molecular biomarkers from NIRST.^2,8 Among the NIRST absorber concentrations, total hemoglobin (HbT) and water have shown the greatest potential to distinguish malignant from benign breast abnormalities.⁹ Yet, from a biophysical standpoint, water values are not expected to be highly related to DCE contrast, but rather DW contrast, since the latter corresponds to the free diffusion of water as measured by MRI.

In this letter, a reconstruction algorithm is presented to directly encode the spatial information derived from DCE and DW MRI into multiple regularization matrices for recovering NIRST HbT and water contents. The motivation for pursuing this algorithm is to improve the accuracy of estimating HbT and water simultaneously. Simulation studies were used to test the approach, and the first image reconstructions from an patient case are presented.

Light propagation in breast tissue was modeled by a diffusion equation solved by the finite element method (FEM) using the NIRFAST software package.¹⁰ Spectral constraints were directly imposed in the reconstruction.¹⁰ The image reconstruction is posed as a least squares optimization problem, defined as

where $f(c)$ is the data fidelity term, $\varphi (c)$ is the regularizer, and $c$ is the concentrations of oxyhemoglobin (HbO), deoxyhemoglobin (Hb), water, lipids, scatter amplitude (SA), and scatter power (SP). The conventional regularizer, $\varphi (c)$, used in NIRST is $\varphi (c)=\lambda ·{\Vert Lc\Vert }_2^2$, where $L$ is an identity matrix or a matrix derived from high-resolution imaging such as MRI, and $\lambda$ is the regularization parameter. Generally, a single image volume is utilized to define $\varphi (c)$ and all optical parameters are regularized with the same spatial behavior.^7,8 However, in practice, vascular and water perfusions in tissue are independent.

Considering that DCE and HbT are associated with tissue vasculature and DW images should be correlated with the distribution of NIRST water content, DCE images were used to construct the regularization matrix ($L_1$) to constrain total hemoglobin ($\mathrm{HbT}=\mathrm{Hb}+\mathrm{HbO}$), whereas DW images were used to construct $L_2$ to constrain water. Thus, the regularizer was formulated as

where $\alpha$ and $\beta$ are regularization parameters, $L_i(i=\mathrm{1,2})$ is the regularization matrices derived from DCE ($L_1$) and DW ($L_2$) MRI images, and have the form⁷ where ${\gamma }_i$ is the grayscale value mapped to node $i$ on the FEM mesh, $r_i$ is the position of node $i$ on the FEM mesh, ${\sigma }_g$ is the characteristic grayscale difference for regularization, and $M_i$ is a normalization factor chosen for each row so that $\sum _i^N{L}_{ij}=0$ $\forall i=1,\dots ,N$, where $N$ is the number of finite element nodes. The function $\theta (·)$ is the Heaviside step function, which determines the local weight assigned to the $i$’th node of the finite element mesh.¹¹ ${\sigma }_d$ is a factor corresponding to the distance of influence of the regularization effect. Generally, the empirical value of ${\sigma }_d$ is 0.4. The operator $L_i$ encodes uniformity by penalizing similarly gray locations in the finite element nodes to have significantly similar optical concentration update values at each iteration.

For choice of regularization parameters $\alpha$ and $\beta$, they are usually obtained from $\rho *\max (J^T*J)$, where $\rho$ is the initial regularization parameter, $J$ is the Jacobin matrix of unknown optical concentrations, and $T$ is the transpose.

To minimize the above optimization problem with respect to unknown optical properties, derivatives of Eq. (1) are set to zero. The reconstruction is stopped when the change between successive iterations is less than 2% or the maximum iteration number is reached. The algorithm scheme is shown in Fig. 1.

To validate the performance of the proposed algorithm, an absolute bias error in the tumor was defined⁷

where $C_i^{\text{true}}$ and $C_i^{\mathrm{recon}}$ represent the true and reconstructed images, respectively, and $M$ is the number of finite element nodes in the tumor.

Numerical simulations were used to test the performance of the proposed algorithm. A three-dimensional (3-D) breast simulation was created from MRI slices of an actual patient with a tumor size of $22\times 20\times 32{\mathrm{mm}}^3$. The corresponding DCE (b) and DW (c) images used to construct the $L_1$ and $L_2$ matrices are shown in Fig. 2. The MRI slices were manually segmented into tumor, adipose, and fibroglandular regions based on the preconstrast T1 and DCE images. After segmentation, typical HbT, water, SA, and SP values for tumor ($34.2\mu \mathrm{M}$, 70.5%, 1.23, and 0.26), adipose ($17.1\mu \mathrm{M}$, 47%, 1.17, and 0.2), and fibrogandular ($25.4\mu \mathrm{M}$, 60%, 1.2, and 0.23) tissues from Ref. 12 were assigned to each region, as shown in Fig. 2(d). For simulation, 16 colocated source-detector positions were placed around the breast. For each source illumination, data were collected at 15 detector locations. Thus, the total number of measurements was 2160 ($240\times 9$) involving nine wavelengths (from 661 to 948 nm). Amplitude-dependent Gaussian noise with a variance of 1% was added to the simulated measurements. To obtain the initial estimates of the optical properties, the simulated data set was calibrated with a known homogeneous circular phantom. Images were reconstructed using a pixel basis of $25\times 25\times 25$. The same ${\sigma }_g$ and ${\sigma }_d$ of 0.01 and 0.4 were used to construct $L_1$ and $L_2$, respectively. The initial parameters $\rho$ used in different methods were same, and their value was 100.

The no-prior images [Fig. 2(e)] yielded very limited tumor contrast relative to the surrounding normal tissue (contrast values of HbT and water: 1.17 and 1.07), whereas the reconstructed images of HbT ($\mu \mathrm{M}$) and water content (%) using the prior information from DCE or/and DW images [Figs. 2(f)–2(h)] produced much better tumor contrast (contrast values of HbT and water: 1.58 and 1.12, 1.44 and 1.11, 1.34 and 1.14). Cross-section profiles of HbT and water through the slice [along white lines in Fig. 2(d)] are presented in Fig. 3. As shown in Fig. 3(a), the reconstructed HbT concentrations overestimated the true HbT in the tumor, especially in the DCE case, while the values from the DCE-DW guided reconstruction are closer to the true results. The cross-sectional results [Fig. 3(b)] also show that the recovered water content by combining DCE and DW images is more accurate than the results only guided by DCE or DW images. Table 1 summarizes the bias errors in HbT and water images in the tumor reconstructed by no-prior, DCE, DW, and DCE-DW image guided methods. The bias errors of HbT and water in the DCE-DW guided image reconstruction were reduced 22% and 18%, 21% and 6%, and 10% and 11%, compared to that in the no-prior, DCE- or DW-guided reconstructed images, respectively. This is aligned with a recent report on systematic validation of prior-guided diffuse optical tomography reconstructions.¹³

Table 1

Absolute bias errors in recovered HbT and water in the tumor.

Finally, a patient image volume with actual optical measurements was reconstructed to test the proposed algorithm with in vivo data. The patient had a pathologically confirmed $25\times 22\times 26{\mathrm{mm}}^3$ invasive ductal carcinoma (IDC) in her right breast. However, her BIRADS score diagnosed by two experienced radiologists (J. Xu, 15 years and K. Wang, 6 years) according to MRI images was 2, which was false-negative. Six wavelengths of frequency domain data and three wavelengths of continuous-wavelength data were collected with an MRI-guided NIRST system.^1,10 Optical data and MR images (T1, T2W, DW, and DCE) were acquired simultaneously.

A uniform finite element mesh was created from the precontrast T1 images without segmentation. Initial estimates of optical property values were obtained by direct calibration method. Images were reconstructed based on a pixel basis of $25\times 25\times 25$. The same ${\sigma }_g$ and ${\sigma }_d$ of 0.02 and 0.4 were used to construct $L_1$ and $L_2$, respectively; the initial regularization parameters used in all methods were the same, and the value was 10.

Figure 4 shows the MR-T1, MR-DCE, MR-DW, and optical images overlaid on T1 images in the plane $z=-0.25\mathrm{mm}$. The optical images reconstructed by the no-prior approach (d), guided by DCE (e), DW (f), combined DCE-DW (g) are also presented. In the no-prior images shown in Fig. 4(d), the tumor is hard to distinguish from the normal surrounding tissue in either the HbT or water images while it can be recognized in the HbT and water images in Figs. 4(e)–4(g). The tumor volume and location were defined by a radiologist based on the DCE MRI scan. The contrast values of HbT and water reconstructed by DCE and DW guided image reconstruction are 1.4 and 1.4, and 1.8 and 1.4, respectively. Comparing to the above contrast values, a better contrast in HbT and water are recovered at 4.6 and 1.6 when using DCE-DW guided image reconstruction. These higher values are more in line with what would be expected biologically, and match the higher contrasts seen in the MRI scans. It is not known exactly what the true contrast is, because of this being a clinical exam, but the benchmarks of expected biological contrast and obsered MRI contrast can be taken as a guide. Optical data correctly identified the lesion as malignant. The contrasts of HbT and water were highest in the DCE-DW case. The improvement in both parameters together is consistent with observations from the simulation studies.

The goal of this work has been to develop a reconstruction algorithm that directly encodes spatial information derived from both DCE and DW MRI images into separate regularization matrices to improve the accuracy of NIRST HbT and water content. The simulated and patient image reconstruction results suggest that the combined approach of DCE-DW priors enhances tumor contrast in both HbT and water images, relative to using only DCE or DW data as priors. In addition, the structural information from DCE and DW MRI is directly encoded without the need for segmentation prior to reconstruction.

The manual selection of regions-of-interest is always needed since the DCE and DW MRI are very complicate (especially the low resolution of MR-DW images), and it is very difficult to adapt a computer-operable definition for the segmentation. In contrast, our direct regularization is a method that automatically combines the structure information from DCE and DW to the optical reconstruction. By this method, the subjective bias and error due to the operator’s training and experience during the segmentation process can be eliminated.⁵

Comparing to the algorithm reported in Ref. 14, which focused on using different regularization coefficients for unknowns from different spatial domains, it is different from our algorithm. Our multiobjective function is similar to that used in Ref. 15. However, the novelty of the proposed algorithm focused on independently regularizing (i) the vascular property values ($\mathrm{HbT}=\mathrm{Hb}+\mathrm{HbO}$) by the DCE MRI images and (ii) the water values by DW MRI images, in order to recover NIRST HbT and water contents without a user-defined segmentation process. Moreover, the weighting factor $0<w<1$ used in Ref. 15, which controls the contribution of each individual regularization matrix to the global regularization matrix, was excluded from our algorithm.

As shown in Fig. 4, the false-negative diagnosis by MRI (DCE + DW) has been corrected by NIRST imaging outcome. The malignance has been confirmed by the surgical pathologic results. In addition, as we reported in our previous paper, the sensitivity and specificity of using MRI (combining T2W, DCE, and DW) to differentiate malignant breast tumor to benign are 95% and 67%, respectively.² We expect the specificity can be improved significantly by adding NIRST outcome, especially, the NIRST outcomes estimated from the images reconstructed by this new approach.

In our experience, ${\sigma }_g$ plays an important role in determining the image quality. Smaller ${\sigma }_g$ may overestimate the peak value, and a larger ${\sigma }_g$ may over-smooth. After numerous tests, we found that ${\sigma }_g$ usually varied between the interval of [0.001, 0.1], and the default value is 0.01. According to the physiologic normal range of HbT, we also introduced an automatic ${\sigma }_g$ selection algorithm. If HbT peak value under certain ${\sigma }_g$ is larger than $1000\mu \mathrm{M}$, ${\sigma }_g$ will be increased by 0.01; if the peak value of HbT is smaller than $10\mu \mathrm{M}$, the ${\sigma }_g$ will be decreased to 0.001. Additionally, if the breast density is extremely dense, ${\sigma }_g=0.001$ is better to construct $L_2$ based on MR-DW images.