An Automated Preprocessing Method for Diffuse Optical Tomography to Improve Breast Cancer Diagnosis

Abstract

The ultrasound-guided diffuse optical tomography is a noninvasive imaging technique for breast cancer diagnosis and treatment monitoring. The technique uses a handheld probe capable of providing measurements of multiple wavelengths in a few seconds. These measurements are used to estimate optical absorptions of lesions and calculate the total hemoglobin concentration. Any measurement errors caused by low signal to noise ratio data and/or movements during data acquisition would reduce the accuracy of reconstructed total hemoglobin concentration. In this article, we introduce an automated preprocessing method that combines data collected from multiple sets of lesion measurements of 4 optical wavelengths to detect and correct outliers in the perturbation. Two new measures of correlation between each pair of wavelength measurements and a wavelength consistency index of all reconstructed absorption maps are introduced. For phantom and patients’ data without evidence of measurement errors, the correlation coefficient between each pair of wavelength measurements was above 0.6. However, for patients with measurement errors, the correlation coefficient was much lower. After applying the correction method to 18 patients’ data with measurement errors, the correlation has improved and the wavelength consistency index is in the same range as the cases without wavelength-dependent measurement errors. The results show an improvement in classification of malignant and benign lesions.

Introduction

Diffuse optical tomography (DOT) is a noninvasive medical imaging technique capable of providing functional information of the tissue. This technology utilizes light in the near-infrared (NIR) window to assess the interior optical parameters of the tissue, which are closely related to the oxygenated and deoxygenated hemoglobin content within the tissue. ^{1 –5} Because of its high sensitivity and low cost, DOT applications are growing in many areas, such as functional brain imaging, breast cancer detection and treatment monitoring, and many others. ^{6 –11}

The intense NIR light scattering in tissue imposes a significant challenge to DOT on its spatial resolution and localization accuracy. Coregistration approaches of DOT with high-resolution imaging modalities, such as X-ray, magnetic resonance imaging, or ultrasound (US), have been investigated. ^{9,12 –15} In these approaches, a lesion is localized and the information of lesion size, shape, and depth obtained from a high-resolution imaging modality can be used to guide the DOT image reconstruction. ^{16 –20} The US-guided DOT approach developed by our group has improved lesion location uncertainty and reconstruction accuracy of optical parameters. ^{21 –23} In recent years, a field of shape-based reconstruction techniques has emerged where a priori information about the topology, the approximate location, and shapes of the unknown subdomains and their optical properties can be incorporated in the inversion directly. ²⁴ In the US-guided DOT approach, a semiautomated level set method that extracts tumor information from US images for DOT image reconstruction was introduced. ²⁵

Previous studies have shown the capacity of US-guided DOT in classifying malignant lesions from benign masses. ^26,27 However, there are number of challenges for the wide clinical acceptance of this technology. First, the contralateral reference breast measurements used to compute weight matrix and perturbation, which is the normalized difference between lesion and reference measurements, may not be homogeneous. This can reduce the accuracy of fitted background optical properties and then computed weight matrix, which in turn causes inaccuracy in reconstructed lesion optical properties. Second, measurement errors due to low signal to noise ratio (SNR) data at longer source and detector distances and bad coupling between skin and breast tissue are wavelength dependent. For example, SNR of wavelength 740 nm at longer source and detector distances is much lower than longer wavelengths due to higher light absorption of darker skin and skin pigment. ^{28 –30} Additionally, our DOT system uses a handheld probe which is placed on top of patient’s breast while the patient is in a supine position. Movements of operator’s hand or patient may lead to a bad coupling between the light guides and the breast, which may result in outliers in measurements and then in the calculated perturbations. These outlier measurements are random and significantly affect the accuracy of the reconstructed lesion optical properties.

