基于旋转不变非线性局部模糊编码的膝骨关节炎辅助诊断

任会峰，伏建雄，鄢锋，李富

（1.无锡学院物联网工程学院，江苏无锡 214105；2.湖南中自工业智能研究院有限公司，湖南长沙 410000；3.长沙有色冶金设计研究院有限公司，湖南长沙 410011）

0 引言

膝骨关节炎（Knee Osteoarthritis，KOA）是世界上致残率第四高的疾病，其病因和发病机制尚不明确，性别、职业、年龄、遗传、运动过度、体重超重等多种因素都可以诱发或加重KOA的进展。控制膝骨关节炎的最佳策略是早发现、早诊断、早治疗。研究表明，磁共振成像（Magnetic Resonance Imaging，MRI）对众多疾病的诊断表现出良好的敏感性、特异性和可重复性。毛家骥等依据MRI加权像分析手腕部关节炎表现，确定有关病证分布及发生规律；王紫仪等探讨磁共振成像加权像高信号比值的病因推测价值及其与听力预后的关系，证明信号强度越高，提示出血可能性越大，听力预后越差。尤其在器官磁共振评分系统中，通过定量测量T2值，可以量化结构成分和生物力学指标变化，相较于关节镜等检查方法，磁共振检查方式准确率更高、漏诊率更低，组间差异显著。因此，T2加权像表现出与膝骨关节炎较强的相关性。依据磁共振加权像推测膝骨关节炎病况的传统方式主要依赖于医生的主观经验，对医生接诊的病例类型数目、积累的医疗诊断经验等要求较高，受专业水平、情绪波动等多种因素影响难免存在误判，容易导致医疗纠纷，而且诊断效率降低。智慧医疗技术的发展为众多疾病提供了客观的诊断依据和高效的治疗手段。相关领域广泛采用图像处理技术辅助疾病诊断，受疾病医学图像样本数量少，尤其是难以获得数据完备的样本，以及图像标注缺乏领域知识等因素影响，一些基于大量数据样本的图像处理技术受到限制。

纹理是一种无处不在的视觉特征，已广泛应用于多种疾病的辅助诊断。刘欣等应用Logistic回归分析行多特征联合回归分析相关纹理特征ROC曲线的AUC，论证了MR T2WI图像纹理对于鉴别类风湿关节炎和痛风性关节炎具有一定价值。研究表明，纹理分析能充分挖掘医学影像中肉眼无法观察到的隐藏信息，能反映组织的潜在病理及生理学特征，在鉴别诊断方面潜力巨大。基于图像纹理的计算机辅助诊断技术用于膝骨关节炎检测，可以减少对诊断的主观性影响，降低对操作者熟练程度的依赖性，避免因疲劳等因素造成误诊。纹理特征提取为纹理分析提供数据基础，从实用角度讲，要求纹理特征数量少、区分度高、鲁棒性好且易于理解。目前主流的纹理提取方法有基于统计的方法、基于建模的方法和基于频域变换的方法。然而，由于运动模糊、斑点噪声等伪影造成膝关节T2加权像的不确定性或不完备性，目前流行的灰度共生矩阵、分形计算、Gabor滤波等方法性能不佳。而纹理直方图计算简单、鲁棒性强、易于理解，是描述局部区域纹理的有效方法。汤卫霞等利用MR图像的纹理直方图特征揭示肿瘤异质性，表明MR图像纹理直方图在肿瘤的早期诊断、精准治疗和预后评估等方面具有巨大潜力。但是，纹理直方图的传统局部硬编码基于邻域像素和种子像素的二值比较提取纹理特征，其前提是不同纹理单元的纹理单元向量不同，这一前提假设与实际不符，而且局部硬编码忽略了像素分布的不确定性和灰度差导致纹理单元数量减少随之图像区分能力下降。在此基础上，Barcelo等提出采用局部模糊编码研究骨髓显微图像中白细胞和红细胞图像变化和内在空间结构敏感性，其采用线性隶属度函数，提高了图像分类准确率，但是利用线性隶属度函数拟合像素灰度差的包络线存在较大偏差，而且得到6 561维的高维特征向量描述局部模糊编码，特征维数偏高、冗余噪声较大且对旋转图像分类效果不佳。

