本文已参加「新人创造礼」活动,一起敞开掘金创造之路

图画纹路特征

 本文首要介绍医学影像范畴常用到的根据图画灰度值改变所衍生出的各项纹路特征,理论部分首要参阅了文献1,一起介绍了相关函数调用办法。

灰度共生矩阵

 灰度共生矩阵2,Gray-Level co-occurance matrix (GLCM)用于描绘两个空间上契合必定散布规律的灰度值对呈现的概率,矩阵(x,y)(x,y)处的值即表明了这样的灰度值对在图画中呈现的频次,用公式描绘为

p(x,y)=∑I{(i,j)∣img(i,j)=x,img(f(i,j))=y,i=1,2…h,j=1,2,..,w}p(x,y)=\sum I\{(i,j)|img(i,j)=x,img(f(i,j))=y, i=1,2…h,j=1,2,..,w\}

其间II为指示函数,w,hw,h为原始图片的宽度和高度,ff为界说空间散布的函数。

一般来说f(i,j)f(i,j)会设置为(i,j+1),(i−1,j+1),(i−1,j),(i−1,j−1)(i,j+1),(i-1,j+1),(i-1,j),(i-1,j-1)其间一个,也就分别对应于原图中(i,j)(i,j)元素的0、45、90、135方向,更进一步的的,有时候也认为这一空间散布应当是双向的,也就是类似f(i,j)=(i,j+1)or(i,j−1)f(i,j)=(i,j+1) or (i,j-1)。 给出当f(i,j)=(i,j+1)or(i,j−1)f(i,j)=(i,j+1)or (i,j-1)时的示例图如下,PP中(1,2)方位的为4,对应原图中灰度值对(1,2)、(2,1)的呈现次数为4次。

图像纹理特征(灰度共生矩阵等)解析和编程调用
图像纹理特征(灰度共生矩阵等)解析和编程调用

细节:核算灰度值对的时候并不考虑循环填充的状况,即并不认为最右边一行和左面一行是衔接的。矩阵的巨细为:灰度阶数*灰度阶数。

灰度游程矩阵

 灰度游程矩阵,Gray-Level run-length matrix (GLRLM)用于描绘沿某一方向的连续灰度值序列呈现的频率,矩阵(x,y)(x,y)处的值即表明原图中沿某一方向连续呈现yy个灰度值为xx像素点频次,用公式描绘为:

p(x,y)=∑I{(i,j)∣∀k∈[0,1,..,y−1],img(f(i,j,k))=x,i=1,2…h,j=1,2,..,w}p(x,y)=\sum I\{(i,j)\ |\ \forall k\in[0,1,..,y-1],img(f(i,j,k))=x ,i=1,2…h,j=1,2,..,w\}

其间II为指示函数,w,hw,h为原始图片的宽度和高度,ff为界说方向的函数,值得注意的是参加了一次较长游程核算的元素不会再参加下一次较短游程的核算,如连续的4个1值不会再被认定为存在2组连续2个1值。

一般而言f(i,j,k)f(i,j,k)会设置为(i,j+k),(i−k,j+k),(i−k,j),(i−k,j−k)(i,j+k),(i-k,j+k),(i-k,j),(i-k,j-k),相同是表征0、45、90、135方向,给出0方向事例如下:

图像纹理特征(灰度共生矩阵等)解析和编程调用
图像纹理特征(灰度共生矩阵等)解析和编程调用

原图II中灰阶3在水平方向上最长的连续长度为3个,呈现过1次;连续长度为2的呈现过1次,连续长度为1的呈现过4次,对应PP中第三行为[4,1,1,0,0]。

细节:相同不认为是循环填充的。矩阵巨细为:灰度阶数*max(w,h)max(w,h)

灰度尺度区域矩阵

 灰度尺度区域矩阵,Gray-Levle size zone martrix (GLSZM),用于描绘平面上相连连续灰度值序列的呈现频次。这个界说十分类似于灰度游程矩阵,但灰度游程矩阵要求连续序列是在某一方向上排布的,是一种一维的概念;而灰度尺度区域矩阵则只要求连续序列见构成连通域即可,是一种二维的概念。其数学公式相同为:

p(x,y)=∑I{(i,j)∣∀k∈[0,1,..,y−1],img(f(i,j,k))=x,i=1,2…h,j=1,2,..,w}p(x,y)=\sum I\{(i,j)\ |\ \forall k\in[0,1,..,y-1],img(f(i,j,k))=x ,i=1,2…h,j=1,2,..,w\}

只不过f(i,j,k)f(i,j,k)现在描绘了一个面积为yy的二维连通域而非一维的线段(连通域的概念可参见opencvn中的漫水填充算法),给出事例如下:

图像纹理特征(灰度共生矩阵等)解析和编程调用
图像纹理特征(灰度共生矩阵等)解析和编程调用

原图中II中由灰阶3组成的连通域有3个,巨细分别为3,5,1(这里的连通域认定是8连通域,因而右侧斜向和竖向的3共同构成巨细为5的连通域),对应PP中第三行为[1,0,1,0,1]。

细节:矩阵的巨细为:灰度阶数*最大连通域巨细

灰度依靠矩阵

 灰度依靠矩阵,Gray-Level dependence matrix (GLDM) ,用于描绘某一灰阶的周边区域某类灰度值散布的呈现频次,矩阵(x,y)(x,y)处的值表明原图中灰度值为xx的点满意:以该点为中心必定规模内灰度值与xx的差小于阈值\delta的点有yy个,这一条件的呈现频次。用数学公式表明如下:

p(x,y)=∑I{(i,j)∣≥∣img(f(i,j))−img(i,j)∣,i=1,2…h,j=1,2,..,w}p(x,y)=\sum I\{(i,j)\ |\ \delta\geq |img(f(i,j))-img(i,j)| ,i=1,2…h,j=1,2,..,w\}

其间f(i,j)f(i,j)用于界说这一规模,一般就是八范畴,而\delta则是能够容忍的差异阈值。给出=0\delta=0时事例如下:

图像纹理特征(灰度共生矩阵等)解析和编程调用
图像纹理特征(灰度共生矩阵等)解析和编程调用

原图中灰度值为3的点共有9个,八邻域规模内具有1个灰度值为3的街坊的点有4个,具有2个灰度值为3的街坊的点有4个,没有灰度值为3的街坊的点有1个,则PP中第三行为[1,4,4,0]。

细节:矩阵巨细为:灰度阶数*最大街坊数+1。加1是因为存在街坊数为0的状况。

相邻灰度差分矩阵

 相邻灰度差分矩阵,Neighboring gray tone difference matrix (NGTDM),用于描绘某灰阶和它周边规模内平均灰度值差异,和前面几个矩阵不同,相邻灰度差分矩阵的每列所代表的含义并不相同,先给出事例再分别叙说:

图像纹理特征(灰度共生矩阵等)解析和编程调用
图像纹理特征(灰度共生矩阵等)解析和编程调用

榜首张为原图,第二张给出对应NGTDM。NGTDM有三列,榜首列给出了不同灰阶在原图中呈现频次,第二列给出不同灰阶在原图中的呈现频率,第三列最为重要,给出某个灰阶在原图一切方位与其周边规模内平均灰度差异之和。一般周边规模就界说为八邻域,以s5s_5为例,原图中共有4个5,(0,2)(0,2)处的5八邻域规模内均值为2+5+1+3+25=2.6\frac{2+5+1+3+2}{5}=2.6,与5的差异为∣5−2.6∣=2.4|5-2.6|=2.4,同理得到其他方位处的差异为2.375、2.5、2.8,合计10.075。

灰度特征矩阵分析概括

 统一起来看的的话,其实灰度尺度区域矩阵能够看作是灰度游程矩阵的整合延申,将多个一维层面的的核算量扩张到了二位层面,都是在描绘某个灰阶值的连续状况;而相邻灰度差分矩阵和灰度依靠矩阵则是一体两面得刻画了某个灰阶与它范畴的关系,前者描绘相似性,后者描绘相异性;灰度共生矩阵则意图描绘图中灰阶的散布状况。这些值都是在描绘图画的纹路特征,如果原始图画的纹路越粗糙则越趋近于呈现大块的连续灰度值,从而灰度共生矩阵非零值呈现在对角线附近;灰度游程矩阵和灰度尺度区域矩阵的非零元素大都集中在右侧;灰度依靠矩阵的非零元素首要散布在矩阵右侧;相邻灰度差分矩阵的最后一列的值会较小。

 值得注意的是,一切的灰度特征矩阵的巨细和灰度阶数有关,一般咱们运用的图片巨细都是256阶的,再加上矩阵核算时的复杂程度于输入图片的巨细成正比,直接在原始图画上核算将是极端消耗时间,贮存结果也是极端消耗资源的。

  • 针对核算耗时问题,咱们能够调整原图的量化等级,下降灰度阶数,经过np.digitize()函数实现3
  • 针对核算结果过大的问题,咱们一般不直接运用核算出的灰度特征矩阵当作终究的特征,而是依此核算出一系列核算量来描绘,此处不再打开,详见pyrandiomics4的文档。