A variety of experimental and modeling approaches have been developed for system calibration of optical source strengths and detector gains, source and detector (optodes) position errors, and coupling errors between skin and optodes. ^{31 –39} In general, source strengths and detector gains can be accurately estimated from homogeneous phantom measurements at each experiment and used to compensate tissue measurements. The compensation can be performed independently from image reconstruction ^31,33 or as part of the inverse problem for image reconstruction. ^32,34 These methods compensate system-related parameters and do not compensate tissue-caused wavelength-dependent measurement errors. Model-based approaches to compensate coupling errors between skin and optical optodes as well as optode position errors have been proposed in several studies. ^{35 –37} These methods model coupling errors as unknowns and include them in the reconstruction of tissue optical properties either sequentially or simultaneously. The advantage of these methods is the adaptive estimation of the coupling errors; however, the methods require the unknown coupling coefficients minimally vary from their initial estimates. Other approaches include the use of the differences between measurements at 2 separate wavelengths to reduce the coupling errors. ³⁸ A nonlinear approach for difference imaging was studied on how this approach tolerates modeling errors like domain truncation, source and detector coupling errors, and domain shape errors. ³⁹ The challenges we are facing are wavelength-dependent measurement errors caused by low SNR and wavelength-dependent skin and optode coupling errors, which can vary depending on patients’ bulk tissue optical properties, skin conditions, and operator’s hand or patient motion during data acquisition. An effective and robust data processing method is needed to remove these measurement errors while preserving the wavelength-dependent lesion optical properties. The method also needs to be automated to minimize the user interface and facilitate clinical translation. In our early approach, the individual source strengths and detector gains were estimated using a least-square method and compensated before fitting each patient’s background optical properties from the contralateral normal breast measurements. ³¹ The perturbation approach, which is the difference between lesion and contralateral normal tissue measurements normalized to the contralateral measurements, was used in imaging reconstruction. This approach has canceled out the unknown source strengths and detector gains. In a recent attempt by our group, Vavadi and Zhu et al ⁴⁰ have introduced a statistical method to automatically remove outliers from contralateral normal breast measurements. This method utilizes multiple sets of reference measurements to produce a robust set of reference. However, in many clinical cases, the outliers due to measurement errors are present in the lesion measurements. Lesion measurements contain wavelength-dependent information and are expected to be more heterogeneous than the reference measurements. To separate the measurement errors from lesion heterogeneity, more information from multiple wavelength measurements can be incorporated in the preprocessing before image reconstruction.

In this study, we introduce a new automated approach for data filtering based on multiple wavelength measurements collected at lesion site. The method combines data collected from multiple sets of lesion measurements to detect and correct outliers caused by wavelength-dependent measurement errors in the perturbation. This method represents an important step toward a fully automated DOT system for its clinical translation. The hypothesis is that there is a strong correlation between multiple wavelength measurements in NIR spectrum collected from the same lesion site and that this correlation has significantly decreased with the presence of wavelength-dependent measurement errors. The hypothesis was tested with phantom and patients’ data.

Method

Ultrasound-Guided DOT System and Data Acquisition

Phantom experiments and clinical studies were performed using our US-guided DOT system. This system consists of 4 laser diodes of wavelength 740, 780, 808, and 830 nm and 10 parallel photomultiplier (PMT) detectors. Laser diodes were modulated at 140 MHz, and the light at each wavelength was sequentially delivered to 9 positions on a handheld probe through optical fibers. Ten light guides couple the reflected light from tissue to 10 parallel PMT detectors simultaneously. As a result, each wavelength data set has 90 measurements. The details of the NIR system can be found elsewhere. ¹⁹ To implement statistical tests for outlier detection, we repeated the data acquisition 3 times for each wavelength at the same lesion location, which results in 12 data sets of 90 measurements each.

Automated Preprocessing Approach

Using the US-guided DOT system, multiple sets of measurements are collected from a lesion site and a contralateral normal breast (a reference site) to calculate normalized perturbation. For each wavelength, the normalized perturbation (U_λ) is given as:

${\text{U}}_{\lambda }\left(\text{m}\right)=\frac{{\text{A}}_1{\text{e}}^{j{\varnothing }_1}\left(\text{m}\right)-{\text{A}}_{\text{r}}{\text{e}}^{j{\varnothing }_{\text{r}}}\left(\text{m}\right)}{{\text{A}}_{\text{r}}{\text{e}}^{j{\varnothing }_{\text{r}}}\left(\text{m}\right)},$

1

where A₁ (m) and ∅₁ (m) are, respectively, measured amplitude and phase from the mth source–detector pair at lesion site, while $A_{\text{r}}\left(m\right)$ and ${\varnothing }_{\text{r}}\left(m\right)$ , respectively, are the measured amplitude and phase from the mth source–detector pair at reference breast. This approach can cancel out the unknown source strengths and detector gains. ³¹

