2.1.

Label-Free Tissue Scanner

Let $T={\mathrm{e}}^{i\varphi }$ be the transmission of the biopsy and $\varphi$ be the optical phase shift introduced by the sample. We have $\varphi =h\mathrm{\Delta }n(2\pi /\lambda )$. Here, $\lambda$ is the central wavelength of $552\mu \mathrm{m}$, $\mathrm{\Delta }n$ is the refractive index difference, and $h$ is the sample thickness. Our system measures $\varphi$ using a commercial inverted phase contrast microscope (Axio Observer, Z1, Zeiss Inc.) connected to an external SLIM module (CellVista SLIM Pro, Phi Optics, Inc.). A schematic of our optical setup is shown in Fig. 1(a). Under uniformly coherent illumination, the imaging field at the output port of the microscope is a magnified version of the transmission, i.e., $U_t=T$. This total field is Fourier transformed into ${\tilde{U}}_t$ by the lens ${\mathrm{L}}_1$ onto the surface of a spatial light modulator (SLM) (Boulder Nonlinear Inc.). On this surface, the spatial spectrum ${\tilde{U}}_t$ is spatially separated into two different components, the nonscattering component, ${\tilde{U}}_o$, and the scattering one, ${\tilde{U}}_s={\tilde{U}}_t-{\tilde{U}}_o$. The nonscattering component ${\tilde{U}}_o$ matches the support of the condenser phase annulus and the phase ring of the phase contrast objective. Meanwhile, the scattering component ${\tilde{U}}_s$ covers the rest of the aperture of the objective. These Fourier components are inversely Fourier transformed by a secondary lens ${\mathrm{L}}_2$ into two fields, the nonscattering field, $U_o$, and the scattering field, $U_s$. The field $U_o$ spatially averages the total field $U_t$ into $U_o={⟨U_t⟩}_{\mathbf{r}}$. The scattering field, $U_s$, on the other hand, is given as $U_s=U_t-U_o$. The SLM further retards $U_o$ by $n\pi /2$, with $n=\mathrm{0,1},\mathrm{2,3}$, resulting in four interference intensities on the camera plane of

Fig. 1

Optical setup and working principle: (a) SLIM optical setup, (b) four intensity images and the phase map computed by combining four frames, and (c) the resulting SLIM image.

Figure 1(b) shows four intensity images for a small section of a prostate biopsy. The first inset is the original phase contrast image at zero external phase modulation while the last inset is the bright-field image with a phase modulation of $3\pi /2$. Figure 1(c) is the extracted phase map for $\varphi$ from these intensities. All details in the intensity images are still visible in the SLIM image with very good contrast. More importantly, the SLIM image is not susceptible to variation in the illumination, as it only captures intrinsic information of the sample.

2.2.

Imaging the Tissue Microarray

The TMA is provided by the Co-operative Prostate Cancer Tissue Resource from the College of Medicine of University of Illinois at Chicago. It consists of 368 prostate cores (one core per patient) with several diagnosis results, including normal, BPH, high-grade prostatic intraepithelial neoplasia (HGPIN), and Gleason scores varying from $2+2$ to $5+5$. After being deparaffinized, unstained cores were first imaged using SLIM under a $40\times$ magnification. The phase images are stitched together to generate one high-resolution image per core. Each such image has $\mathrm{10,000}\times \mathrm{10,000}\text{pixels}$ with a pixel ratio of 14 pixels per micron. Afterward, the cores were stained with H&E and scanned by a bright-field tissue scanner. The results of all cores can be found in Fig. 2(a). Figures 2(b) and 2(c) are a zoomed-in H&E image of a core and its corresponding SLIM image, respectively.

Fig. 2

H&E image of the whole TMA with diagnosis results. (a) H&E image of the whole TMA slide consisting of 368 cores. (b) A zoomed-in H&E image of a prostate core with annotations. The region highlighted in green represents normal glands, region in blue is Gleason grade 3 prostate cancer glands, and region in red corresponds HGPIN. (c) A zoomed-in SLIM image of the same core as in (b), obtained prior to staining. Morphological features in the H&E image are recapitulated by SLIM.