本文在分析传统局部硬编码的基础上，采用指数函数代替线性隶属度函数拟合邻域像素与中间像素灰度差的包络线，通过局部模糊编码频率直方图统计描述局部单元的灰度变化，借鉴串移位操作将旋转后局部模糊编码相同的局部单元归为一类，构建一种旋转不变的非线性局部模糊编码纹理描述方法，从众多局部单元中选取最具代表性的局部单位向量描述图像纹理，从而大幅降低纹理特征的维数。在此基础上，通过构建一个局部模糊编码描述符进一步压缩纹理特征的维数，并用于纹理图像分类和膝骨关节炎辅助诊断。

1 传统局部硬编码

传统局部硬编码（Local Grid Coding，LGC）是基于图像块邻域像素灰度值大小关系的纹理描述算法，其基本思想是：任意图像

I

均由多个3×3共9个像素大小的局部单元组成，将局部单元中间像素

V

灰度值记为

V

，周围的8个领域像素

V

到

V

灰度值记为

V

，

i

=1…8。假设

V

为局部单元左上角像素，从

V

开始按照顺时针方向将每个邻域像素与中间像素的灰度值进行比较，得到一个8维的局部单元向量（Local Unit Vector，LUV），

LUV

=（

E

，…，

E

），其中：

根据式（1）可知，局部单元可能有3～3种

LUV

，不同局部单元的

LUV

数用局部单元值

N

区分，其中：

所有局部单元值的频率直方图即为传统局部硬编码。可见，局部硬编码需要高达6 561维的单元向量描述，而且，从局部硬编码的计算过程可知，该算法忽略了像素分布的不确定性和像素灰度值差异程度。为了描述邻域像素灰度值差异程度，文献［24］提出线性局部模糊编码（Linear Local Fuzzy Coding，LLFC），但该算法采用的模糊隶属度函数为线性，利用线性隶属度函数拟合邻域像素灰度差异大小的包络线存在较大偏差，因此本文提出非线性局部模糊编码。

2 非线性局部模糊编码

为便于描述，将局部硬编码中

E

隶属于0、1、2的第

i

个邻域像素的隶属度值记为

μ

（

d

）、

μ

（

d

）和

μ

（

d

），构成局部模糊单元向量

FUV

。

由式（3）可知，邻域像素

V

和中间像素

V

之间的灰度差

d

越小，邻域像素和中间像素灰度值越相似，反之亦然，则图像

I

邻域像素与中间像素灰度差直方图的包络线如图1中实线所示，采用指数函数对该包络线进行拟合，图1中虚线为指数拟合曲线。

Fig.1 Envelope and fitting curve图1 包络线与拟合曲线

从图1可知，|

d

|＞40的概率小于0.1，故可设置非线性模糊隶属度函数如下：

由式（1）-式（6）所示，对任意图像

I

，从图像

I

的第2行第2列像素作为中间像素开始到倒数第2行第2列像素为止，以3*3的滑动窗口扫描整幅图像以计算所有局部单元值出现的频率，图像

I

的非线性局部模糊编码就是所有局部单元值的频率直方图。可见，非线性局部模糊编码很好地消除了局部硬编码采用绝对差值的缺陷，但仍然需要很高的维数向量以描述局部模糊编码，最高可达6 561维，而且该算法依然没有考虑图像角度旋转问题。文献［25］指出，对称不变的纹理模式应该视为同一纹理模式，对一幅图像旋转一定角度后，图像的纹理并没有发生变化，但是按照非线性局部模糊编码的计算过程却会出现完全不一样的非线性局部模糊编码，因此传统局部硬编码以及非线性局部模糊编码均不具备旋转不变性。此外，高达6 561维的向量不仅导致较大的计算量，而且会为图像的解读增加干扰噪声，影响机器的判断。为此，对非线性局部模糊编码加以改进，得到旋转不变非线性局部模糊编码（Rotation Invariant Nonlinear Local Fuzzy Coding，RLFC）。

3 旋转不变非线性局部模糊编码

图2为任意图像中都可能出现的12种局部硬编码的局部单元，大致分为图2（a）和图2（b）两大类。借鉴字符串中的移位操作可知，在顺时针右移位之后，图2（a）中的局部单元向量可以转换为下一个局部单元向量。在图2（b）中也有相同的结构。

Fig.2 Examples of local unit图2 局部单元示例

在此，将右移后局部单元向量相同的局部单元归为一类，则局部单元维数从6 561维缩减到834维，大幅削减描述图像纹理的特征维数，而且具备旋转不变性，描述图像纹理的能力增强，但描述图像纹理的向量维数依然较高。

因此，

f/f

比值越小，图像越粗糙。当

f

＜0.01*

m*n

时，图像极度粗糙。