As noted earlier, outliers in the normalized perturbation are mainly caused by measurement errors due to wavelength-dependent low SNR and skin–fiber coupling. In this method, perturbation outliers are detected and corrected by combining measurements of multiple wavelengths. A block diagram of the method is given in Figure 1. This figure illustrates that multiple data sets are acquired from a lesion site and a contralateral site of the normal breast. For each data set, our system is capable of acquiring measurements at 4 different wavelengths: 740, 780, 808, and 830 nm. The perturbation is calculated for each wavelength measurement, which is the normalized difference between the examined lesion and the contralateral normal breast measurements. Next, the correlation coefficients between all wavelengths are calculated to find the weakest correlated wavelength measurement. An outlier detection method is applied to the weakest correlated wavelength measurement to detect the outliers and correct them. This process is followed until all outliers are corrected. The following text is a detailed explanation of the method.

Figure 1.

Block diagram of the proposed automated data preprocessing method.

In the first step, the normalized perturbations are used to calculate the correlation coefficients between each pair of wavelengths as ^41,42 :

$\text{Corr}\left({\text{λ}}_i,{\text{λ}}_j\right)=\frac{1}{\text{M}-1}\sum _{\text{m}=1}^{\text{M}}\left(\frac{\overline{{\text{U}}_{\text{λ}}{}_i\left(m\right)-{\text{μ}}_{{\text{λ}}_i}}}{{\text{σ}}_{{\text{λ}}_i}}\right)\left(\frac{\overline{{\text{U}}_{\text{λ}}{}_j\left(m\right)-{\text{μ}}_{{\text{λ}}_j}}}{{\text{σ}}_{{\text{λ}}_j}}\right),$

2

where ${\text{μ}}_{{\text{λ}}_{\text{i}}}$ and ${\text{σ}}_{{\text{λ}}_{\text{i}}}$ are the mean and standard deviation of the normalized perturbations of the wavelength ${\text{λ}}_{\text{i}}$ , respectively, and ${\text{μ}}_{{\text{λ}}_{\text{j}}}$ and ${\text{σ}}_{{\text{λ}}_{\text{j}}}$ are corresponding values of wavelength ${\text{λ}}_{\text{j}}$ . M is the total number of measurements for each data set. The correlation coefficients are calculated for both the real and imaginary parts of perturbation. Based on the average correlation coefficients of both real and imaginary parts, the measurement set with the weakest correlated wavelength with correlation value below 0.6 determined from phantom experiments is selected for outlier detection. If all correlation coefficients are above 0.6, the reconstruction will be performed for each wavelength and then total hemoglobin concentration (tHb) will be calculated.

In the second step, the method examines each source–detector pair of both real and imaginary perturbation in the selected wavelength to determine outliers. This is achieved by combining 4 wavelength data of 3 lesion data sets, so that for each source–detector pair we have 12 measurements. Using multiple data sets measured at each source–detector pair, we are able to apply statistical test to detect outlier measurements for each source–detector pair and eliminate them.

We use the maximum normed residual (MNR) method for outlier detection. The MNR is a statistical test used to detect outliers based on the largest absolute deviation from the sample mean. By calculating the t distribution with k − 2 degrees of freedom, a threshold value for each source–detector pair i can be obtained as:

${\text{G}}_{\text{Thershold}}\left(i\right)=\frac{\text{k}-1}{\sqrt{\text{k}}}\sqrt{\frac{{\text{t}}_{\frac{\alpha }{2\text{k}},\text{k}-2}^2\left(\text{i}\right)}{\text{k}-2+{\text{t}}_{\frac{\alpha }{2\text{k}},\text{k}-2}^2\left(\text{i}\right)}},$

3

where ${\text{G}}_{\text{Thershold}}\left(i\right)$ is the outlier threshold for ith source–detector pair, ${\text{t}}_{\frac{\alpha }{2\text{k}},\text{k}-2}\left(i\right)$ denotes the upper critical value of the t distribution with k − 2 degrees of freedom, and α represents the level of significance that determines the strictness of outlier removal procedure. By changing this value between 0 and 1, the total number of the outliers and the significance of these outliers removed from the database can be changed. To find the optimal value of α, the outlier removal process is performed for different significance level ranging from .01 to .5 and the optimal value is set to .05 based on visual examination of the removed outliers. This optimal value is selected in a way that the test only removes the significant outlier data. A G value is determined as an absolute deviation of the data point from mean value of the measurements and normalized by standard deviation. The data point corresponding to the maximum G value which has absolute deviation higher than the threshold is considered as an outlier and removed from the data set. The test is iterated until no further outliers are detected beyond the threshold. This test is done for real and imaginary parts of perturbations separately. The details of this method can be found elsewhere. ^{40,42 –44} If an outlier is detected for any source–detector pair, its value is corrected with the average value of other wavelengths of the same source–detector pair. Both the real and imaginary parts are updated.

