1. Introduction

According to the American Cancer Society, breast cancer is estimated to be the most common type of cancer among women with more than 246000 new cases and approximately 40,000 deaths projected in the United States in 2016 [1]. Mammography, ultrasound (US) and magnetic resonance imaging (MRI), are widely used to detect and diagnose breast lesions. Mammography is the most clinically used imaging modality for breast cancer screening, but it suffers from significant limitations such as relatively low sensitivity in early cancer diagnosis, reduced sensitivity in women with dense breasts, and low specificity that results in a large number of unnecessary biopsies [2]. A new study demonstrates comparable results for breast cancer detection rate between US and mammography. However, US still lacks sensitivity and specificity with a large number of false positive cases reported [3]. Although MRI has a higher sensitivity in breast cancer detection, the low cancer detection yields in the general screening population and high costs of MRI prohibit its applicability as a screening tool [2].

Diffuse optical tomography (DOT) is a noninvasive functional imaging modality that utilizes near-infrared (NIR) light to probe tissue optical properties. Minimal light absorption in the NIR wavelength range allows for several centimeters of light penetration in soft tissue such as breast and brain and facilitates probing lesions deep in the tissue. Image reconstruction is performed using reflected or transmitted light at the tissue surface [4–7]. NIR diffuse optical tomography has emerged as a promising modality to detect and monitor functional changes related to blood flow and tumor angiogenesis. Using multi-wavelength data acquisition, it is possible to quantify tissue characteristics such as oxygenated, deoxygenated, and total hemoglobin (HbO2, HbR, and HbT) concentrations as well as hemoglobin oxygen saturation (SO2) and lipid and water concentrations [8,9]. Clinical studies of female breasts have demonstrated higher HbT in the malignant lesions compared to both healthy tissue and benign lesions [8,10,12,14,15].

High sensitivity and cost effectiveness of DOT make it an ideal alternative for conventional breast imaging modalities. Currently DOT is not widely used in clinical settings mainly due to the low spatial resolution and lesion location uncertainty caused by intense light scattering in soft tissue. However, DOT guided by mammography, US, and MRI [10–13] have overcome lesion location uncertainty and improved the light quantification accuracy, and have been demonstrated as promising complementary methods to the existing imaging modalities. Ultrasound-guided DOT system and technique has been developed by our group and clinical results have demonstrated its potential in differentiating malignant and benign breast lesions and reducing the need for unnecessary benign biopsies [10,15]. In this approach, DOT is used to characterize the US detected lesions and to improve the diagnostic accuracy of US by providing complementary optical contrast. In a recent study, we have demonstrated that when US detection and diagnosis used together with optical hemoglobin contrast, the sensitivity of breast cancer detection and diagnosis has researched 96.6-100% [10].

The path toward commercialization of imaging guided DOT techniques depends upon the improvement of robustness and user-friendliness of this technique in hardware and software. Several groups have investigated various approaches to improve the robustness of the hardware, system calibration, data selection, target localization, and image reconstruction [16–20]. Reducing user interaction via automated processing of calibration and data selection is an important step toward clinical applicability of DOT.

In US-guided DOT or other hand-held DOT operation, a hand-held probe is placed on top of a patient’s breast while the patient is in a supine position [21]. Movements of patient or operator’s hand could cause a bad coupling between the light guides and the breast which may result in some outliers in measurements. In other imaging guided DOT approaches, the probe may be fixed, however, the bad coupling between tissue and light guides can occur which may result in outlier measurements. Additionally, tissue heterogeneity can cause measurement errors at some source-detector pairs. Recovered background and lesion optical properties highly depend on the boundary measurements of light propagating through the tissue underneath and any errors in these measurements could cause inaccuracy in fitted background and reconstructed lesion optical properties [22,23]. In this study, we introduce an automated outlier removal, data selection, and perturbation filtering method for imaging-guided DOT to improve its robustness in breast imaging. This method utilizes multiple sets of reference measurements acquired at a normal contralateral breast in one imaging session to produce a robust set of reference measurements. Then, the lesion measurement set is subtracted and scaled by this robust reference to form the normalized perturbation of the scattered field. Next, the perturbation was filtered for outliers and measurement errors based on analysis from simulations. Finally, the filtered perturbation is used for reconstruction of absorption maps and total hemoglobin distributions to characterize malignant and benign breast masses. The method of combining multiple measurements to perform statistical test for noise and outlier removal and saturation data rejection can be applied to other methods of imaging-guided DOT, such as mammogram or MRI guided diffuse optical tomography with some modifications.

