光学相干断层扫描技术(OCT)是一种非接触高分辨率层析成像技术，其利用近红外光线及光学干涉原理，通过对活体组织扫描可得到微米分辨率下的生物组织断层影像，已广泛用于眼科医疗成像的应用^[1-4]。OCT的微米级分辨率使多个视网膜细胞层可视化，可观察到细微而有价值的病理结构，诸如老年性黄斑变性(AMD)和糖尿病性黄斑水肿(DME)等许多黄斑和眼部疾病可以在早期阶段确定。在临床诊断中，对于由一系列2维OCT B-scan构成的3维OCT数据即OCT体数据，医生需要手动识别OCT体数据中每个B-scan的各种黄斑病变然后确定疾病类型。这种人工手动分析费时且对医生的专业性要求很高。因此，迫切需要开发一种有效的计算机辅助OCT图像分析技术。

近年来，视网膜OCT图像分类获得越来越多的研究。2011年，Liu等人^[5]利用局部二值模式(LBP)直方图提取OCT视网膜图像的局部特征，然后通过主成分降维和多尺度空间金字塔(SP)形成图像的全局特征，最后送入分类器对OCT视网膜图像进行分类。2014年，Srinivasan等人^[6]利用多尺度方向梯度直方图(HOG)提取OCT图像的特征并送入分类器完成分类。这些方法均取得不错的分类效果，但计算量大且耗时，仅适用于小数据集的分类，难以进行大规模的图像训练和分类任务。2017年，Sun等人^[7]利用稀疏表示的稀疏降维特点，先提取OCT视网膜图像中每片(Patch)的尺度平移不变性(SIFT)特征，再对这些SIFT特征稀疏编码，然后通过多尺度空间金字塔形成图像的全局特征，最后送入线性SVM分类器中完成分类，降低了分类计算量且提高了分类效率。这些方法的共同点是所采用的LBP、HOG和SIFT特征都是人工设计特征，而临床上获得的OCT图像通常包含非常复杂的病理结构，因此直接从OCT图像中学习特征更为可取，例如最近大受关注的卷积神经网络(CNN)可以学习图像更深层的特征，对类内差异的波动有很好的鲁棒性，克服了人工特征的局限性。Karri等人^[8]就用卷积神经网络(CNN)^[9-11]来学习OCT视网膜图像特征，分类结果很好。但一般来说，CNN有许多网络层，通常包含至少数百万层参数，因此需要大量的训练数据集和计算成本来训练网络。

针对CNN这一问题，Chan等人^[12]提出了一种结构简单、参数少、训练高效的PCANet算法。该网络利用一系列已学的PCA滤波器对图像多层卷积，然后进行二值化和块直方图最大池化处理，最后送入线性支撑向量机(SVM)分类器中完成图像分类。实验证明，PCANet在MNIST、FERET和CIFAR10等数据集上都有很好的图像分类和识别效果。2017年，Fang等人^[13]基于PCANet提出了PCANet-CK算法应用于3D视网膜OCT图像的分类。这种方法的主要优点是能够利用OCT体数据中相邻B-Scan图像之间的相关性。

在文献[12]中，Chan等人还提出了一种有类标签监督训练的LDANet用于图像分类，该网络是将PCANet中的PCA滤波器替换为一种通过线性判别分析学习的LDA滤波器，但实验结果显示其分类效果不如PCANet。本文深入分析了PCANet和LDANet的优点，提出了对PCANet获得的图像特征向量进行LDA后再利用SVM分类的算法即PCANet-LDA。与PCANet算法相比，该算法能够提高图像特征的类别区分性，从而提升图像的分类准确率。本文将PCANet-LDA算法用于OCT视网膜图像的分类，在两个数据集上与PCANet算法^[12]及先进的ScSPM算法^[7]进行对比实验，验证所提出算法的有效性和先进性。

PCANet-LDA分类算法是一种类标签监督的多层卷积网络，该网络主要分为3部分：1)OCT图像预处理；2)PCANet特征提取；3)LDA监督层。网络最后连接一个线性SVM分类器，框架如图 1所示。对于给定的$N$张带类标签(AMD、DME、NOR)的视网膜OCT图像训练集，先进行预处理获得感兴趣的视网膜区域，然后利用PCANet-LDA网络在PCA卷积层学习PCA滤波器，在LDA监督层利用提取的PCANet特征和类标签学习LDA矩阵，最后由LDA投影特征训练线性的SVM分类器。

