1.

Introduction

Ovarian cancer has the highest mortality of all gynecologic cancers because it is predominantly diagnosed in the late stage due to the nonspecific early stage symptoms and lack of efficacious screening techniques.^1,2 Basic examinations including pelvic examination, blood test for CA-125, and transvaginal ultrasound have a low positive predictive value.^3–8 X-ray computed tomography (CT) is also poor in detecting small metastases of less than 2 cm in diameter.⁶ Furthermore, magnetic resonance imaging (MRI) is expensive to use and is only reserved as a secondary imaging method when transvaginal ultrasound is not deterministic, though it is more sensitive than x-ray CT.⁹ Positron emission tomography (PET) is useful for detection of recurrent or residual ovarian cancer and for monitoring response to therapy, but it may yield false-negative results in patients with small, necrotic, mucinous, cystic, or low-grade tumors.¹⁰

Based on ultrasound and CA-125 levels, the risk of malignancy index (RMI) was developed to estimate the risk of ovarian cancer. An RMI2 (RMI of method 2) has been estimated to have a sensitivity of 74% to 80%, a specificity of 89% to 92%, and a positive predictive value of approximately 80% for ovarian cancer detection.^11,12 Unfortunately, the subjective scoring with different physicians may result in different diagnostic results for the same patient. Other complex models, such as artificial neural network models and relevance vector machines, have been used to improve ovarian cancer diagnosis, with the relevance vector machine showing the best performance.^13,14

Recently, photoacoustic tomography (PAT) has been widely investigated as a noninvasive tool for cancer detection and diagnosis.^15–18 PAT utilizes high optical contrast and high resolution ultrasound to image tumor vasculature distribution. In PAT, biological tissue is illuminated by short laser pulses with wavelengths in visible red or near-infrared range, to achieve better penetration depth. Laser energy absorbed by tissue produces transient thermoelastic expansion and generates wideband acoustic waves, which can be detected by ultrasound transducer arrays for reconstruction of the optical absorption distribution of the tissue. Optical absorption distributions at multiple wavelengths can be used to map tumor hemoglobin distribution and oxygen saturation which are directly related to tumor proliferation, growth, and metastasis.¹⁹

Our research group has developed several co-registered ultrasound and photoacoustic imaging systems for small animal imaging and ovarian cancer detection and characterization; the ultrasound provides anatomical information, while the PAT provides the corresponding light absorption distribution within the anatomical structure.^18,20–23 The study reported in Ref. 18 evaluated 33 ex vivo ovaries and showed that malignant ovaries on average exhibited a much higher total absorption than normal ovaries. The quantitative parameter used to evaluate the total absorption was the measured “average maximum radio-frequency” (AMRF) PAT signal for each ovary, with the maximum taken from the RF signal at each ultrasound array element and the average taken across all the array elements. This single parameter method achieved a sensitivity and specificity of 83% between malignant and normal ovaries ($p=0.0237$) in the postmenopausal group. The study reported in Ref. 23 focused on the spatial frequency and statistical properties of the PAT images. The support vector machine (SVM) classifier in this study achieved 100% sensitivity and specificity for the training group of 33 ovaries, and 76.92% sensitivity and 95.12% specificity for the testing group of 37 ovaries. These results suggest that photoacoustic imaging is a promising modality for improving ultrasound diagnosis of ovarian cancer. However, the above two methods have their own drawbacks. For the single parameter testing used in Ref. 18, the measurement of the PAT signal can be affected by imaging system parameters and signal-to-noise ratio (SNR). When the PAT signal from an ovary has a low SNR, these RFs data-based method may not be sensitive enough.