2. Data set and system

2.1 Patients

US-guided DOT is used to characterize 10 malignant and 10 benign lesions of 20 patients. The study protocol was approved by the local internal review board (IRB). Signed informed consent was obtained from all patients and the study was compliant with the Health Insurance Portability and Accountability Act. For each patient, US images and optical measurements were simultaneously acquired before biopsy at the lesion site and at the normal contralateral breast of the same quadrant as the lesion which is used as reference. Mammograms, US images before the NIR scan were reviewed by attending radiologists and all contralateral measurements were taken at normal areas based on the available information [10].

2.2 System

The DOT system used in this study consisted of four laser diodes of wavelengths 740, 780, 808 and 830 nm. The outputs of the diodes were multiplexed by two optical switches to nine positions and were then coupled to a hand-held probe through optical fibers. Ten optical light guides coupled diffused light reflected from tissue to 10 photomultiplier tubes (PMT). The outputs of the laser diodes were amplitude modulated at 140 MHz and the detector outputs were demodulated to 20 kHz. The entire data acquisition took approximately 3 to 4 seconds, which was fast enough to repeat several times during one imaging session and obtain multiple measurements from each patient at the lesion site and the contralateral normal breast at the same quadrant as the lesion, referenced as reference measurement in the following text. For each patient, we typically collect 4 to 6 sets of measurements from the contralateral breast at slightly different locations. Each set consists of 3 consecutive measurements which takes around 4 to 5 seconds. The separation between each set of measurements varies from few seconds to tens of seconds because of hand-held operation. Therefore, a total 12 to 18 measurements may take about one to two minutes. A commercially available US transducer located in the middle of the hand-held probe of 9 cm diameter was used to acquire US images and optical source fibers and light guides distributed at the periphery were used to collect the diffused light from tissue [10].

3. Methods

Fast data acquisition of US-guided DOT system allows the collection of several data sets during one imaging session. The data contains multiple data sets from both lesion and contralateral reference breast. A novel automated method for outlier removal and data selection is introduced to eliminate the effect of inaccurate measurements. The block diagram of the procedures is given in Fig. 1.

 

Fig. 1 Block diagram of the automated outlier removal and data selection method. In the block diagram MNR stands for maximum normed residual and IRLS stands for iterative reweighted least square.

Download Full Size | PDF

Multiple sets of reference measurements were used to form a single high quality data set. First, an outlier removal procedure is incorporated to eliminate the highly inaccurate measurements with a criterion based on the distribution of data collected at each source-detector pair. Second, a piece-linear fitting is used to reject the source-detector pair measurements obtained from the saturated PMTs. Third, an iterative fitting of residue of the remaining data is calculated to further eliminate inaccurate measurements based on the linearity of fitted results of the reference measurements of all source-detector pairs. Fourth, a least-square error method is used to form the most accurate reference data set from the remaining measurements. Finally, a perturbation filtering based on analysis obtained from the semi-infinite analytical solution of light propagation in tissue [24] is utilized to form accurate perturbation set that is more robust to outlier and inaccurate measurements. The aforementioned steps are described in details below.

3.1 Outlier rejection in reference measurements

