1.

Introduction

Magnetic resonance imaging (MRI)-guided near-infrared spectral tomography (NIRST) has been developed at just a few institutions, but has the potential to add molecular information to the spatial maps of MR imaging of the breast, and thereby increases the specificity of contrast-enhanced MRI exams.^1,2 In this imaging modality, the NIRST data measured at nine distinct wavelengths (ranging from 660 to 950 nm) and clinical MRI images are simultaneously acquired. Through model-based image recovery, the intrinsic biophysical composition of tissue can be estimated in terms of the concentrations of total hemoglobin (Hbt) and oxyhemoglobin, and the percentages of water and lipids.^3–5 Additionally, the ultrastructural cellular density and size ensemble associated with the extracellular matrix and subcellular constituents of breast tissue can be interrogated from the NIRST scattering spectrum. These parameters are spatially recovered according to anatomical regions identified from segmentation of the clinical breast MRI exams into regions of adipose, fibroglandular, and suspicious tissues.^6,7

NIRST image recovery⁸ is nonlinear and ill-posed and has been the subject of many years of research.^9–11 Reconstruction of the MRI/NIRST images requires the inversion of the Jacobian matrix, which encodes the spatial sensitivity to relative changes in the absorbing and scattering parameters of all locations for each source and detector pair. A widely used approach to solve the inverse problem is the Newton-Raphson technique regularized by a modified Levenberg–Marquardt (LM) algorithm.^12,13 Using segmented regions from MRI dramatically reduces the ill-conditioning of the problem to the point where it is overdetermined if the optical properties of only three regions are recovered (adipose, fibroglandular, and lesion). Nonetheless, the diffuse propagation of NIRS light in tissue still generates a poorly conditioned matrix that requires inversion. By regularizing the matrix to be inverted, amplification of the effects of measurement noise on the inversion is reduced, although the resulting solution is modified. Thus, the choice of regularization value is of great importance and has been studied by many researchers.^10,14 High regularization values smooth the spatial modes contained in the resulting images and may yield inadequate contrast recovery of the real NIRS parameters, whereas low regularization values do not sufficiently inhibit the noise effects from being amplified in the inversion process. Implementation of a method to automate the choice of regularization is critical for patient imaging.^15–18 While methods to choose this parameter are automatically well established in computational studies, little investigation of how the selection influences the diagnostic performance of the imaging method in actual practice has been reported. In this study, optimization of the regularization parameter based on L-curve analysis was pursued in clinical MRI/NIRST imaging with the specific goal of retrospectively maximizing the discrimination between known malignant and benign breast scans.¹⁴

A classically defined L-curve optimization algorithm was developed for regularization parameter selection during image reconstruction of data from 25 patients collected in a clinical study of women with breast abnormalities (BI-RADS category 4-5) of unknown diagnosis at the time of the imaging exam. The methodology was also used to study NIRST reconstruction with either amplitude data or both amplitude and phase data acquired with a multichannel 100 MHz frequency-domain NIR spectroscopy system. The performance of the optimal regularization was compared to fixed regularization in this setting of differentiating malignant from benign breast abnormalities. To the best of our knowledge, this study represents the first time the effects of regularization have been investigated in MRI-guided NIRST based on a relatively large dataset of clinical exams with the goal of optimizing task-based discrimination.

2.

Methods

2.4.

Optimization of Regularization using L-Curve

A common approach to selecting the regularization parameter, $\lambda$, is to use a fixed empirical number based on prior tissue-phantom studies, and to apply this value to data acquired from patient exams. In this study, an L-curve approach was used to find the optimal regularization parameter.^15,16 The L-curve is a parametric plot of the L-2 norm of the data-model mismatch $\rho (\lambda )$ versus the difference between optical properties of two iterations $\eta (\lambda )$, where

Eq. (5)

$$\rho (\lambda )=J\delta \mu -\delta \mathrm{\Phi },$$

Eq. (6)

$$\eta (\lambda )=\sum _{i=1}^{NN}{({\mu }_j-{\mu }_0)}^2.$$

In other words, $\rho (\lambda )$ and $\eta (\lambda )$ represent the model-data error and spatial prior error, respectively, each of which is being minimized during the model inversion. With a relatively small $\lambda$, $\rho (\lambda )$ dominates the objective function, and a lower model error is expected at the cost of a larger prior error, and vice-versa for the case of large $\lambda$. Plotting $\rho (\lambda )$ versus $\eta (\lambda )$ for a range of $\lambda$ illustrates the trade-off between these two types of errors, which typically exhibit an L-shaped curve. The corner of the L-curve is commonly regarded as the optimal regularization¹⁶ because it minimizes the two error terms. In this study, the L-curve method was applied to determine the optimal regularization at each iteration.