Figure 3 displays three H&E [Fig. 3(a)] and SLIM [Fig. 3(b)] images of three cores in the TMA. To provide ground truth for automatic diagnosis, different regions of interest (ROIs) in each core are studied and color-coded by certified pathologists for the diagnosis results using H&E images. In Fig. 3, green color indicates normal areas, red color indicates HGPIN, and blue color indicates tumors. Based on these annotations, corresponding regions in the SLIM images are extracted and further used for the automatic diagnosis.

Fig. 3

H&E images versus SLIM images. (a) H&E images of three cores in the TMA. Each core in the training data set includes an annotation of diagnosis results. (b) Corresponding SLIM images of those in (a).

2.3.

Automatic Diagnosis Framework

To obtain the automatic diagnosis from SLIM images, we use an approach summarized in Fig. 4. First, texture-based features are extracted for each pixel in the SLIM image and passed into a classifier to do automatic segmentation based on pixel classification. Each pixel is assigned into one of three following classes (gland, stroma, and lumen). Second, using the label map obtained from the previous step, morphological and phase-based features are evaluated for all glands in the current field of view and its surrounding stroma. These features are later passed into a subsequent classifier to produce diagnosis results, as described below.

Fig. 4

Automatic diagnosis scheme of steps from an SLIM image to a diagnosis result.

2.4.

Feature Extraction from Phase Images

To capture the texture expression of the biopsies, we use the “texton” framework proposed by Julesz³⁷ and later expanded by Leung and Malik.^38,39 The framework has demonstrated great success in solving several computer vision problems, e.g., material classification and characterization,^38–41 due to its ability to accurately imitate mechanism of human’s textural perception.^37,42,43 In this work, the framework is used to train a texton dictionary [see Fig. 5(a)] from a set of training images and extract a feature vector for each pixel as follows.

Fig. 5

Feature extraction from SLIM images: (a) extracting the texton dictionary from a training set of SLIM images, (b) feature extraction using vector quantization and the trained dictionary of textons from (a), (c) RF training using descriptors obtained from (b), and (d) pixel classification using a trained RF from (c).

First, each SLIM image $\varphi$ is convolved with a Leung–Malik filter bank,⁴³ consisting of $L$ filters $h_1,\dots ,h_L$ to generate $L$ filter responses $\varphi *h_1,\dots ,\varphi *h_L$. The bank consists of 90 filters, $L=90$, with 10 symmetric ones and 80 directional ones, oriented at eight different angles over five scales. The directional filters are generated from the first- and second-order derivative of Gaussian kernels with an elongation factor of 3. After this step, each pixel has a feature vector of dimension 90. A phase image in the training set generates totally $3072\times 3072$ feature vectors. To reduce the size of the training set, a subset is formed from four million feature vectors randomly selected from a larger pool of all feature vectors, accounting for 0.18% of the total number of feature vectors. Finally, the K-mean clustering algorithm with $K=50$ is applied on this subset to divide it into $K$ clusters whose centers are chosen as textons. Here, the value of $K$ is chosen to balance between the complexity of the model and the estimation error, i.e., avoiding cases where there is not enough textons to capture texture variation or those when some textons come from clusters only associated to a very small number of filter responses.