Each data set contains measurements from s sources and d detectors with the total number of m = s × d measurements. The system used in this study provides 90 source-detector measurements per data set. A total of k data sets of total k × s × d measurements collected at the reference site were used for selecting a best reference data set. In general, k is in the range of 12 to 18. Since we are using a frequency domain DOT system, each data set consists of amplitude and phase data. The maximum normed residual (MNR) is a widely used statistical test to address the problem of outlier rejection [25] that has shown outstanding performance in both linear and nonlinear data [26]. The MNR test is based on the largest absolute deviation from the sample mean in units of the sample standard deviation. This test is applied to remove outliers for each source-detector pair of total k data sets. The outlier measurement is expunged from each source-detector pair data set based on the following criterion. An upper critical value of the t-distribution with k−2 degrees of freedom is calculated, and a threshold is obtained based on Eq. (1).

In this equation G_Threshold(i) is the outlier threshold for i^th source-detector pair, t_α/(2k),k-2(i) denotes the upper critical value of the t-distribution with k−2 degrees of freedom, and α represents the level of significance which 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 have changed. To find the optimal value of α, the outlier removal process is performed for different significance level ranging from 0.01 to 0.5 and the optimal value is set to 0.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 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 both amplitude and phase measurements separately. The details of this method can be found in [27,28]. If a data point in either amplitude or phase measurements is considered as an outlier, both amplitude and phase measurements are removed from the data sets.

3.2 Saturation and noise data rejection

In addition to outliers in the reference measurements, detector saturation is another common problem in DOT that can happen as a result of higher light intensity detected at a shorter source detector distance. Each PMT may saturate at a different light intensity level. A semi-infinite analytic solution predicts that the logarithm of the detected amplitude for each source-detector pair multiplied by square distance of that specific source-detector pair, referred as logarithmic amplitude, should linearly decrease with the source-detector distance for homogeneous reference measurements. The phase measurement should increase linearly with the source-detector distance. A piece-linear fitting method is implemented for the amplitude measurements of all remaining source-detector pairs in the reference data after outlier rejection. In general, three sections of shorter source-detector distance, mid-range, and longer distance range were used. If a measured logarithmic amplitude at a shorter source-detector distance does not follow the linear profile plotted as a function of the source-detector separation in the mid and longer range, we can assume that the PMT is saturated at this detector distance. The measurements that fit the linear profile are kept for further processing. Additionally, the phase data corresponding to the saturated amplitude data are not reliable and are removed from both reference and lesion data sets.

Besides the saturated measurements which usually happen for the source-detector pairs with shorter distances, there are measurements of source-detector pairs with longer distances which are dominated by the noise of the system and are not reliable. To improve the data set, any measurements of longer source detector pairs with amplitudes below the electronic noise of the system are considered as noisy measurements and expunged from the data. Since the corresponding phase data with amplitudes at noise level are not reliable and these data have also removed from the data.

3.3 Iterative reweighted least square fitting

The MNR test is based on each source detector measurements separately and it removes the outliers at each source-detector pair. At this stage, all remaining data sets after MNR and saturation data removal are inputted to an iterative reweighted least square (IRLS) method to obtain an accurate linear fit with the minimum fitting residue for both log scaled amplitude and phase measurements as functions of source and detector separation. This method has given great results in different p-norm minimization problems ranging from compressive sensing to baseline correction in industrial settings [29,30].

As shown in Eq. (2) IRLS iteratively minimizes the bi-square weighted residual in the least-square sense [29]:

Here i is the index of the measurement, w is the bi-square weight function, y is the measurement value, β includes slope and intercept of the line fitted to the data and f(i,β) is the fitted measurement based on the current β . This method reduces the influence of large residuals in the fitted reference parameters and improves the fitted results. After the fitting for both logarithmic amplitude and phase is completed, the distance of each amplitude and phase measurement from the corresponding value on the fitted line of same source-detector distance is calculated. All the measurements with higher absolute residue compared to the threshold in either amplitude or phase measurements are selected as non-accurate measurements and removed from the data set. The absolute residue threshold of 0.5 has been empirically selected for both amplitude and phase based on trial and error from clinical data. Since the IRLS-based minimization is robust and less sensitive to noise bursts, it improves the robustness of the data set selection [30].

