图解医学影像纹理特征

芷亦止兮

在医学影像分析中，纹理（texture）特征被广泛地应用于定量化描述病灶的特性。相比病灶影像的统计学特征和形状特征，纹理特征的计算往往较为复杂，不易于搞懂其背后的原理。然而，学习过程中建立对概念直观的认识（intuition）对理解问题的本质很有帮助，这也是本文的出发点。
  本文尝试解决如下问题：
  

• 医学影像特征包的官方文档给了纹理特征数学定义，但是理解起来不够直观
  
• 网络上有许多介绍纹理特征的帖子，但是每个帖子介绍的内容有限；
  
• 少有文章关注不同类型纹理特征之间的比较，不易于理解纹理特征背后的逻辑
  

  本文希望通过一个统一的描述，为研究医学影像的入门者提供帮助。如有疏漏还请大佬们多多指正。
  0 简介

  在介绍具体的特征之前，我们先说明一下什么是“致密的”和“粗糙的”和纹理。下图中，左侧的纹理相对致密，右侧的纹理相对粗糙，因为右侧的图中相邻区域的颜色变化相对较小，我们肉眼可见若干个子区域。但是这种对于纹理的描述是定性的，不利于我们建立数学模型来对图像的性质做出判断，因此需要引入定量的纹理特征。
  

  致密的纹理（左）和粗糙的纹理（右）
   
  简单来说，纹理特征是一种“二阶特征”，他们不是从图像上直接得出的，而是先通过某种计算将原始图像的特性提取出来并存在一个中间矩阵中，之后在这个中间矩阵上定义一系列的统计量，作为图像的纹理特征。在医学影像研究中，根据所用的中间矩阵的不同，可以将纹理特征的分为若干类别。这里使用影像特征提取包pyradiomics中对于纹理特征的分类：
  

• Gray-level co-occurrence matrix (GLCM) based features
  
• Gray-level size zone matrix (GLSZM) based features
  
• Gray-level run length matrix (GLRLM) based features
  
• Gray-level dependence matrix (GLDM) based features
  
• Neighboring gray tone difference matrix (NGTDM) based features
  

  我们一个一个来看。
   

  1 GLCM与相关纹理特征

  灰度共生矩阵 (GLCM) 是最常用的计算图像纹理的手段。这里通过MATLAB帮助文档里的一个例子说明GLCM是怎么计算的：
  

  图的左侧是原始图像矩阵，右侧是计算出的GLCM。所谓“灰度共生”，顾名思义是研究两个灰阶同时出现的情况。那么这里就产生了两个问题：1. 要考虑哪两个灰阶；2. 在图上怎么考虑两个灰度。
  第一个问题的答案很简单，所有的灰阶组合我们都考虑。GLCM是一个方阵，矩阵的维度是原始图像的灰阶数。上面例子中原图一共有8个灰度值，因此GLCM的维度是8。要考虑的两个灰度分别作为行和列，比如灰度1和2的共生次数是2，这个值就保存在GLCM(1,2)中。
  第二个问题需要我们设定“共生灰度”对应像素的距离 

（像素个数）和方向  

（0~360度）。在上面的例子中显然距离为1，方向为0度（向右），相当于对于灰度1和2的共生，我们仅考虑1在2的左侧，且对应像素之间的距离为1个像素的情况。
  所以GLCM其实是一个计数矩阵，它保存了图像中所有灰度组合在我们定义的距离和方向条件下“共生”的频数。直观上来说，当GLCM中数值集中在主对角线上时，原图的纹理更为粗糙。这是因为GLCM主对角线附近保存了两个相差不大的灰度值“共生”的频数，主对角线附近元素的值越大，说明原图中存在很多相邻像素灰度值相差不大的情况，肉眼看上去影像中就会存在若干个面积较大的子区域，回想一下前面石子图片的纹理便能够理解。
  实际计算GLCM时候一般距离都设为1，同时会尝试所有的角度（每45度计算一次）并取平均。因此得到的GLCM一般是一个对称的矩阵。
  下面通过一个例子来说明如何在GLCM上计算纹理特征，本文所有具体的纹理特征定义均来自pyradiomics包的文档。
  https://pyradiomics.readthedocs.io/en/latest/
  根据GLCM计算图像对比度（contrast）的公式如下：
  

  式中  

是图像中的灰阶数，  

