基于灰度共生矩阵的纹理提取

基于灰度共生矩阵的纹理提取

 

1 前言

纹理特征在地物光谱特征比较相似的时候常作为一种特征用于图像的分类和信息提取。本文主要介绍基于灰度共生矩阵的图像纹理提取。

2 灰度共生矩阵

纹理石油灰度分布在空间位置上反复出现而形成的，因而图像空间中相隔某距离的两个像素之间存在一定的灰度关系，即图像中灰度的空间相关特性。灰度共生矩阵是一种通过研究灰度的空间相关特性来描述纹理的常用方法。

有一个网站提供了非常详细的关于灰度共生矩阵的介绍，我在这里就不赘述了。网址如下：http://www.fp.ucalgary.ca/mhallbey/texture_calculations.htm

3 纹理提取的算法流程

本文重点说明一下遥感影像纹理提取的算法流程，如下图1所示：

 

图 1 纹理提取算法流程

具体描述如下：

1)        灰度降级，对原始影像进行灰度降级如8,16,32,64等；

纹理计算的灰度降级策略来源于IDL的bytscl函数介绍，具体描述如下：

 

图 2 灰度降级

2)        根据设定好的窗口大小，逐窗口计算灰度共生矩阵；

3)        根据选择的二阶统计量，计算纹理值。

4  纹理算子

协同性（GLCM_HOM）：对应ENVI的Homogeneity

 

反差性（GLCM_CON）：

 

非相似性（GLCM_DIS）：

 

均值GLCM_MEAN：对应ENVI的Mean

 

方差GLCM_VAR：对应ENVI的Variance

 

角二阶矩GLCM_ASM：对应ENVI的Second Moment

 

相关性GLCM_COR：对应ENVI的Correlation

 

GLDV角二阶矩GLDV_ASM：

 

熵GLCM_ENTROPY：对应ENVI的Entropy

 

归一化灰度矢量均值GLDV_MEAN：对应ENVI的Dissimilarity

 

归一化角二阶矩GLDV_CON：对应ENVI的Contrast

 

5  代码实现

     同样实现了CPU和GPU两个版本。但是目前效率较低，要是有哪位大神可以指点下如何提高纹理计算的效率感激不尽~

先上传头文件和源文件

 1 #pragma once
 2 #include <iostream>
 3 #include <string>
 4 #include <vector>
 5 
 6 class CCo_occurrenceTextureExtration
 7 {
 8 public:
 9     CCo_occurrenceTextureExtration(void);
10     ~CCo_occurrenceTextureExtration(void);
11 
12 
13     typedef enum ANGLE_TYPE
14     {
15         DEGREE_0,
16         DEGREE_45,
17         DEGREE_90,
18         DEGREE_135
19     };
20     typedef enum LIST_TYPE
21     {
22         GLCM_HOM,
23         GLCM_CON, 
24         GLCM_DIS,
25         GLCM_MEAN, 
26         GLCM_VAR, 
27         GLCM_ASM,
28         GLCM_COR,
29         GLDV_ASM,  
30         GLCM_ENTROPY,
31         GLDV_MEAN, 
32         GLDV_CON,
33     };
34 
35 
36     //影像灰度降级（包括直方图均衡化）
37     bool Init();
38     bool ReduceGrayLevel();
39     bool ExtraTextureWinEXE();
40     bool ExtraTexturePixEXE();
41     int ExtraTexturePixOpenCLEXE();
42 private:
43     char* LoadProgSource(const char* cFilename, const char* cPreamble, size_t* szFinalLength);
44     bool GetGrayCoocurrenceMatrix(std::vector<std::vector<int>>grayValue,
45         std::vector<std::vector<float>> &grayMatrix,int &m_count);
46     float CalGLCMStatistics(std::vector<std::vector<float>> grayMatrix,
47         const int count);
48     double* cdf(int *h,int length);
49 
50 public:
51     std::string m_strInFileName;
52     std::string m_strOutFileName;
53     std::string m_strReduceLevelFileName;
54     int m_iWidth;
55     int m_iHeigth;
56     int m_iBandCount;
57     int *m_BandMap;
58     int m_AngleMode;
59     int m_TextureMode;
60     int m_dis;
61     int m_grayLevel;
62     int m_winWidth;
63     int m_winHeigth;
64     
65 };

