python 提取灰度共生矩阵纹理特征

一阶统计量特征，或者说灰度直方图特征，主要思想是对整张图像，或者图像中的感兴趣区域进行一些统计学计算，求得相应的统计量，用于在灰度层面描述图像。需要注意的是，一阶统计量特征仅适用于单通道的灰度图像，如果想对彩色图像提取一阶统计量特征，需要先对彩色图像进行灰度化操作。本文的代码展示的是如何对图像的感兴趣区域提取一阶统计量特征。

本文使用的例子是一例骨肿瘤患者的CT图像和它对应的肿瘤区域标记，原始CT图像数据是三维的.nii格式，在Matlab中打开和处理这种格式的数据需要借助NIFTI工具库，NIFTI工具库大家可以在GitHub中搜索下载，或者通过下面的链接下载，下载解压后将其添加到Matlab的路径中即可使用。

链接：https://pan.baidu.com/s/1giKDEtsVs-YVA5p-dtuhyA 
提取码：grcm 

提取三维数据中的一层进行演示，代码如下，使用的软件版本为Matlab 2018b，首先介绍数据的读取和预处理。

% 读取三维原始图像image和对应的感兴趣区域标记文件label
% image和label的尺寸为512×512×65，65是层数
image = load_nii('C:\Users\Administrator\Desktop\image.nii');
image = image.img;
label = load_nii('C:\Users\Administrator\Desktop\label.nii');
label = label.img;
% 选取原始图像和标记文件的第35层，并对标记文件进行二值化操作
% 使得标记中感兴趣区域为1，其他区域为0
img = image(:,:,35);
mask = label(:,:,35);
mask = imbinarize(mask);
% 医学影像灰度值跨度很大，所以将img的像素值归一化到0-255范围内，便于后续分析
img = double(img);
ymin = 0;
ymax = 255;
xmin = min(min(img));
xmax = max(max(img));
img = round((ymax-ymin)*(img-xmin)/(xmax-xmin) + ymin);
img = uint8(img);
% 将第35层的原图和标记展示出来
imshow(img, []);
figure
imshow(mask);
% 求原图在感兴趣区域上的灰度直方图并可视化出来
% 变量level存放的是图像的灰度级，这里为256×1的向量，存放0-255之间的整数
% 变量count存放的是感兴趣区域中每个像素值的数量
[count, level] = imhist(img(mask == 1));
figure
imhist(img(mask == 1));

上面代码得到的数据读取和预处理结果如下所示，第一幅图是像素值归一化到0-255后的第35层原始数据，第二幅图是它对应的二值化后的感兴趣区域标记图像，白色区域为手工勾画的肿瘤区域，第三幅图是将感兴趣区域投射回原图的示意图，第四幅图是原始图像感兴趣区域内的灰度直方图。

下面的代码根据感兴趣区域和灰度直方图信息对原图提取一些统计量，用于在灰度层面上描述图像。

% 变量roi中存放的是对应于感兴趣区域的原始图像中的灰度值
roi = img(mask == 1);
roi = double(roi);
% 变量pixelnum中存放的是感兴趣区域（白色区域）内的像素数
pixelnum = sum(sum(mask == 1));
% 定义第一维特征为感兴趣区域内所有像素值的平均值
histofeats(1) = sum(roi) / pixelnum;
% 定义第二维特征为感兴趣区域内所有像素值的中位数
histofeats(2) = median(roi);
% 定义第三维特征为感兴趣区域内所有像素值的标准差
histofeats(3) = std(roi);
% 定义第四维特征为感兴趣区域内所有像素值的偏度
histofeats(4) = skewness(roi);
% 定义第五维特征为感兴趣区域内所有像素值的峰度
histofeats(5) = kurtosis(roi);
% 定义第六维特征为感兴趣区域内所有像素值的最小值
histofeats(6) = min(roi);
% 定义第七维特征为感兴趣区域内所有像素值的最大值
histofeats(7) = max(roi);
% 定义第八维特征为感兴趣区域内所有像素值的熵
histofeats(8) = entropy(roi);
% 定义第九维特征为感兴趣区域内所有像素值的能量
p = count/ pixelnum;
histofeats(9) = p' * p;
% 定义第十维特征为感兴趣区域内所有像素值的范围
histofeats(10) = max(roi) - min(roi);
% 定义第十一维特征为感兴趣区域内所有像素值的平均绝对离差
minu = roi - sum(roi) / pixelnum;
histofeats(11) = sum(abs(minu)) / pixelnum;
% 定义第十二维特征为感兴趣区域内所有像素值的均方根
histofeats(12) = sqrt(sum(roi.^2) / pixelnum);

代码中的特征仅仅是例子，大家可以发挥自己的想象力，设计出更多特征用于实际任务。