Figure 1(a) outlines the sequence for optimization of regularization based on L-curve analysis at each iteration. To begin, prior error, which is the difference in optical property solutions between two iterations, is plotted versus model error over a range of discrete regularization values from 0.001 to 100. Figure 1(b) shows a typical L-curve at the first iteration of image reconstruction for a patient with a malignant breast abnormality. Here, the prior error, $\eta (\lambda )$, is defined as the difference of the sum of squared reconstructed optical parameters between the first iteration and initial estimate. The model-data mismatch, $\rho (\lambda )$, depicts the model error in the first iteration. When regularization increases, the model error increases and prior error decreases. Next, constrained regions are selected. In the range of regularization values from 0.001 to 0.02, the behavior of the prior error versus model error does not show a linear relationship because a small regularization results in unstable solutions in the ill-conditioned inverse problem. Once the constrained regions are selected, the L-curve can be fitted through least squares to slopes of $m_v$ and $m_h$ for vertical and horizontal regions, respectively, and a point on the L-curve with a maximal second derivative can be obtained. If this L-metric is larger than a threshold, an optimal regularization is assigned as the point which is closest to the intersection of the two fitting lines. Otherwise, the fitting at this iteration is discarded, and the optimal regularization value from the previous iteration is retained. The algorithm is repeated for each iteration until the stopping criterion for image reconstruction is satisfied.

Fig. 1

(a) Flowchart outlining the sequence for the optimization algorithm. (b) The L2 norm of the prior property error versus L2 norm of the model error creates the L-curve of values for each regularization parameter values from 0.001 to 100. The optimization algorithm is implemented at each iteration, and optimal regularization parameter from previous iteration is used once the L metric falls behind the threshold.

3.

Results

Posthoc analysis performed using Tukey’s HSD test revealed significant differences between optical property contrast when using a fixed regularization of 1 compared to using a fixed regularization of 0.1 ($p<0.001$), and using a fixed regularization of 1 compared to optimal regularization ($p<0.05$).

Figure 2 shows the recovered images of HbT, ${\mathrm{StO}}_2$, water, lipid, SA, and SP of a patient with a biopsy-confirmed breast malignancy. The tumor-to-adipose contrasts of HbT, ${\mathrm{StO}}_2$, water, lipid, SA, and SP were found to be 1.2, 1.1, 1.3, 1.6, 1.4, and 1.02, respectively. A fixed regularization parameter of 0.1 was used through the image reconstruction procedure for this patient.

Fig. 2

Images from a 33-year-old patient with a $11\times 21\times 14\text{-}{\mathrm{mm}}^3$ biopsy-confirmed invasive ductal carcinoma (IDC) in her right breast. (a) Noncontrast T1 MRI with tumor location indicated (arrow); (b) Reconstructed images for total hemoglobin (HbT), (c) oxygen saturation (${\mathrm{StO}}_2$), (d) water, (e) lipid, (f) scattering amplitude (SA), and (g) scattering power (SP) is overlaid on the MR scan. The value of each parameter in the adipose region is suppressed for clarity of visualization.

Figure 3 shows the tumor-to-adipose contrast in HbT versus regularization for benign and malignant patients when images were estimated from amplitude data only, and box plots of HbT contrast in the two diagnostic groups for fixed regularization parameters of either 0.1 or 1, respectively. Compared to the malignant group, the HbT contrast was less sensitive to variation in regularization in the benign group, especially over the range of values from 0.01 to 0.2. Larger separation or difference in mean HbT contrast, between the malignant and benign groups resulted from the lower regularization of 0.1, even though less variation occurred within the same groups with a higher regularization of 1, although a statistically significant difference existed between the means of the malignant and benign groups with $p$ $\text{values}<0.05$ in either case.

Fig. 3

Tumor-to-adipose contrast in HbT versus regularization for (a) benign, and (b) malignant cases. Circles have a regularization of 0.1, and asterisks have a regularization of 1. Box plots of contrast for the two pathologies are shown with fixed regularizations of (c) 0.1 and (d) 1 for all patients.