In feature extraction, given an input image ($I$), one wants to obtain a set of descriptors for each pixel in the image. Figure 5(b) shows how to calculate these descriptors. First, filter responses of the image $I$ to the Leung–Malik filter bank are computed. Then, we apply vector quantization to associate the filter response at each pixel to the closest texton in the dictionary obtained in the previous step. For each pixel $i$, let us use $t_i$ to denote the index of the closest texton. By definition, $t_i$ can take one of the $K$ following indices {$\mathrm{1,2},\dots ,K$}. The output of this step is an indexing map where each pixel is assigned to a number, the index of the closest texton. Using this indexing map, for each pixel $i$, we further obtain a histogram of texton indices evaluated over pixels in its neighborhood. To control the trade-off between the richness of texture information and locality of the descriptor, we apply a Gaussian weight to the histogram calculation where larger histogram contribution is given to pixels closer to the center of the neighborhood. By trial and error, we determine that a neighborhood radius of $\sigma =45\text{pixels}$ is suitable to characterize the pixel. This radius corresponds to $\sim 10\mu \mathrm{m}$ in the downsampled SLIM image.

2.5.

Random Forest for Automatic Segmentation

To perform image segmentation, we use a random forest (RF) classifier, a method introduced by Breiman^44,45 and Shotton et al.⁴⁶ This classifier has shown success in several problems, such as object segmentation,^40,47,48 human-pose estimation,^46,49 and medical image analysis using magnetic resonant imaging,^50–52 due to its ability to reduce the tree dependence with “feature bagging” and “bootstrap” sampling, i.e., random sampling with replacement. In this work, we use the RF to classify each pixel in the image into one of three classes, i.e., epithelial gland, connecting stroma, and lumen. Lumen pixels are classified first based on the proximity of their phase values to that of the background. Then, remaining pixels are classified into either gland or stroma. Here, we train an “extremely randomized” forest⁵³ of 50 trees. Our implementation is written in MATLAB^® with the MexOpenCV wrapper that allows us to call OpenCV routines. The wrapper is obtained from Ref. 54. The training set consists of 4.92 million histograms with two possible labels (“gland” or “stroma”). These histograms of texton are randomly sampled from a larger pool of 2.2 billion histograms. Each tree in the forests is trained on 11 randomly selected features (out of 50, the total number of textons). For each feature of interest, 100 possible thresholds are considered between the minimum and maximum of the feature values for splitting. Let us use $S_n$ to denote the training data set of samples reaching node $n$’th. At this note, this set is partitioned into a left set, $S_{\mathrm{l}}$, and a right one, $S_{\mathrm{r}}$, based on an optimal selected feature ${\nu }_n$ and its optimal threshold $t_n$, namely $S_{\mathrm{l}}=\{s\in S_n|v_n<t_n\}$ and $S_{\mathrm{r}}=S_n\backslash S_{\mathrm{l}}$. Here, (${\nu }_n,t_n$) are chosen to maximize the expected gain of information on category, i.e.,

After the RF has been trained, automatic segmentation results can be obtained by classifying each pixel in the query image to either gland or stroma using its pixel descriptor [Fig. 5(d)]. The optimal class $g^{*}$ for each pixel is determined by summing the class likelihood values over all the trees and picking the class that maximizes the combined score, i.e.,

2.6.

Postprocessing of Segmentation Results

The output of the pixel classifier is usually noisy. We use the following procedure to obtain good segmentation results and resolve the glandular fusion, which is more frequent in high-grade prostate carcinoma.

Fig. 6

Steps of image postprocessing for automatic segmentation: (a) the stroma likelihood map produced by the RF classifier, (b) the raw binary gland map, (c) the raw binary nongland map (in white), (d) the distance map from each gland pixel to the nearest nongland pixel, (e) the oversegmentation map produced by watershed segmentation, (f) the map of gland separating lines, (g) the refined segmentation map where some similar regions in (e) have been merged, and (h) the final segmentation map.

3.1.

Automatic Segmentation

Figure 7 shows automatic segmentation results overlaid with the SLIM images. It can be seen that the label map has a good correlation with the H&E images. Figure 8 shows other segmentation examples with H&E and SLIM images of increasing Gleason grade from 3 to 5. Their automatic segmentation results from SLIM images are further shown in Fig. 9.

Fig. 7

H&E versus automatic segmentation: (a)–(c) H&E images of three cores and (d)–(f) corresponding automatic segmentation results overlaid on the top of the SLIM images.

