2.1.

PCANet Feature Extraction

As a modified version of the typical deep CNN, the PCANet method³³ adopts a series of PCA convolution filters to extract the features from an input image. A commonly used PCANet model consists of two PCA filter convolution stages and one output stage. Figure 2 shows the typical architecture of the PCANet method, as also described in the following.

Fig. 2

Illustration of the PCANet architecture with two convolution stages and one output stage.

In the first convolution stage, given $H$ training images $X_h\in {\mathbf{R}}^{M\times N}$, $h\in \{\mathrm{1,2},\dots ,H\}$, we convolve each image with $L_1$ PCA filters in a patch-based pattern to obtain $L_1\times H$ feature maps $X_h^{1,l}$, $l\in \{\mathrm{1,2},\dots ,L_1\}$. Specifically, we first extract the patches $x_h^i$ (of size $n_1\times n_2$) for each image $X_h$ and remove the mean for each patch. Then, the extracted patches of each image $X_h$ are vectorized and can be combined into a matrix $P_h=[{\overline{P}}_h^1,\dots ,{\overline{P}}_h^i,\dots {\overline{P}}_h^{M^{*}{N}^{*}}]$, where ${\overline{P}}_h^i\in {\mathbf{R}}^{n_1{n}_2\times 1}$ is a vector for the related patch, $M^{*}=M-\lceil n_1/2\rceil$, $N^{*}=N-\lceil n_2/2\rceil$, and $\lceil n\rceil$ gives the smallest integer greater than or equal to $n$. The matrices of all the training images are also constructed as a matrix:

After that, we compute the eigenvectors of ${PP}^T$ and select the first $L_1$ principal eigenvectors as the $L_1$ PCA filters $W_l^1$, $l=\mathrm{1,2},\dots ,L_1$. Finally, $L_1$ filters are separately applied to each training image $X_h$:

For each feature map $X_h^{1,l}$ from the first stage, the second stage convolves this map with $L_2$ PCA filters to obtain $L_2$ feature maps $X_h^{2,l}$, $l\in \{\mathrm{1,2},\dots ,L_2\}$. Similar to the first stage, we also first extract the mean-removed patches from all $L_1\times H$ maps $X_h^{1,l}$ and use the vector version of patches to construct the matrix $U$. Then, we compute the eigenvectors of ${UU}^T$, and select the first $L_2$ principal eigenvectors as the $L_2$ PCA filters $W_l^2$, $l=\mathrm{1,2},\dots ,L_2$. As in Eq. (2), by applying the $L_2$ PCA filters on each feature map $X_h^{1,l}$ from the first stage, we can obtain $L_1\times L_2\times H$ feature maps. Given an input image $X_h$, the first stage has $L_1$ feature maps, and the second stage can create $L_1\times L_2$ feature maps.

In the output stage, we first binarize the $L_1\times L_2$ feature maps of the second stage using a Heaviside step function and apply the hashing encoding to convert them into $L_1$ integer-valued images $Z_h^l$, $l\in \{\mathrm{1,2},\dots ,L_1\}$.³³ Then, each integer-valued image $Z_h^l$ is partitioned into $B$ blocks (of size $n_{B1}\times n_{B2}$), in which the local histogram (overlapping or nonoverlapping) is calculated. After that, all $B$ histograms are concatenated into a vector ${\mathrm{Bhist}(Z}_h^l)$. Finally, the feature vector of the one input image $X_h$ can be extracted as

The extracted feature vector can be fed into a classifier (e.g., SVM or ELM) for the classification.

2.2.

ELM Classifier

ELM is a very efficient supervised learning model, and its objective is to find a decision rule for the classification.^35,36 To be specific, let ${{\{x}_i,y_i\}}_{i=1}^H$ be the training set. $\{x_1,x_2,\dots ,x_H\}\in {\mathbf{R}}^{m\times 1}$ denote the $H$ input training samples and $y_i={[0,\dots ,1,\dots ,0]}^{\mathrm{T}}\in {\mathbf{R}}^{G\times 1}$ represent the corresponding class labels, which are a vector with one for the true label and zero entries for the others. $G$ denotes the number of training classes. In general, the ELM aims to simultaneously minimize the training error and the norm of output weights of the objective function

