Reducing image artifact in diffuse optical tomography by iterative perturbation correction based on multiwavelength measurements.

Ultrasound (US) guided diffuse optical tomography has demonstrated great potential for breast cancer diagnosis, treatment monitoring, and chemotherapy response prediction. Optical measurements of four different wavelengths are used to reconstruct unknown optical absorption maps, which are then used to calculate the hemoglobin concentration distribution of the US visible lesion. Reconstructed absorption maps are prone to image artifacts from outliers in measurement data from tissue heterogeneity, bad coupling between tissue and light guides, and motion by patient or operator. We propose an automated iterative perturbation correction algorithm to reduce image artifacts based on the structural similarity index (SSIM) of absorption maps of four optical wavelengths. The initial image is estimated from the truncated pseudoinverse solution. The SSIM was calculated for each wavelength to assess its similarity with other wavelengths. An absorption map is repeatedly reconstructed and projected back into measurement space to quantify projection error. Outlier measurements with highest projection errors are iteratively removed until all wavelength images are structurally similar with SSIM values greater than a threshold. Clinical data demonstrate statistically significant improvement in image artifact reduction.

Introduction

Breast cancer is the second leading cause of death for women in the United States. Diffuse optical tomography (DOT) and diffuse optical spectroscopy are functional imaging modalities being explored for breast cancer diagnosis and treatment prediction. DOT provides quantitative estimations of oxygenated and deoxygenated hemoglobin concentrations calculated from reconstructed optical absorption maps of breast lesions. Hemoglobin concentrations are directly related to tumor angiogenesis, a hallmark of cancer. Since diffused light suffers from poor lesion localization and inaccurate target quantification due to intense scattering, image reconstruction is often guided by another high-resolution imaging modality, such as ultrasound (US), magnetic resonance imaging,, or x-ray, which provide structural information about the tumor. While DOT guided by other modalities can significantly reduce the number of voxels with unknown optical properties for more accurate reconstruction, image artifacts caused by outlier measurements due to wavelength-dependent tissue heterogeneity, bad coupling between tissue and sources and detectors, and patient motion or operator hand motion can distort the reconstructed images. Several experimental and modeling approaches have been developed for system calibration of optical source strengths, detection channel gains and phase shifts, source and detector (optodes) position errors, and coupling errors between skin and optodes., For correcting motion artifacts, different algorithms have been proposed, including cubic spline interpolation, adaptive Kalman filtering, and wavelet-based motion artifact removal. To improve the imaging quality, projection error-based adaptive regularization techniques have been employed, which outperform standard Tikhonov regularization. However, these approaches do not compensate for wavelength-dependent problems. For example, the 740-nm wavelength is prone to measurement errors from tissue heterogeneity caused by dark skin and skin pigment. As another example, the 830-nm wavelength has lower signal to noise ratio at longer source and detector distances due to the reduced sensitivity of photomultiplier tube (PMT) detector beyond 800 nm.

In our US-guided DOT approach to assessing breast cancer, perturbation, which is the normalized difference between the lesion breast and the contralateral normal breast (reference) measurements, is used for mapping lesion absorption at each wavelength. The total hemoglobin map is computed from the absorption maps of four optical wavelengths. The tissue heterogeneity of the reference measurements contributes to outliers in the perturbation measurements. In a recent investigation by our group, Vavadi et al. introduced a statistical method based on the semi-infinite tissue model to automatically remove outliers from contralateral normal breast measurements. However, this method cannot be used for perturbation measurements because lesion measurements 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. Recently, Althobaiti et al. introduced an approach for data filtering based on multiple wavelength measurements collected at the 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. However, this approach requires that two to three wavelength perturbation datasets must be correlated, and then the rest of the wavelength-dependent distortion can be compensated for.

In this paper, we propose an iterative perturbation correction algorithm by using structural similarity index (SSIM) as an image quality assessment criterion. The initial estimate of the absorption map is obtained from the truncated pseudoinverse solution. In subsequent iterations, the average SSIM for each wavelength and errors between the measurement and the projected data are computed and outliers are removed from the measurements based on the errors. This procedure is iterated until the SSIM reaches a preset threshold. We demonstrate the effectiveness of this approach in phantoms and clinical data.

Methods

US-guided DOT System, Data Acquisition, and Calibration

US-guided DOT system is a frequency domain imaging system consisting of a hand-held DOT imaging probe with an US transducer located in the middle. Four laser diodes with wavelengths 740, 780, 808, and 830 nm are sequentially switched by a and a optical switch to deliver light modulated at 140 MHz to each of the nine source positions on the probe. Fourteen parallel PMTs detect reflected light via light guides from the tissue. A custom A/D board samples detected signals from all channels and stores data in a PC. Using data collected from a homogenous intralipid solution, the system is calibrated to compensate for system gains and phase delays. Multiple datasets acquired from the contralateral normal breast are used to compute a robust reference. The selected reference is considered as a homogeneous reference. The fitted optical absorption () and reduced scattering () are used to calculate a weight matrix, , for image reconstruction using a dual mesh scheme.