f+f

值越大，图像越粗糙；

f/f

比值越小，图像越平滑，当

f

＜0.01m*n，图像极度平滑。

f+f

值越大，图像越平滑，其中，

m*n

为图像的分辨率。记

f

=

f

+

f

，

f

=

f

+

f

，可知，对于任意图像，若

f

越大，

f

越小，则图像越粗糙；

f

越大，

f

越小，图像越平滑。从而，可用一个2维向量

X=

（

f

，

f

）表示局部模糊编码直方图的复杂度，大幅降低了描述图像纹理的特征维数，且削减了干扰信息。

4 实验结果及分析

4.1 标准图像集测试

为了评价本文所提方法的有效性，将所提算法RLFC与传统局部硬编码LGC、线性局部模糊编码LLFC、非线性局部模糊编码NLFC进行比较，从Kylberg纹理图像库和Brodatz纹理图像库中抽取测试图像集，这两个图像库常用于评价和比较纹理度量。其中，Kylberg图像库包括28个类，每个类有160个图像，每幅图*像的大小为576×576像素，图像数据集可从http：//www.cb.uu.se/～gustaf/texture下载，从每个类中随机抽取一幅图像均分为4个大小为288×288像素的子图像，将每一幅子图像以10°为间隔，旋转10°到360°，从旋转后的子图像中心提取大小为128×128像素的孙图像，共得到28×4×36=4 036幅图像，作为本次测试的Kylberg图像集；Brodatz图像库可从http：//www.ee.oulu.fi/mvg/page下载，包括112类纹理图像，从每一类中随机抽取一张纹理图像以10°为间隔，旋转10°到360°，从旋转后的纹理图像中心提取大小为128×128像素的子图像，共得到112×36=4 032幅图像，作为本次测试的Brodatz图像集。采用式（8）基于

L

-距离的

K

-means聚类：

其中，

X

、

X

是两幅图像的局部模糊编码描述向量，

x

、

x

是

X

和

X

的分量。将聚类结果与原始标签一一对比，正确归类的样本数与总样本数之比定义为精度，以此作为评价指标。随机初始化聚类中心，并分别执行10次聚类计算。表1为在构建的Kylberg图像集和Brodatz图像集上分别采用传统局部硬编码LGC、线性局部模糊编码LLFC、非线性局部模糊编码NLFC和本文提出的旋转不变非线性局部模糊编码RLFC聚类精度对比。

从表1可以看出，局部硬编码可以在一定程度上反映邻域像素和种子像素的关系，但灰度值的硬划分导致聚类精度较低。此外，对于不同的初始聚类中心，该算法表现出较弱的鲁棒性。然而，无论是否对纹理单元向量进行三进制编码，局部模糊编码、旋转不变非线性局部模糊编码都能较大程度地提高准确率，在经过旋转预处理所构建的图像数据集上，本文提出的方法相比局部模糊编码表现出更高的准确率。10次测试中，最大准确率分别提高9%和1%，最小准确率分别提高6%和3%，平均准确率分别提高14%和7%。这是因为本文提出的方法虽然减少了纹理单元数，但通过局部单元平移，提高了图像的旋转不变性和对不同初始聚类中心的鲁棒性。

Table 1 Comparison of Clustering accuracy of different algorithm on texture ionage dataset表1 不同算法在纹理图像集上的聚类精度对比

为了进一步体现所提方法的优越性，在上一步构建的Kylberg和Brodatz图像集上添加均值为0、信噪比（Signal to Noise Ratio，SNR）由20变化到10的高斯白噪声，以评价各算法抗噪声干扰能力，依然采用式（8）的K-means聚类，以准确率作为评价指标，图3（a）和图3（b）分别为各算法在增加白噪声后的Kylberg和Brodatz图像集的实验结果。

从图3可知，随着SNR的减小，噪声水平逐渐增加，各算法在两个图像集上的变化趋势相同，聚类的准确率都呈下降趋势。其中，传统局部硬编码LGC的聚类精度下降最快，表明该方法的抗噪声能力最差；线性局部模糊编码和非线性局部模糊编码在两个数据集上随着噪声水平的增加，聚类精度略有下降，但总体变化较小，所提方法RLFC在不同高斯白噪声影响下的性能均优于其它方法，而且在增加高斯白噪声的两个纹理图像集上表现都比较稳定。综合表1和图3可见，所提方法无论在聚类精度还是在抗噪声干扰能力上都较好。