As described in Refs. 34 and 37, kernel can further transform the features into a higher dimensional feature space and thus improve the discriminative capacity for the features. Since the mapping $\varphi (•)$ in the ELM learning is represented by the inner function, a kernel function $K$ can be defined by

Finally, by incorporating Eq. (6) into Eq. (5), the decision rule of the ELM for any test sample $x$ is determined by

3.1.

Feature Extraction with PCANet

Since OCT images have large variations in intensity ranges, we first conduct the intensity normalization on each B-scan of the training and testing OCT volumes, which linearly rescales the intensity values to [0, 1]. Then, as in Ref. 23, for training and testing OCT images, we adopt the BM3D algorithm³⁸ to remove the noise of each B-scan and use the flattening to reduce the variations for imaged retinas in their angles of inclination, positions, and natural curvatures among B-scans. After the above preprocessing steps, the PCANet feature extraction consists of training and testing phases. Specifically, in the training phase, we first input $T$ 3-D OCT volumes $V_t^{\text{train}}\in {\mathbf{R}}^{M\times N\times H_t}$, $t=\mathrm{1,2},\dots ,T$ into the PCANet training model. Each volume contains $H_t$ B-scans $X_h^{\text{train}}\in {\mathbf{R}}^{M\times N}$, $h=\{\mathrm{1,2},\dots ,H_t\}$. As described in Sec. 2.1, the PCANet will take all B-scans of all training OCT volumes together as input. Then, we decompose the B-scans into patches and use the patches to train $L_1$ PCA filters $W_l^1$, $l=\mathrm{1,2},\dots ,L_1$ and $L_2$ PCA filters $W_l^2$, $l=\mathrm{1,2},\dots ,L_2$, in the first and second stages, respectively. Finally, after the binary quantization and calculating the block-wise histogram in the output stage, for each B-scan $X_h^{\text{train}}$, one feature vector $f_h^{\text{train}}$ is obtained. For each volume $V_t^{\text{train}}$, its extracted features can be constructed as a matrix $F_t^{\text{train}}=[f_1^{\text{train}},f_2^{\text{train}},\dots ,f_{H_t}^{\text{train}}]$.

In the testing phase, given an input 3-D OCT volume $V_t^{\text{test}}\in {\mathbf{R}}^{M\times N\times H_t}$, we apply the PCA filters $W_l^1$, $l=\mathrm{1,2},\dots ,L_1$ and $W_l^2$, $l=\mathrm{1,2},\dots ,L_2$ from the training phase on each B-scan of this OCT volume. Then, after the binary quantization and computation of the block-wise histogram, we obtain one feature matrix $F_t^{\text{test}}=[f_1^{\text{test}},f_2^{\text{test}},\dots ,f_{H_t}^{\text{test}}]$.

Figure 4 shows two examples of the extracted feature maps from different layers on the AMD and ME OCT images, respectively. As can be observed, for the input AMD image, the drusen regions will have high responses in the different feature maps of different layers [see the zoomed rectangle regions in Fig. 4(a)], while the edema is also very prominent in feature maps extracted from the ME image [see the zoomed rectangle regions in Fig. 4(b)]. Therefore, the learned PCA filters tend to capture meaningful pathology structure information, and the extracted feature maps can be used to achieve an effective classification, even without quantifying the size of the lesions.

Fig. 4

Two examples of feature maps extracted by PCANet from OCT images on subjects (a) AMD and (b) ME.

3.2.

Classification with Composite Kernel