是频率形式的GLCM矩阵（即GLCM每个元素除以GLCM所有元素之和，得到每个元素的频率）。
  如果原图中相邻两个像素灰度值差别很大（即  i-j 的绝对值很大）的情况很多（即 p(i,j)很大），图像的纹理就越致密，contrast的取值就会很大。
   
  2 GLSZM与相关特征

  

  灰度尺寸区域矩阵（GLSZM）的计算过程如上图所示（底图来自Wikipedia）。左侧是原始图像，右侧是计算出的GLSZM矩阵，上下分别是两个例子。GLSZM矩阵的每行代表一个灰阶，每列代表连通域大小。对于每个灰阶（这里有1-4共4个不同的灰阶），我们逐个考察它在图像上的连通域（2D图像一般采用8连通域）情况。比如对于灰阶3，图像中它只有一个连通域，该连通域的大小为3，于是我们在GLSZM(3,3)的位置填入1，同时在GLSZM(3,1)和(3,2)填入0。同理，对于灰阶4，图像中有三个连通域，两个大小为1，一个大小为3，所以GLSZM(4,1)为2，GLSZM(4,3)为1。
  GLSZM也是一个计数矩阵，它保存了图像中所有灰阶的连通域大小和个数信息。直观上来说，当GLSZM较宽，非零元素集中在右侧的时候，原图的纹理更为粗糙。这是因为GLSZM的宽度是由最大连通域的大小决定的，同时GLSZM右侧的元素对应面积较大的连通域的个数。GLSZM越宽，说明原图中存在较大的连通域；GLSZM右侧非零元素越多，说明原图中连通域的面积都较大，肉眼看上去影像中就会存在若干个面积较大的子区域。
  看一个由GLSZM计算得到的纹理特征的例子：Large Area Emphasis
   

  式中  

是图像中的灰阶数（行数）,  

是图像中的最大连通域面积（列数） ，

 是频率形式的GLSZM矩阵（即GLSZM每个元素除以GLSZM所有元素之和，得到每个元素的频率）。
  对于原图中面积较大的区域，  j 的取值也较大，因此LAE给予面积大的区域更大的权重，起到了“emphasis”的作用。
   
  3 GLRLM与相关特征

  

  灰度游程矩阵（GLRLM）的计算过程与GLSZM十分相似，如上图所示（底图来自Xu et al., Applied Sciences, 2019）。左侧是原始图像，右侧是计算出的GLRLM矩阵。GLRLM矩阵的每行代表一个灰阶，每列代表Run-Length大小。这里的Run-Length表示沿着某个方向的连通域（称为一个“Run”）长度。对于每个灰阶（这里有1-4共4个不同的灰阶），我们逐个考察它在图像上的Run-Length情况。比如对于灰阶2，如果设定方向为向右，图像中它有3个连通域，其中两个的Run-Length为1，一个的Run-Length为2，于是我们在GLRLM(2,1)的位置填入2，同时在GLRLM(2,2)填入2。
  GLRLM同样是一个计数矩阵，它保存了图像中所有灰阶的“Run”的长度和个数信息。直观上来说，当GLRLM较宽，非零元素集中在右侧的时候，原图的纹理更为粗糙。这是因为GLRLM的宽度是由最长的Run的长度决定的，同时GLRLM右侧的元素对应较长的Run的个数。GLRLM越宽，说明原图中存在较长的Run；GLRLM右侧非零元素越多，说明原图中的Run都比较长，肉眼看上去影像中就会存在若干个面积较大的长条状子区域。
  实际计算GLRLM时候一般会尝试所有的角度（每45度计算一次）并取平均，以考虑所有方向上的连通信息。
  看一个由GLRLM计算得到的纹理特征的例子：Large Run Emphasis
   

  式中

 是图像中的灰阶数（行数）,  

 是图像中的最大的Run-Length（列数） ，  

是频率形式的GLRLM矩阵（即GLRLM每个元素除以GLRLM所有元素之和，得到每个元素的频率）。
  对于原图中较长的Run，  j 的取值也较大，因此LRE给予Run-Length大的区域更大的权重，起到了“emphasis”的作用。
   

  4 GLDM与相关特征

  

  灰度依赖矩阵（GLDM）的计算过程与GLSZM也十分类似，如上图所示（底图来自pyradiomics官方文档）。上边是原始图像，下边是计算出的GLDM矩阵。GLDM矩阵的每行代表一个灰阶，每列代表“依赖”像素的个数。这里对于“依赖”的定义如下：设定一个距离

