2.1.

Patients and Specimen Preparation

Breast tissue was studied from 51 random, nonconsecutive patients (patient age range, 39 to 80 years; mean age, 59 years) undergoing either wide local excision (WLE) or mastectomy at Addenbrooke’s Hospital in Cambridge and Guy’s Hospital in London.⁶ Approval for the study was granted by the respective Research Ethics Committees; signed informed consent, agreeing to research on tissue removed at the time of surgery, was obtained from all patients. All samples were fresh at the time of measurement.

Samples were first inked according to a standard pathology protocol so that margins could be identified during examination of the slides. They were then sliced to expose any palpable lesion; slices were typically $>1\mathrm{cm}$ thick. If there was no palpable tumor the specimen was sliced through suspected regions of tumor identified from the radiographs of the excised sample.

2.2.

Data Acquisition

Data from all breast specimens were collected using TPIscan (TeraView Ltd., Cambridge, UK). The system uses photoconductive methods to generate and detect terahertz pulses in reflection mode.³¹ A detailed description of the system is given by Wallace et al.³² The data were collected over an area that was typically 20 mm by 20 mm, consisting of 75 × 75 pixels, and took less than 5 min to acquire. At each pixel position an entire THz waveform was acquired.

The orientation of the sample on the imaging window and location of the inked boundaries were recorded and a photograph was taken by a camera within the system so that the THz image could later be correlated with the histology section. After imaging, all samples were prepared routinely by being formalin fixed, wax embedded, and then sliced and mounted onto histology slides and stained.

A terahertz impulse function was obtained from each raw terahertz waveform by deconvolving the system response as described in Ref. 33. Each impulse function, referred to as a THz pulse, contained 512 time-domain points which covered a time range of 33.8 ps.

2.3.

Correlation of Terahertz Images with Histology

The histology was read by a pathologist and all areas of pathology (tumor, normal tissue, adipose, fibrosis, etc.) were identified along with inked margins and other landmarks. The histology was then oriented to the THz images in an image-processing package so that the areas of pathology could be mapped directly onto the THz images. An example to demonstrate the orientation of the terahertz image with histology is given in Fig. 1.

Fig. 1

Correlation of terahertz image with histology. The scale facilitates size correction, and inks and photographs enable orientation.

This mapping enabled the regions of the terahertz image to be identified and grouped as normal, tumor, or alternative pathologies. Regions of interest within each image were then applied to extract terahertz pulses with each type of pathology. These extracted pulses were grouped into a database containing pulses grouped as normal, tumor, and other pathologies. The database contained 31,236 pulses from normal tissue and 30,405 pulses from tissue with tumor. These pulses from the database were used in the subsequent data reduction and classification algorithms.

2.4.

Data Dimensionality Reduction and Classification

Given the large numbers of pulses and number of values for each pulse (512 data points) it was important for the application of classification techniques to reduce the dimensionality of the data and the size of the entire data set. Data dimensionality reduction was explored in two ways. The first method used parameters, often applied in THz imaging, that were heuristic and based on characteristics of the terahertz pulse in the time and frequency domains.²¹ The second method made use of the well-known unsupervised method, PCA. Classification, for outputs of all data reduction techniques was performed using the SVM classification algorithm.

2.5.

Data Reduction Method 1: THz Image Parameters

THz data can be represented in either the time-domain or, through Fourier transform, in the spectral domain. The mean time-domain pulses for the dataset for normal and tumor tissue, together with their power spectra, are shown in Fig. 2 to illustrate the often subtle differences between the two.

Fig. 2

The mean pulses for normal fibrous tissue (black line) and tumor (grey line) in the (a) time-domain and (b) frequency domain.

Ashworth and colleagues³⁴ performed THz spectroscopy measurements on normal and tumor breast tissue which they entered into a computer model³⁵ to simulate pulses which showed very similar profiles to those in Fig. 2. It appears that physical and biochemical changes in the tissue lead to changes in the dielectric properties that affect the absorption coefficient and refractive index across the entire frequency range in such a way that the pulse shape is altered with broadening and retardation after reflection from tumor tissue, when compared with normal tissue. This is similar to what was found with THz data from basal cell carcinoma.³⁶

Parameters can be derived from these pulses that can be used to characterize the response of the underlying tissue. Hence parameters that are sensitive to these changes in pulse shape will be useful in enabling differentiation of tumor from normal tissue with THz imaging.