As described above, given an OCT volume, the PCANet can extract one feature vector for each single B-scan in both the training and testing phases. Since the nearby B-scans of the 3-D OCT images are very similar, high correlations should also exist in their extracted features, which can be utilized to enhance the classification.³⁹ One possible way for exploiting the correlations among features within one 3-D OCT volume is to fuse them together. As described in Ref. 40, using the kernel to map the features into the higher dimensional feature space and then combining them as a composite kernel can better utilize the correlations and differences among different features. In addition, the parameters of different features are integrated into the composite kernel function and jointly optimized during the training process,³⁷ which can also effectively combine the information of different features for classification. Specifically, for each feature $f_h$ of a 3-D OCT volume, we use the linear kernel function in Eq. (7) to compute the corresponding kernel $K_h$. As introduced in Refs. 41 and 42, directly stacking these kernels or a linear combination of them can exploit the correlations among them. Since stacking these kernels will greatly increase the dimension and computational cost, we create the composite kernel by a linearly weighted combination of these kernels

Note that, since PCANet extracts many feature vectors (e.g., more than 30 in our test) from each volume, a number of corresponding weights ${\mu }_h$, $h=\mathrm{1,2},\dots ,H$ in Eq. (9) are required to be set, and searching the optimal values for so many weights would be very hard. In addition, different OCT volumes may consist of different numbers $H_t$ of B-scans; thus, the number of kernel weights for different volumes is also varied. To address this issue, before the construction of multiple kernels, we employ the PCA transform⁴³ to reduce the number of features and create the same number of kernel weights for each volume. Specifically, the PCA transform maps the feature matrix $F_t=[f_{1,},f_{2,},\dots ,f_{H_t}]$ to a new set of principal components (PCs): $F_t^{\mathrm{PCA}}=[f_{\mathrm{PC}1},f_{\mathrm{PC}2},f_{\mathrm{PC}3}]$ using the SVD technique on the covariance matrix of $F_t$, where several eigenvectors are calculated as the optimal projection axes. In this paper, we only retain the first three PCs, which account for the most information in the feature matrix $F_t$. After such a feature selection strategy, only three kernels (see Fig. 3) are generated, which are then fused as a composite kernel.

4.1.

Datasets

To evaluate the effectiveness of the proposed PCANet-CK method, we tested it on two real OCT datasets acquired from Duke University and First Affiliated Hospital of Hunan University of Chinese Medicine (HUCM), respectively. The Duke dataset was acquired using the spectral domain (SD)-OCT imaging system (Heidelberg Engineering Inc., Heidelberg, Germany) by Duke University, Harvard University, and the University of Michigan, which was publicly available in Ref. 44. The Duke dataset consists of 45 SD-OCT volumes from 45 patients (15 AMD, 15 DME, and 15 normal), and each volume contains many B-scans (ranging from 31 to 97). The original axial and lateral resolutions of the B-scans are $3.87\mu \mathrm{m}/\text{pixel}$ and 6 to $12\mu \mathrm{m}/\text{pixel}$, respectively. More details of the scanning protocols for this dataset can be found in Ref. 23. Three examples from NM, ME, and AMD subjects of the Duke dataset are shown in Fig. 5(a).

Fig. 5

Examples of SD-OCT images from (a) Duke dataset and (b) HUCM dataset.

The HUCM dataset was captured using the spectral SD-OCT imaging system (Heidelberg Engineering, Heidelberg, Germany) at the First Affiliated Hospital HUCM. The HUCM dataset is composed of 54 SD-OCT volumes (18 AMD, 18 DME, and 18 normal) from 48 patients. All volumes are of size $768\times 496\times 31\text{voxels}$, covering $8.8\times 8.8\times 2.0{\mathrm{mm}}^3$. Figure 5(b) shows three B-scans from NM, ME, and AMD subjects of the HUCM dataset, respectively. Note that only small pathological structures (e.g., drusen) exist in some B-scans of HUCM dataset, which are very challenging for recognition. The above two datasets used in our study were approved by the local Investigational Review Board and were performed in accordance with the tenets set forth in the Declaration of Helsinki. Written informed consent was obtained before enrolling patients in EUGENDA.