3.4 Compounded reference

Even though accurately fitted lines for both log scaled amplitude and phase are obtained from the previous steps, there may still be more than one measurement for each remaining source-detector separation. Now one amplitude and one phase measurement per source-detector pair need to be selected to form a single robust reference data set. Here, a least square method is utilized to select the measurements with minimum distance from the center of the distribution of remaining measurements for each source-detector pair. This processed is done separately on the remaining amplitude and phase data. Therefore, a final reference data set with high similarity to the fitted slope and intercept of outlier removed combination of all the reference data sets is achieved. This reference data set consists of the selected amplitude and phase measurements for the remaining source-detector pairs. This robust set of reference measurements is named as the compounded reference and is less sensitive to outliers, PMT saturation and noise.

To visualize the effect of the proposed process, the amplitudes and phase profiles of one clinical case before and after preprocessing, as well as the final compound reference have been illustrated in Fig. 2. The preprocessing includes, outlier removal, saturation and noise rejection and iterative reweighted fitting with higher residue removal.

 

Fig. 2 Log scaled amplitude and phase profiles of reference data sets before (first column) and after preprocessing (second column) as well as final compound reference (third column). Saturated source-detector pairs have been marked with red rectangle in the amplitude part of the first column.

Download Full Size | PDF

3.5 Perturbation filtering

The procedures described above provide a robust set of reference measurement. For reconstruction of lesion absorption map at each wavelength, perturbation is calculated by subtracting the compound reference from the lesion data as shown in Eq. (3), in which A_l and A_r are amplitude and φ_l and φ_r are phase of each source-detector pair obtained at the lesion and compound reference, respectively.

the outliers may present in the perturbation data due to 1) measurements errors caused by movements of patient or operator’s hand as well as bad coupling between the light guides and the breast; and 2) heterogeneity of the background tissue and the lesion. Lesion measurements are expected to be more heterogeneous than the reference measurements because the heterogeneity is partially caused by the lesion and partially by the background tissue heterogeneity. Therefore, outlier removal procedures applied to reference measurements cannot be implemented to the lesion measurements. Instead, a filtering method is applied to the perturbation based on the constraints imposed on the phase difference between the lesion data and the reference data given in Eq. (3). These conditions are determined based on the predictions obtained from the semi-infinite analytical solution derived from diffusion approximation [24, 31].

Simulations were performed using different background optical properties for both reference and lesion breasts as well as different optical properties for lesions of different sizes located at different depths [32]. The simulations used the same probe geometry and the same number of the sources and detectors as the experiments. Table 1 shows the range of the parameters used for simulations. The results show that the maximum phase difference between lesion and reference measurements of all source detector pairs for most scenarios listed in Table 1 are in the range of few degrees. For a 4 cm larger lesion of optical contrast of 10 times higher than the background, the maximum phase difference is only 22 degrees. This implies that the maximum phase difference of none of the source detector pairs can exceed 90 degrees even in extreme cases. This indicates that the cosine term in Eq. (3) should always be positive. Therefore, the real part of the perturbation cannot be less than −1 assuming that the amplitude measured from lesion $A_l$ is smaller than that of the amplitude measured at reference $A_r$ due to a higher absorption of the lesion in general. Furthermore, the mean and the standard deviation of the imaginary part of the perturbation are calculated. If the imaginary part of a data point is farther than three standard deviations from the mean, it is considered an outlier. Any data points in the perturbation that do not meet these two criteria are rejected. This step removes outliers in the perturbation likely caused by measurement errors rather than heterogeneity of the lesion data.

Table 1. Range of parameters used for analytical model and obtained maximum phase delay.

View Table | View all tables in this article

