計算灰度共生矩陣GLCM

http://blog.youkuaiyun.com/cxf7394373/article/details/6988229

應用

    由上面的敘述知道，可以根據各種間距和角度計算灰度共生矩陣，下面程序中給定了間距，根據傳入的參數計算：

[cpp]  view plain copy

1. #define GLCM_DIS 3  //灰度共生矩陣的統計距離  
2. #define GLCM_CLASS 16 //計算灰度共生矩陣的圖像灰度值等級化  
3. #define GLCM_ANGLE_HORIZATION 0  //水平  
4. #define GLCM_ANGLE_VERTICAL   1  //垂直  
5. #define GLCM_ANGLE_DIGONAL    2  //對角  
6. int calGLCM(IplImage* bWavelet,int angleDirection,double* featureVector)  
7. {  
8.     int i,j;  
9.     int width,height;  
10.   
11.     if(NULL == bWavelet)  
12.         return 1;  
13.   
14.     width = bWavelet->width;  
15.     height = bWavelet->height;  
16.   
17.     int * glcm = new int[GLCM_CLASS * GLCM_CLASS];  
18.     int * histImage = new int[width * height];  
19.   
20.     if(NULL == glcm || NULL == histImage)  
21.         return 2;  
22.   
23.     //灰度等級化---分GLCM_CLASS個等級  
24.     uchar *data =(uchar*) bWavelet->imageData;  
25.     for(i = 0;i < height;i++){  
26.         for(j = 0;j < width;j++){  
27.             histImage[i * width + j] = (int)(data[bWavelet->widthStep * i + j] * GLCM_CLASS / 256);  
28.         }  
29.     }  
30.   
31.     //初始化共生矩陣  
32.     for (i = 0;i < GLCM_CLASS;i++)  
33.         for (j = 0;j < GLCM_CLASS;j++)  
34.             glcm[i * GLCM_CLASS + j] = 0;  
35.   
36.     //計算灰度共生矩陣  
37.     int w,k,l;  
38.     //水平方向  
39.     if(angleDirection == GLCM_ANGLE_HORIZATION)  
40.     {  
41.         for (i = 0;i < height;i++)  
42.         {  
43.             for (j = 0;j < width;j++)  
44.             {  
45.                 l = histImage[i * width + j];  
46.                 if(j + GLCM_DIS >= 0 && j + GLCM_DIS < width)  
47.                 {  
48.                     k = histImage[i * width + j + GLCM_DIS];  
49.                     glcm[l * GLCM_CLASS + k]++;  
50.                 }  
51.                 if(j - GLCM_DIS >= 0 && j - GLCM_DIS < width)  
52.                 {  
53.                     k = histImage[i * width + j - GLCM_DIS];  
54.                     glcm[l * GLCM_CLASS + k]++;  
55.                 }  
56.             }  
57.         }  
58.     }  
59.     //垂直方向  
60.     else if(angleDirection == GLCM_ANGLE_VERTICAL)  
61.     {  
62.         for (i = 0;i < height;i++)  
63.         {  
64.             for (j = 0;j < width;j++)  
65.             {  
66.                 l = histImage[i * width + j];  
67.                 if(i + GLCM_DIS >= 0 && i + GLCM_DIS < height)   
68.                 {  
69.                     k = histImage[(i + GLCM_DIS) * width + j];  
70.                     glcm[l * GLCM_CLASS + k]++;  
71.                 }  
72.                 if(i - GLCM_DIS >= 0 && i - GLCM_DIS < height)   
73.                 {  
74.                     k = histImage[(i - GLCM_DIS) * width + j];  
75.                     glcm[l * GLCM_CLASS + k]++;  
76.                 }  
77.             }  
78.         }  
79.     }  
80.     //對角方向  
81.     else if(angleDirection == GLCM_ANGLE_DIGONAL)  
82.     {  
83.         for (i = 0;i < height;i++)  
84.         {  
85.             for (j = 0;j < width;j++)  
86.             {  
87.                 l = histImage[i * width + j];  
88.   
89.                 if(j + GLCM_DIS >= 0 && j + GLCM_DIS < width && i + GLCM_DIS >= 0 && i + GLCM_DIS < height)  
90.                 {  
91.                     k = histImage[(i + GLCM_DIS) * width + j + GLCM_DIS];  
92.                     glcm[l * GLCM_CLASS + k]++;  
93.                 }  
94.                 if(j - GLCM_DIS >= 0 && j - GLCM_DIS < width && i - GLCM_DIS >= 0 && i - GLCM_DIS < height)  
95.                 {  
96.                     k = histImage[(i - GLCM_DIS) * width + j - GLCM_DIS];  
97.                     glcm[l * GLCM_CLASS + k]++;  
98.                 }  
99.             }  
100.         }  
101.     }  
102.   
103.     //計算特征值  
104.     double entropy = 0,energy = 0,contrast = 0,homogenity = 0;  
105.     for (i = 0;i < GLCM_CLASS;i++)  
106.     {  
107.         for (j = 0;j < GLCM_CLASS;j++)  
108.         {  
109.             //熵  
110.             if(glcm[i * GLCM_CLASS + j] > 0)  
111.                 entropy -= glcm[i * GLCM_CLASS + j] * log10(double(glcm[i * GLCM_CLASS + j]));  
112.             //能量  
113.             energy += glcm[i * GLCM_CLASS + j] * glcm[i * GLCM_CLASS + j];  
114.             //對比度  
115.             contrast += (i - j) * (i - j) * glcm[i * GLCM_CLASS + j];  
116.             //一致性  
117.             homogenity += 1.0 / (1 + (i - j) * (i - j)) * glcm[i * GLCM_CLASS + j];  
118.         }  
119.     }  
120.     //返回特征值  
121.     i = 0;  
122.     featureVector[i++] = entropy;  
123.     featureVector[i++] = energy;  
124.     featureVector[i++] = contrast;  
125.     featureVector[i++] = homogenity;  
126.   
127.     delete[] glcm;  
128.     delete[] histImage;  
129.     return 0;  
130. }