On the other hand, the method developed in Ref. 23 focused on PAT images and their spatial patterns. PAT images represent the amplitude of light absorption but cannot provide details such as diagnostic information offered by the spectra of PAT beamformed data or beams, which is the delay-and-sum of all channel RF signals. Studies from other groups have revealed that the spectrum parameters from a PAT channel signal^24–27 or beam,^28,29 such as the spectral slope and 0-MHz intercept, are affected by the photoacoustic target dimension and chromospheres’ concentration;^24–29 these parameters could provide valuable information for cancer characterization and classification. Furthermore, in the formation of PAT images, the spatial patterns of PAT images are affected by the manually set dynamic range and the mapping method between beam envelope and pixel values. Improper setting of the dynamic range or nonlinear mapping of beam envelope and pixel values may lower the sensitivities and specificities for image-based classifiers.²³ In this paper, we extract new features from photoacoustic beams and their envelopes to explore the spectral information and spatial homogeneity of light absorption in ovarian tissues and to combine them with statistical and spatial features obtained from PAT images to improve the performance of classifiers for ovarian cancer diagnosis.

2.

Methodology

The details of the experimental setup and imaging system can be found in Ref. 30. Ovaries were freshly excised and obtained within 90 min after the oophorectomy at the University of Connecticut Health Center. All patients signed the informed consent.

The typical size of human ovaries is in the range of 2 to 7 cm with approximately elliptical shapes. In ex vivo studies, each ovary was placed inside intralipid solution and supported by a holder made of thin optical fibers. The distance between the sample and the imaging probe surface was adjusted to approximately 1 cm, which is the typical vaginal muscle wall thickness. The intralipid had a reduced scattering coefficient of $4{\mathrm{cm}}^{-1}$ and absorption coefficient of $0.02{\mathrm{cm}}^{-1}$ to emulate the optical properties of tissue. Based on the measurements of pig vaginal muscles, the reduced scattering coefficient is in the range of 3 to $5{\mathrm{cm}}^{-1}$ and the absorption coefficient is in the range of 0.02 to $0.04{\mathrm{cm}}^{-1}$. The details of the co-registered ultrasound and PAT probe used for the reported experiments were given in Ref. 30. Briefly, the probe consists of a transvaginal ultrasound probe (6-MHz central frequency and 80% bandwidth, Gore and Associates, Inc., Newark, Delaware) and 36 fibers of $200\text{-}\mu \text{m}$ core diameter mounted around the ultrasound transducer to uniformly deliver the light. Both ultrasound probe and the fibers were shielded inside a protective sheath internally lined with a highly reflecting aluminum sheet, with which uniform fluence underneath the probe surface can be achieved. The laser energy density on the tissue surface was $\sim 17\mathrm{mJ}/{\mathrm{cm}}^2$ at 750-nm wavelength, which is below the American National Standards Institute safety limit of $24\mathrm{mJ}/{\mathrm{cm}}^2$ at this wavelength.

Signals from 128 channels of the transvaginal probe were used to form the co-registered ultrasound and PAT images. Both the pulse-echo (PE) ultrasound image and PAT image were formed by using conventional delay-and-sum beamforming techniques.^31,32 In PE imaging, one angular beam is formed by applying a delay-and-sum algorithm on the received RF signals from all channels or a certain number of channels (65 in our study) centered in the direction of the transmission. A total of 93 evenly distributed PE beams over 120 degrees are used to form PE images. In PAT, a sector scan is also used to co-register with the corresponding PE image. Each angular beam is formed by applying a delay-and-sum algorithm on the received RF photoacoustic signals synchronized with each laser pulse. A total of 151 PAT beams were formed and evenly distributed across the 120-deg image plane. The PE and PAT beam envelopes were subsequently obtained by base-band demodulation and amplitude detection. Spectral information was lost in the procedure of demodulation and amplitude detection of beam envelope. However, it was preserved in the PE and PAT beams. In the preprocessing stage, the mapping between the PAT beam envelopes and PAT image pixel values was set to be linear instead of the conventional logarithmic mapping in ultrasound imaging. The dynamic range of the image was set by mapping the system noise threshold of the PAT beam envelope to a pixel value of 0 and mapping the maximum of the envelope to a pixel value of 255. The system noise threshold for each image frame was estimated from the averaged envelope peaks of nonregion-of-interests (nonROIs) identified from the co-registered PE image. NonROIs are regions away from ovarian samples in the imaging plane where the light absorption was negligible. Both RF channel data from the transducer array and PAT beams contained spectral information. However, the channel data were the integration of signals from the arc with the center at the channel element and its SNR was low in general. In beamforming, the coherent RF signals from the same absorber detected at all channels were appropriately delayed and summed and the SNR was improved. In the reported study, we have used the averaged spectral parameter (see Sec. 2.2.1 for details) of PAT beams in the suspicious area to extract features for classification.