4.2.

Experimental Setting

In the proposed PCANet method, the main parameters include the PCA filter size, the block size for local histograms, and the filter numbers $L_1$ and $L_2$. In our experiment, we set the PCA filter size $n_1{=n}_2=11$, the block size $n_{B1}{=n}_{B2}=11$, and the block overlap ratio as 0.1. The effect of the PCA filter size and block size was analyzed in Sec. 4.4. As described in Ref. 33, $L_1$ and $L_2$ were set to be 8. In addition, for the composite kernel creation, the kernel weights ${\mu }_1$, ${\mu }_2$, and ${\mu }_3$ for the three kernels were set as 0.85, 0.10, and 0.05, respectively. The regularization parameter $C$ in Eq. (5) for the ELM classifier training was set to be 5.

In our experiments, we utilized the cross validation to evaluate the classification performance of the proposed PCANet-CK method. The cross validations were repeated with different random seeds to avoid the dataset splitting bias. The accuracy, specificity, and sensitivity were adopted to evaluate the classification performances, which were defined in a binary classification problem:

In our three-class classification problem, the sensitivity of a class label is also the prediction accuracy, and the specificity is defined in the same way for each class label, where the negative samples are the samples not in the considered class. Therefore, the overall sensitivity (Ov-Se), overall specificity (Ov-Sp), and overall accuracy (Ov-Acc) are averaged over the three-class labels. If the number of samples in each class is equal, the Ov-Se is also defined as the ratio: number of correctly predicted samples/total number of samples. In addition, the mean and standard deviation values of the above metrics are calculated.

4.3.

Results Comparisons

The performance of the proposed PCANet-CK method was first compared with the well-known OCT classification method: HOG-SVM.²³ The 15-fold (Duke dataset) and 18-fold (HUCM dataset) cross validations was repeated 10 times with different random seeds. The HOG-SVM method utilizes the HOG descriptor to extract the feature vectors and trains three binary SVMs⁴⁵ for the classification. The leave-three-out cross-validation strategy was applied on the two test datasets; this was achieved by choosing three volumes once from each class as the test dataset and the remaining 42 (in Duke dataset) or 51 (in HUCM dataset) volumes as the training dataset. For different experiments, three different volumes were chosen as test volumes, and 15 or 18 experiments were conducted to cover all the volumes in the two datasets.

The quantitative results on the Duke and HUCM datasets are tabulated in Tables 1 and 2, respectively. As can be observed, in the two test datasets, the proposed PCANet-CK method consistently performs better than the HOG-SVM method in terms of all quantitative metrics. Specifically, in the Duke dataset, the PCANet-CK can accurately classify all the test volumes of the three classes, whereas the Ov-Acc for the HOG-SVM method is about 96.6%. In the HUCM dataset, the Ov-Acc for the PCANet-CK method is 96.9%, whereas the Ov-Acc for the HOG-SVM method is about 90.7%. Moreover, in the classes of AMD and NM on HUCM dataset, the gain of the mean sensitivity of the proposed method over the HOG-SVM method is more than 9%, which demonstrates the effectiveness of the PCANet feature extraction and the composite kernel for exploiting the 3-D information.

Table 1

Classification results (%) on Duke dataset.

Table 2

Classification results (%) on HUCM dataset.

