提取共生矩阵特征

　　 共生矩阵用两个位置的象素的联合概率密度来定义，它不仅反映亮度的分布特性，也反映具有同样亮度或接近亮度的象素之间的位置分布特性，是有关图象亮度变化的二阶统计特征。它是定义一组纹理特征的基础。

　　 灰度共生矩阵能反映出图象灰度关于方向、相邻间隔、变化幅度的综合信息。设f(x,y)为一幅二维数字图象，其大小为M×N，灰度级别为Ng,则满足一定空间关系的灰度共生矩阵为：

　　 P(i,j)=#｛(x1,y1),(x2,y2)∈M×N｜f(x1,y1)=i,f(x2,y2)=j｝

　　 其中#(x)表示集合x中的元素个数，显然P为Ng×Ng的矩阵，若(x1,y1)与(x2,y2)间距离为d,两者与坐标横轴的夹角为θ，则可以得到各种间距及角度的灰度共生矩阵P(i,j,d,θ)。

1. class CGlcm  
2. {  
3. public:  
4.     CGlcm(void);  
5.     virtual ~CGlcm(void);  
6. public:  
7.     int       *m_pMat1;  
8.     int       *m_pMat2;  
9.     int       *m_pMat3;  
10.     int       *m_pMat4;  
11. protected:  
12.     int       *m_pLine[4][256];  
13.     int        m_nMin;  
14.     int        m_nMax;  
15.     int        m_nSum;  
16. public:  
17.     // 计算共生矩阵   
18.     // 参数:   
19.     // 1. pImageData: 图像数据指针，单通道，8位。   
20.     // 2. nLeft,nTop,nWidth,nHeight: 计算共生矩阵的区域。   
21.     // 3. nWidthStep: 行偏移量。   
22.     // 4. nScale: 尺度。   
23.     // 5. nReduction: 灰度级压缩。   
24.     bool CalGlcm(unsigned char *pImageData, int nLeft, int nTop, int nWidth, int nHeight, int nWidthStep, int nScale, int nReduction);  
25. };  
26. // 灰度共生矩阵   
27. CGlcm::CGlcm(void)  
28. : m_pMat1(NULL)  
29. , m_pMat2(NULL)  
30. , m_pMat3(NULL)  
31. , m_pMat4(NULL)  
32. , m_nMin(0)  
33. , m_nMax(0)  
34. , m_nSum(0)  
35. {  
36.     int i, j;  
37.     unsigned int nSize = sizeof(int) * 256 * 256;  
38.     // 创建共生矩阵   
39.     m_pMat1 = (int*) malloc(nSize);  
40.     m_pMat2 = (int*) malloc(nSize);  
41.     m_pMat3 = (int*) malloc(nSize);  
42.     m_pMat4 = (int*) malloc(nSize);  
43.     if (m_pMat1 && m_pMat2 && m_pMat3 && m_pMat4)  
44.     {  
45.         for (i = 0, j = 0; i < 256; i++, j += 256)  
46.         {  
47.             m_pLine[0][i] = m_pMat1 + j;  
48.             m_pLine[1][i] = m_pMat2 + j;  
49.             m_pLine[2][i] = m_pMat3 + j;  
50.             m_pLine[3][i] = m_pMat4 + j;  
51.         }  
52.     }  
53. }  
54. CGlcm::~CGlcm(void)  
55. {  
56.     // 释放共生矩阵   
57.     if (m_pMat1) free(m_pMat1);  
58.     if (m_pMat2) free(m_pMat2);  
59.     if (m_pMat3) free(m_pMat3);  
60.     if (m_pMat4) free(m_pMat4);  
61. }  
62. // 计算共生矩阵   
63. bool CGlcm::CalGlcm(unsigned char *pImageData, int nLeft, int nTop, int nWidth, int nHeight, int nWidthStep, int nScale, int nReduction)  
64. {  
65.     bool bResult = false;  
66.     int x0, x1, x2;  
67.     int y0, y1, y2;  
68.     int nGray;  
69.     unsigned int nSize;  
70.     unsigned char *pLine[3];  
71.     if (pImageData)  
72.     {  
73.         // 灰度最值   
74.         m_nMin = 0xFF;  
75.         m_nMax = 0;  
76.         pLine[1] = pImageData + nWidthStep * nTop + nLeft;  
77.         for (y1 = 0; y1 < nHeight; y1++)  
78.         {  
79.             for (x1 = 0; x1 < nWidth; x1++)  
80.             {  
81.                 // 灰度级压缩   
82.                 if (nReduction > 0)  
83.                 {  
84.                     pLine[1][x1] = pLine[1][x1] >> nReduction;  
85.                 }  
86.                 nGray = pLine[1][x1];  
87.                 if (nGray < m_nMin)  
88.                 {  
89.                     m_nMin = nGray;  
90.                 }  
91.                 else if (nGray > m_nMax)  
92.                 {  
93.                     m_nMax = nGray;  
94.                 }  
95.             }  
96.             pLine[1] += nWidthStep;  
97.         }  
98.         // 累加和   
99.         m_nSum = nWidth * nHeight * 2;  
100.         if (m_nMax >= m_nMin)  
101.         {  
102.             // 清空内存   
103.             nSize = sizeof(int) * (m_nMax - m_nMin + 1);  
104.             for (y1 = m_nMin; y1 <= m_nMax; y1++)  
105.             {  
106.                 memset(&m_pLine[0][y1][m_nMin], 0, nSize);  
107.                 memset(&m_pLine[1][y1][m_nMin], 0, nSize);  
108.                 memset(&m_pLine[2][y1][m_nMin], 0, nSize);  
109.                 memset(&m_pLine[3][y1][m_nMin], 0, nSize);  
110.             }  
111.             // 计算共生矩阵   
112.             pLine[1] = pImageData + nWidthStep * nTop + nLeft;  
113.             for (y0 = -nScale, y1 = 0, y2 = nScale; y1 < nHeight; y0++, y1++, y2++)  
114.             {  
115.                 pLine[0] = pLine[1] - nWidthStep;  
116.                 pLine[2] = pLine[1] + nWidthStep;  
117.                 for (x0 = -nScale, x1 = 0, x2 = nScale; x1 < nWidth; x0++, x1++, x2++)  
118.                 {  
119.                     nGray = pLine[1][x1];  
120.                     // 0 度   
121.                     if (x2 < nWidth)  
122.                     {  
123.                         m_pLine[0][nGray][pLine[1][x2]]++;  
124.                     }  
125.                     // 45 度   
126.                     if (y0 > 0 && x2 < nWidth)  
127.                     {  
128.                         m_pLine[1][nGray][pLine[0][x2]]++;  
129.                     }  
130.                     // 90 度   
131.                     if (y0 > 0)  
132.                     {  
133.                         m_pLine[2][nGray][pLine[0][x1]]++;  
134.                     }  
135.                     // 135 度   
136.                     if (y0 > 0 && x0 > 0)  
137.                     {  
138.                         m_pLine[3][nGray][pLine[0][x0]]++;  
139.                     }  
140.                     // 180 度   
141.                     if (x0 > 0)  
142.                     {  
143.                         m_pLine[0][nGray][pLine[1][x0]]++;  
144.                     }  
145.                     // 225 度   
146.                     if (x0 > 0 && y2 < nHeight)  
147.                     {  
148.                         m_pLine[1][nGray][pLine[2][x0]]++;  
149.                     }  
150.                     // 270 度   
151.                     if (y2 < nHeight)  
152.                     {  
153.                         m_pLine[2][nGray][pLine[2][x1]]++;  
154.                     }  
155.                     // 315 度   
156.                     if (y2 < nHeight && x2 < nWidth)  
157.                     {  
158.                         m_pLine[3][nGray][pLine[2][x2]]++;  
159.                     }  
160.                 }  
161.                 pLine[1] += nWidthStep;  
162.             }  
163.             bResult = true;  
164.         }  
165.     }  
166.     return bResult;  
167. }  