Once outliers are corrected for the selected wavelength, the correlation coefficients are recalculated between wavelengths with the updated values. If there are still any other weakly correlated wavelengths measurements (correlation <0.6), the above procedure is iteratively followed until all wavelengths measurements are corrected. The corrected normalized perturbation is used for DOT image reconstruction.

Diffuse Optical Tomography Image Reconstruction

The US-guided dual-zone mesh method introduced by our group ^9,19 was utilized to perform the DOT image reconstruction. In this framework, we use the coregistered US image to divide the DOT imaging volume into a lesion region as a region of interest and a background region. A finer mesh is selected for the lesion and a coarse mesh for the background region. Thus, the method employs smaller voxel size for lesion and larger coarse voxel size for the background region. This technique notably reduces the total number of voxels with unknown optical properties and keeps it at the same scale as of the total measurements. As a result, the ill-posed DOT reconstruction problem has significantly improved. Born approximation of light propagation in the tissue is utilized for computing the weight matrix and an iterative optimization based on conjugate gradient method is used to compute the lesion distribution at each wavelength.

The reconstructed absorption maps of the 4 wavelengths (740, 780, 808, and 830 nm) are used to compute the tHb map directly using the summation of oxygenated hemoglobin and deoxygenated hemoglobin maps with extinction coefficients given in reference. ⁴⁵

Computation of Wavelength Consistency Index

To assess the improvement in consistency between reconstructed results obtained from different wavelength measurements, we calculate the absorption coefficient ${\text{μ}}_{\text{a}}$ consistency index ( ${\text{μ}}_{\text{a}}\text{CI}$ ) as:

${\text{μ}}_{\text{a}}\text{CI}=\frac{1}{\text{V}}\sum _{\text{p}=1}^{\text{V}}\left\{\frac{1}{\text{W}\left(\text{W}-1\right)/2}\left[\sum _{i=1}^{\text{W}}\sum _{j=i+1}^{\text{W}}|{\text{μ}}_{\text{a}}^{\text{λ}i}\left(\text{p}\right)-{\text{μ}}_{\text{a}}^{\text{λ}j}\left(\text{p}\right)|\right]\right\},$

4

where “V” is number of voxels and “W” is number of wavelengths. This index is calculated voxel by voxel for absorption coefficient µ_a(p) of each pair of wavelengths. The value of this index is used as an indicator to compare the improvement in the consistency between wavelengths before and after the application of the method. Since this index calculates the absolute difference between absorption maps at different wavelengths, a low index value is an indicator of high absorption consistency between wavelengths. Likewise, a high index value is an indicator of low absorption consistency between wavelengths.

To evaluate the range of ${\text{μ}}_{\text{a}}\text{CIs}$ obtained from breast lesions without evidence of wavelength-dependent measurement errors, we have applied the correction method to 12 patients’ data. For each patient’s data, the correlation coefficient between each pair of wavelength measurements was calculated and found to be larger than 0.6, and ${\text{μ}}_{\text{a}}\text{CIs}$ before and after correction were calculated. Result was used as a reference to characterize breast lesion wavelength-dependent absorption heterogeneity.

Results

Phantom Study

We performed phantom study to evaluate the correlation coefficients between measurements of multiple wavelengths at the same target site. In phantom experiments, we used Intralipid solution as a reference to calculate the normalized perturbation. The Intralipid solution is a homogenous medium, which eliminates any reference heterogeneity effect on the calculated perturbation. Moreover, the DOT handheld probe was fixed during data acquisition in the phantom experiments. Outliers caused by patient’s or user’s movements of the handheld DOT probe during data acquisition did not exist in the phantom data.

