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] = '