A large number of parameters may be derived from the THz pulses and spectra based on time or frequency domain characteristics. To reduce this number we used two steps. First we determined the optimal time and frequency index points at which these parameters could be evaluated, and secondly we removed any highly correlated parameters.

To aid analysis and parameter calculation, all pulses were baseline corrected, and the peaks aligned in time. For the purposes of parameter definition, the term $E_{\min }$ is defined as the amplitude at which the pulse is a minimum. All pulses are aligned in time so that $E_{\min }$ occurs at the same time in the time range. The time at which $E_{\min }$ occurs is defined as $T_{E\min }$ and is identically set to zero as the reference point on the time-domain axis. All time indices relating to the parameters are relative to this point unless otherwise specified.

For many of the characteristics, such as the value of the amplitude of the pulse at a given time, there was a large range of time or frequency data points at which the parameter could be evaluated. For example, the amplitude of the pulse at a given time point could be evaluated at 512 possible time points over a time range of 33.8 ps. To refine the selection of optimal time and frequency points for classification, a receiver operator characteristics (ROC) analysis^37–39 was performed on the time and frequency range for each dependent parameter. ROC is a fast and efficient method for determining the differentiation ability of each parameter, so it was suitable for determining the optimal parameters from the large number of time and frequency data.

ROC determines discrimination ability by adjusting the discrimination threshold and plotting, for the range of threshold values, sensitivity, or true-positive rate (normals classified as normal), against 1-specificity, which is the false positive rate (tumors classified as normal). This produces a graph, and the area under the curve (AUROC), is a measure of the predictive power of each parameter. The higher the AUROC value the better the classification strength of that parameter. To demonstrate this procedure, Fig. 3 shows the ROC and AUROC for the parameter P7; the amplitude of the pulse parameter over a range of time indices. In this way, the best performing time index and frequency index points were chosen for parameters related to the values of the pulse or the spectra.

Fig. 3

ROC analysis to evaluate the discrimination ability of parameter P7, the amplitude of the pulse at a given time. (a) AUROC values, which are a measure of the discrimination ability, for the time range around the time at which the pulse is minimum ($T_{E\min }=0$). From this analysis it can be seen that the time that gives the largest AUROC is $-0.26\mathrm{ps}$, hence the amplitude at this time was used as a parameter feature. (b) The ROC curve for parameter P7 with the time index of $-0.26\mathrm{ps}$.

To determine independent parameters, the linear correlation coefficient was calculated between all parameters and the correlated parameters removed, leaving the better performing parameter, according to the AUROC. The ten remaining optimal and uncorrelated parameters were then used for classification with the SVM technique. Table 1 contains the list of parameters and their descriptions, whereas Fig. 4 illustrates an example terahertz time-domain pulse and power spectrum, with sections of the pulse that were curve-fitted by various functions to give parameters related to aspects of the pulse shape.

Table 1

Time-domain and frequency domain terahertz parameters for the ten parameter space.

Fig. 4

(a) The sections of the pulse in the time-domain that were fitted with linear (section 1), quadratic (section 2) and exponential (section 3) functions in the time-domain. Section (1) was fitted with linear regression for parameter P4. Section (2) was fitted with a quadratic for parameter P5, whereas section (3) was fitted with an exponential for parameter P6. (b) For parameter P10, the log power spectrum was fitted with a linear fit.

The ten parameters (Table 1) were assessed for their individual ability to discriminate between normal and tumor pulse responses by applying the SVM classifier (Table 2). The values obtained from the SVM, given in Table 2, were then used to select the order in which parameters were added to the combinations of parameters that could give better classification than the individual parameters alone. It is not immediately obvious in which order parameters should be added in the subsequent combinations to give the optimal results. To investigate whether order was important we compared a number of different combinations. We compared combining two individual parameters that were the best performing, according to the individual SVM results, against combining two that were the worst performing. The combinations that used the best performing individual parameters always equaled or outperformed those using combinations of lower accuracy parameters (results not shown). With this justification we use the combinations formed by adding the individual parameters in order according to their classification accuracy from the SVM, starting with the best performing parameter, and adding one parameter at a time until all ten parameters were included in the combination. Results for the combinations of parameters are given in Table 3.

Table 2

Performance of the SVM on the parameter feature space.