A phantom target with calibrated value of μ_a = 0.23 cm⁻¹ and diameter 3 cm was used for experiments. The target was merged in the Intralipid solution (which calibrated as μ_a =.03 cm⁻¹) at different depths of 2.0 to 3.5 cm in 0.5 cm increments. The same Intralipid solution was used to acquire the reference measurements. Table 1 shows the correlation coefficients between measurements of multiple wavelengths for the target. Here, the presented correlation coefficients are the average correlation values of the real and imaginary parts of the perturbation. The table also shows the reconstructed target absorption coefficients at different depths. The phantom experiments show strong correlation coefficients between measurements of different wavelengths. The correlation coefficients decrease with depths due to lower SNR of the measurements, but all coefficients are above 0.6. Also, we see that reconstructed absorption coefficients are consistent among all wavelengths (Table 1). However, consistency between wavelengths also decreases with depth; the ${\text{μ}}_{\text{a}}\text{CI}$ does not exceed 0.06 in the reconstructed absorption profiles obtained at different depth. Note that the lower sensitivity at the shallower depth 2 cm (target center depth) is related to the lack of shorter source and detector pairs across the center of the probe where the US transducer is housed. This low sensitivity can be calibrated.

Table 1.

The Correlation Coefficients, Maximum Reconstructed Absorption, and ${\text{μ}}_{\text{a}}\text{CI}$ Values for the Phantom Target at Different Depths.

	Depth (cm)
	2	2.5	3	3.5
Correlation coefficients
Corr (740, 780)	0.9561	0.9223	0.7706	0.7130
Corr (740, 808)	0.9004	0.9564	0.9545	0.7063
Corr (740, 830)	0.9884	0.9728	0.9699	0.8797
Corr (780, 808)	0.9494	0.9119	0.7803	0.7409
Corr (780, 830)	0.9252	0.9738	0.8411	0.7711
Corr (808, 830)	0.9798	0.9660	0.9654	0.8720
Maximum reconstructed μa (cm−1)
λ = 740 nm	0.1093	0.1744	0.2435	0.2230
λ = 780 nm	0.0931	0.1565	0.2280	0.1994
λ = 808 nm	0.0832	0.1628	0.2363	0.2124
λ = 830 nm	0.0894	0.1767	0.2548	0.2283
${\text{μ}}_{\text{a}}\text{CI}$	0.0137	0.0159	0.0269	0.0476

Clinical Study

In this study, we have examined a total of 28 patients. These patients are divided into 2 groups based on the wavelength correlation analysis and ${\text{μ}}_{\text{a}}\text{CI}$ values. Sixteen patients’ data were wavelength inconsistent cases, where the wavelength correlation coefficients were low and ${\text{μ}}_{\text{a}}\text{Cis}$ were high. Based on biopsy results, 8 of these patients had malignant lesions and 8 had benign lesions. The other group includes 12 wavelength consistent cases, where the wavelength correlation coefficients were high and ${\text{μ}}_{\text{a}}\text{CI}$ was low. The study was approved by institutional review boards and all patients signed the informed consent. The data used for this study were deidentified.

In Figure 2, we present the real and imaginary plots versus source–detector distance of normalized perturbation profiles for all 4 wavelengths of a wavelength consistent case with no preprocessing. Both real and imaginary parts of each wavelength show similar trend with other wavelengths, with no visible outliers. The range of the real and imaginary parts of the perturbation is similar across 4 wavelengths.

Figure 2.

Real and imaginary plots of normalized perturbation (equation 1) of a wavelength consistent case with no data preprocessing. x-Axis is the source and detector distance in centimeter.

Table 2 shows the correlation coefficients, the maximum reconstructed absorption coefficients of all wavelengths, and ${\text{μ}}_{\text{a}}\text{CI}$ for the wavelength consistent case presented in Figure 2. In this case, we see higher correlation between all wavelengths and the maximum reconstructed absorption coefficients of all wavelengths are in the same range. In addition, ${\text{μ}}_{\text{a}}\text{CI}$ shows low value, which indicates high consistency in ${\text{μ}}_{\text{a}}$ between wavelengths.

Table 2.

The Correlation Coefficients, Maximum Reconstructed Absorption, and ${\text{μ}}_{\text{a}}\text{CI}$ of a Wavelength Consistent Case.

	No Preprocessing
Correlation coefficients
Corr(740, 780)	0.91
Corr (740, 808)	0.89
Corr (740, 830)	0.79
Corr (780, 808)	0.92
Corr (780, 830)	0.89
Corr (808, 830)	0.84
Maximum reconstructed μa (cm−1)
λ = 740 nm	0.225
λ = 780 nm	0.216
λ = 808 nm	0.237
λ = 830 nm	0.226
${\text{μ}}_{\text{a}}\text{CI}$	0.023