图 1 PCANet-LDA算法流程图

Fig. 1 Flow chart of PCANet-LDA algorithm

OCT图像通常充满斑点，并且视网膜的位置在扫描中变化很大，这使得将所有视网膜区域对准到相对统一的位置并不容易。因此，需要对OCT图像进行预处理，提取视网膜OCT图像的感兴趣区域，提高图像分类准确率。本文采用Sun等人^[7]提出的自动对齐和剪切视网膜区域技术对视网膜OCT图像预处理，该方法分为感知、拟合、归一化3个阶段，具体描述如图 2所示。

图 2 图像预处理过程

Fig. 2 Image preprocessing((a) original image; (b) BM3D denoising; (c) binarizing; (d) median filtering; (e) morphological closing; (f) morphological opening; (g) polynomial fitting; (h) retina aligning)

在感知阶段，需要检测视网膜的整体形态，感知视网膜而不估计它们的视网膜色素上皮边界。具体方法是：1)利用BM3D去噪和阈值滤波降低原图像的噪声来感知视网膜的结构；2)通过中值滤波去除视网膜内部分离的黑点；3)用形态学闭运算和形态学开运算分别去除视网膜内的大黑斑和视网膜外因OCT扫描而存在的大噪声斑点。在拟合阶段，通过对数据点的自动选择与拟合来模拟去噪图像中视网膜的轮廓与形态：先从图像中提取两组数据点，然后自动选择其中一组数据点并用一阶或二阶多项式来拟合。在归一化阶段：将拟合后的视网膜区域对齐到一个相对统一的形态，并裁剪图像以修剪无关紧要的空间。

经过这3个阶段，可以得到预处理后的训练集 $\left\{ {{\mathit{\boldsymbol{I}}_i}} \right\}_{i = 1}^N$，${{\mathit{\boldsymbol{I}}_i}}$∈${{\bf{R}}^{{m_i} \times {n_i}}}$($m_i$和$n_i$为图像${{\mathit{\boldsymbol{I}}_i}}$的长宽尺寸)，接下来提取$\left\{ {{\mathit{\boldsymbol{I}}_i}} \right\}_{i = 1}^N$的PCANet特征。

PCANet网络^[12]是CNN的一种简单变形，一般采用两层PCA卷积层和一个非线性输出层来提取图像的特征，具体过程如下：

1) PCA卷积层第1层。将$N$幅预处理后的训练图像$\left\{ {{\mathit{\boldsymbol{I}}_i}} \right\}_{i = 1}^N$作为第1层卷积层的输入图像，对每幅都进行分块、去均值和矩阵化处理，那么每一幅图像都可以用一个矩阵表示；然后计算这$N$个矩阵的前$L_1$个主成分并作为该层的$L_1$个PCA滤波器 $\left\{ {{\mathit{\boldsymbol{W}}_l^1}} \right\}_{l = 1}^{{L_1}}$；最后，用这$L_1$个滤波器对$N$幅输入图像作卷积，得到$L_{1}N$个特征图。因此，一幅图像${{\mathit{\boldsymbol{I}}_i}}$通过第1层PCA卷积后可以得到$L_1$个特征图 $\left\{ {{\mathit{\boldsymbol{I}}_i} * \mathit{\boldsymbol{W}}_l^1} \right\}_{l = 1}^{{L_1}}$。

2) PCA卷积层第2层。将第1层卷积得到的$L_{1}N$幅特征图作为第2层卷积层的输入图像，与第1层类似，每一幅输入图像都可以用一个矩阵表示，然后计算这$L_{1}N$个矩阵的前$L_2$个主成分并作为第2层的PCA滤波器 $\left\{ {\mathit{\boldsymbol{W}}_l^2} \right\}_{l = 1}^{{L_2}}$；最后，用这$L_2$个滤波器对$L_{1}N$幅输入图像作卷积，得到$L_2$$L_{1}N$个特征图。因此，一幅图像${{\mathit{\boldsymbol{I}}_i}}$通过两层PCA卷积后可以得到$L_1$$L_2$个特征图。当然，此处可以根据应用需求增添PCA卷积层数，学习出更深层次的特征图。