Table 3

Performance of the SVM on combinations of parameter features for method 1.

2.6.

Data Reduction Method 2: PCA on Database Pulses

PCA is a well-known unsupervised linear feature extraction algorithm traditionally used to reduce data dimensionality. It involves an orthogonal transformation of the feature space to form a set of uncorrelated values with each successive principal component having the highest variance possible under the orthogonality constraint.

For the purposes of this study the Matlab implementation of PCA was used (MATLAB 7.7, The Mathworks Inc., Natick, MA, 2008). PCA was applied to the database of tumor and normal pulses. Combinations of the projections in the principal components were then used for classification with the SVM. The first combination used only the first principal component, the second combination used the first two principal components, and so on, up to the tenth combination which used the highest ten principal components. Results for the SVM classification on these ten combinations of principal components are given in Table 4.

Table 4

Performance of the SVM with increasing numbers of principal components from the time-domain THz pulses for method 2.

2.7.

Data Reduction Method 3: PCA on Parameter Feature Space

With method 1, although we combined parameters in the order according to the individual accuracy from the SVM, there was no obvious optimal order in which to combine them. To investigate if there is a better grouping based on variation under the orthogonality constraint, we applied the PCA to the ten parameter space described in method 1. There is evidence in the literature to suggest that PCA performs better on discrete data, such as the parameter feature space, compared to time-domain waveforms such as the pulses used in method 2.²⁶ Combinations of the principal components (PCs) from the highest 10 PCs were used for the SVM classification in the same manner as method 2. Results for the SVM classification on the PCs of the parameter space are given in Table 5.

Table 5

Performance of the SVM on the PCA parameter space for method 3.

2.8.

Classification of TPI Data by SVM

All classification was performed using the SVM algorithm to find decision boundaries between the classes in the reduced dimension feature spaces. This nonlinear classification function uses an iterative method to maximize the margin between the classes with an estimated optimal boundary by selecting a minimum number of support vectors.⁴⁰

For SVM classification in this study we used the Matlab implementation developed by the researchers at the National Taiwan University⁴¹ with the radial basis function. The kernels and coefficients were optimized by grid search for each of the methods, giving optimized cost and gamma coefficients of 32 and 128, respectively for methods 1 and 3, the parameter space and PCA of the feature space; and 8 and 8 for method 2, the PCA on the pulses.

The SVM was trained on 75% of the normal and tumor database pulses, randomly selected, and tested on the other 25%. From the SVM, we obtained the total overall classification percentage accuracy, as well as the accuracy rates for classifying tumor and normal pulses.

The precision or error on the classification accuracy, which is dependent on sample size and the desired confidence limit, is given by Fisher and van Belle⁴² as:

3.1.

Method 1: SVM Classification of Parameter Features

The ten parameters used for classification are compared in Table 2 for their ability to classify tumor and normal tissue from THz pulses. It is apparent that certain parameters more accurately identify tumor tissue than normal tissue, for example parameters 1 and 2, and vice versa (e.g., parameters 8 and 9). Combinations of parameters clearly perform better at classifying tumor and normal tissue than single parameters alone, and the more parameters used in a combination, the better classification accuracy as shown in Fig. 5.

Fig. 5

Comparison of SVM total accuracy for method 1, parameter space (Δ, grey line), method 2, PCA on the pulses (□, dashed line), and method 3 PCA of the parameter space (○, full line).

3.2.

Method 2: Classification of PCA for THz Pulses

The first principal component of the database pulses was compared with the difference between the mean database normal and tumor THz pulses. As might be expected, the highest variance principal component had a very similar profile to the difference between the mean tumor and normal THz pulses (Fig. 6). The first three principal components are illustrated in Fig. 7, showing the importance of the region around the pulse minimum and the decay, accompanied with pulse broadening on either side as observed in the parameter space. Table 4 shows that the SVM classification accuracy increases for greater numbers of PCs. Fig. 5 shows that the increase in the classification accuracy is largest for the first three PCs, which contribute the most to the overall variance, and then increases steadily with the addition of the next seven PCs. This continued increase could be due to the contribution of variation on the pulse being spread widely across all ten PCs.

Fig. 6

The first PC of the pulses (grey line) compared to difference of mean normal and tumor waveforms (black lines). Both have been normalized to their absolute maximum value for comparison.