We have applied the automated preprocessing method to all cases. For wavelength inconsistent cases, measurements with outliers are seen to have weak correlation coefficients with others. An example of a normalized perturbation of a malignant case is illustrated in Figure 3. This figure shows the calculated perturbation profile for each wavelength versus source–detector distance. The wavelength 740 nm is weakly correlated with the other wavelengths, with correlation coefficients of less than 0.6, as summarized in Table 3. This lower correlation value is due to outliers observed in both real and imaginary parts of the 740 nm wavelength shown in Figure 3A. When the automated method is applied, outliers are corrected (Figure 3B) for the 740 nm wavelength and the correlation coefficients are improved (Table 3).

Figure 3.

Real and imaginary plots of normalized perturbation (equation 1) of a malignant lesion with (A) no data preprocessing and (B) after applying the proposed automated method. x-Axis is the source and detector distance in centimeter.

Table 3.

The Correlation Coefficients, Maximum Reconstructed Absorption, and ${\text{μ}}_{\text{a}}\text{CI}$ Values Before and After Automated Method for a Malignant Case.

	No Preprocessing	Automated Preprocessing Method
Correlation coefficients
Corr (740, 780)	0.50	0.81
Corr (740, 808)	0.39	0.75
Corr (740, 830)	0.79	0.82
Corr (780, 808)	0.91	0.91
Corr (780, 830)	0.74	0.74
Corr (808, 830)	0.73	0.73
Maximum reconstructed μa (cm−1)
λ = 740 nm	0.52	0.33
λ = 780 nm	0.36	0.36
λ = 808 nm	0.33	0.33
λ = 830 nm	0.35	0.35
${\text{μ}}_{\text{a}}\text{CI}$	0.10	0.04

The reconstructed absorption maps of this case are illustrated in Figure 4. Both the reconstructed lesion shape and the maximum reconstructed absorption value for the wavelength 740 nm are not consistent with the other wavelengths (Figure 4A and Table 3). After applying the automated method, we see the shape of this wavelength becomes more consistent with other wavelengths (Figure 4B), and ${\text{μ}}_{\text{a}}\text{CI}$ reduced from .10 to .04. Likewise, the correlation coefficients between all wavelength measurements are improved (Table 3).

Figure 4.

The reconstructed absorption maps of a malignant lesion (A and B) along with the ultrasound B scan of the lesion (C). Reconstructed absorption coefficients are shown for 4 wavelengths with (A) no data preprocessing and (B) proposed automated method. Each image marked with wavelength represents a spatial image in x and y dimensions of 9 cm scale at the target depth of 2 cm, other depths are not shown.

A benign case is also presented in Figure 5. The wavelength correlation coefficients and the reconstructed absorption coefficients are summarized in Table 4. Based on the wavelength correlation analysis, the wavelength 808 nm is weakly correlated with the other wavelengths. Table 4 summarizes the average correlation coefficients of real and imaginary parts of the perturbation between multiple wavelengths. With no perturbation preprocessing, the reconstructed absorption map at 808 nm is too high and distorted as compared to other wavelengths. When the automated method is used, the ${\text{the μ}}_{\text{a}}\text{CI}$ has decreased from 0.14 to 0.05. The maximum reconstructed absorption coefficient of the wavelength 808 nm became more consistent with the rest of the wavelengths and decreases from 0.27 to 0.13 cm⁻¹. Besides that, the tHb decreased from 68 to 55 μmol/L when the automated method was used.

Figure 5.

The reconstructed absorption map of a benign lesion (A and B) along with the ultrasound B scan of the lesion (C). Reconstructed absorption coefficients are shown for 4 wavelengths with (A) no data preprocessing and (B) proposed automated method.

Table 4.

The Correlation Coefficients, Maximum Reconstructed Absorption, and ${\text{μ}}_{\text{a}}\text{CI}$ Values Before and After Automated Method for a Benign Case.

	No Preprocessing	Automated Preprocessing Method
Correlation coefficients
Corr (740, 780)	0.72	0.72
Corr (740, 808)	0.33	0.77
Corr (740, 830)	0.76	0.76
Corr (780, 808)	0.19	0.67
Corr (780, 830)	0.79	0.79
Corr (808, 830)	0.20	0.81
Reconstructed μa (cm−1)
λ = 740 nm	0.148	0.148
λ = 780 nm	0.154	0.154
λ = 808 nm	0.276	0.137
λ = 830 nm	0.152	0.152
${\text{μ}}_{\text{a}}\text{CI}$	0.14	0.05