3) 非线性输出层。非线性输出层将训练图像${{\mathit{\boldsymbol{I}}_i}}$的$L_1$$L_2$个特征图融合成一个特征。首先，将$L_1$$L_2$个特征图分成$L_1$组；然后采用二值哈希、加权求和、块直方图和空间金字塔级联(SPM)，将每组的$L_2$个特征图融合成一个特征$\mathit{\boldsymbol{Z}}_i^l \in {{\bf{R}}^{2^{{L_2}} \times 21}}$；最后将$L_1$个特征组合起来作为图像${{\mathit{\boldsymbol{I}}_i}}$的PCANet特征 $ {\mathit{\boldsymbol{f}}_i} = \left[{\mathit{\boldsymbol{Z}}_i^1, \mathit{\boldsymbol{Z}}_i^2, \cdots, \mathit{\boldsymbol{Z}}_i^{{L_1}}} \right] \in {{\bf{R}}^{2^{{L_2}} \times 21{L_1}}}$。

LDA监督层是根据线性判别分析思想，将带有AMD、DME、NOR类标签的PCANet特征 $\left\{ {{\mathit{\boldsymbol{f}}_i}} \right\}_{i = 1}^N $投影到一个低维空间，使投影后同类图像的特征相距更近，异类图像的特征相距更远，即投影后特征的总类间散布矩阵的迹与总类内散布矩阵的迹之比最大。

设投影矩阵为LDA矩阵 ${\mathit{\boldsymbol{Q}}^{\rm{T}}} \in {{\bf{R}}^{q \times {2^{{L_2}}}}}$，$q \in \left[{1, {2^{{L_2}}}} \right]$，经过投影${\mathit{\boldsymbol{g}}_i} = {\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{f}}_i} $后，$\left\{ {{\mathit{\boldsymbol{g}}_i}} \right\}_{i = 1}^N$的总类间散布矩阵和总类内散布矩阵分别为${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_b} = {\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_b}\mathit{\boldsymbol{Q}}$，${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_w} = {\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_w}\mathit{\boldsymbol{Q}}$，其中，${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_b}$和${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_w}$分别为$\left\{ {{\mathit{\boldsymbol{f}}_i}} \right\}_{i = 1}^N$的总类间散布矩阵和总类内散布矩阵。那么投影矩阵$\mathit{\boldsymbol{Q}}$满足^[12]

式中，${\mathit{\boldsymbol{E}}_q}$表示 $q \times q$的单位矩阵。注意，这里要求$\mathit{\boldsymbol{Q}}$是正交矩阵。上述问题的解即 $\mathit{\boldsymbol{ \boldsymbol{\varPhi} = \boldsymbol{\varSigma} }}_w^ + {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_b}$的前$q$个主成分特征向量组成的矩阵，其中 $q \in \left[{1, {2^{{L_2}}}} \right]$，表示伪逆，用来求解非全满秩矩阵的逆。

最后，将监督层输出的LDA特征 $\left\{ {{\mathit{\boldsymbol{g}}_i}} \right\}_{i = 1}^N$送入线性SVM分类器中训练分类器。

可以看出，PCANet-LDA算法的核心是对PCANet提取的图像特征利用LDA进行降维。由于PCANet能提取图像的纹理特征，线性判别分析能对特征降维的同时提高特征的类别区分性，因此，PCANet-LDA算法结合了PCANet和LDA两种方法的优点。

本节分别在医院临床数据集和杜克数据集上对PCANet-LDA算法进行实验，并与图像基线分类PCANet^[12]及当前先进的ScSPM算法^[7]进行比较。由于PCANet-CK^[13]处理的对象是OCT体数据，本文处理的是2维OCT图像数据，因此，本文的对比算法不包括PCANet-CK。所有实验程序均基于MATLAB R2015b版本编写，实验平台为一台装载Intel(R) Core(TM) i7-7700k CPU @4.5 GHz处理器和32 GB内存的组装台式机，系统版本为Windows 10 Professional (10586)，使用固态硬盘提升IO性能。其中，PCANet和ScSPM的源代码分别来自文献[12, 14]。

首先介绍实验的两个数据集:

1) DUKE数据集^[16]。用SD-OCT频域在杜克大学、哈佛大学和密歇根大学成像，且获得机构审查委员会批准的议定书。该数据集包含45个OCT体数据，其中AMD、DME和正常体数据均为15个，而每个体数据包含36~97幅OCT B-scans图像。本文将每类15个体数据中的所有B-scans混合在一起，得到AMD、DME、NOR 3类视网膜的OCT B-scans图像集，每类依次包含722、1 101、1 194幅OCT图像，从每类中取700幅图像用于分类实验。

