2018, Vol. 39
形状、纹理和颜色是机器视觉中三大底层特征，在诸多图像识别、理解等任务中得到了广泛的应用.但在医疗影像的应用中，大量图像由于获取方法、分辨率和色彩通道等方面的局限性无法对颜色与纹理特征进行有效地识别和分类.因此如何合理地提取形状特征并完成形状建模[1]以供后续算法进行识别与分类就变得尤为重要.
本文选择计算机断层扫描数据进行形状建模研究，对不同病理、类型、分期的肺结节的形状进行研究，并着重对其中恶性肿瘤进行了病灶影像形状建模的研究.关于形状描述和识别的方法可以被简单分为三类：基于区域的方法、基于边界的方法和基于骨架的方法.
其中，基于区域的方法[2]考虑形状边界内部的全部像素计算一组形状特征的描述符，通过该描述符可以较好地处理复杂形状图像.基于边界的方法[3]把形状边界看作一组有序的物理或数学意义的坐标点，用向量表示坐标点.通过对每个向量表示的形状进行统计分析，得到一个描述形状的统计形状模型.基于骨架的方法[4]是提取图像形状的骨架图，图像形状的识别是利用匹配骨架图来实现.
现有方法在不同场景下表现出多种不足，此外三类方法在物体发生非刚性形变时其形状描述符鲁棒性较差，对于恶性肺结节形状描述存在着较高的难度.针对上述问题，本文提出了一种利用复杂网络理论[5]和分形维数[6]的形状描述新方法.
1 相关理论基础 1.1 局部二进制模式(LBP)LBP(local binary patterns)是由Ojala等[7]首先提出的一种用来描述局部纹理特征的描述子，采用一种圆形邻域(如图 1所示)进行LBP计算，对于在多尺度上运用LBP，则可以利用大小不一的圆形邻域(如图 2所示)来实现.
|
图 1 LBP特征描述过程 Fig.1 LBP feature description process |
|
图 2 大小不同的圆形邻域结构 Fig.2 Circular neighborhood structures of different sizes (a)—(1, 8)；(b)—(2, 8); (c)—(2，16). |
如图 1所示，该圆形邻域的中心是像素xc，半径为r，P个邻域像素点等角间距均匀分布在圆周上，所有像素点的灰度值为
|
(1) |
定义在圆形邻域系统上的LBP模式计算为
|
(2) |
|
(3) |
其中:s(x)是一个标记函数; xr, p, n是某个像素点的像素值，圆形邻域的半径为r，P个邻域像素点等角间距均匀分布在圆周上，n是圆形邻域系统上的像素点的相对位置，没有位于图像像素中心的邻域采样点采取双线性差值方式获得.
上述LBPr, p描述子是灰度不变的，但不具有旋转不变性，旋转图像会得到不同的LBP描述子.为了使LBP描述子具有旋转不变性，Ojala等又提出了一种旋转不变LBPr, pri描述子，通过一定的规则旋转圆形邻域得到一系列原始定义的LBP值，则旋转不变
|
(4) |
LBPr, pri描述子具有旋转不变性，但是，将旋转不变描述子LBPr, pri推广到多尺度上时，只是单独地计算各个尺度上的LBPr, pri，而没有考虑各个尺度之间的关系.
针对上述问题，本文将采用文献[8]提出的pi-LBP，该方法考虑了不同尺度上像素点之间的相关性(如图 3所示)，将更好地描述形状纹理.具体定义为
|
(5) |
|
图 3 pi-LBP的例子 Fig.3 Example of pi-LBP |
其中:G0=(g0, 1, …, g0, k)是从中心像素xc出发的一条路径; g0, 1为中心像素点的像素值，顺时针旋转该路径可得到新的路径Gp=(gp, 1, …, gp, k)；f=(f(1), …, f(k))是满足条件
|
(6) |
在本文的应用中，计算每个形状轮廓上的像素点的LBP值时，pi-LBPriP, G0, f中的两个圆形邻域系统分别为(r, P)=(1, 8), (2, 8)，滤波器为f=(1, -2, 1).
1.2 盒覆盖算法可以就复杂性来描述复杂网络.文献[9]将分形对象的复杂性度量定义为分形维数, 它也可以被理解为这些对象中的一个自相似性的表征.用一个在度量空间中量化分形密度的非整数来表示分形维数.
计算复杂网络的分形维数的盒覆盖算法由Song等[10]在2007年首次给出.借助一个N×lBmax的二维矩阵cil就可以实现贪心着色盒覆盖算法，N是组成网络的总的节点个数，网络直径再加上1即为最大盒子尺寸lBmax.NB(lB)为盒子总数，其算法流程图如表 1所示.
|
表 1 盒覆盖算法流程图 Table 1 Flow chart of box covering algorithm |
盒覆盖算法需要整体覆盖网络，且其是等概率地覆盖网络中任意部分，即网络中的节点是无差别的.该算法借助贪心着色的算法思路很巧妙地解决了确定所需的最少盒子数问题，并且该算法很容易在计算机上实现，对于用不同尺寸的盒子覆盖网络，其所需的最少盒子数可以通过一个矩阵实现.因此，盒覆盖算法是计算复杂网络分形维数最常用的方法.
2 形状网络化建模 2.1 网络建模方法本文把形状轮廓看作是相连的像素点，因此，一个形状的轮廓信息可以用一个向量S=[s1, s2, …, sN]来描述，N为形状轮廓上的像素点的个数，si=(xi, yi)是离散数值，表示轮廓上点i的坐标.
为了把复杂网络理论应用到形状建模的问题上，用图来表示形状轮廓(图 4).建立一个图G=(V, E)，轮廓中的每个点s∈S对应图G中每个顶点v∈V.如果两个点si和sj之间的欧式距离小于阈值r则将它们连起来，权值为它们的LBP值之差的绝对值.因此，用N×N的加权矩阵W来表示网络：
|
(7) |
|
图 4 形状边界的图表示法 Fig.4 Graphic representation of shape boundary |
式中:LBP(si)为点si的LBP值；dist(si, sj)为si和sj之间的欧式距离; wij为加权矩阵W的第i行第j列的值.
归一化到[0, 1]区间：
|
(8) |
在这个阶段，此网络的拓扑结构特征不明显，不存在本文需要的特征，因此需要在此网络基础上进行动态演化以便得到具有本文感兴趣的属性的复杂网络.
2.2 网络动态演化为了得到一个新的网络G*=(V, E*)，可以为前一小节提出的形状表示法给定一个变换δTi.应用阈值Ti选择边的集合E*，E*中每条边都有一个小于或等于阈值Ti的值.将表示为A=δTi (W′)的操作应用到加权矩阵W′的每个元素生成一个无权矩阵A，即
|
(9) |
这个变换δTi允许在演化的中间过程研究网络的属性，即子网络的属性是随着边的最大权重的增加而产生的.通过使用不同的阈值可能实现一组更丰富的描述网络动态的测量值.所以，用各种δTi转换描述网络特征，其中Tmin是Ti的初始阈值，其值以定值Tinc递增直到达到最大阈值Tmax(图 5).
|
图 5 随着阈值Ti变化的网络的动态演化 Fig.5 Network dynamic evolution by a threshold Ti (a)—Ti=0.1; (b)—Ti=0.2; (c)—Ti=0.3. |
本文的复杂网络不在自然的欧式空间中，因此复杂网络中两个节点之间的距离不是欧式距离，两个节点之间的最短路径包含边的数目为它们之间的距离.因此，在本文中，用盒覆盖算法计算复杂网络的分形维数时，本文将上述提到的盒覆盖算法加以改进.首先，定义盒子尺寸LB为连接两个节点最短路径包含的边的数目，复杂网络被尺寸为LB的盒子覆盖，LB大于盒子中任意两个节点之间的距离，一组相互之间不重叠的盒子可以覆盖一个网络，网络中的每个节点只能属于一个盒子；其次，本文根据节点度和LBP值对节点进行排序，节点先按节点度递减排序，当节点度相等时，再按LBP值递减排序.
原始的盒覆盖算法比较依赖于网络节点的初始编号，而初始编号是随机的，后续给节点着色的过程中也存在随机性，往往需要进行多次计算进行拟合.而改进后的盒覆盖算法利用节点度和LBP值对节点进行排序，排序之后的随机性主要来源于度值和LBP值相同时节点的排序问题，这种随机性比原始的盒覆盖算法的随机性小得多，可忽略不计.因此改进后的算法不仅不需要进行多次计算，还使得到的结果更精确，又因为可以并行处理节点入盒，所以提高了算法的计算速度.
基于以上工作，复杂网络估计分形维数d的常用方法是基于覆盖复杂网络所需的最少盒子数NB(LB)和盒子边长LB之间存在幂函数关系：
|
(10) |
因此，复杂网络的盒维数计算公式为
|
(11) |
通过拟合lg (NB(LB))和lg(LB)直线的斜率可以得到分形维数d，分形维数d为斜率的相反数，如图 6所示.
|
图 6 分形维数d Fig.6 Fractal dimension d |
对于本文的应用，对通过阈值Ti演化得到的每个复杂子网络可得到其分形维数di.因此，对于每个阈值Ti都会有一个对应的分形维数di，最后通过如非线性回归计算可得到阈值Ti和分形维数di的曲线，用此曲线的系数来表示复杂网络代表的形状，如图 7所示.
|
图 7 阈值Ti与分形维数d之间的关系 Fig.7 Relationship between the threshold Ti and the fractal dimension d |
通过复杂网络的分形维数d可以研究形状的物理性质.根据所用的阈值，凸起或其他轮廓特征会使形状表现出不同的复杂性，从而得到不同的分形维数.关于形状特性的附加信息加入到曲线中使得该方法对网络结构的细微变化很敏感，因此对形状非常敏感.
改进的盒覆盖算法流程图如表 2所示.
|
|
表 2 改进的盒覆盖算法流程图 Table 2 Flow chart of improved box covering algorithm |
本文中，用一个向量表示分形维数描述符，该向量包含通过回归计算得到的拟合多项式曲线的系数.
3 实验和结果为了验证本文提出的方法，使用LIDC-IDRI数据库和沈阳盛京医院的CT资料来进行实验.实验所用的CT图片的格式均为DICOM、像素均为512×512，尺寸大小也一致(为了让一些实验结果看得更清楚进行了适当放大)，肺结节在CT图片中的形状坐标已被医生标注.从肺癌数据库中选择500个肺结节作为实验样本，随机选择50%作为训练样本，剩下的作为测试样本.表 3为实验所用的部分数据.
|
|
表 3 实验所用部分数据 Table 3 Part data used in the experiment |
在所有实验中，本文选取形状轮廓上两点之间欧式距离最大值的四分之一为阈值r，初始阈值Tmin=0.05，按照相同间隔Tinc=0.025递增, 直到达到最终的阈值Tmax=0.85的δTi转换生成不同表现形式的复杂网络，再通过改进的盒覆盖算法计算分形维数d.这样对于每个阈值Ti都对应着一个分形维数di.图 8为肺结节对应的复杂子网络的分形维数d和分形维数di与阈值Ti之间的关系.
|
图 8 肺结节对应的复杂子网络的分形维数d和分形维数di与阈值Ti之间的关系 Fig.8 Fractal dimension d of complex sub-networks corresponding to pulmonary nodules and the relationship between the fractal dimension di and the threshold value Ti (a)—肺结节; (b)—复杂子网络分形维数d; (c)—阈值Ti与分形维数di之间的关系. |
本文用TP，FN，FP和TN来评价本文提出的形状统计模型的准确性，再和4种基于描述符的方法用线性判别分析(LDA)[11]进行分类性能的比较：傅里叶描述符、Zernike矩、曲率描述符、Bouligand-Minkowski多尺度维和复杂网络[12].
|
(12) |
|
(13) |
|
(14) |
其中：SEN是本文算法用于肺结节分类的敏感度，表示使用本文算法得到的结果中真实肺结节的发现率；SPE是算法的特异性，表示本文算法得到的结果中非肺结节的发现率；ACC是本文算法的准确率，表示算法用于分类的判定能力.表 4为本文算法对各种肺结节的分类结果统计，从中可看出用本文算法建立的肺结节形状模型对各类型肺结节进行分类，其平均敏感度SEN为86.00%，平均特异度SPE为69.26%和平均正确率ACC为78.80%.
|
|
表 4 本文算法对各种肺结节的分类结果 Table 4 Classification results of various pulmonary nodules by proposed algorithm |
图 9为本文算法与其他4种算法分类性能对比的ROC曲线图.表 5为本文算法与其他4种算法分类性能的对比表.从图 9和表 5中能看出，本文算法的分类性能均要好于其他算法，这是因为本文在分析形状的时候考虑了形状的局部纹理特征，且用pi-LBP提取形状的局部纹理特征，这种方法提取的LBP值更准确；本文把形状轮廓看作是由相连的像素点组成的，不需要手动选择采样点来表示形状，减少人工操作带来的误差；还不需要将所有的形状在同一坐标系下进行对齐，以便采样点的选择；本文将形状建模为一个复杂网络，这对形状的刚性形变(如旋转、平移和缩放)具有更好的鲁棒性.因此，提高了对形状建模的准确率.
|
图 9 ROC曲线 Fig.9 ROC curves |
|
|
表 5 各种算法的分类性能的对比 Table 5 Comparison of classification performance of various algorithms |
本文提出了一种基于复杂网络理论和分形维数的肺结节形状模型.首先利用像素点之间的欧式距离建立基础网络，然后利用在基于阈值Ti的动态演化过程中得到的复杂网络有效地表示一个形状轮廓，并通过分形维数d来估计网络的复杂性来描述其特征.对提出的方法进行实验之后，结果表明本文方法能够建立表现良好的肺结节形状模型.
本文还有一些问题有待进一步研究，如研究新的复杂网络建立方法和新的复杂网络演化方式，以便得到描述肺结节形状特征的多尺度描述符，不再使用单一的分形维数对网络特征进行描述.
| [1] |
Zhang M, Golland P.
Statistical shape analysis:from landmarks to diffeomorphisms[J]. Medical Image Analysis, 2016, 33: 155–158.
DOI:10.1016/j.media.2016.06.025 |
| [2] |
Ma Z, Kang B, Lyu K, et al.
Nonlinear radon transform using Zernike moment for shape analysis[J]. Computational and Mathematical Methods in Medicine, 2013, 2013(3): 208402.
|
| [3] |
Zhao L, Wang J H.
Research on wood cell shape analysis methods based on Fourier descriptors[J]. International Journal of Multimedia & Ubiquitous Engineering, 2015, 10(2): 61–74.
|
| [4] |
Shen W, Wang Y, Bai X, et al.
Shape clustering:common structure discovery[J]. Pattern Recognition, 2013, 46(2): 539–550.
DOI:10.1016/j.patcog.2012.07.023 |
| [5] |
Albert R.
Barab'asi:statistical mechanics of complex networks[J]. Reviews of Modern Physics, 2002, 74(1): 47–97.
DOI:10.1103/RevModPhys.74.47 |
| [6] |
Florindo J B, Backes A R, De Castro M, et al.
A comparative study on multiscale fractal dimension descriptors[J]. Pattern Recognition Letters, 2012, 33(6): 798–806.
DOI:10.1016/j.patrec.2011.12.016 |
| [7] |
Ojala T, Pietikainen M, Maenpaa T.
Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2002, 24(7): 971–987.
|
| [8] |
Lin Q, Qi W.Multi-scale local binary patterns based on path integral for texture classification[C]// IEEE International Conference on Image Processing.London, 2015: 26-30.
|
| [9] |
Tricot C.Curves and fractal dimension[M].Berlin: Springer, 1995.
|
| [10] |
Song C, Gallos L K, Havlin S, et al.
How to calculate the fractal dimension of a complex network:the box covering algorithm[J]. Journal of Statistical Mechanics Theory & Experiment, 2007(3): 297–316.
|
| [11] |
Liu Z P.Linear discriminant analysis[M].New York: Springer, 2013.
|
| [12] |
Backes A R, Casanova D, Bruno O M.
Texture analysis and classification:a complex network-based approach[J]. Information Sciences, 2013, 219: 168–180.
DOI:10.1016/j.ins.2012.07.003 |