编程实现

 现在包括了这些纹路特征的python库较少,skimage中包括了核算灰度共生矩阵和相关特征的函数,但没有其他矩阵的核算办法;pyradiomics包括核算上述一切特征的办法,而且注释中具有具体的理论解说,但它只针对于医学图画,常规输入图片的特征无法进行求取(也可能是我没找到办法)。以下给出灰度游程矩阵代码,改写自5

def getGrayLevelRumatrix( array, theta):
    '''
    核算给定图画的灰度游程矩阵
    参数:
    array: 输入,需求核算的图画
    theta: 输入,核算灰度游程矩阵时采用的角度,list类型,可包括字段:['deg0', 'deg45', 'deg90', 'deg135']
    glrlm: 输出,灰度游程矩阵的核算结果
    '''
    P = array
    x, y = P.shape
    min_pixels = np.min(P)   # 图画中最小的像素值
    run_length = max(x, y)   # 像素的最大游行长度
    num_level = np.max(P) - np.min(P) + 1   # 图画的灰度级数
    deg0 = [val.tolist() for sublist in np.vsplit(P, x) for val in sublist]   # 原图分化按水平方向分化为多个序列
    deg90 = [val.tolist() for sublist in np.split(np.transpose(P), y) for val in sublist]   # 原图分化按竖直方向分化为多个序列
    diags = [P[::-1, :].diagonal(i) for i in range(-P.shape[0]+1, P.shape[1])]   # 运用diag函数获取45方向序列
    deg45 = [n.tolist() for n in diags]
    Pt = np.rot90(P, 3)   # 旋转后重复前一操作获取135序列
    diags = [Pt[::-1, :].diagonal(i) for i in range(-Pt.shape[0]+1, Pt.shape[1])]
    deg135 = [n.tolist() for n in diags]
    seqs={'deg0':deg0,'deg45':deg45,'deg90':deg90,'deg135':deg135}
    glrlm = np.zeros((num_level, run_length, len(theta)))  
    for angle in theta:
        for splitvec in range(0, len(seqs)):
            flattened =seqs[angle][splitvec] # 获取某一方向上的某个序列
            answer = []
            for key, iter in groupby(flattened):   # 运用itertools.groupby()获取连续灰度值长度
                answer.append((key, len(list(iter))))
            for ansIndex in range(0, len(answer)):
                glrlm[int(answer[ansIndex][0]-min_pixels), int(answer[ansIndex][1]-1), theta.index(angle)] += 1   # 每次将核算像素值减去最小值就能够填入GLRLM矩阵中
    return glrlm

更新:在github上找到一个十分强大的包叫做pyfeat6,里面包括了大部分根据上述灰度特征矩阵核算的核算量的函数,如果有需求的话也能够直接导入源文件核算原始的灰度特征矩阵。以灰度共生矩阵为例,调用pyfeats.glcm_features(img)就能够输出原图在0、45、90、135方向上的灰度共生矩阵核算出的14个核算量的均值。但pyfeat默许是不核算灰度值为0对应矩阵元素,而且有时候这一设置不能在features核算时显式修正,比方核算根据GLRLM特征时就没有给出可选项,需求自行去源文件pyfeats\textural\glrm.py中修正.

参阅

Footnotes

  1. 图解医学影像纹路特征 ↩

  2. 根据Python探求灰度共生矩阵(GLCM)那点事儿 ↩

  3. 灰度共生矩阵 ↩

  4. Radiomic Features ↩

  5. 特征提取之灰度游程(行程)矩阵-GLRLM ↩

  6. pyfeats 0.0.11 ↩