2) 医院临床数据集。从临床获得的OCT视网膜数据集，使用Cirrus TM检查装置。该数据集包含AMD、DME、NOR 3类视网膜的OCT B-scans图像，每类均有560幅，所有的SD-OCT图像读取和评估均由医学专家完成。由于医院临床中获得的OCT图像噪声大，且存在不同病变程度，处理起来非常具有挑战性。例如图 3中第1行是从临床数据集中选取的3幅OCT图像，分别对应AMD、DME和NOR，由图可见，图像噪声很大，视网膜在整幅图像中的位置和姿态各不一样，且周围有大片对分类没有意义的非视网膜区域。图 3第2行是对这3幅OCT图像预处理后的结果，从图 3中可以看到，预处理过程将OCT图像中非视网膜区域裁剪掉，只留下视网膜这些有用的区域；同时视网膜被旋转到一个统一的水平状态，减少了方向不一致对分类的影响。

图 3 AMD、DME和NOR视网膜OCT图像预处理示例

Fig. 3 Three examples of OCT images preprocessing on different subjects((a) AMD; (b) DME; (c) NOR)

为了说明PCANet特征图提取的有效性，图 4展示了图 3中预处理后的AMD和DME OCT图像在不同层的特征图提取结果。其中，PCANet的两层滤波器大小均设置为$k_1$×$k_2$=15×15，两层滤波器个数均为8($L_1$=$L_2$=8)。从图 4(a)中可见，对于AMD图像，在不同层的不同特征图中，玻璃膜疣都具有很高的响应(如图 4(a)中方框内所示)；而在DME图像的各层各个特征图中(图 4(b))，水肿也比较突出。因此，PCANet训练出来的PCA滤波器倾向于捕获有意义的病理结构信息，有助于视网膜OCT图像的分类。

图 4 AMD和DME视网膜OCT图像的PCANet特征图提取示例

Fig. 4 Two examples of feature maps extracted by PCANet from OCT images on different subjects((a) AMD; (b) DME)

在PCANet中两层滤波器大小和滤波器个数参数如上选定情况下，设置线性输出层的直方图块大小为$b_1$×$b_2$=15×15，块重叠率$ratio$=0.5。考虑到LDA矩阵的主成分个数$q$对PCANet-LDA算法的分类效果有一定影响，因此用临床数据集，随机选定一半图像作为训练集，另一半图像作为测试集，固定以上参数，对$q$∈[1, 256]进行调优，实验结果如图 5所示。可以看出，当LDA矩阵的主成分个数$q$为33~35时，总体分类正确率最高；当$q$低于33或高于35时，总体分类正确率都较低，虽然略有起伏但都低于最高正确率，有逐渐降低的趋势，因此本文设置$q$=35。

图 5 LDA矩阵的主成分个数($q$)对PCANet-LDA在临床数据集上的分类正确率的影响

Fig. 5 The influence of principal component number ($q$) of LDA matrix on correct classification rate of PCANet-LDA on clinical dataset

分别在医院临床数据集和杜克数据集上进行分类实验比较PCANet、ScSPM和PCANet-LDA 3种算法的分类效果。

对于临床数据集，每次随机取一半的图像作为训练集，另一半作为测试集，分别计算PCANet、ScSPM和PCANet-LDA算法关于AMD、DME和NOR图像的分类正确率和总体分类正确率；将上述过程重复10次，计算每种算法10次分类结果的均值和标准差。实验中，ScSPM算法的图像尺寸设置为16×16像素，字典原子个数为512，字典迭代次数为5，最近邻个数为40。实验结果如表 1所示。其中，分类正确率的计算方法为

表 1 3种算法在临床数据集上的分类正确率对比
Table 1 Correct classification rate comparison of three different algorithms on the clinic dataset

下载CSV

/%
	PCANet	ScSPM	PCANet-LDA
AMD	91.86±1.98	97.00±0.93	96.57±1.51
DME	90.89±1.58	95.64±0.95	96.54±1.38
NOR	97.54±1.54	98.32±0.81	98.50±1.19
总体	93.43±0.82	96.99±0.37	97.20±0.99