Fig. 7

First three PCs of the pulses.

3.3.

Method 3: Classification of PCA for THz Parameters

Table 5 shows the SVM classification accuracy for the PCA on the parameter space. Again, increasing numbers of PCs leads to improved classification accuracy. Also, as before, the first three PCs contribute the most to the overall variance (80%), and therefore contribute to the classification accuracy. However, after these first three PCs, the subsequent PCs appear to contribute less to the overall variance, and hence classification.

4.1.

Classification of THz Image Data Using SVM

As a visual demonstration of the classification technique we apply the data reduction and classification to several example images. During breast surgery the intraoperative THz probe will be used to classify pulses and not form images as such. However, the exercise of forming and classifying pixels in images can inform the practical development of portable clinical devices based on THz technology, for example the intraoperative breast probe.

Because data reduction method 3, the PCA on the parameter feature space using 10 PCs, leads to best classification accuracy, this method was used in the images. Training of the SVM model was done using the PCA results from all pulses in the database. The SVM output was a probability that each pixel in the image was tumor or normal. The values near 1 represented the probability that the pixel was classified by the SVM model as tumor, whereas a value near zero represented tissue classified as normal.

Prior to data reduction and classification, there were two structures in the images that were excluded from the analysis. These were the adipose tissue, and gaps in the tissue due to the internal structure or damage due to the slicing of the tissue. These regions were determined by an algorithm that identified them through their pulse shapes, which are very different to normal fibrous and tumor pulses. To represent the results, a classification probability map was plotted, with a neutral background for all pixels that were calculated to be adipose or gaps; these results are presented in Fig. 8.

Fig. 8

Demonstration of SVM classifications on example THz images of breast tissue. Probability maps (top row), obtained from the SVM classification using the PCA on the parameter space are displayed with red representing the probability that the pixel is tumor, whereas green represents the probability from the SVM that the pixel is normal. Adipose and gaps in the tissue are removed from the analysis and shown as white background. The bottom row shows the corresponding histology for the various cases (a) normal tissue, (b) normal tissue, (c) mixed normal tissue and multifocal tumor (invasive ductal and DCIS), (d) mixed normal and tumor (invasive lobular), (e) tumor (invasive lobular), (f) tumor (invasive ductal), and (g) tumor (invasive ductal and necrosis). The square box around the classification maps and histology represents 20 mm by 20 mm.

Adipose tissue can be easily identified and removed in images with a simple algorithm, because the pulse shape due to its refractive index and absorption³⁴ is so different to normal or tumor fibrous tissue. However, at the border of adipose and fibrous tissue, an edge effect occurs that is more complicated to deal with in classification. At this border, the pulse shape undergoes a transition where the finite dimension of the THz pulse profile samples a portion of both adipose and fibrous tissues. This edge has been observed with other THz imaging applications^45,46 and is frequency dependent because the beam diameter, and therefore resolution, is smaller at higher frequencies.⁴⁷ This distortion of the recorded pulse shape leads to features not recognized in the classification database, and hence not well classified by the algorithms. It may be possible to train the SVM classifier for edges to deal with these features; however, in practice this may not be necessary as it does not appear to affect the classification significantly. In the few cases where it does appear to be an issue, for example in case 2, the edge is mostly misclassified on the top or side edges, suggesting there is a directional component, probably due to the 30-deg angle of incidence of the THz pulse.⁴⁸ In practice with the concept of the intraoperative probe, the angle of incidence is much lower, approximately 5 deg, so the effect is likely to be less profound. However, if the effect remains and causes ambiguity, the surgeon would simply be able to change the orientation and determine if it is an edge effect or not. There is also a second type of edge effect present in case 1, whereby the edge of the tissue presents a discontinuous boundary. It occurs in this case because the tissue has been removed by incision and has edges that are thin compared to the thicker central portion. This type of edge effect would not normally occur during surgery in situ as the surrounding breast tissue would still be intact. The combination of the discontinuous edge and thin profile in the third dimension produced an altered profile not accounted for in the database.

In general the image pixels appear to be well classified. In most cases there appear to be small numbers of pixels, around 10% to 15%, misclassified, which is similar to the accuracy rate seen in the classification algorithm on the database pulses. In some images however, e.g., case 6, nearly all the pixels are misclassified, highlighting the importance of considering interpatient variability in this approach.