灰度共生矩阵GLCM、图像均值、方差、熵、Percentile mean and SD、Kurtosis、Skewness

Halcom

灰度共生矩阵GLCM、图像均值、方差、熵、Percentile mean and SD、Kurtosis、Skewness等特征：

1. function feature = feature_extracted2( im2,x,y )
   
2. [h,w] = size(im2);
   
3. % Gray-level histogram features
   
4. % [n_count,x_value]  = imhist(im2);
   
5. Nt = h*w;
   
6. % prob = n_count./Nt;     % Pi概率值
   
7. % x_value = prob;
   
8. % x_value = im2(:);
   
9. 
   
10. for i=1:length(x)
    
11.     x_value(i) = im2(x(i),y(i));
    
12. end
    
13. 
    
14. mean1 = mean(x_value);  % 均值
    
15. SD1 = std(x_value);     % 方差
    
16. % Percentile mean and SD
    
17. x_value1 = sort(x_value,'ascend');
    
18. mean_beta_10 = mean( x_value1(floor(0.90*length(x_value1)):end ) );
    
19. mean_beta_25 = mean( x_value1(floor(0.75*length(x_value1)):end ) );
    
20. mean_beta_50 = mean( x_value1(floor(0.50*length(x_value1)):end ) );
    
21. SD_beta_10 = std( x_value1(floor(0.90*length(x_value1)):end ) );
    
22. SD_beta_25 = std( x_value1(floor(0.75*length(x_value1)):end ) );
    
23. SD_beta_50 = std( x_value1(floor(0.50*length(x_value1)):end ) );
    
24. % Kurtosis
    
25. % kurtosis1 = sum((x_value-mean1).^4)/Nt ./ sqrt( sum((x_value-mean1).^2)./Nt ).^4;
    
26. % kurtosis1 = sum((x_value-mean1).^4) ./ sqrt( sum((x_value-mean1).^2) ).^4;
    
27. % Skewness
    
28. % skewness1 = sum((x_value-mean1).^3)/Nt ./ sqrt( sum((x_value-mean1).^2)./Nt ).^3;
    
29. % skewness1 = sum((x_value-mean1).^3) ./ sqrt( sum((x_value-mean1).^2) ).^3;
    
30. kurtosis1 = kurtosis(double(x_value));
    
31. skewness1 = skewness(double(x_value));
    
32. 
    
33. % Gray-Level Co-Occurrence Matrix(GLCM)
    
34. %                    Angle     OFFSET
    
35. %                    -----     ------  
    
36. %                    0         [0 D]   
    
37. %                    45        [-D D]
    
38. %                    90        [-D 0]
    
39. %                    135       [-D -D]  
    
40. glcms_0 = graycomatrix(im2, 'offset', [0 1]);       % Angle = 0
    
41. glcms_45 = graycomatrix(im2, 'offset', [-1 1]);     % Angle = 45
    
42. glcms_90 = graycomatrix(im2, 'offset', [-1 0]);     % Angle = 90
    
43. glcms_135 = graycomatrix(im2, 'offset', [-1 -1]);   % Angle = 135
    
44. stats_0 = graycoprops(glcms_0, 'all');
    
45. stats_45 = graycoprops(glcms_45, 'all');
    
46. stats_90 = graycoprops(glcms_90, 'all');
    
47. stats_135 = graycoprops(glcms_135, 'all');
    
48. 
    
49. Contrast_0 = stats_0.Contrast;
    
50. Correlation_0 = stats_0.Correlation;
    
51. Energy_0 = stats_0.Energy;
    
52. Homogeneity_0 = stats_0.Homogeneity;
    
53. % 求解熵
    
54. tGLCM = glcms_0./sum(glcms_0(:));
    
55. entropy_0 = -sum(sum(tGLCM.*log(tGLCM+eps)));
    
56. 
    
57. Contrast_45 = stats_45.Contrast;
    
58. Correlation_45 = stats_45.Correlation;
    
59. Energy_45 = stats_45.Energy;
    
60. Homogeneity_45 = stats_45.Homogeneity;
    
61. % 求解熵
    
62. tGLCM = glcms_45./sum(glcms_45(:));
    
63. entropy_45 = -sum(sum(tGLCM.*log(tGLCM+eps)));
    
64. 
    
65. Contrast_90 = stats_90.Contrast;
    
66. Correlation_90 = stats_90.Correlation;
    
67. Energy_90 = stats_90.Energy;
    
68. Homogeneity_90 = stats_90.Homogeneity;
    
69. % 求解熵
    
70. tGLCM = glcms_90./sum(glcms_90(:));
    