Fig.3 Experimental results after adding Gaussian white noise图3 增加高斯白噪声后的实验结果

4.2 膝骨关节炎辅助诊断

图4依次为KOA逐渐严重的4幅T2加权像，与图像对应的病症都具有明确的病理诊断结果，而且相关患者在MRI检查前未接受穿刺取样及手术治疗，将图像进行灰度均衡化处理，以降低不同设备不同参数的图像对比度和亮度对图像灰度值的影响，并由具有多年MRI诊断经验的医师进行阅片分析及标注KOA分期，对应KOA的初期、中期、后期和晚期，4幅图像中存在不同程度的关节软骨损伤，但对病变严重程度肉眼难以区分。

Fig.4 Images of KOA图4 KOA样本图像

图5为基于本文提出方法的4幅KOA图像的局部模糊编码直方图。从图5可以看出，不同的KOA图像在局部模糊编码描述符上有比较明显的差异。为了进一步直观展示不同样本图像之间的纹理差异，按照式（7）对各图像分别计算4个典型局部单元类及出现的频率。

Fig.5 Local fuzzy coding histogram of T2 weighted image图5 T2加权像的局部模糊编码直方图

表2是各幅图像局部模糊编码直方图复杂度的二维向量（

f

，

f

）对照表。从表2可见，KOA1的

f

最大，

f

最小，表明KOA1对应的图像纹理最平滑，这是因为初期血管破坏较少，软骨纤维化较轻。从初期到晚期，

f

逐渐减小，

f

逐渐增大，KOA4对应的图像纹理最粗糙，这是因为晚期滑膜破坏超过全面积的75%，软骨纤维化严重，甚至形成骨质增生骨赘，血管破坏严重到血管闭塞。

Table2 RLFC descriptor of KOA images表2 KOA样本图像的RLFC描述符

图6为

f

/

f

比值和KOA严重程度之间的关系曲线。由图6可知，从初期到晚期，

f

/

f

比值逐渐减小，进一步说明随着膝骨关节炎的不断恶化，T2加权像的纹理越来越粗糙。从曲线的斜率变化趋势看，KOA发展的4个阶段中，从初期到中期变化的过程相对缓慢，这是由于KOA属于多因素导致的相应部位病理产物对人体刺激的综合结果，这与初期人体生理组织相对较好，发生病变的组织相对较少有关，一般由于膝关节后侧半膜肌腱的薄弱部位发生病变有关；从中期到后期恶化最为迅速，这是因为中期伴随外侧半月板前外缘游离，容易在扭转时发生横形破裂等不可逆的继发损伤，而这种损伤会逐渐加剧；后期一般发生于前侧膑腱劳损，从后期到晚期发展也比较迅速，在后期发生髌骨软化后膝关节内侧半月板由于失去保护而加剧磨损形成KOA晚期。

Fig.6 Relationship between f T/f H ratio and KOA severity图6 f T/f H比值与KOA严重程度的关系曲线

由表2和图6可知，

f

越小，则

f

越大。KOA图像的平滑度随

f

/

f

的减小而增加，对应关节炎更为严重，测试结果与实际情况相吻合，表明所提出的方法对膝骨关节炎的分期和诊断有较好的辅助作用。

5 结语

采用MRI的纹理分析有助于术前评估膝骨关节炎的病理分期，是一种KOA辅助诊断的有益探索。局部模糊编码摒弃了传统纹理谱中的刚性比较，能够更准确地描述图像的局部特征。通过引入移位操作和三进制编码，降低了纹理向量的维数，减少了计算量，且具有更好的旋转不变性和抗噪性，在标准图像集和膝骨关节炎T2加权像的实验结果表明，提出的旋转不变非线性局部模糊编码的纹理分析方法是一种无创性KOA辅助分期手段，可以为KOA早期发现提供快速的客观依据，有利于早期指导患者生活避免病程快速发展。

本研究为单中心的回顾性研究，虽然对获取的实验图像进行了预处理等，一定程度降低了图像数据的差异性，但仅仅基于MRI图像纹理分析，并未结合临床特征，如性别、年龄等，有待后续对影像数据进一步完善后结合更深入的影像组学方法与大数据关联挖掘进行综合分析。作为一个相对新兴的领域，基于纹理分析的KOA辅助诊断处于起步阶段，仍面临很多挑战，需加强实验和数据处理技术以增加实验样本数，研究小样本机器学习的辅助诊断技术和采用半自动或全自动方法实现病变位置的自动定位和病理解释将是下一步工作的重点。