Figure 4 shows the L-curves for the first three iterations during optical image reconstruction of a patient with a malignant breast abnormality. The same range of regularization values (0.001 to 100) was used to plot the L-curve at each iteration. Both the range of model error and prior error decrease with an increasing number of iterations. Optimization of regularization was implemented using the algorithm outlined in Fig. 1(a); the optimal point is marked by a red asterisk in the first two iterations and had values of 0.18 and 0.22, respectively. At the third iteration, the L-curve metric was below the threshold, and an optimal regularization could not be found using the same process. Instead, the regularization value for the third iteration was set to be the same as used in the second iteration, namely 0.22. The L-curves of the optical data acquired from most patient exams exhibited similar behaviors. The number of patient datasets which had an optimal regularization based on L-curve analysis was 22, 12, and 0 for the first three iterations. The average optimal regularization for the first and second iterations was 0.21 and 0.24, respectively. After three iterations, none of the exam data had an L-curve with a metric higher than the threshold.

Fig. 4

L-curves are shown for the first three iterations of a single case with a malignant tumor. The regularization was varied from 0.001 to 100 at each iteration, and the optimal regularization was (a) 0.18, (b) 0.22, and (c) undefined for the three iterations, respectively.

The number of iterations used in the reconstruction was also an important factor. In Fig. 5(a), the projection error, which represents the residual of the inversion equation, is plotted on a log scale as a function of the number of iterations and decreases with an increasing number of iterations. The largest decrease in the projection error occurred at the first iteration, which is typical behavior for a Newton-type iterative method. The change in the projection error by the ninth iteration was relatively small, resulting in convergence of the tumor-to-adipose contrast as shown in Fig. 5(b). In this study, nine iterations were chosen as the stopping criteria for image reconstruction of all patients.

Fig. 5

(a) Log scale of projection error versus number of iterations. (b) Tumor-to-adipose contrast in HbT versus number of iterations.

Table 1 summarizes the statistical and diagnostic results in terms of HbT and TOI contrast using fixed regularization values of 0.1 and 1, and optimal regularization based on L-curve analysis. Data from 22 patient exams which had an optimal regularization in the first iteration were included in the analysis. A significant difference ($p<0.05$) was found in HbT contrast for malignant versus benign patients for all three regularization selections.

Table 1

Comparison of statistics using fixed regularization of 0.1 and 1, and optimal regularization. Only amplitude data (AMPL) was used for image reconstruction.

Both amplitude and phase data were used in the image reconstruction to see if the involvement of phase data obtained at the six FD domain wavelengths further improved the accuracy of the optical image reconstruction.

Figure 6 is organized the same as Fig. 3, except that the reconstructions utilized both amplitude and phase data, and the three clinical exams that did not have an apparent optimal regularization value in the first iteration were excluded. In Figs. 6(a) and 6(b), the tumor-to-adipose contrast in HbT was plotted as a function of regularization for benign and malignant patients, respectively. When regularization decreased, the contrast increased in most cases. As a result, better separation occurred between the malignant and benign groups, although with a higher standard deviation at the lower regularization of 0.1 as shown in Figs. 6(c) and 6(d). Here, the HbT contrast in the benign group was more sensitive to the change in regularization for low values when compared with the results in Fig. 2(a). Apparently, determining an appropriate regularization for image reconstruction using both amplitude and phase data is even more important than when using only amplitude measurements.

Fig. 6

Tumor-to-adipose contrast of HbT versus regularization for (a) benign conditions and (b) malignant tumors. Circles have a regularization of 0.1, and asterisks have a regularization of 1. Box plots of contrast are shown for regularizations of (c) 0.1, and (d) 1 for all patients. Both amplitude and phase data were used in the image reconstructions.

We also applied the optimization algorithm (Fig. 1) during image reconstruction with both amplitude and phase data; but, L-curve analysis failed to find an optimal regularization. The compromise between model error and prior error was not significant on the L-curve at any point; hence, no optimal regularization parameter could be found that satisfied the criterion of the optimization algorithm. One approach to improve the L-metric might be to use separate regularization parameters for amplitude and phase data. In practice, we found that a fixed regularization of 2 for phase data, and an optimal regularization based on L-curve analysis for amplitude data worked well.