2.1.

Center Estimation for Suspicious Area

The ultrasound imaging field of view covers 120 deg in the lateral dimension and 7.7 cm in depth, while reconstructed PAT images are mainly confined to several centimeters’ area underneath the probe. Therefore, limiting data analysis within the suspicious area will greatly reduce the total computation time. To crop the suspicious areas, the centers of these areas should be estimated. A simple method introduced in Ref. 23 was applied as follows: the radon transforms of a PAT image along 0 ($x$ axis) and 90 deg ($y$ axis) can be simply calculated from the summation of the PAT pixel values along the column and row directions. After normalization to their respective peaks, the radon transforms along 0 and 90 deg were used to estimate the center of the suspicious area using the normalized Gaussian model in the following equation:

Eq. (1)

$$f(x)=\exp [-\frac{{(x-\mu )}^2}{2{\sigma }^2}].$$

In Eq. 1, $\mu$ is the index of the estimated center and $\sigma$ represents the spatial expansion of vasculature along the projection direction. A 3.0 cm by 1.5 cm window (large enough in most cases) with the estimated center was then applied to the PAT image to crop the suspicious area as shown in Fig. 1. Note that the normalized Gaussian model is used to estimate the center of the suspicious area from 0 and 90 deg or $x$ and $y$ projections.

Fig. 1

(a) A representative co-registered ultrasound/PAT image of malignant ovary and its radon transform along $x$ axis and $y$ axis; (b) a representative co-registered ultrasound/PAT image of normal ovary and its radon transform along $x$ axis and $y$ axis. The red rectangles are the cropped suspicious areas.

2.2.

Feature Extraction

Seventeen features were selected from three categories: PAT beams, PAT beam envelopes, and PAT images. Spectral parameters from the PAT beams revealed the distribution of the frequency components within the frequency range of the transducer array which could be related to vasculature dimension and concentration.^24–29 Features from PAT beam envelopes, such as the total summation of envelopes (PAT summation)³³ and the peak envelope, described the light absorption distribution within the suspicious area. Statistical analysis was performed on the PAT images to obtain their mean values and variances. By separately applying several spatial filters constructed from the joint frequency spectrum of the malignant or normal cases, the common features which were not visibly distinguishable but were embedded in the spatial spectrum of the PAT image could be extracted.²³

2.2.1.

Spectrum parameters of photoacoustic tomography beams

Power spectra were calculated using fast Fourier transform (FFT) on the PAT beams inside the cropped window. The cropped rectangular window defined the angular sector and the depth range of PAT beams inside the window, as shown in Fig. 2. To ensure a high quality spectral estimation, any beams with SNR below 15 dB were excluded from spectral information extraction. The spectrum of PAT beam normalized to the spectrum of an approximate point-like target (a $100\text{-}\mu \mathrm{m}$ black thread orthogonal to PAT imaging plane) was used to estimate the spectral slope and 0-MHz intercept (Fig. 3).²⁹ A linear fitting was performed on the normalized spectrum in logarithmic scale within the 80% bandwidth of transducer’s 6-MHz central frequency. The spectrum of an approximate point-like target characterizes the frequency response of the transducer and the electrical receiving system.²⁹ The PAT beam and its spectrum obtained from the approximate point target-like were shown in Fig. 4. After taking the average of the spectral slopes and 0-MHz intercept obtained from all high quality beams, the mean spectral slope and 0-MHz intercept were chosen as spectral features.

Fig. 2

Illustration of photoacoustic tomography (PAT) beams and the selection of beams for extracting spectral information of malignant [Fig. 1(a)] and normal [Fig. 1(b)] ovarian samples. PAT beams fan out from the transducer surface (arc LM) and form an angular scan range of 120 deg with the center at the origin O. The sector OLM is physically occupied by transvaginal probe and the PAT beam length is 7.7 cm. Enclosure DEGF defines the minimal annular segment covering the whole cropped window. PAT beam segment JK contains noise only and PAT beam segment HI covers the absorbing area. The waveform, spectra, and linear spectral fitting of PAT segment HI are shown in Fig. 3. The rectangles confined by two horizontal dashed lines in (a) and (b) are shown in Fig. 1.