Very recently, the deep CNN method was also tested on the Duke dataset.²⁷ The deep CNN method first uses a pretrained CNN model (GoogleNet) and then fine-tunes it on the Duke dataset to extract the features for the identification of AMD, DME, and NM. In the experimental setting of Ref. 27, the cross validation was also utilized by dividing the whole dataset into 15 folds, with each fold containing three volumes (one from each class). However, different from the above validation, each experiment involved eight folds (24 volumes) for training and seven folds (21 volumes) for testing. Folds were sequentially rather than randomly chosen. Here, the proposed PCANet-CK and HOG-SVM methods were tested under the same experimental setting as in Ref. 27. Since only the mean sensitivity was reported in Ref. 27, the value of the deep CNN method was compared with those of the PCANet-CK and HOG-SVM methods, as reported in Table 3. As can be seen, the proposed PCANet-CK method generally performs much better than the HOG-SVM and deep CNN methods for classifying the images of AMD and ME subjects. In addition, for the classification of images from NM subject, the PCANet-CK delivers better performance than the HOG-SVM, while is very close to the deep CNN. These results show the superiority of the feature learning strategy used in PCANet-CK and deep CNN methods over the hand-crafted HOG feature extraction strategy adopted in the HOG-SVM method.

Table 3

Classification performances (sensitivities in %) of different methods on Duke dataset.

In our experiments, the proposed PCANet-CK method was implemented on a desktop computer with an Intel (R) Core i7-6700K CPU and 64 GB of RAM under the environment of MATLAB R2016b. The average time for testing one volume requires about 7.2 and 3.1 s on the Duke and HUCM datasets, respectively. Note that, since the training phase of the PCANet-CK method can be an offline process, it does not need to be considered in the testing phase. In addition, our code is not optimized for speed. The processing time is expected to be reduced significantly by more efficient coding coupled with a general purpose graphics processing unit.

4.4.

Effect of the PCA Filter Size and Block Size

In this section, the effect of the PCA filter size and the block size on the proposed PCANet model was analyzed on the Duke and HUCM datasets. The leave-three-out cross validation and the Ov-Se were also used here. For the evaluation of the PCA filter size, we varied the PCA filter size ($n_1{=n}_2$) in the first two convolution stages from 3 to 19 and kept other parameters (e.g., block size, overlap ratio, and kernel weights) the same as that in Sec. 4.2. The results with different PCA filter sizes are shown in Fig. 6(a). As can be observed, the performance of the proposed PCANet-CK method will generally improve as the PCA filter size increases from 3 to 11. When the PCA filter size further increases, our performance will become stable or even worse. Since utilizing a larger PCA filter will create higher computational cost, the PCA filter size is set to 11. For the evaluation of block size, we varied the block size ($n_{B1}{=n}_{B2}$) in the output stage in a range from 3 to 19. Figure 6(b) shows the classification results with different block sizes on two test datasets. As can be seen, our PCANet-CK method is comparatively stable and achieves the best performance with the block sizes from of $3\times 3$, $5\times 5$, and $11\times 11$ on two datasets. In this paper, we set the block size to be $11\times 11$.

Fig. 6

Overall sensitivities of the PCANet-CK method on Duke and HUCM datasets for different (a) PCA filter sizes in the first two stages and (b) block sizes in the output stage.

4.5.

Comparisons of the Composite Kernel with Single Kernel

We also conducted additional experiments using a single-kernel ELM classifier to validate the superiority of the composite kernel classifier. Specifically, after obtaining three PCs for one volume, we use them to separately create three different kernels. However, instead of fusing these kernels together, only one of them is utilized to train the ELM classifier and test one volume. The method using the first kernel is called the PCANet-${\mathrm{K}}_1$. The method using the second kernel is denoted the PCANet-${\mathrm{K}}_2$. The method using the third kernel is termed the PCANet-${\mathrm{K}}_3$. The results for the Duke and HUCM datasets are tabulated in Tables 4 and 5. As can be observed, in both the Duke and HUCM datasets, the composite kernel method consistently performs better than the single kernel method in terms of all quantitative metrics, which demonstrates the effectiveness of the composite kernel for exploiting the correlations among adjacent B-scans within one OCT volume. The results also show that the composite kernel ELM can utilize complementary information of different single kernels to further enhance classification.

Table 4

Classification results (%) on Duke dataset using single kernel.

Table 5

Classification results (%) on HUCM dataset using single kernel.