PET/CT图像的纹理特征提取

---

目的

检验纹理特征对3d-PET/CT图像分类的效果。

简介

在使用传统分类器的时候，和深度学习不一样，我们需要人为地定义图像特征，其实CNN的卷积过程就是一个个的滤波器的作用，目的也是为了提取特征，而这种特征可视化之后往往就是纹理、边缘特征了。因此，在人为定义特征的时候，我们也会去定义一些纹理特征。在这次实验中，我们用数学的方法定义图像的纹理特征，分别计算出来后就可以放入四个经典的传统分类器（随机森林，支持向量机，AdaBoost，BP-人工神经网络）中分类啦。

工具

我使用的工具是MATLAB 2014b，建议版本高一点好，因为里面会更新很多的函数库。实验过程尽量简化，本实验的重点是检验纹理特征对PET/CT图像分类的效果，因此，有些常规的代码我们就用标准的函数库足够啦。

参考文档

PORTS 3D Image Texture Metric Calculation Package

---

1. 直方图-histogram

直方图描述的是一幅图像中各个像素的分布情况，也就是一个对像素做的统计图。
对于一幅灰度图像 I，它每个像素值的范围是0-255，我们对这些像素点做一个统计，遍历整幅图像，统计像素值0,1,2,3,…,255分别出现的次数。统计完以后相当于我们有了256个频数（次数），再把它们转化成频率，也就是每个频数除以总频数：

p(i) = P(i) / ∑P

以像素值作为横坐标，对应的频率作为纵坐标，就可以得到这个灰度图像 I 的直方图啦。

1.1 举栗子：CT图像的直方图

1. CT图像的像素值范围是-1000~1000。相当于我们需要统计2000个像素值的频数，这样划分的粒度有点太细了，因此

2. 将这-1000~1000的区间20等分，每个像素值投射到20个值。直接导致的结果是图像看上去不那么丰富了，但是这样有利于计算。

3. 分别统计这20个像素值出现的频数，除以总频数转化成频率。这样频率介于[0,1]，并且加和为1.

4. 以20个像素值为横坐标，对应的频率为纵坐标，即可画出这个CT图像的直方图。

The end of this 栗子.

1.2 直方图的代码实现

%%%%%
%%%%% Histogram-based computations:
%%%%%

% Compute the histogram of the ROI and probability of each voxel value:
vox_val_hist = zeros(num_img_values,1);
for this_vox_value = 1:num_img_values
    vox_val_hist(this_vox_value) = length(find((img_vol_subvol == this_vox_value) & (mask_vol_subvol == 1) ));
end

% Compute the relative probabilities from the histogram:
vox_val_probs = vox_val_hist / num_ROI_voxels;


% Compute the histogram_based metrics:
texture_metrics(1:6) = compute_histogram_metrics(vox_val_probs,num_img_values);

1.3 基于直方图的PET/CT纹理特征

包括六个值，分别是：

(1) Mean

(2) Variance

(3) Skewness – set to 0 when σ=0

(4) Kurtosis – set to 0 when σ=0 (NOTE: “Kurtosis” and “Excess Kurtosis” differ in that Excess Kurtosis = Kurtosis – 3).

(5) Energy

(6) Entropy (NOTE: We will differentiate between the various entropy calculations in this document, specifying the distribution from which the entropy is computed)

1.4 纹理特征计算实现

%%% Overhead:

% The numerical values of each histogram bin:
vox_val_indices = (1:num_img_values)';

% The indices of non-empty histogram bins:
hist_nz_bin_indices = find(vox_val_probs);


%%% (1) Mean 
metrics_vect(1) = sum(vox_val_indices .* vox_val_probs);

%%% (2) Variance
metrics_vect(2) = sum( ((vox_val_indices - metrics_vect(1)).^2) .* vox_val_probs );

%%%%% IF standard variance is zero, so are skewness and kurtosis:
if metrics_vect(2) > 0

    %%% (3) Skewness
    metrics_vect(3) = sum( ((vox_val_indices - metrics_vect(1)).^3) .* vox_val_probs ) / (metrics_vect(2)^(3/2));

    %%% (4) Kurtosis
    metrics_vect(4) = sum( ((vox_val_indices - metrics_vect(1)).^4) .* vox_val_probs ) / (metrics_vect(2)^2);
    metrics_vect(4) = metrics_vect(4) - 3;

else

    %%%