Fig. 3

Waveform of representative PAT beams within cropped windows, their spectra and linear fittings. Subfigures (a) and (c) are the waveform, spectrum (blue), and linear fitting (red) of PAT beam HI of malignant ovarian tissue shown in Fig. 2(a). Subfigures (b) and (d) are the waveform, spectrum (blue), and linear fitting (red) of PAT beam HI of normal ovarian tissue shown in Fig. 2(b). The malignant sample has a smaller absolute value of fitted slope than that of the normal sample, indicating more high frequency components in its beams.

Fig. 4

Waveform of PAT beam of a point-like target (a) and its spectrum (b).

2.2.2.

Features of photoacoustic tomography beam envelope

PAT beam envelopes contain the amplitude information of PAT beams. The system noise threshold was estimated first from the average peaks of the PAT beam envelopes in the nonROIs, where almost no photoacoustic signals were generated. The system noise threshold was used as a cutoff value for PAT image formation. The count of image pixels with light absorption higher than the system noise threshold is defined as the PAT area and the summation of the envelopes of these pixels is defined as PAT summation.³³ PAT summation describes the total light absorption within the suspicious area. Peak envelope is the maximum of PAT beam envelope within the suspicious area, representing the maximum light absorption. It is a replacement for the single parameter AMRF used in Ref. 18. Compared with AMRF, peak envelope is less affected by noise since the delay-and-sum beamforming can effectively increase the SNR of the PAT beams. We characterized the homogeneity of the PAT beam envelopes by the mean correlation coefficients of adjacent envelopes within the cropped windows since the light fluences along adjacent beam directions were similar. The histogram of PAT beam envelope was used to carry out an estimation of probability density function of Rayleigh distribution.³⁴ The scale parameter of the Rayleigh distribution was also chosen as one of the features.

2.2.3.

Features of photoacoustic tomography images

The features from the PAT images were mainly adopted from Ref. 23. The mean and variance of the PAT images describes the average light absorption intensity and the spatial variation of the absorption. From our experimental observations, the PAT images of malignant cases usually showed clustered distributions due to the abundant and localized smaller vessels, whereas the distribution was more diffused, scattered, and spatially spread out in normal cases.²³ The center estimation method outlined in Sec. 2.1 was extended to quantify this observation as follows: the normalized radon transforms from 0 deg to 90 deg with an interval of 1 deg were calculated and fed to the Gaussian model after averaging. The parameter sigma ($\sigma$) estimated from the mean radon transform was used as a feature to describe the overall spatial spread of photoacoustic intensity (Fig. 5). A high value of Gaussian fitting error meant a high irregularity of PAT spatial distribution.²³ This error was defined as the norm of the difference between the mean radon transform and the fitting value of the Gaussian model.

Fig. 5

(a) Mean radon transform and its Gaussian fitting of the malignant ovary shown in Fig. 1(a); (b) mean radon transform and its Gaussian fitting of the normal ovary shown in Fig. 1(b); malignant ovaries often have a clustered absorption distribution while normal ovaries are more diffused. More than one spike may present in the mean radon transform of malignant PAT images, resulting larger fitting errors than those of normal ovaries.

Two sets of optimum spatial filters were constructed from the existing malignant²³ and normal ovarian tissue PAT images, which served as image templates for malignant and normal ovaries, respectively (see Ref. 23 for more details). When the order parameter $k$ of the spatial filter was set to 1, the filter was linear. Additional nonlinear versions (cubic root and binary) of this filter were achieved by setting $k$ to $1/3$ and 0.01, respectively, to make them more tolerant to distortions and obtain better recognition SNR.³⁵ The peak outputs of the two sets of spatial filters were chosen as features for two classifiers. High peak outputs of spatial filters constructed from malignant PAT image templates and low outputs of spatial filters constructed from normal PAT image templates were expected to be an indication of cancer. On the other hand, high peak outputs of spatial filters constructed from normal PAT image templates and low outputs of spatial filters constructed from malignant PAT image templates were more likely to be an indication of normal cases.