Table 2 summarizes the statistical results of 22 patient exams, in which the optimal regularization (on amplitude only, phase fixed at 2) was compared with fixed regularization values of 0.1 and 1, when both amplitude and phase data were used in the image reconstruction. Reconstruction using a fixed regularization of 0.1 for amplitude and 2 for phase data was compared as well. Significant differences ($p<0.05$) exist between the malignant and benign groups in terms of both HbT and TOI for either a fixed regularization of 1 or optimal regularization. Optimal regularization provides the highest AUC for TOI (0.94) among all the regularization choices. Meanwhile, the optimization method still gives a relatively high AUC for HbT, without a significant difference compared with the fourth approach. To summarize, the optimization technique provides a systematic and automated way to find the optimal regularization parameter in each iteration, which gives the best separation between malignant and benign groups in terms of recovered optical parameter TOI. The L-curve-based optimization technique utilized in this paper aims at finding the tradeoff between prior error and model error in terms of Hb and deoxy-Hb, water, and lipids. As a result, the increase of AUC for TOI is more obvious than HbT compared with fixed regularization, since TOI ($\mathrm{HbT}\times \mathrm{Water}/\mathrm{Lipid}$) is defined such that it represents all the recovered chromophore concentrations. Interestingly, the addition of phase information appears to degrade the statistical performance of the image reconstructions in the fixed regularization cases.

Table 2

Comparison of statistics using fixed regularization of 0.1 and 1, optimal regularization, and fixed regularization of 0.1 for amplitude and 2 for phase. Both amplitude and phase data (AMPL/PH) were used for image reconstruction.

Instead of using a single predictor, either HbT or TOI, we evaluated the combination of HbT and TOI as an indicator of malignant versus benign contrast-enhancing MRI regions of interest. Specifically, we applied both HbT and TOI as predictors, and fit a logistic regression using pathology-confirmed malignancy to obtain the outcome. The predicted score from the logistic regression model was used to construct the ROC curve.²⁶ We repeated this procedure for three different regularization parameters and compared their AUCs with a bootstrapping method.²⁷ As shown in Fig. 7, under these conditions, the optimal regularization improved the AUC (94.4%), relative to the fixed regularizations of 0.1 (88.2%) or 1 (84.4%).

Fig. 7

Receiver-operating characteristic (ROC) curve for fixed regularization of 1 (black) and 0.1 (blue), and optimal regularization (red), when HbT and TOI are combined. Both amplitude and phase data were used for image reconstruction.

Figure 8 presents the boxplots of tumor-to-adipose contrast for (a) HbT, (b) ${\mathrm{StO}}_2$, (c) TOI, (d) SA, and (e) SP by applying the optimization algorithm during image reconstruction with amplitude and phase data. HbT was the most significant indicator for differentiating the malignant and benign groups, and provided the largest mean difference in tumor-to-adipose contrast in the malignant and benign groups ($1.55\times$ versus $0.89\times$). A significant difference in TOI and SP contrast was also observed. The average contrast in both HbT and TOI was significantly higher in the malignant group than in the benign group. No significant difference in the ${\mathrm{StO}}_2$ contrast was found.

Fig. 8

Box plots of the contrast for: (a) HbT, (b) ${\mathrm{StO}}_2$, (c) tissue optical index (TOI), (d) scattering power (SA), and (e) scattering power (SP), as recovered using the optimal regularization and amplitude and phase data. ^aSignificant difference.

4.

Discussion

The results presented in this paper show that the contrasts recovered in both HbT and TOI were diagnostically significant, but also depended on the choice of regularization parameter. Thus, an objective methodology to select/identify the regularization parameter for reconstructing data from individual patient exams is critical for practical application of the combined MRI/NIRST imaging approach in the clinical diagnostic setting. Instead of applying an empirical value for every subject exam, the optimization algorithm sought an exam-specific regularization parameter (or series of parameters) for image reconstruction that was derived from the measured optical data. The optimization algorithm developed and applied here was based on L-curve analysis, which is widely used to stabilize ill-conditioned problems by balancing the model error with the prior error as constrained by the prior information. However, previous approaches of finding the “elbow” on the L-curve have been primarily limited to theoretical or simulation studies with few efforts being based on actual patient data. To deal with an L-curve generated from clinical exam data, we defined a metric describing the characteristics of the resulting L-shape, and developed an optimization approach based on least-square fitting of the response. The L-metric determined whether the L-curve was sufficient for selecting an optimal regularization parameter, and if so, the optimal value from the curve was utilized. Otherwise, the value from the previous iteration remained. The optimization algorithm appeared to be robust and effective in the clinical dataset we applied. We found that the average L-metric of all exams decreased with increasing number of iterations (although the number of iterations was set to 9 in all cases, i.e., we did not apply the L-metric or L-curve as a stopping criterion). Within four iterations, the projection error and measurement noise amplification were balanced and the contrast converged. Table 1 showed the comparisons of the diagnostic statistics resulting from three choices of regularization parameters. Moreover, although ${\mathrm{StO}}_2$ or scattering properties were not found to be significant indicators when taken alone, their values when used to form TOI contributed to the optimal indicator for differentiation of malignant from benign lesions.