This normalized perturbation is used for reconstructing the absorption map at each wavelength. The total hemoglobin map was calculated from the four wavelength absorption data. The dual-zone mesh scheme is used for inversion [32]. Briefly, the imaging volume is segmented into two regions consisting of the lesion and the background regions. These two regions are identified by the co-registered ultrasound images. This method reduces the total number of voxels with unknown optical properties by using smaller mesh size for lesion region and a larger coarse mesh size for background. The conjugate gradient method is utilized for iterative optimization of the inverse problem. Patient results are calculated from the selected data based on this automated outlier removal and data selection procedures.

4. Results

To quantify the effectiveness of the proposed method, the results obtained from manually selected and processed data set were compared with the results obtained from the automated data processing method. The malignant versus benign classification was based on biopsy result which is the “gold” standard in current clinical practice. The manual selection was performed by a trained user to obtain the best outlier removal result. The absorption maps obtained at four optical wavelengths of benign and malignant examples are shown in Figs. 3 and 4, respectively. The effect of outlier and bad measurements on the reconstructed image with no preprocessing is clearly seen from the 740 nm image of first row in Fig. 4. Both manual and automated methods were able to improve wavelength consistency by removing outlier measurements from the data set. The automated method provides similar reconstruction results as the manually processed ones with only minor differences.

 

Fig. 3 Reconstructed absorption maps of a benign breast lesion obtained at 740, 780, 808 and 830 nm with no preprocessing (1st row), with manual data selection (2nd row) and automated data selection (3rd row). Each slide is 9 cm by 9 cm reconstructed at the center depth of the lesion. Images from other depths were not shown. Vertical bars are absoption in cm⁻¹ unit.

Download Full Size | PDF

 

Fig. 4 Reconstructed absorption maps of a malignant breast lesion at 740, 780, 808 and 830 nm with no preprocessing (1st row), with manual data selection (2nd row) and automated data selection (3rd row). Each slide is 9 cm by 9 cm reconstructed at the center depth of the lesion. Images from other depths were not shown. Vertical bars are absorption in cm⁻¹ unit.

Download Full Size | PDF

The corresponding ultrasound images and the calculated total hemoglobin maps obtained by both manual and automatic methods along with calculated total hemoglobin maps with no outlier removal and data selection are presented in Fig. 5 for the benign and malignant cases. The lower reconstructed absorption at 740 nm of the malignant case with no data selection caused wavelength inconsistency and reduced the calculated total hemoglobin concentration. It can be seen that both manual and automated methods have improved wavelength consistency and better malignant to benign separation by reducing the effect of saturated, noisy and outlier data. The HbT concentrations are also comparable between manual and automated methods.

 

Fig. 5 Ultrasound image (1st column) and total hemoglobin map calculated with no data selection (2nd column), with manual data selection (3rd column) and automatic data selection (4th column) of the benign (top row) and the malignant breast lesion (bottom row) reconstructed with automated data selection method. Each slide is 9 cm by 9 cm and vertical bars are hemoglobin concentrations in uM unit.

Download Full Size | PDF

To compare manual and automated methods statistically, both have applied to the data sets collected from 20 patients (10 malignant and 10 benign). The maximum reconstructed absorption coefficient with each method was calculated. The mean value and the standard deviation of the maximum reconstructed absorption coefficient for malignant and benign groups along with the ratio of the mean values are presented in Table 2. The proposed automated data selection and outlier removal method shows slightly higher malignant to benign ratio of all wavelengths due to improved data selection.

Table 2. Comparison of maximum reconstructed absorption coefficients of malignant and benign groups using manual and automated data selection method

View Table | View all tables in this article

The total hemoglobin concentrations for all 20 patients are calculated using both manual and automated data selection methods. Table 3 shows the comparison of the mean and the standard deviation of the maximum total hemoglobin concentration reconstructed by both methods as well as the ratio of the average HbT concentration between malignant and benign cases. Figure (6) shows the box plots of the reconstructed HbT concentrations of malignant and benign cases. The results indicate a slight increase of the malignant to benign ratio of HbT from 2.01 to 2.55 by using the automated data selection method.