3.

Results and Discussions

Table 1 lists 17 extracted features and indices of the 13 features after feature correlation analysis and ranking. Among the 17 features, features from the first two categories and the maximum output of the spatial filters constructed from normal ovaries are new features in addition to those reported in Refs. 18 and 23. The new features provide additional information from the spectrum of PAT beams and the quantitative description of light absorption in suspicious areas. Figure 6 shows the mean MCE of the two classifiers versus the number of features after ranking. The mean MCE decreases quickly with the feature number until 7. After 7, it slightly increases for both logistic and SVM classifiers. The logistic and SVM classifiers achieve minimum MCEs of $7.06\pm 3.29\%$ and $6.05\pm 3.67\%$, respectively, when the first seven features are used. Figure 7 compares the sensitivity, specificity, PPV, and NPV of the two classifiers. Though not every one of these four parameters obtains its maximum value when the first seven ranked features are used, the improvement of these parameters beyond the use of the first seven features is relatively small. Figure 8 shows the mean values, standard deviations, and $p$-values of $t$-tests of the first seven ranked features for all malignant and normal samples. They are ranked in the following order: the peak outputs of spatial filters (with cancerous template, $k=1$), PAT envelope homogeneity, peak outputs of spatial filters (with normal template, $k=1/3$), mean spectral slope of PAT beams, PAT area, Gaussian fitting error of mean radon transform, and maximum PAT envelope. The $p$-values for the first seven features are all smaller than 0.005, which is statistically significant. The $p$-values for the eighth and ninth features are around 0.01. They are statistically significant but do not improve the overall performance of these two classifiers. The $p$-values for the other four features are larger than 0.14 and the inclusion of these features increases the MCE of the two classifiers. Smaller absolute values of spectral slopes and lower values of the PAT envelope homogeneity for malignant ovaries in Fig. 8 demonstrate that the malignant cases have relatively more high frequency components in the PAT beams and abrupt changes in the PAT beam envelopes because of rich microvasculatures in the tumor region. Malignant ovaries show higher values in the maximum output when the corresponding cropped PAT images were input into the spatial filter constructed from the malignant template, and lower values when fed into the spatial filter constructed from the normal template. At the same time, malignant ovaries have a larger PAT area, higher peak envelope, and higher (irregularity) fitting error of mean radon transform of PAT images. The following performance evaluation of classifiers is based on these seven features.

Table 1

Seventeen extracted photoacoustic tomography (PAT) features and indices of the 13 features after feature correlation analysis and ranking.

Fig. 6

Mean misclassification error (MCE) of both logistic and SVM classifiers for the testing group. Both logistic and SVM classifiers obtain minimum mean misclassification error (7.06% and 6.05%, respectively) when the first seven ranked features are used.

Fig. 7

Sensitivity, specificity, positive predictive value, and negative predictive value with respect to different numbers of ranked features.

Fig. 8

Means, standard deviations, and $p$-values of the seven selected features for all 124 malignant samples and 45 normal samples.

We also obtained an average feature value of all samples from each patient and calculated the average feature values of four patients with malignant ovaries and six patients with normal ovaries. The results of these two groups are shown in Fig. 9 and the average features between these two groups are statistically significant.

Fig. 9

Means, standard deviations, and $p$-values of the average values of seven selected features for four patients of malignant ovaries and six patients of normal ovaries.

Table 2 lists the sensitivity, specificity, PPV, and NPV for the training and testing data. For the training group, both classifiers have achieved high accuracy. For the testing group, the sensitivity and NPV of the SVM classifier are similar to those of the logistic classifier while the specificity and PPV of the SVM classifier are 7% to 11% better.

Table 2

Performance of two classifiers on training and testing data.

Figure 10 shows the ROC curve of logistic classifier for one set of testing data with an AUC of 0.98. The average AUC for the total 50 sets of testing data is $0.92\pm 0.05$.

Fig. 10

Receiver operating characteristic (ROC) curve for one set of testing data of logistic classifier with an area under ROC of 0.98; the mean area under ROC curve for the 50 sets of testing data is $0.92\pm 0.05$.