An example of 2 inconsistent wavelength measurements is presented in Table 5. In this example, the wavelengths 740 and 830 nm are weakly correlated with 2 other wavelengths (low correlation coefficients). The correlation between the wavelengths 780 and 808 nm is high. With no preprocessing, the reconstructed absorption coefficient is too high at 740 nm and too low at 830 nm. We applied the automated method to this case and correlation coefficients have improved for the wavelengths 740 and 830 nm. The maximum reconstructed absorption coefficient for 740 nm decreased from 0.35 to 0.16, and the maximum reconstructed absorption for 830 nm has increased to 0.14 and became more consistent with the rest of the wavelengths. Also, the wavelength consistency index ( ${\text{μ}}_{\text{a}}\text{CI}$ ) reduced from 0.12 to 0.03.

Table 5.

The Correlation Coefficients, Maximum Reconstructed Absorption, and ${\text{μ}}_{\text{a}}\text{CI}$ Values Before and After Automated Method for a Case With 2 Inconsistent Wavelengths.

	No Preprocessing	Automated Preprocessing Method
Correlation coefficients
Corr (740, 780)	0.21	0.74
Corr (740, 808)	0.35	0.79
Corr (740, 830)	0.36	0.71
Corr (780, 808)	0.87	0.87
Corr (780, 830)	0.47	0.85
Corr (808, 830)	0.45	0.89
Reconstructed μa (cm−1)
λ = 740 nm	0.354	0.167
λ = 780 nm	0.155	0.155
λ = 808 nm	0.156	0.156
λ = 830 nm	0.125	0.141
${\text{μ}}_{\text{a}}\text{CI}$	0.12	0.03

The average value of the ${\text{μ}}_{\text{a}}\text{CI}$ for all wavelength inconsistent clinical cases has reduced from 0.12 to 0.04 when the automated method is applied, as shown in Figure 6. For the wavelength consistent clinical cases, the average ${\text{μ}}_{\text{a}}\text{CI}$ with no preprocessing is 0.03 and has remained the same after applying the automated method (Figure 6). This indicates an improvement in the consistency of the reconstructed absorption coefficients between all wavelengths for the wavelength inconsistent cases.

Figure 6.

The calculated μ_aCI values for the 2 groups of clinical cases: the wavelength inconsistent cases and the wavelength consistent cases. The figure shows the index values of each group before and after the automated method is used.

The average wavelength correlation coefficients calculated for all 28 clinical cases are summarized in Table 6. Here, the wavelength correlation coefficients of the wavelength inconsistent cases are compared to those of the wavelength consistent cases. The results show the correlation coefficients of the wavelength inconsistent cases after applying the automated method are similar to the correlation coefficients of the wavelength consistent cases.

Table 6.

Average (Standard Deviation) of the Correlation Coefficients of the Clinical Cases.

	Wavelength Inconsistent Cases: No Preprocessing	Wavelength Inconsistent Cases: Automated Method	Wavelength Consistent Cases: No Preprocessing
Wavelength correlation coefficients	0.35 (0.22)	0.89 (0.09)	0.83 (0.22)

The reconstructed absorption coefficients of the 4 wavelengths (740, 780, 808, and 830 nm) are used to compute the tHb. The computed tHb for the wavelength inconsistent cases is presented in Figure 7, which shows the calculated tHb for malignant and benign cases before and after applying the automated method. A 2-sided t test was performed between malignant and benign groups before and after applying the automated method. The proposed method shows the statistical significance, where the P value improved from.006 to.001 using proposed automated method. The tHb ratio of malignant to benign group is improved from 2.01 to 2.42 with the automated preprocessing method.

Figure 7.

Total hemoglobin concentration of wavelength inconsistent clinical cases before and after the automated data preprocessing method. Red and blue boxes indicate the malignant and benign cases, respectively.

Discussion and Summary