Fig. 8

H&E versus SLIM: (a)–(c) H&E images of three cores in the TMA that have Gleason grade of 3, 4, and 5. (d)–(f) Corresponding SLIM images of these cores. Grayscale bar represents phase shift in radians.

Fig. 9

H&E versus automatic segmentation: (a)–(c) H&E images of three cores of different Gleason grades in Fig. 7, as indicated. (d)–(f) Automatic segmentation results of the cores overlaid on the SLIM images.

To quantify the performance of our segmentation, we summarize in Fig. 10(a) receiver operating characteristic (ROC) curves for the segmentation of different diagnosis groups using 10-fold cross-validation. Since these ROC curves are very close to each other, it is difficult to compare their performance. To address this issue, we show in Fig. 10(b) a zoomed-in image to the upper-left corner of Fig. 10(a), where the true positive rate and the false negative rate range from [0.75, 1.0] and [0, 0.5], respectively. Figure 10(c) shows the corresponding AUC values for each group. The ground truth for the ROC evaluation is created by manually labeling glandular and stromal regions directly on the SLIM images after validating them with H&E images. In the noncancer cases, the AUC is at least 0.97, which indicates that gland and stroma pixels are classified with high accuracy. Meanwhile, in the malignant cases, as cancer progresses to higher grades, e.g., from $3+3$ to $5+5$, the AUC reduces from 0.98 to 0.87. This result can be explained by the fact that more glandular distortions and deformations are observed at higher grades, which leads to a reduction in discrimination between stroma and glands. Therefore, it is not surprising that the classifier has the smallest AUC with a Gleason score of $5+5$. However, these high grades are very easily diagnosed by the pathologist and, thus, do not represent our main focus. Using these segmentation results, we solve the automatic Gleason grading problem, with particular emphasis on discriminating between grades 3 and 4.

Fig. 10

Automatic segmentation performance. (a) ROC curves for classifying gland versus nongland pixels for all diagnosis groups. The accuracy reduces in the case of higher Gleason grade due to the less defined glands. (b) Zoomed-in plot for a red-dashed region of (a). (c) Bar plot of AUC values of the ROC curves for all diagnosis groups.

3.2.

Automatic Gleason Grading

Next, we demonstrate the use of our technique in clinical applications with automatic Gleason grading. To generate the ground truth for the training, all cores were reviewed, manually marked, and graded unanimously by two trained pathologists. Figure 2 shows an example of markup results for all cores: 129 regions with Gleason grade 3, 92 regions with grade 4, and 75 regions with grade 5. Since Gleason grades 2 and 5 are rarely diagnosed, we study the automatic diagnosis problem of differentiating Gleason grade 3 versus 4. It has been shown by Allsbrook et al.⁵⁶ that grading 3 versus 4 has a reproducibility problem due to interobserver variation. Also, there is the crucial turning point from active surveillance to aggressive treatment when moving from Gleason grade 3 to 4. Distinguishing between 2 and 5 is not a problem of clinical relevance since it is quite trivial. For example, the authors report an experiment where 38 biopsies with known “consensus” Gleason grade were sent to 41 pathologists to measure interobserver variability. The result was that Gleason grade 4 was undergraded by 21%. Furthermore, there was consistent undergrading of Gleason scores of 5 to 6 (47%), 7 (47%), and 8 to 10 (25%). Clearly, a computer-driven, unbiased procedure for grading is a potential way to tackle this challenge.

Figure 11 shows five different types of features extracted from each region. These features include the mean of glandular distortion, the fusing ratio of glands, the mean number of lumen areas per gland, the coefficient of variation for gland variation, and the mean of stroma anisotropy. Other features include the maximum number of lumen areas, the median of the glandular distortion, median of stroma anisotropy, and the mean circularity. Some features are explained here in detail.

Fig. 11