源文件如下：

  1 #include "Co_occurrenceTextureExtration.h"
  2 #include <gdal_alg_priv.h>
  3 #include <gdal_priv.h>
  4 #include <math.h>
  5 #include <omp.h>
  6 #include <CL/cl.h>
  7 #include <stdlib.h>
  8 #include <malloc.h>
  9 #pragma comment(lib,"OpenCL.lib")
 10 
 11 CCo_occurrenceTextureExtration::CCo_occurrenceTextureExtration( )
 12 {
 13     m_BandMap = NULL;
 14 }
 15 CCo_occurrenceTextureExtration::~CCo_occurrenceTextureExtration(void)
 16 {
 17     if(!m_BandMap){
 18         delete []m_BandMap;
 19     }
 20 }
 21 
 22 bool CCo_occurrenceTextureExtration::Init()
 23 {
 24     GDALAllRegister();
 25     GDALDataset *poDS = (GDALDataset*)GDALOpen(m_strInFileName.c_str(),GA_ReadOnly);
 26     m_iWidth = poDS->GetRasterXSize();
 27     m_iHeigth = poDS->GetRasterYSize();
 28     m_iBandCount = poDS->GetRasterCount();
 29 
 30     m_BandMap = new int[m_iBandCount];
 31     int i = 0;
 32     while(i < m_iBandCount){
 33         m_BandMap[i] = i+1;
 34         i++;
 35     }
 36     GDALClose((GDALDatasetH)poDS);
 37     return TRUE;
 38 }
 39 
 40 bool CCo_occurrenceTextureExtration::ReduceGrayLevel()
 41 {
 42     //直方图均衡化
 43     ///// 统计直方图
 44     GDALAllRegister();
 45     GDALDataset *poDS = (GDALDataset*)GDALOpen(m_strInFileName.c_str(),GA_ReadOnly);
 46     double adfMSSGeoTransform[6] = {
      
      0};
 47     poDS->GetGeoTransform(adfMSSGeoTransform);
 48     GDALDataType eDT = poDS->GetRasterBand(1)->GetRasterDataType();
 49     //创建文件
 50     GDALDriver *poOutDriver = (GDALDriver*)GDALGetDriverByName("GTIFF");
 51     GDALDataset *poOutDS  = poOutDriver->Create(m_strReduceLevelFileName.c_str(),
 52         m_iWidth,m_iHeigth,m_iBandCount,GDT_Byte,NULL);
 53     poOutDS->SetGeoTransform(adfMSSGeoTransform);
 54     poOutDS->SetProjection(poDS->GetProjectionRef());
 55     for(int i = 0;i<m_iBandCount;i++){
 56         double *bandMINMAX = new double[2];
 57         GDALRasterBand *poBand = poDS->GetRasterBand(m_BandMap[i]);
 58         GDALRasterBand *poNewBand = poOutDS->GetRasterBand(i+1);
 59         poBand->ComputeRasterMinMax(FALSE,bandMINMAX);
 60         for(int j = 0;j<m_iHeigth;j++){
 61             double *readBuff = new double[m_iWidth];
 62             unsigned char*writeBuff = new unsigned char[m_iWidth];
 63             poBand->RasterIO(GF_Read,0,j,m_iWidth,1,readBuff,m_iWidth,1,GDT_Float64,0,0);
 64             int k = 0;
 65             while(k<m_iWidth){
 66                 int tmp = 0;
 67                 if(eDT == GDT_Float32 || eDT == GDT_Float64){
 68                     tmp = int((m_grayLevel-1+0.99999)*(readBuff[k] - bandMINMAX[0])/(bandMINMAX[1] - bandMINMAX[0]));
 69                 }else{
 70                     tmp = int((m_grayLevel)*(readBuff[k] - bandMINMAX[0] - 1)/(bandMINMAX[1] - bandMINMAX[0]));
 71                 }
 72                 writeBuff[k] = unsigned char(tmp);
 73                 k++;
 74             }
 75             CPLErr str = poNewBand->RasterIO(GF_Write,0,j,m_iWidth,1,writeBuff,m_iWidth,1,GDT_Byte,0,0);
 76             delete []readBuff;readBuff = NULL;
 77             delete []writeBuff;writeBuff = NULL;
 78         }
 79         delete []bandMINMAX;bandMINMAX = NULL;
 80     }
 81     GDALClose((GDALDatasetH) poDS);
 82     GDALClose((GDALDatasetH) poOutDS);
 83     return TRUE;
 84 }
 85 
 86 double* CCo_occurrenceTextureExtration::cdf( int *h,int length )
 87 {
 88     int n = 0;  
 89     for( int i = 0; i < length - 1; i++ )  
 90     {  
 91         n += h[i];  
 92     }  
 93     double* p = new double[length];  
 94     int c = h[0];  
 95     p[0] = ( double )c / n;  
 96     for( int i = 1; i < length - 1; i++ )  
 97     {  
 98         c += h[i];  
 99         p[i] = ( double )c / n;  
100     }  
101     return p; 
102 }
103 
104 char* CCo_occurrenceTextureExtration::LoadProgSource( const char* cFilename, const char* cPreamble, size_t* szFinalLength )
105 {
106     FILE* pFileStream = NULL;
107     size_t szSourceLength;
108 
109     // open the OpenCL source code file
110     pFileStream = fopen(cFilename, "rb");
111     if(pFileStream == 0) 
112     {     
113         return NULL;
114     }
115 
116     size_t szPreambleLength = strlen(cPreamble);
117 
118     // get the length of the source code
119     fseek(pFileStream, 0, SEEK_END); 
120     szSourceLength = ftell(pFileStream);
121     fseek(pFileStream, 0, SEEK_SET); 
122 
123     // allocate a buffer for the source code string and read it in
124     char* cSourceString = (char *)malloc(szSourceLength + szPreambleLength + 1); 
125     memcpy(cSourceString, cPreamble, szPreambleLength);
126     if (fread((cSourceString) + szPreambleLength, szSourceLength, 1, pFileStream) != 1)
127     {
128         fclose(pFileStream);
129         free(cSourceString);
130         return 0;
131     }
132 
133     // close the file and return the total length of the combined (preamble + source) string
134     fclose(pFileStream);
135     if(szFinalLength != 0)
136     {
137         *szFinalLength = szSourceLength + szPreambleLength;
138     }
139     cSourceString[szSourceLength + szPreambleLength] = '\0';
140 
141     return cSourceString;
142 }
143 
144 bool CCo_occurrenceTextureExtration::GetGrayCoocurrenceMatrix(std::vector<std::vector<int>>grayValue,
145                                                               std::vector<std::vector<float>> &m_grayMatrix,
146                                                               int &m_count)
147 {
148     int i,j,r,c;
149     m_count = 0;
150     int ii = 0;
151 
152     switch(m_AngleMode)
153     {
154     case DEGREE_0:
155         i = 0;
156         while(i<m_winHeigth){
157             j = 0;
158             while(j<m_winWidth){
159                 int endR = i;
160                 int endC = j + m_dis;
161                 if(endC < m_winWidth){
162                     r = grayValue[i][j];
163                     c = grayValue[endR][endC];
164                     m_grayMatrix[r][c] += 1;
165                     m_grayMatrix[c][r] += 1;
166                     m_count = m_count+2;
167                 }
168                 j++;
169             }
170             i++;
171         }
172         break;
173     case DEGREE_45:
174         i = 0;
175         while(i<m_winHeigth){
176             j = 0;
177             while(j<m_winWidth){
178                 int endR = i - m_dis;
179                 int endC = j + m_dis;
180                 if(endR>=0 && endR < m_winHeigth && endC >=0 && endC < m_winWidth){
181                     r = grayValue[i][j];
182                     c = grayValue[endR][endC];
183                     m_grayMatrix[r][c] += 1;
184                     m_grayMatrix[c][r] += 1;
185                     m_count = m_count+2;
186                 }
187                 j++;
188             }
189             i++;
190         }
191         break;
192     case DEGREE_90:
193         i= 0;
194