71. entropy_90 = -sum(sum(tGLCM.*log(tGLCM+eps)));
    
72. 
    
73. Contrast_135 = stats_135.Contrast;
    
74. Correlation_135 = stats_135.Correlation;
    
75. Energy_135 = stats_135.Energy;
    
76. Homogeneity_135 = stats_135.Homogeneity;
    
77. % 求解熵
    
78. tGLCM = glcms_135./sum(glcms_135(:));
    
79. entropy_135 = -sum(sum(tGLCM.*log(tGLCM+eps)));
    
80. 
    
81. feature = [mean1,SD1,mean_beta_10,mean_beta_25,mean_beta_50,SD_beta_10,SD_beta_25,SD_beta_50,...
    
82.     kurtosis1,skewness1,...
    
83.     Contrast_0,Correlation_0,Energy_0,Homogeneity_0,entropy_0...
    
84.     Contrast_45,Correlation_45,Energy_45,Homogeneity_45,entropy_45...
    
85.     Contrast_90,Correlation_90,Energy_90,Homogeneity_90,entropy_90...
    
86.     Contrast_135,Correlation_135,Energy_135,Homogeneity_135,entropy_135];

复制代码

图像滤波函数：

1. function im5 = freqfilter2(im1,delta)
   
2. % 滤波器
   
3. % input:
   
4. %     M,N:频域滤波器的尺寸
   
5. %     delta：带阻滤波器的截止频率
   
6. % output:
   
7. %       H：M x N的矩阵，表示频域滤波器矩阵，数据类型为double，
   
8. if ~isa(im1,'double')
   
9.     im1 = double(im1)/255;
   
10. end
    
11. % im2 = fft2(im1);      % 傅里叶变换
    
12. % im3 = fftshift(im2);  % 中心化
    
13. 
    
14. % [N1, N2] = size(im3);
    
15. for i = 1:3
    
16.     for j = 1:3
    
17.         h(i,j) = -1/3.1415926/(delta.^4).*(1-(i.^2+j.^2)/2/delta/delta).*exp( -((i.^2+j.^2)./2/delta/delta) );
    
18.     end
    
19. end
    
20. im5 = conv2(im1,h,'same');
    
21. 
    
22. % result = ifftshift(result);    % 反中心化
    
23. % im4 = ifft2(result);           % 反变换
    
24. % % im5 = im2uint8(real(im4));     % 滤波图像
    
25. % im5 = abs(real(im4));     % 滤波图像
    
26. 
    
27. minV = min(min(im5));
    
28. maxV = max(max(im5));
    
29. im5 = 255.0*(im5-minV)./(maxV-minV);
    
30. 
    
31. end

复制代码

图像特征求解，主程序：

1. clc,clear,close all
   
2. warning off
   
3. im = dicomread('例图2.dcm');
   
4. im = mat2gray(im);
   
5. % figure(1),imshow(im)
   
6. % 扣出感兴趣区域
   
7. bw = im>0;
   
8. % figure(2),imshow(bw)
   
9. [x,y] = find(bw==1);
   
10. im1 = im(min(x):max(x),min(y):max(y));
    
11. bw = im2bw(im1);
    
12. [x,y] = find(bw==1);
    
13. minV = min(min(im1));
    
14. maxV = max(max(im1));
    
15. im1 = 255.0*(im1-minV)./(maxV-minV);
    
16. % im1 = im2uint8(im1);
    
17. figure(3),imshow(im1,[])
    
18. [h,w] = size(im1);
    
19. 
    
20. delta = [0.0,1.0,1.5,2.0,2.5];
    
21. feature_total = [];
    
22. for i = 1:length(delta)
    
23.    
    
24.     if isequal(delta(i),0)
    
25.         im2 = im1;
    
26.     else
    
27.         % 图像滤波
    
28.         im2 = freqfilter2(im1,delta(i));
    
29.     end
    
30.     % 特征求解
    
31.     feature = feature_extracted2( im2,x,y );
    
32.    
    
33.     figure,imshow(im2,[])
    
34.     feature2(i,:) = feature;
    
35.     feature_total = [feature_total,feature];
    
36.    
    
37. end
    
38. %% 特征绘图
    
39. figure(10),plot(feature_total,'ro-','linewidth',2);
    
40. grid on

复制代码

参考资料：【1】结肠癌 淋巴转移 nomogram JCO：http://pan.baidu.com/s/1cgYzO6
【2】Development and validation of a radiomics nomogram for preoperative predicting lymph node metastasis in colorectal cancer ：http://pan.baidu.com/s/1jHDBSEU
【3】图像：http://pan.baidu.com/s/1nuXmiXr