We have introduced a new automated data preprocessing method based on multiple wavelength measurements. In the absence of wavelength-dependent measurement errors in phantom experiments, the data show strong correlations between measurements of multiple wavelengths in the NIR range. In addition, the phantom experiments show a high consistency between reconstructed images obtained from measurements of different wavelengths. Even though both the correlation and the consistency between wavelengths are reduced with target depth, correlation coefficients have maintained beyond 0.6. Similarly, the calculated ${\text{μ}}_{\text{a}}\text{CI}$ has not exceeded .06 for all target depths in the phantom experiments.

The automated preprocessing method was used to analyze 16 malignant and benign wavelength inconsistent cases. The results of this method were compared with those obtained with no data preprocessing. The wavelength correlation coefficients and the consistency indexes, ${\text{μ}}_{\text{a}}\text{CI},$ were further compared with those values obtained from 12 wavelength consistent cases. Both the correlation coefficients and the ${\text{μ}}_{\text{a}}\text{CI}$ of the wavelength consistent cases and the wavelength inconsistent cases after applying the automated method are in the same range (Table 6). Because the correlation coefficients for wavelength consistent cases were above 0.6 and correlation coefficients at all depths in the phantom study were also above 0.6, we have chosen 0.6 as a threshold value in this method. The automated method helped to improve the consistency between reconstructed absorption maps of different wavelengths for all malignant and benign cases (Figures 4 and 5). The high ${\text{μ}}_{\text{a}}\text{CI}$ values of wavelength inconsistent cases have been improved after applying the automated preprocessing method. These values are all decreased to the similar values of wavelength consistent cases and to phantom data ${\text{μ}}_{\text{a}}\text{CI}$ values of less than.06. This improvement in the ${\text{μ}}_{\text{a}}\text{CI}$ shows that the wavelength-dependent measurement errors have minimized using other wavelength measurements to improve the data correlation before performing reconstruction.

The lesion shape is also improved and become more consistent with the results of other wavelengths as we see in Figures 4 and 5. With no preprocessing, the reconstructed absorption coefficients of the benign case (Figure 5) at 808 nm are high and distorted as compared to the other wavelengths. The automated method demonstrated the ability to improve the consistency for this case. The shape of the lesion at 808 nm wavelength became more consistent as well.

These improvements were reflected in the tHb. The average value of the tHb for malignant cases increased by 10 μmol/L, while benign group showed similar improvement with tHb reduced by 12 μmol/L. As a result, the malignant to benign tHb ratio was improved by 20% after applying the automated method.

The automated method used in this study focuses only on the cases with one poorly correlated wavelength with other wavelengths, and all analyzed inconsistent cases are having problems in one wavelength. The results show strong improvement in terms of improving the wavelength consistency and robustness of the image reconstruction method. In cases that more than one wavelength measurements are distorted, the method can still be applied using the remaining consistent wavelength measurements to predict the distorted wavelengths and correct the distortion. However, our method requires that the wavelength with measurement outliers should be at most 2 if 4 optical wavelengths are used.

There are 4 major steps for a fully automated US-guided DOT system. First, the system calibration should be done automatically and the code has been implemented in our imaging algorithm. Second, the estimation of background optical properties of the contralateral breast and therefore the computation of the weight matrix should be done automatically after filtering out outliers in the measurements. Third, the perturbation, which is the normalized difference between the examined lesion and the background, should be automated after correcting wavelength-dependent errors and also eliminating other noise due to low SNR. Fourth, the US lesion measurements should be done automatically and these parameters can be inputted into the DOT imaging reconstruction for near-real-time US-guided DOT imaging. In this study, the automated approach for data filtering based on multiple wavelength measurements collected at lesion site is intended to automate the third step data processing, while the method introduced in Vavadi and Zhu ⁴⁰ was used to automate the second step data processing. The semiautomated method proposed by Mostafa et al ²⁵ is aimed to solve problems of the fourth step. The future work will combine these 4 steps of data preprocessing and validate and refine the combined approach in a large prospective patient study.

In conclusion, this study has introduced an automated data preprocessing method based on measurements of multiple wavelengths. The method eliminates the effect of wavelength-dependent measurement errors in the DOT perturbation, which helps to achieve more accurate reconstruction of optical properties of breast lesions. The automated method also helps to minimize both the user interface and the time for data preprocessing. The average time for the experienced user to manually perform data preprocessing is from 15 to 30 minutes for one patient’s data. The automated method could help to reduce this time to less than a minute. Although the method is demonstrated using US-guided DOT data, it is applicable to any DOT data preprocessing obtained with multiple wavelengths.