Feature extraction for automatic diagnosis. Each subfigure shows a feature, how it is calculated and the distributions of the feature values for G3 (blue) and G4 (orange). (a) Mean distortion feature: $D$, mean distortion of a gland; $P$, perimeter of a gland; and $A$, area of a gland. (b) Average number of lumens: ANL, average number of lumens; NL, number of lumen; and NG, number of glands. (c) Average stroma anisotropy. (d) Fusing ratio: FR, fusing ratio; TG2, total areas of glands with at least two lumens; and TG, total area of all glands in current field of view. (e) Coefficient of variation: CV, coefficient of variation; STDA, standard deviation of areas of all glands; and MA, mean of the areas of all glands.

These features are designed to measure the distortion, area homogeneity of the glands, and the amount of gland fusion. They are often used by pathologists for Gleason grading.^5,6,60 Each ROI is characterized by a feature vector of 11 elements. Diagnosis grades from pathologists are used as ground truth for automatic diagnosis. After computing the feature vectors for ROIs, we use logistic regression to classify the grade 3 ROIs versus grade 4. This classifier models the conditional probability of the Gleason grade ($G$) given the feature vector $\mathbf{X}$ as $p(G=G_3|\mathbf{X})=\exp ({\mathbf{\beta }}^T\mathbf{X})/[1+\exp ({\mathbf{\beta }}^T\mathbf{X})]$, and $p(G=G_4|\mathbf{X})=1/[1+\exp ({\mathbf{\beta }}^T\mathbf{X})]$. Here, $\mathbf{\beta }$ is a coefficient vector estimated from the data by maximizing the log-likelihood function of the observations of the training set; see Ref. 61 for more details.

The performance of the logistic regression classifier is shown in Fig. 12 in terms of the ROC using leave-one-out cross-validation. This type of validation trains classifiers using all samples in the dataset except one. A remaining sample is used for testing. The classifier provides two conditional probabilities $p(G_{\mathrm{t}}=G_4|{\mathbf{X}}_{\mathbf{t}})$ and $p(G_{\mathrm{t}}=G_3|{\mathbf{X}}_{\mathbf{t}})$ for the testing sample $X_{\mathrm{t}}$. Here, $G_{\mathrm{t}}$ is the Gleason grade of the test sample. Then, we alternate the selection of the testing sample to make sure every sample in our dataset becomes a testing sample once. Finally, the conditional probabilities $p(G_{\mathrm{t}}=G_3|{\mathbf{X}}_{\mathbf{t}})$ and $p(G_{\mathrm{t}}=G_4|{\mathbf{X}}_{\mathbf{t}})$ of all samples in the dataset are combined to produce the ROC curve for Gleason grade 3 and 4 classification. This curve has an AUC value of 0.82. Note that this error is well within that for interobserver variability reported by Allsbrook et al.⁴⁴ The reason why this ROC curve has a stair-case response while the ROC curves shown in Figs. 10(a) and 10(b) are smooth is because of the difference in the number of samples used. To obtain the ROC curve in Fig. 12, we used a dataset consisting of 141 samples. This number is much smaller than that of texture descriptors, about 3.9 millions, used for segmentation in Fig. 10. Since the number of samples is large, the steps of these ROC curve are very small. This fact makes them look smooth. It can be also seen from the curve that, to detect Gleason grade 4 at an accuracy of 90%, the false positive rate will be $\sim 60\%$. The inset presents a horizontal bar plot of AUC values when the classifier is trained separately on each individual feature, also using leave-one-out cross-validation. The two largest AUC values are obtained on the coefficient of variation for the areas of the glands and the fusing ratio. Our results demonstrate that label-free imaging and machine learning can provide an objective alternative to pathology, even in the case of difficult tasks, such as classifying Gleason grade 3 and 4.

Fig. 12

ROC curves for classifying regions with Gleason G3 versus G4 under leave-one-out cross-validation. These features are extracted from all glands within the ROIs. Diagnosis results from pathologists are used as ground truths. The inset shows AUCs of classification using individual features.