（采用8邻域）和阈值

 ，对于原图中的每个像素（称为目标像素），考察与其距离在

 内的所有像素点，如果这些像素的灰度值与目标像素灰度值的差距小于 

 ，则认为它们和目标像素是“依赖”的。对于每个灰阶（这里有1-4共4个不同的灰阶），我们逐个考察该灰阶所有像素点在图像上的"依赖"像素个数情况。比如对于灰阶1，如果设定距离为1，阈值为0，图像中(3,2)位置的像素点有1个依赖像素，(2,4)和(3,4)位置的像素点各有2个依赖像素，(3,3)位置的像素点有1个依赖像素，因此GLDM(1,2)填入1，(1,3)填入2，(1,4)填入1。
  GLDM同样是一个计数矩阵，它保存了图像中所有像素点“依赖”像素的数目和出现次数的信息。直观上来说，当GLDM较宽，非零元素集中在右侧的时候，原图的纹理更为粗糙。这是因为GLDM的宽度是由“依赖”元素数量的最大值决定的，同时GLDM右侧的元素对应灰阶有较多数目的“依赖”元素。GLDM越宽，右侧非零元素越多，说明原图中有限距离的邻域内连通域的面积越大。
  实际计算GLDM时候一般设定距离  

为1，阈值  

 为0。当然也可以根据需要加以调整。
  看一个由GLDM计算得到的纹理特征的例子：Large Dependence Emphasis
   

  式中  

是图像中的灰阶数（行数）,  

是图像中的最大的“依赖”元素数目（列数） ，  

是频率形式的GLDM矩阵（即GLDM每个元素除以GLDM所有元素之和，得到每个元素的频率）。
  对于原图中“依赖”像素数目较多的情况，  j 的取值也较大，因此LDE给予“依赖”像素多的区域更大的权重，起到了“emphasis”的作用。
   

  5 NGTDM与相关特征

  

  相邻灰度差分矩阵（NGTDM）的计算过程如上图所示（底图来自pyradiomics官方文档，有修改）。上方是原始图像，中间是计算出的NGTDM矩阵，矩阵中的  i 是灰阶的值，  

是该灰阶对应像素的计数， 

是该灰阶出现的频率，即  

， 

是灰度差的绝对值。
  计算示例如下方子图所示。对于灰阶1，原图中共有6个像素，

；出现的频率  

；对于每个灰度为1的像素，计算其与邻域内像素灰度差的绝对值，比如上图中红色框内的像素都属于左上角灰度值为1像素的相邻像素，这3个像素值的和为10，因此对应灰度差的绝对值是  

，其他像素计算同理，  

就是每个灰阶所有像素与其相邻像素灰度差绝对值的求和。
  由NGTDM计算得到的纹理特征的例子：粗糙度（Coarseness）
  

  式中

 是图像中的灰阶数（行数）。不难理解，  

是灰度差绝对值的加权平均，这个加权平均值越小，说明灰阶对应像素与其相邻像素之间的差异越小，图像上的空间变化也就越小，图像的纹理就更为粗糙。
   

  总结

  总结一下几类影像特征计算过程中的共同点：
  

• 都是针对原图中的每个灰阶，或者灰阶之间的相互关系计算出一个中间矩阵，再从中间矩阵中计算获得一个统计量作为影像特征
  
• 除了NGTDM，其他中间矩阵均为计数矩阵
  

  从上面的分析不难看出，GLSZM, GLRLM 和 GLDM 这三种矩阵十分相似，实际上基于这三种中间矩阵所定义的影像特征的数量和意义都是相近的。我们把他们放到一起来看，尝试寻找其中的异同点：
  

  

• GLSZM与相关特征衡量的是图像中2D/3D连通域的大小和数目
  
• GLRLM与相关特征衡量的是图像中1D连通域（这里称为Run）的大小和数目
  
• GLDM与相关特征衡量的是图像中2D/3D连通域的大小和数目，但是限制了关注邻域的面积
  

  同时这正因为这种相似性，我们得到的影像特征之间并非完全独立，因此在应用到下游分析时，可以考虑做一下PCA或者特征选择，避免特征之间的correlation对结果的影像。
  最后，
  

• 医学影像研究的发展过程中，人们提出了多种类型的影像特征。但也由此产生了一些问题——大家各自用自己的一套特征，不利于人们将现有的研究统一起来。不过好在现在学界已经有了共识，Radiology杂志今年5月发文提出了IBSI（其实2016年就已经发布在arxiv上了）。使用IBSI中提供的影像特征有利于不同研究结果的对比分析，而且有助于医学影像研究工作的可重复性。
  
• 虽然深度学习在医学影像研究中的应用越来越广泛，但是影像特征能够提取影像的全局特征，同时在可解释性等方面仍然具有优势。因此影像特征（尤其是纹理特征）在医学影像中的应用仍然呈快速增长的趋势。
  
• 随着技术的发展，很可能会有新的特征被提出、验证和应用。因此在深入认识目前已有影像特征的同时，也要保持对新特征的关注。