式中，${\eta _i}$和$\eta $分别表示第$i$类分类正确率和总体分类正确率，$A_i$和$B_i$分别表示第$i$类测试集的数量以及该类测试集中分类正确的数量。第$i$类分类正确率是用来衡量该类中分类正确的百分比，总体分类正确率是用来衡量数据集(包含AMD、DME及NOR)中分类正确的百分比。

由表 1可知，在临床数据集上，PCANet-LDA算法的分类正确率高出PCANet算法4%，性能明显提高；稍高于ScSPM算法的分类正确率，达到了先进水平。

在杜克数据集上，采用上述类似的实验方法及评估标准，实验结果如表 2所示。其中，在实现细节上与上述实验稍有不同，这里每次只用数据集中20%的图像作为训练集，其余80%图像用于测试集。主要原因是，杜克数据集中的图像总体上比临床数据集中的图像简单，用较少数量的图像作为训练集能够更好地区分不同算法的分类性能。

表 2 3种算法在杜克数据集上的分类正确率对比
Table 2 Correct classification rate comparison of three different algorithms on DUKE dataset

下载CSV

/%
	PCANet	ScSPM	PCANet-LDA
AMD	98.29±0.53	99.38±0.57	99.82±0.21
DME	98.14±0.74	99.59±0.60	99.77±0.32
NOR	97.21±1.56	99.30±0.55	98.98±0.88
总体	97.88±0.55	99.42±0.25	99.52±0.30

由表 2可知，在杜克数据集上，3种算法的分类精度都高于表 1中的结果。这里3种算法的精度差距不大，但总体结论一致，即PCANet-LDA算法的分类正确率与PCANet算法相比显著提高，同时优于ScSPM算法的结果。

文献[12]中，LDANet是从图像中学习LDA矩阵，并将其每一列作为一个LDA滤波器对图像作卷积，如此进行两层卷积后再连接一个非线性输出层和一个线性SVM分类器。与LDANet不同，本文PCANet-LDA是从图像的PCANet特征中学习LDA矩阵，并用矩阵对特征作投影，再连接线性SVM分类器。

由上文已知, 由PCA滤波器卷积获得的特征图包含了OCT图像中重要的病理结构信息，再经过非线性输出层就是PCANet特征 ${\mathit{\boldsymbol{f}}_i}$。当两层PCA卷积层的滤波器个数均为8时，特征${\mathit{\boldsymbol{f}}_i}$的维度为43 008，这样高维度的特征不仅冗余度高不利于分类，而且直接送入线性SVM分类器中训练会很耗时。加入LDA监督层后，PCANet特征${\mathit{\boldsymbol{f}}_i}$可以从43 008维降至5 880维($q$=35时)，维度仅为原来的13.7%。由表 1和表 2可知，添加LDA监督层后，在医院临床数据集和杜克数据集上，OCT图像的分类正确率都显著提高。可见，图像特征通过LDA降维后，特征更有类别区分性。

ScSPM算法是从OCT图像的片中提取人工特征(SIFT特征)后，通过稀疏编码进一步提取图像的局部特征；PCANet是通过两个卷积层学习得到OCT图像的PCA特征。从实验结果可以看出，ScSPM提取的图像特征优于PCANet提取的图像特征。而PCANet-LDA算法的分类结果稍优于ScSPM算法，也进一步说明了LDA监督层对PCANet特征的降维有助于OCT图像的分类。

由于临床数据集比杜克数据集的图像更复杂，更具挑战性，因此实验结果中PCANet、ScSPM和PCANet-LDA 3种算法在杜克数据集上的分类效果都优于在临床数据集上的分类效果。

PCANet作为测试图像分类效果的基准(baseline)算法非常有效。本文在PCANet算法的基础上提出了一种新的PCANet-LDA算法并应用于视网膜OCT图像分类。该算法在PCANet提取图像特征的基础上，进一步通过LDA方法得到类别区分性更强的特征。在两个视网膜OCT图像数据集上的测试结果表明，与PCANet算法相比，PCANet-LDA算法的分类正确率显著提高，达到或优于当前具有代表性的ScSPM OCT图像分类算法。可见，PCANet-LDA可作为以后视网膜OCT图像分类的基准算法，来比较不同算法的有效性。未来的工作是在更大的OCT数据集以及其他典型的图像数据集上验证算法的有效性。