1 前⾔
纹理特征在地物光谱特征⽐较相似的时候常作为⼀种特征⽤于图像的分类和信息提取。本⽂主要介绍基于灰度共⽣矩阵的图像纹理提取。
2 灰度共⽣矩阵
纹理⽯油灰度分布在空间位置上反复出现⽽形成的，因⽽图像空间中相隔某距离的两个像素之间存在⼀定的灰度关系，即图像中灰度的空间相关特性。灰度共⽣矩阵是⼀种通过研究灰度的空间相关特性来描述纹理的常⽤⽅法。 有⼀个⽹站提供了⾮常详细的关于灰度共⽣矩阵的介绍，我在这⾥就不赘述了。⽹址如下：
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
6class CCo_occurrenceTextureExtration
7 {
8public:
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//影像灰度降级（包括直⽅图均衡化）
37bool Init();
38bool ReduceGrayLevel();
39bool ExtraTextureWinEXE();
40bool ExtraTexturePixEXE();
41int ExtraTexturePixOpenCLEXE();
42private:
43char* LoadProgSource(const char* cFilename, const char* cPreamble, size_t* szFinalLength);
44bool GetGrayCoocurrenceMatrix(std::vector<std::vector<int>>grayValue,
45 std::vector<std::vector<float>> &grayMatrix,int &m_count);
46float CalGLCMStatistics(std::vector<std::vector<float>> grayMatrix,
47const int count);
48double* cdf(int *h,int length);
49
50public:
51 std::string m_strInFileName;
52 std::string m_strOutFileName;
53 std::string m_strReduceLevelFileName;
54int m_iWidth;
55int m_iHeigth;
56int m_iBandCount;
57int *m_BandMap;
58int m_AngleMode;
59int m_TextureMode;
60int m_dis;
61int m_grayLevel;
62int m_winWidth;
63int 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 {
17if(!m_BandMap){
18delete []m_BandMap;
19 }
20 }
21
22bool 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];
31int i = 0;
32while(i < m_iBandCount){
33 m_BandMap[i] = i+1;
34 i++;
35 }
36 GDALClose((GDALDatasetH)poDS);
37return TRUE;
38 }
39
40bool CCo_occurrenceTextureExtration::ReduceGrayLevel()
41 {
42//直⽅图均衡化
43///// 统计直⽅图
44 GDALAllRegister();
45 GDALDataset *poDS = (GDALDataset*)GDALOpen(m_strInFileName.c_str(),GA_ReadOnly);
46double 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());
55for(int i = 0;i<m_iBandCount;i++){
56double *bandMINMAX = new double[2];
57 GDALRasterBand *poBand = poDS->GetRasterBand(m_BandMap[i]);
58 GDALRasterBand *poNewBand = poOutDS->GetRasterBand(i+1);
59 poBand->ComputeRasterMinMax(FALSE,bandMINMAX);
60for(int j = 0;j<m_iHeigth;j++){
61double *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);
64int k = 0;
65while(k<m_iWidth){
66int tmp = 0;
67if(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);
76delete []readBuff;readBuff = NULL;
77delete []writeBuff;writeBuff = NULL;
78 }
79delete []bandMINMAX;bandMINMAX = NULL;
80 }
81 GDALClose((GDALDatasetH) poDS);
82 GDALClose((GDALDatasetH) poOutDS);
83return TRUE;
84 }
85
86double* CCo_occurrenceTextureExtration::cdf( int *h,int length )
我是电影通87 {
88int n = 0;
89for( int i = 0; i < length - 1; i++ )
90 {
91 n += h[i];
92 }
93double* p = new double[length];
94int c = h[0];
95 p[0] = ( double )c / n;
96for( int i = 1; i < length - 1; i++ )
97 {
98 c += h[i];
99 p[i] = ( double )c / n;
100 }
101return p;
102 }
103
104char* 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");
111if(pFileStream == 0)
112 {
113return 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
124char* cSourceString = (char *)malloc(szSourceLength + szPreambleLength + 1);
125 memcpy(cSourceString, cPreamble, szPreambleLength);
126if (fread((cSourceString) + szPreambleLength, szSourceLength, 1, pFileStream) != 1)
127 {
128 fclose(pFileStream);
129free(cSourceString);
130return0;
131 }
132
133// close the file and return the total length of the combined (preamble + source) string
134 fclose(pFileStream);
135if(szFinalLength != 0)
136 {
137 *szFinalLength = szSourceLength + szPreambleLength;
138 }
139 cSourceString[szSourceLength + szPreambleLength] = '\0';
140
141return cSourceString;
142 }
143
144bool CCo_occurrenceTextureExtration::GetGrayCoocurrenceMatrix(std::vector<std::vector<int>>grayValue,
145 std::vector<std::vector<float>> &m_grayMatrix,
146int &m_count)
147 {
148int i,j,r,c;
149 m_count = 0;
150int ii = 0;
151六盘水马拉松
152switch(m_AngleMode)
153 {
154case DEGREE_0:
155 i = 0;
156while(i<m_winHeigth){
157 j = 0;
158while(j<m_winWidth){
159int endR = i;
160int endC = j + m_dis;
161if(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 }
172break;
173case DEGREE_45:
174 i = 0;
175while(i<m_winHeigth){
176 j = 0;
177while(j<m_winWidth){
178int endR = i - m_dis;
179int endC = j + m_dis;
180if(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 }
191break;
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议