1. #include <math.h>   
2. // 共生矩阵特征   
3. struct GlcmFeature  
4. {  
5.     double ASM[4];  // 角二阶矩/能量   
6.     double ENT[4];  // 熵   
7.     double CON[4];  // 对比度   
8.     double HOM[4];  // 逆差矩/同质性   
9.     double COR[4];  // 相关性   
10. };  
11. // 灰度共生矩阵   
12. class CGlcm  
13. {  
14. public:  
15.     CGlcm(void);  
16.     virtual ~CGlcm(void);  
17. public:  
18.     int       *m_pMat1;  
19.     int       *m_pMat2;  
20.     int       *m_pMat3;  
21.     int       *m_pMat4;  
22. protected:  
23.     int       *m_pLine[4][256];  
24.     int        m_nMin;  
25.     int        m_nMax;  
26.     int        m_nSum;  
27. public:  
28.     // 计算共生矩阵   
29.     // 参数:   
30.     // 1. pImageData: 图像数据指针，单通道，8位。   
31.     // 2. nLeft,nTop,nWidth,nHeight: 计算共生矩阵的区域。   
32.     // 3. nWidthStep: 行偏移量。   
33.     // 4. nScale: 尺度。   
34.     // 5. nReduction: 灰度级压缩。   
35.     bool CalGlcm(unsigned char *pImageData, int nLeft, int nTop, int nWidth, int nHeight, int nWidthStep, int nScale, int nReduction);  
36.     // 共生矩阵特征   
37.     // 参数:   
38.     // 1. Feature: 输出共生矩阵特征。   
39.     void GetFeature(GlcmFeature& Feature);  
40. };  
41. // 共生矩阵特征   
42. void CGlcm::GetFeature(GlcmFeature& Feature)  
43. {  
44.     int x, y;  
45.     int nTheta;  
46.     int nValue;  
47.     int nTemp;  
48.     double dValue;  
49.     double dMean, dStdDev;  
50.     double dSum[256];  
51.     // 清空内存   
52.     memset(&Feature, 0, sizeof(Feature));  
53.     // 方向循环   
54.     for (nTheta = 0; nTheta < 4; nTheta++)  
55.     {  
56.         dMean = 0;  
57.         dStdDev = 0;  
58.         // 清空内存   
59.         memset(dSum, 0, sizeof(dSum));  
60.         for (y = m_nMin; y <= m_nMax; y++)  
61.         {  
62.             for (x = m_nMin; x <= m_nMax; x++)  
63.             {  
64.                 nValue = m_pLine[nTheta][y][x];  
65.                 if (nValue != 0)  
66.                 {  
67.                     // 归一化共生矩阵   
68.                     dValue = (double) nValue / (double) m_nSum;  
69.                     nTemp = (x - y) * (x - y);  
70.                     // 角二阶矩/能量   
71.                     Feature.ASM[nTheta] += (dValue * dValue);  
72.                     // 熵   
73.                     Feature.ENT[nTheta] -= (dValue * log(dValue));  
74.                     // 对比度   
75.                     Feature.CON[nTheta] += (nTemp * dValue);  
76.                     // 逆差矩/同质性   
77.                     Feature.HOM[nTheta] += (dValue / (1 + nTemp));  
78.                     // 相关性   
79.                     Feature.COR[nTheta] += (x * y * dValue);  
80.                     dSum[y] += dValue;  
81.                 }  
82.             }  
83.         }  
84.         for (y = m_nMin; y <= m_nMax; y++)  
85.         {  
86.             dMean += (y * dSum[y]);  
87.         }  
88.         for (y = m_nMin; y <= m_nMax; y++)  
89.         {  
90.             dStdDev += ((y - dMean) * (y - dMean) * dSum[y]);  
91.         }  
92.         // 相关性   
93.         if (abs(dStdDev) > 1e-15)  
94.         {  
95.             Feature.COR[nTheta] = (Feature.COR[nTheta] - dMean * dMean) / dStdDev;  
96.         }  
97.         else  
98.         {  
99.             Feature.COR[nTheta] = 0;  
100.         }  
101.     }