Data Preprocessing

DOT data acquisition is performed on both a lesion breast and a contralateral normal breast, referred to as the reference breast. Amplitude and phase measurements are extracted from the detected radio frequency signal using the Hilbert transform. For the ’th source–detector pair, reference measurement is given as , and the lesion measurement is , where , and is the total number of source–detector pairs or the total number of measurements. Perturbation, , is defined as the normalized difference between the reference and target measurements:

The first term is the real part of the perturbation, and the second term is the imaginary part of the perturbation. A typical 2-D representation of perturbation data for a phantom is shown in Fig. 1, with real perturbation on the axis and imaginary perturbation along the axis. The unit circle represents the expected boundary that perturbation data should lie within. From simulations of different target contrasts and locations in depths, it was shown that maximum phase difference for any source–detector pair should not exceed 90 deg, even in extreme cases. Since for for all , real perturbation should be greater than . For a high-contrast target, , so perturbation is more skewed toward the negative real axis. For a low-contrast target, perturbation may be small and clustered around the origin, or distributed more toward the positive real axis, or distributed evenly across both positive and negative sides around the origin. Note that the imaginary part of the perturbation can be positive or negative because for .

Fig. 1

Phantom perturbation data. (a) Data measured from a high-contrast phantom target imbedded in intralipid solution. (b) Data measured from a low-contrast phantom target imbedded in intralipid solution.

Figure 2 shows one set of perturbation data of a highly absorbing malignant breast lesion and a low absorbing benign breast lesion. As evident from the figure, clinical data is more scattered because of tissue heterogeneity, bad coupling of the tissue and source and detector fibers, and patient movement or operator hand motions.

Fig. 2

Clinical perturbation data. (a) A malignant breast lesion and (b) a benign breast lesion.

Multiple perturbation datasets are compiled together, and a multivariate Gaussian is fitted, and data points are removed if there are any with a Mahalanobis distance greater than the threshold computed from the inverse chi-square distribution with a cumulative probability of 99%. Based on chest wall matching of the reference and lesion breast, a single measurement dataset is selected from multiple measurements. The structural similarity-based perturbation correction algorithm depicted in Fig. 3 is applied to perturbation data to obtain corrected perturbation and artifact-free images.

Fig. 3

Data preprocessing and iterative perturbation correction algorithms.

DOT Image Reconstruction

The DOT inverse problem is typically linearized by Born approximation. By digitizing the imaging space into voxels, the resulting integral equations are formulated as follows: where is the measured scattered photon density wave, is the number of measurements, and denotes the unknown differences in the absorption coefficients of each voxel. The weight matrix, , describes the distribution of the diffused wave in the homogenous medium and characterizes the measurement’s sensitivity to the absorption and scattering changes. In the end, the reconstruction problem can be formulated as a regularized optimization problem:

A two-step imaging reconstruction was proposed earlier to solve this regularized optimization problem after obtaining an initial image estimate, , by truncated pseudoinverse in the first step and refining the solution using the regularized conjugate gradient (CG) algorithm. The implementation of CG is adapted from Ref. 26.

Image Quality Assessment

In DOT reconstruction, quantitative assessment of imaging quality poses a challenge. Previously, image distortion and inconsistent images for different wavelengths were visually inspected, and perturbation was manually corrected by an experienced operator. Such manual data processing is operator dependent and time consuming. In this manuscript, we propose to use the SSIM to quantitatively evaluate imaging quality by taking all four wavelength images into account. The SSIM measure is a function of the image’s luminance, contrast, and structure., The SSIM between two images and is defined as in Ref. 27: where , , and are the luminance, contrast, and structure similarity, respectively, and , , are three parameters used to adjust relative importance of the three components of the similarity measure. The luminance, contrast, and structure of an image are computed from mean, standard deviation, and normalized images as follows:

Here, , , , , and are the means of pixel values of image and image , the standard deviation of image and image , and the covariance of image and , respectively. , , and are constants.

For each wavelength, , the other three wavelength images are used as references to compute SSIMs for three image pairs. An average of the three SSIMs is the quantitative image quality index, , used to evaluate the reconstructed image quality of wavelength as given below:

Iterative Perturbation Correction

Iterative perturbation correction is performed based on for each wavelength. The wavelength with the minimum is corrected first. The initial estimate is from the truncated pseudoinverse. If is lower than a preset threshold (0.9), perturbation from wavelength is corrected based on the original perturbation and projected perturbation. The reconstructed image, , for is projected into measurement space by multiplying the weight matrix, , to obtain projected data:

Based on the Euclidean distance of original perturbation data, , and projected data, , projection error, , is calculated as follows:

The data point with maximum projection error is removed from . Modified perturbation is again used to reconstruct the absorption map for wavelength using regularized CG. is recomputed and compared with the threshold. This process is repeated until the lowest is greater or equal to the threshold. This iterative correction procedure is performed for each wavelength until the values for all four wavelengths are above the threshold.

Note that there is a wavelength-dependent variation in absorption values of different wavelengths for different oxygen conditions. However, the same lesion should have similar absorption distributions over the narrow wavelength window of 730–830 nm. Thus, we use SSIM to characterize the similarity between wavelengths. The algorithm does not enforce a perfect or 100% image similarity but somewhat similar not totally different as determined by this preset threshold.

Results

Phantom Experiments

In phantom experiments, intralipid solution was used to simulate a homogeneous background medium. The experiment was repeated for solid spherical balls with different contrasts and sizes simulating different types of tumors in intralipid solution in different depths. The average image similarity index was computed for all phantom absorption map images. The reconstructed absorption map for a ball of 2-cm diameter at 2-cm depth is shown in Fig. 4. The average structural similarity indices for four wavelengths —740, 780, 808, and 830 nm—are 0.98, 0.97,0.99, and 0.96, respectively.

Fig. 4

Reconstructed image similarity for phantom data: (a) US image and (b) reconstructed absorption maps (two layers at and ) for all four wavelengths. Each 2-D layer is . Average SSIMs are 0.98, 0.97, 0.99, and 0.96 for 740,780, 808, and 830 nm, respectively.

Pairwise SSIMs () are presented in Table 1. Large similarity indices indicate strong structural similarity among different wavelengths, which is visually apparent in Fig. 4.

Table 1

SSIM (mean ± standard deviation) for phantom data.

	740 nm	780 nm	808 nm	830 nm
740 nm	—	0.976±0.004	0.988±0.003	0.942±0.015
780 nm	0.976±0.004	—	0.985±0.006	0.943±0.023
808 nm	0.988±0.003	0.985±0.006	—	0.947±0.020
830 nm	0.942±0.015	0.943±0.023	0.947±0.020	—
Average	0.954±0.019	0.969±0.008	0.957±0.018	0.974±0.009

Clinical Study

A clinical study was approved by the local Institutional Review Board and was compliant with the Health Insurance Portability and Accountability Act. Informed consent was given by each patient. Data used in this study have been deidentified. A total of 40 patients were studied including 13 malignant and 27 benign lesions, based on biopsy results. All patients were categorized into two categories: patients with image artifact present in one or more wavelength absorption maps (17 patients) and patients with no image artifact (23 patients). This categorization was based on a preset cutoff value of 0.9 for structural similarity among the four wavelengths reconstructed absorption map images. An example of benign fibroadenoma is shown in Fig. 5. Figure 5(a) shows the US image with the lesion marked by a white ellipse. Figure 5(b) shows reconstructed absorption maps for four wavelengths. Each wavelength absorption map has one 2-D layer at depth, . The mean SSIMs for the four wavelengths 740, 780, 808, and 830 nm are 0.87, 0.91, 0.87, and 0.82. The reconstructed maximum absorption coefficients are 0.2463, 0.2069, 0.1326, and , respectively. Image SSIM indicates that there is an image artifact at wavelength 830 nm, and visual inspection confirms this. Figure 5(c) shows reconstructed absorption maps after perturbation correction. The mean SSIMs for the four wavelengths change to 0.95, 0.97, 0.97, and 0.96 while maximum absorption coefficients change to 0.1582, 0.1470, 0.1345, and , respectively.

Fig. 5

Image artifact reduction for a benign case: (a) US image, 1 cm lesion depth, (b) absorption maps for original data before perturbation correction, and (c) absorption maps after perturbation correction.

Iterative changes in the perturbation and absorption map for this case at 830 nm are shown in Fig. 6. For iteration 0, we see the original dataset and the absorption map: the map is similar to that in Fig. 5(b). In successive iterations, we removed perturbation points denoted by red dots. We continue to remove perturbations until absorption map is structurally similar to the maps of other wavelengths.

Fig. 6

Iterative changes in absorption map and perturbation filtering for 830 nm for the benign case. Red dots denote removed data points.

An example of malignant breast cancer with mixed ductal and lobular features is shown in Fig. 7. Figure 7(a) shows an US image with a lesion marked by a white ellipse. Figure 7(b) shows reconstructed absorption maps for the four wavelengths. Each wavelength absorption map shows two layers at depths, and 2 cm. Mean image similarity indices for the four wavelengths 740, 780, 808, and 830 nm are 0.83, 0.82, 0.77, and 0.79, and reconstructed maximum absorption coefficients are 0.254, 0.237, 0.070, and , respectively. Image artifact is present in 808- and 830-nm absorption maps. Figure 7(c) shows reconstructed absorption maps after perturbation correction. Mean SSIMs for four wavelengths improve to 0.94, 0.94, 0.93, and 0.91 while reconstructed maximum absorption coefficients changed to 0.254, 0.237, 0.228, and , respectively.