Table 3. Comparison of maximum reconstructed absorption coefficients of malignant and benign groups using manual and automated data selection method

View Table | View all tables in this article

 

Fig. 6 Comparison of maximum total hemoglobin concentration reconstructed with both manual and automated data selection methods for 10 malignant and 10 benign cases. Vertical axis is the total hemoglobin concentration in μmol/Liter.

Download Full Size | PDF

5. Discussion and conclusion

In this study an automated pre-processing method for imaging-guided DOT is introduced. This method combines multiple data sets obtained from a normal breast which is used as a reference and implements outlier and saturation data removal, and linear fitting to form a robust compound set of reference measurements. Because the reference site does not have lesions inside, the location of the probe relative to the imaged area is not as critical as lesion site. Therefore, multiple reference data sets obtained from slightly different positions can be combined and used as the inputs of the automated compound reference selection method. On the other hand, in the affected breast, due to the presence of the lesion, the described outlier removal and compound data selection method is not appropriate. Instead, a perturbation filtering method, based on simulations, is introduced to eliminate the outliers and inaccurate data due to measurement errors.

The clinical results using the automated compound reference selection method were compared with the manually processed results. An experienced user performed the manual outlier removal and data selection after the imaging session and the best manually selected data sets have been used for the image reconstruction to compare with the automated method. Although the results calculated by our experienced user shows high significance for malignant and benign classification, the proposed automated compound outlier removal and data selection method outperforms the experienced user without the need for user interactions. On average, an expert needs about 15 to 20 minutes to perform the outlier removal and data selection manually, while the automated method takes about 20 seconds using a typical dual core computer of 2.3 to 2.8 GHz speed. T-test is performed on the total hemoglobin results obtained from each method and the test shows p-values of less than 0.001 in both cases which indicates strong statistical significance between malignant and benign groups. Although both methods provide strong statistical significance between malignant and benign groups, but the ratio for malignant to benign has been improved with the automated method. In particular, the ratio of malignant to benign lesions has been increased from 2.01 using manual data selection to 2.55 using automated data selection method with about 27% improvement.

Although our automated compound data selection method is introduced in the US-guided DOT approach, it can be readily applied to other imaging guided DOT approaches. In general, diffuse optical tomography relies on reference measurements, either reference data from a normal area of the same breast with lesions or a contralateral normal breast. The reference selected from a normal area of the same patient is considered the best healthy control because each patient has different background tissue optical properties depending on age, tissue composition, and menopausal status as reported by several research groups in the past. Some groups used average background tissue by excluding the tumor area in mammography guided or MRI guided diffuse optical tomography [12]. Some groups used a localized area of healthy tissue in the lesion breast [14] and some including our group used the mirror region in the contralateral breast [8,10,15]. Depending on the heterogeneity of the healthy tissue, the choice of reference area and the data quality of the area impacts the level of tissue contrast. The proposed automated robust reference selection method removes data outliers, rejects saturated and noise data, and performs a linear fitting if a linear model is used such as approaches using hand-held probes where semi-infinite boundary condition is applicable. Thus, the proposed automated reference selection method is applicable to these approaches to improve the reference data quality and reduce data processing time. This is an important step toward clinical translation and commercialization of diffuse optical tomography techniques.

In conclusion, an automated compound data selection method is introduced. This method forms a compound set of measurements from multiple data sets and outperformed the manual method in two different aspects. First, this method is able to eliminate the effect of inaccurate and noisy measurements by combining multiple data sets and improves the differentiation between malignant and benign cases. The results show improved differentiation between malignant and benign lesions compared to the manual data selection. Second, this method makes outlier removal, data selection, and processing automated by removing the need for experienced users to perform the manual data selection and outlier removal. Utilizing the automated and compound data selection method, the data processing speed is significantly improved and the user dependency of data processing is reduced. Future work includes automatic measurements of lesion boundaries in US images and input the measurements into image reconstruction.