Selection of the optimal set of features is critically important for achieving the best performance for classification of malignant versus normal/benign ovarian tissue. Though all selected seven features are of low $p$-values, the classification result from single feature is not as good as that of the combined seven optimal features. Figure 11(a) shows the MCEs of the logistic and SVM classifiers for each individual feature from the 7 optimal features set, 6 features from PAT images, all 13 features after ranking, and the optimal 7 features. The inclusion of features from PAT beams and PAT beam envelopes results in lower MCE when all 13 features are used; however, the MCE achieved the lowest value when the seven optimal features are used for classification. For 50 resampled datasets, the output from the logistic function can be used to estimate the AUCs and their 95% confidence intervals with 1000 stratified bootstrap replicates.⁴³ The averaged AUC and averaged 95% confidence intervals for different feature combinations are shown in Fig. 11(b). The optimal seven feature set generates the highest average AUC (0.92) and narrowest average 95% confidence interval (0.86, 0.97), while feature # 3 gives the lowest average AUC (0.58) and widest average 95% confidence interval (0.42, 0.73).

Fig. 11

(a) The MCE of logistic and SVM with RBF kernel for each individual feature, six PAT image features in Table 1 and all 13 features after ranking and the selected optimal seven features. The optimal seven features result in the lowest MCE comparing with other feature combinations for both classifiers. (b) Average areas under ROC curve (AUC) for logistic classifier and their average 95% confidence interval. The optimal seven features yield the highest average AUC (0.92) and the narrowest average 95% confidence interval (0.86, 0.97).

4.

Summary

We extracted a total of 17 features from 169 datasets of 19 ex vivo ovaries. In addition to features reported in Ref. 23, new features from PAT beams and envelopes were extracted to obtain spectral information and light absorption distribution in the ovarian tissue. New spatial filters were constructed from normal ovary PAT images to characterize normal ovaries. The performances of the logistic and SVM classifiers were evaluated using misclassification error criterion. Both classifiers achieved the lowest MCE when seven optimal features were used. With these seven optimal features, the logistic and SVM classifiers obtained sensitivities of $96.39\pm 3.35\%$ and $97.82\pm 2.26\%$, and specificities of $98.92\pm 1.39\%$ and 100%, respectively, for the training group. For the testing group, the logistic and SVM classifiers achieved sensitivities of $92.71\pm 3.55\%$ and $92.64\pm 3.27\%$, and specificities of $87.52\pm 8.78\%$ and $98.49\pm 2.05\%$, respectively. SVM has a better performance than the logistic classifier because its sensitivity and PPV are higher than that of the logistic classifier while the specificity and NPV of both classifiers are similar.

Among the four patients with malignant ovaries, one patient had an early stage cancer on both ovaries. Three features in the selected seven optimal features for this case were different from those of other cancer cases. These features are PAT envelope homogeneity, PAT spectral slope, and maximum PAT envelope as summarized in the Table 3. As shown in Figs. 8 and 9, malignant ovaries have a lower PAT envelope homogeneity, a smaller absolute value of PAT spectral slope, and a higher maximum PAT envelop than those of normal ovaries. The early stage cancer case has shown 31.31% lower PAT envelope homogeneity and 33.45% lower absolute spectral slope than those of other cancer cases. Thus, these two features for the early stage cancer case show an even larger difference that those of normal cases. However, the maximum PAT envelope of this early stage cancer case is 39% lower than those of other three cancer cases. This example further highlights the importance of using a set of optimal features rather than any individual feature to characterize the malignant versus normal/benign ovarian tissue.

Table 3

Comparison of features between the early stage cancer and other cancer cases.

Ultrasound features have not been included in this study and will be explored in our future studies. We believe this will further improve the performance of our classifiers. Real-time feature extraction and implementation of the SVM-based classifier can also be achieved by using well-established SVM libraries.³⁷ With the help of our real-time co-registered ultrasound and photoacoustic imaging system, the real-time classifier has a great potential for in vivo detection and characterization of ovarian cancer. Recently, we have received an institutional review board approval for a human subject protocol and the results of the in vivo patient study will be reported in the near future.