Fig. 7

Image artifact reduction for a malignant case: (a) US image, 1.5- and 2-cm lesion depths, (b) absorption maps for original data before perturbation correction, and (c) absorption maps after perturbation correction.

Figure 8 shows iterative changes of perturbation and absorption maps for the malignant case at 808 nm. In iteration 0, we have original dataset and absorption map similar to that in Fig. 7(b). In successive iterations, we removed perturbation points denoted by red dots. We continue to remove perturbations until the absorption map is structurally similar to maps of other wavelengths.

Fig. 8

Iterative changes of absorption map and perturbation filtered at 808 nm for the malignant case. Red dots denote removed data points.

Perturbation correction statistically improves the SSIM among different wavelengths, as depicted in Fig. 9. A two-tailed paired -test was done for images with artifacts both before and after perturbation correction, and the SSIM is statistically higher after perturbation correction, with a -value less than 0.001. Student’s -test on images with no artifacts shows no significant change in terms of structural similarity (-value 0.52), which is expected.

Fig. 9

Comparison of SSIMs of reconstructed images before perturbation correction (blue box) and after perturbation correction (red box).

Discussion and Summary

In summary, an iterative perturbation correction algorithm based on image similarity is introduced and its performance in image artifact reduction is demonstrated using clinical data. This algorithm follows two simple assumptions. First, absorption map images for all four wavelengths are assumed to be structurally similar. Since we are imaging the same tissue region with closely spaced wavelengths, the image structure should be similar, even though the local absorption coefficients might differ due to wavelength-dependent absorption variations. Second, image artifacts in all four wavelengths are assumed to be dissimilar. Data acquisition is done sequentially, from one wavelength after another, in a few seconds. Motion or experimental errors can affect one or more wavelengths, but these effects are random and are unlikely to generate structurally similar artifacts at all four wavelengths. Additionally, certain tissue heterogeneity caused artifacts may not be present at all wavelengths, for example, 740 nm is very sensitive to dark skin pigment than other wavelengths. The total hemoglobin distribution, which is calculated by linear weighting of multiwavelength absorption maps based on extinction coefficients, is significantly improved due to artifacts reduction. However, the average maximum total hemoglobin levels which we have used to classify malignant versus benign lesions remain statistically the same as compared with no perturbation correction. This is because one or two absorption distributions are often significantly improved on artifacts, but the maximum total hemoglobin level may not change much.

Other approaches, such as the weighted least-square (WLS) approach, can compensate measurement differences between different wavelengths (or source/detector channels). Such a method depends on the accurate modeling of the system noise. However, the noise in the DOT measurements includes coherent and incoherent components as well as random motion produced noise. We have evaluated the WLS approach to compensate measurement differences between different wavelengths by downweighting the noisy data instead of completely removing them. However, we have found that WLS may not be suitable for the measurement noise we experienced in the patient data.

In the past, we have performed simultaneously the reconstruction of oxygenated hemoglobin, , and deoxygenated hemoglobin, Hb. However, because the is about twice as large as Hb, the Hb is over-reconstructed and is under-reconstructed due to the weight distribution in weight matrix for Hb and in simultaneous reconstruction. When each wavelength data is used for reconstructing each absorption map and then all absorption maps are used to compute and Hb, there is no cross-coupling of and Hb.

In our earlier studies, we have attempted to simultaneously reconstruct both absorption and scattering distributions of breast lesions. We had success in phantom data but these algorithms are not robust for patient data when they applied to a large patient database. In the simultaneous reconstruction, the distribution of lesion diffusion coefficient, , is reconstructed with the distributin of lesion absorption, , where . However, is one order of magnitude smaller than and cannot be reconstructed reliably for all patient data. Additionally, simultaneously reconstruction of both absorption and scattering distributions doubles the unknowns in the image reconstruction. Since is directly related to tumor angiogenesis and we have focused on this important parameter in our algorithm development

In summary, the proposed iterative artifact reduction algorithm significantly reduces the effect of wavelength-dependent measurement errors in DOT perturbation, which helps to achieve a more accurate reconstruction of the optical properties of breast lesions. This automated method also helps to minimize both the user interface and the time for data preprocessing. The average time for an experienced user to manually perform data preprocessing for one patient’s data is from 15 to 30 min. The automated method could reduce this time to less than a minute and facilitate the clinical translation of US-guided DOT technology. Although the method is demonstrated using US-guided DOT data, it is applicable to any DOT data preprocessing obtained with multiple wavelengths.