When the regularization parameter decreased, less dependence in the recovered optical solution occurred with the prior information. The prior estimate was obtained under the assumption that the whole breast had the same optical properties. As a result, higher regularization will lower the contrast between tumor and background leading to a smaller separation in the mean HbT and TOI contrast between the malignant and benign diagnostic groups at a regularization of 1. On the other hand, a high regularization also suppresses the high spatial frequency variation in the recovered optical property distribution resulting from measurement noise, and produced lower variation in HbT and TOI within the malignant and benign groupings. A compromise between separation of the mean value and noise driven variability was obtained by applying the optimal regularization approach during image reconstruction.

We also investigated the effects of optimal regularization on the recovered optical image when both amplitude and phase data were used (relative to only amplitude measurements). In this case, the regularization algorithm failed to find an optimal value based on L-curve analysis, apparently because the relative contribution of phase data noise was higher than amplitude data noise, requiring much higher regularization of the phase data. Here, separate regularization values for phase and amplitude data improved the resulting image outcomes. Specifically, improvement occurred when regularization of the amplitude data was obtained through L-curve analysis, and regularization of the phase data was fixed at 2 ($\sim 1$ order of magnitude higher than the value used in the amplitude regularization). With this approach, improvement in diagnostic performance in terms of higher AUC for both HbT and TOI occurred relative to using fixed regularization. When the phase data were added to the Jacobian matrix for image reconstruction, the optimal regularization (with separate but fixed phase regularization) achieved maximal separation of the malignant from benign diagnoses (highest AUC). Finally, the optimal regularization generated the best AUC value (0.94) relative to the other regularization choices considered when HbT and TOI were combined as the diagnostic indicator (Fig. 7).

Although the best diagnostic performance occurred using both amplitude and phase data, the three regularization approaches, either fixed or optimal, were still able to separate the malignant and benign groups in terms of HbT ($p<0.05$), when applying only amplitude data for image reconstruction. These results suggest the possibility of simplifying the MRI/NIRS system into one with only CW channels, which would significantly reduce costs while not sacrificing much in terms of diagnostic performance.

In the results presented here, data acquired in 22 out of 25 patient exams had an optimal regularization parameter identified by L-curve analysis in the first iteration and were included in the statistical analyses. For completeness, we compared the outcomes of the optimization algorithm in all 25 patients in Table 3 in terms of the $p$ value and AUC for HbT and TOI contrast when applying the three regularization approaches. Both amplitude and phase data were used for image reconstruction in these results. For the subset of selected patients (22), optimal regularization provided the best AUC for both TOI and HbT relative to fixed regularization, but only for HbT when all exam data were evaluated. Thus, the L-curve approach may be a pragmatic way to identify which clinically acquired datasets are reliable. NIRS data can be compromised by a number of factors during their acquisition, such as fiber-tissue coupling or reflections; hence, the ability to objectively and conclusively determine which datasets will not lead to accurate spectroscopic parameter recovery could be important in clinical practice.

Table 3

Comparison of the three regularization approaches in all clinical exams, relative to selected exams when the optimal regularization parameter was found. The group of all exams included 16 malignant and 9 benign pathology-confirmed diagnoses, whereas the selected exams included 15 malignant and 7 benign cases (three exams from the former did not have optimal regularization at the first iteration). Both amplitude and phase data were used during image reconstruction.

5.

Conclusion

In this work, a robust optimization algorithm for selection of the regularization to be applied during image reconstruction based on exam-specific data acquired during MRI/NIRST examination of the breast was developed. The optical contrast values for HbT, ${\mathrm{StO}}_2$, TOI, and scattering parameters were estimated by applying the optimization algorithm when amplitude only and both amplitude and phase data. A statistical difference ($p<0.05$) occurred between malignant and benign groups for both absorption-derived contrasts (HbT and TOI) as well as reduced-scattering derived contrast (SP).

To the best of our knowledge, these results represent the first time an extensive study of regularization has been conducted on a relatively large amount of clinical breast exam data with the MRI/NIRST multimodality imaging approach. The optimization algorithm better differentiated malignant from benign cases compared to a fixed regularization parameter. The best diagnostic performance occurred with optimal regularization values selected from the individual’s exam data, and when combining HbT and TOI estimated from both amplitude and phase data as the diagnostic indicator.