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         while(i<m_winHeigth){
195             j = 0;
196             while(j < m_winWidth){
197                 int endR = i - m_dis;
198                 int endC = j;
199                 if(endR >= 0){
200                     r = grayValue[i][j];
201                     c = grayValue[endR][endC];
202                     m_grayMatrix[r][c] += 1;
203                     m_grayMatrix[c][r] += 1;
204                     m_count = m_count+2;
205                 }
206                 j++;
207             }
208             i++;
209         }
210         break;
211     case DEGREE_135:
212         i = 0;
213         while(i<m_winHeigth){
214             j = 0;
215             while(j<m_winWidth){
216                 int endR = i - m_dis;
217                 int endC = j - m_dis;
218                 if(endR >= 0 && endC >= 0){
219                     r = grayValue[i][j];
220                     c = grayValue[endR][endC];
221                     m_grayMatrix[r][c] += 1;
222                     m_grayMatrix[c][r] += 1;
223                     m_count = m_count+2;
224                 }
225                 j++;
226             }
227             i++;
228         }
229         break;
230     default:
231         return FALSE;
232     }
233     return TRUE;
234 }
235 
236 float CCo_occurrenceTextureExtration::CalGLCMStatistics(std::vector<std::vector<float>> m_grayMatrix,
237                                                         const int m_count)
238 {
239     int i,j,k;
240     float imean = 0.0f;
241     float jmean = 0.0f;
242     float ivar = 0.0f;
243     float jvar = 0.0f;;
244     float res = 0.0f;
245     switch(m_TextureMode)
246     {
247     case GLCM_HOM:  //同质性
248         i = 0;
249         while(i < m_grayLevel){
250             j = 0;
251             while(j < m_grayLevel){
252                 if(m_grayMatrix[i][j] != 0){
253                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
254                     res += Pij/(1 + (i-j)*(i-j));
255                 }
256                 j++;
257             }
258             i++;
259         }
260         break;
261     case GLCM_CON:  //对比度
262         i = 0;
263         while(i < m_grayLevel){
264             j = 0;
265             while(j < m_grayLevel){
266                 if(m_grayMatrix[i][j] != 0){
267                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
268                     res += Pij*(i-j)*(i-j);
269                 }
270                 j++;
271             }
272             i++;
273         }
274         break;
275     case GLCM_DIS:  //非相似性
276         i = 0;
277         while(i < m_grayLevel){
278             j = 0;
279             while(j < m_grayLevel){
280                 if(m_grayMatrix[i][j] != 0){
281                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
282                     res += Pij*std::abs(i-j);
283                 }
284                 j++;
285             }
286             i++;
287         }
288         break;
289     case GLCM_COR:  // 相关性
290         i = 0;
291         while(i < m_grayLevel){
292             j = 0;
293             while(j < m_grayLevel){
294                 if(m_grayMatrix[i][j] != 0){
295                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
296                     imean += Pij * i;
297                     jmean += Pij * j;
298                 }
299                 j++;
300             }
301             i++;
302         }
303         //  var
304         i = 0;
305         while(i < m_grayLevel){
306             j = 0;
307             while(j < m_grayLevel){
308                 if(m_grayMatrix[i][j] != 0){
309                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
310                     ivar += Pij * (i - imean)*(i - imean);
311                     jvar += Pij * (j - imean)*(j - imean);
312                 }
313                 j++;
314             }
315             i++;
316         }
317         // xiangguanx
318         i = 0;
319         while(i < m_grayLevel){
320             j = 0;
321             while(j < m_grayLevel){
322                 if(m_grayMatrix[i][j] != 0){
323                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
324                     res += (Pij *(i - imean)*(j - jmean))/sqrt(ivar*jvar);
325                 }
326                 j++;
327             }
328             i++;
329         }
330         break;
331     case GLCM_MEAN:  // 均值
332         i = 0;
333         while(i < m_grayLevel){
334             j = 0;
335             while(j < m_grayLevel){
336                 if(m_grayMatrix[i][j] != 0){
337                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
338                     imean += Pij * i;
339                     //jmean += Pij * j;
340                 }
341                 j++;
342             }
343             i++;
344         }
345         res = imean;
346         break;
347     case GLCM_VAR:  // 方差
348         //  mean
349         i = 0;
350         while(i < m_grayLevel){
351             j = 0;
352             while(j < m_grayLevel){
353                 if(m_grayMatrix[i][j] != 0){
354                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
355                     imean += Pij * i;
356                 }
357                 j++;
358             }
359             i++;
360         }
361         //  var
362         i = 0;
363         while(i < m_grayLevel){
364             j = 0;
365             while(j < m_grayLevel){
366                 if(m_grayMatrix[i][j] != 0){
367                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
368                     ivar += Pij * (i - imean)*(i - imean);
369                 }
370                 j++;
371             }
372             i++;
373         }
374         //res = (ivar + jvar)/2;
375         res = ivar;
376         break;
377     case GLCM_ASM:  // 角二阶矩
378         i = 0;
379         while(i < m_grayLevel){
380             j = 0;
381             while(j < m_grayLevel){
382                 if(m_grayMatrix[i][j] != 0){
383                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
384                     res += Pij*Pij;
385                 }
386                 j++;
387             }
388             i++;
389         }
390         break;
391     case GLCM_ENTROPY:  // 熵
392         i = 0;
393         while(i < m_grayLevel){
394             j = 0;
395             while(j < m_grayLevel){
396                 if(m_grayMatrix[i][j] != 0){
397                     float Pij = m_grayMatrix[i][j]*1.0f/m_count;
398                     res += -1.0f*Pij * log(Pij);
399                 }
400                 j++;
401             }
402             i++;
403         }
404         break;
405     case GLDV_MEAN:
406         k = 0;
407         while(k<m_grayLevel){
408             i = 0;
409             while(i < m_grayLevel){
410                 j = 0;
411                 while(j < m_grayLevel){
412                     if(m_grayMatrix[i][j] != 0 && std::abs(i-j)== k){
413                         float Pij = m_grayMatrix[i][j]*1.0f/m_count;
414                         res += Pij*k;
415                     }
416                     j++;
417                 }
418                 i++;
419             }
420             k++;
421         }
422         break;
423     case GLDV_CON:
424         k = 0;
425         while(k<m_grayLevel){
426             i = 0;
427             while(i < m_grayLevel){
428                 j = 0;
429                 while(j < m_grayLevel){
430                     if(m_grayMatrix[i][j] != 0 && std::abs(i-j)== k){
431                         float Pij = m_grayMatrix[i][j]*1.0f/m_count;
432                         res += Pij*k*k;
433                     }
434                     j++;
435                 }
436                 i++;
437             }
438             k++;
439         }
440         break;
441     case GLDV_ASM:
442         k = 0;
443         while(k<m_grayLevel){
444             i = 0;
445             while(i < m_grayLevel){
446                 j = 0;
447                 while(j < m_grayLevel){
448                     if(m_grayMatrix[i][j] != 0 && std::abs(i-j)== k){
449                         float Pij = m_grayMatrix[i][j]*1.0f/m_count;
450                         res += Pij*Pij;
451                     }
452                     j++;
453                 }
454                 i++;
455             }
456             k++;
457         }
458         break;
459     }
460     return res;
461 }
462 
463 bool CCo_occurrenceTextureExtration::ExtraTextureWinEXE()
464 {
465     //读影像
466     GDALAllRegister();
467     GDALDataset *poDS = (GDALDataset*)GDALOpen(m_strReduceLevelFileName.c_str(),GA_ReadOnly);
468     if(!poDS){
469         return FALSE;
470     }
471     //GDALRasterBand *poBand = poDS->GetRasterBand(1);
472     double adfMSSGeoTransform[6] = {0};
473     poDS->GetGeoTransform(adfMSSGeoTransform);
474     //创建文件
475     GDALDriver *poOutDriver = (GDALDriver*)GDALGetDriverByName("GTIFF");
476     GDALDataset *poOutDS  = poOutDriver->Create(m_strOutFileName.c_str(),
477         m_iWidth,m_iHeigth,m_iBandCount,GDT_Float32,NULL);
478     poOutDS->SetGeoTransform(adfMSSGeoTransform);
479     poOutDS->SetProjection(poDS->GetProjectionRef());
480     int b = 0;
481     while(b<m_iBandCount){
482         GDALRasterBand *poBand = poDS->GetRasterBand(b+1);
483         GDALRasterBand *poNewBand = poOutDS->GetRasterBand(b+1);
484         int i = 0;
485         while(i<m_iHeigth){
486             int startR = i;    
487             int endR = startR + m_winHeigth - 1;
488             if(endR > m_iHeigth-1){
489                 endR = m_iHeigth - 1;
490                 startR = endR - m_winHeigth + 1;
491             }
492             unsigned char *readBuff = new unsigned char[m_iWidth * m_winHeigth];
493             float *writeBuff = new float[m_iWidth* m_winHeigth];
494             poBand->RasterIO(GF_Read,0,startR,m_iWidth,m_winHeigth,readBuff,m_iWidth,
495                 m_winHeigth,GDT_Byte,0,0);
496 
497             int j = 0;
498             while(j < m_iWidth){
499                 int startC = j;            
500                 int endC = startC + m_winWidth - 1;
501                 if(endC > m_iWidth - 1){
502                     endC = m_iWidth - 1;
503                     startC = endC - m_winWidth + 1;
504                 }
505                 std::vector<std::vector<int>> grayValue;
506                 grayValue = std::vector<std::vector<int>>(m_winHeigth,std::vector<int>(m_winWidth,0));
507                 std::vector<std::vector<float>> grayMatrix;
508                 grayMatrix = std::vector<std::vector<float>>(m_grayLevel,std::vector<float>(m_grayLevel,0.0f));
509                 int k = 0;
510                 while(k < m_winHeigth){
511                     int l = 0 ;
512                     while(l < m_winWidth){
513                         grayValue[k][l] = int(readBuff[k*m_iWidth +startC+l]);
514                         l++;
515                     }
516                     k++;
517                 }
518                 int count = 0;
519                 if(!GetGrayCoocurrenceMatrix(grayValue,grayMatrix,count)){
520                     return FALSE;
521                 }
522                 float tmp = CalGLCMStatistics(grayMatrix,count);            
523                 k = 0;
524                 while(k < m_winHeigth){
525                     int l = 0 ;
526                     while(l < m_winWidth){
527                         writeBuff[k*m_iWidth +startC+l] = tmp;
528                         l++;
529                     }
530                     k++;
531                 }
532                 j = j + m_winWidth;
533             }
534             poNewBand->RasterIO(GF_Write,0,startR,m_iWidth,m_winHeigth,writeBuff,m_iWidth,
535                 m_winHeigth,GDT_Float32,0,0);
536             //更新灰度共生矩阵
537             delete []readBuff;readBuff = NULL;
538             delete []writeBuff;writeBuff = NULL;
539             i = i + m_winHeigth;
540         }
541         b++;
542     }
543     GDALClose((GDALDatasetH)poDS);
544     GDALClose((GDALDatasetH)poOutDS);
545     return TRUE;
546 }
547 
548 bool CCo_occurrenceTextureExtration::ExtraTexturePixEXE()
549 {
550     //读影像
551     GDALAllRegister();
552     GDALDataset *poDS = (GDALDataset*)GDALOpen(m_strReduceLevelFileName.c_str(),GA_ReadOnly);
553     if(!poDS){
554         return FALSE;
555     }
556     double adfMSSGeoTransform[6] = {0};
557     poDS->GetGeoTransform(adfMSSGeoTransform);
558     //创建文件
559     GDALDriver *poOutDriver = (GDALDriver*)GDALGetDriverByName("GTIFF");
560     GDALDataset *poOutDS  = poOutDriver->Create(m_strOutFileName.c_str(),
561         m_iWidth,m_iHeigth,m_iBandCount,GDT_Float32,NULL);
562     poOutDS->SetGeoTransform(adfMSSGeoTransform);
563     poOutDS->SetProjection(poDS->GetProjectionRef());
564     int nLineSpace,nPixSpace,nBandSpace;
565     nLineSpace = sizeof(unsigned char)*m_iWidth*m_iBandCount;
566     nPixSpace = 0;
567     nBandSpace = sizeof(unsigned char)*m_iWidth;
568 
569     int iNumRow = 128;  // 分块读取的行数
570     int overloadPix = m_winHeigth-1;
571     if(m_iHeigth < iNumRow && m_iHeigth > m_winHeigth){
572         iNumRow = m_winHeigth;
573     }
574     if(m_iHeigth < m_winHeigth){
575         return FALSE;
576     }
577     int loopNum = (m_iHeigth + iNumRow - 1)/(iNumRow - overloadPix);
578     float stepProcess = 100.0f/(100.0f*loopNum*(iNumRow - overloadPix));
579 
580     for(int i = 0;i<loopNum;i++){
581         int tmpRowNum = iNumRow;
582         int startR = i*iNumRow - i*(overloadPix/2);
583         if(startR < 0)
584         {
585             startR = 0;
586         }
587         int endR = startR + tmpRowNum - 1;
588         if(endR > m_iHeigth -1 ){
589             endR = m_iHeigth -1;
590             tmpRowNum = endR - startR + 1;
591         }
592         unsigned char *readBuff = new unsigned char[tmpRowNum*m_iWidth*m_iBandCount];
593         float *writeBuff = new float[m_iWidth*m_iBandCount*(tmpRowNum - overloadPix)];
594         poDS->RasterIO(GF_Read,0,startR,m_iWidth,tmpRowNum,readBuff,m_iWidth,
595             tmpRowNum,GDT_Byte,m_iBandCount,m_BandMap,nPixSpace,nLineSpace,nBandSpace);
596 
597         for(int iR = overloadPix/2;iR < tmpRowNum - overloadPix/2;iR++)
598         {
599             // zhuyong 2016 10 17
600 #pragma omp parallel for num_threads(2)
601             for(int j = 0;j < m_iWidth;j++)
602             {
603                 int startC = j - m_winWidth/2;
604                 int endC = j + m_winWidth/2;
605                 if(startC<0){
606                     startC = 0;
607                     endC = startC + m_winWidth -1;
608                 }
609                 if(endC > m_iWidth-1){
610                     endC = m_iWidth -1;
611                     startC = endC - m_winWidth +1;
612                 }
613                 // 波段计算
614                 int b = 0;
615                 while(b < m_iBandCount){
616                     // 获取灰度值
617                     std::vector<std::vector<int>> grayValue;
618                     grayValue = std::vector<std::vector<int>>(m_winHeigth,std::vector<int>(m_winWidth,0));
619                     std::vector<std::vector<float>> grayMatrix;
620                     grayMatrix = std::vector<std::vector<float>>(m_grayLevel,std::vector<float>(m_grayLevel,0.0f));
621                     int k = 0;
622                     while(k < m_winHeigth){
623                         int l = 0 ;
624                         while(l < m_winWidth){
625                             grayValue[k][l] = int(readBuff[(iR  - overloadPix/2 + k)*m_iWidth*m_iBandCount + b*m_iWidth + startC+l]);
626                             l++;
627                         }
628                         k++;
629                     }
630                     int count = 0;
631                     if(!GetGrayCoocurrenceMatrix(grayValue,grayMatrix,count)){
632                     }
633                     float tmp = CalGLCMStatistics(grayMatrix,count);
634                     writeBuff[(iR- overloadPix/2)*m_iWidth*m_iBandCount+ b*m_iWidth + j] = tmp;
635                     b++;
636                 }
637             }
638         }
639         poOutDS->RasterIO(GF_Write,0,startR+overloadPix/2,m_iWidth,tmpRowNum - overloadPix,writeBuff,
640             m_iWidth,tmpRowNum - overloadPix,GDT_Float32,m_iBandCount,m_BandMap,nPixSpace,sizeof(float)*m_iWidth*m_iBandCount,
641             sizeof(float)*m_iWidth);
642         delete []readBuff;readBuff = NULL;
643         delete []writeBuff;writeBuff = NULL;
644     }
645 
646     GDALClose((GDALDatasetH)poDS);
647     GDALClose((GDALDatasetH)poOutDS);
648     return TRUE;
649 }
650 
651 int CCo_occurrenceTextureExtration::ExtraTexturePixOpenCLEXE()
652 {
653     //读影像
654     GDALAllRegister();
655     GDALDataset *poDS = (GDALDataset*)GDALOpen(m_strReduceLevelFileName.c_str(),GA_ReadOnly);
656     if(!poDS){
657         return FALSE;
658     }
659     double adfMSSGeoTransform[6] = {0};
660     poDS->GetGeoTransform(adfMSSGeoTransform);
661     //创建文件
662     GDALDriver *poOutDriver = (GDALDriver*)GDALGetDriverByName("GTIFF");
663     GDALDataset *poOutDS  = poOutDriver->Create(m_strOutFileName.c_str(),
664         m_iWidth,m_iHeigth,m_iBandCount,GDT_Float32,NULL);
665     poOutDS->SetGeoTransform(adfMSSGeoTransform);
666     poOutDS->SetProjection(poDS->GetProjectionRef());
667 
668     int nLineSpace,nPixSpace,nBandSpace;
669     nLineSpace = sizeof(unsigned char)*m_iWidth*m_iBandCount;
670     nPixSpace = 0;
671     nBandSpace = sizeof(unsigned char)*m_iWidth;
672 
673     int iNumRow = 32;  // 分块读取的行数
674     int overloadPix = m_winHeigth-1;
675     if(m_iHeigth < iNumRow && m_iHeigth > m_winHeigth){
676         iNumRow = m_winHeigth;
677     }
678     if(m_iHeigth < m_winHeigth){
679         return FALSE;
680     }
681     int loopNum = ceil((m_iHeigth*1.0f - iNumRow)/(iNumRow - overloadPix)) + 1;
682 
683     // OpenCL部分 =============== 1 创建平台
684     cl_uint num_platforms;
685     cl_int ret = clGetPlatformIDs(0,NULL,&num_platforms);
686     if(ret != CL_SUCCESS || num_platforms < 1){
687         printf("clGetPlatformIDs Error\n");
688         return -1;
689     }
690     cl_platform_id platform_id = NULL;
691     ret = clGetPlatformIDs(1,&platform_id,NULL);
692     if(ret != CL_SUCCESS){
693         printf("clGetPlatformIDs Error2\n");
694         return -1;
695     }
696 
697     // OpenCL部分 =============== 2 获得设备
698     cl_uint num_devices;
699     ret = clGetDeviceIDs(platform_id,CL_DEVICE_TYPE_GPU,0,NULL,
700         &num_devices);
701     if(ret != CL_SUCCESS || num_devices < 1){
702         printf("clGetDeviceIDs Error\n");
703         return -1;
704     }
705     cl_device_id device_id;
706     ret = clGetDeviceIDs(platform_id,CL_DEVICE_TYPE_GPU,1,&device_id,NULL);
707     if(ret != CL_SUCCESS){
708         printf("clGetDeviceIDs Error2\n");
709         return -1;
710     }
711 
712     // OpenCL部分 =============== 3 创建Context
713     cl_context_properties props[] = {CL_CONTEXT_PLATFORM,
714         (cl_context_properties)platform_id,0};
715     cl_context context = NULL;
716     context = clCreateContext(props,1,&device_id,NULL,NULL,&ret);
717     if(ret != CL_SUCCESS || context == NULL){
718         printf("clCreateContext Error\n");
719         return -1;
720     }
721 
722     // OpenCL部分 =============== 4 创建Command Queue
723     cl_command_queue command_queue = NULL;
724     command_queue = clCreateCommandQueue(context,device_id,0,&ret);
725     if(ret != CL_SUCCESS || command_queue == NULL){
726         printf("clCreateCommandQueue Error\n");
727         return -1;
728     }
729 
730     // OpenCL部分 =============== 6 创建编译Program
731     const char *strfile = "D:\\PIE3\\src\\Test\\TextOpecCLResample\\TextOpecCLResample\\Co_occurrenceMatrixKernel.txt";
732     size_t lenSource = 0;
733     char *kernelSource = LoadProgSource(strfile,"",&lenSource);
734     cl_program *programs = (cl_program *)malloc(loopNum*sizeof(cl_program));
735     memset(programs,0,sizeof(cl_program)*loopNum);
736 
737     cl_kernel *kernels = (cl_kernel*)malloc(loopNum*sizeof(cl_kernel));
738     memset(kernels,0,sizeof(cl_kernel)*loopNum);
739 
740     int startR = 0;
741     for(int i = 0;i<loopNum;i++)
742     {
743         int tmpRowNum = iNumRow;
744         if(startR < 0)
745         {
746             startR = 0;
747         }
748         int endR = startR + tmpRowNum - 1;
749         if(endR > m_iHeigth -1 ){
750             endR = m_iHeigth -1;
751             tmpRowNum = endR - startR + 1;
752         }
753         unsigned char *readBuff = new unsigned char[tmpRowNum*m_iWidth*m_iBandCount];
754         float *writeBuff = new float[m_iWidth*m_iBandCount*(tmpRowNum - overloadPix)];
755         poDS->RasterIO(GF_Read,0,startR,m_iWidth,tmpRowNum,readBuff,m_iWidth,
756             tmpRowNum,GDT_Byte,m_iBandCount,m_BandMap,nPixSpace,nLineSpace,nBandSpace);
757 
758         // OpenCL部分 =============== 5 创建Memory Object
759         cl_mem mem_read = NULL;
760         mem_read = clCreateBuffer(context,CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR,
761             sizeof(cl_uchar)*tmpRowNum*m_iWidth*m_iBandCount,readBuff,&ret);
762         if(ret != CL_SUCCESS || NULL == mem_read){
763             printf("clCreateBuffer Error\n");
764             return -1;
765         }
766 
767         cl_mem mem_write = NULL;
768         mem_write = clCreateBuffer(context,CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR,
769             sizeof(cl_float)*m_iWidth*m_iBandCount*(tmpRowNum - overloadPix),writeBuff,&ret);
770         if(ret != CL_SUCCESS || NULL == mem_write){
771             printf("clCreateBuffer Error\n");
772             return -1;
773         }
774 
775         // OpenCL部分 =============== 6 创建编译Program
776         //const char *strfile = "D:\\PIE3\\src\\Test\\TextOpecCLResample\\TextOpecCLResample\\Co_occurrenceMatrixKernel.txt";
777         //size_t lenSource = 0;
778         //char *kernelSource = LoadProgSource(strfile,"",&lenSource);
779         //cl_program program = NULL;
780         programs[i] = clCreateProgramWithSource(context,1,(const char**)&kernelSource,
781             NULL,&ret);
782         if(ret != CL_SUCCESS || NULL == programs[i]){
783             printf("clCreateProgramWithSource Error\n");
784             return -1;
785         }
786         ret = clBuildProgram(programs[i],1,&device_id,NULL,NULL,NULL);
787         if(ret != CL_SUCCESS){
788             char* build_log;
789             size_t log_size;
790             //查询日志的大小
791             clGetProgramBuildInfo(programs[i], device_id, CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
792             build_log = new char[log_size+1];
793             //获得编译日志信息
794             ret = clGetProgramBuildInfo(programs[i], device_id, CL_PROGRAM_BUILD_LOG, log_size, build_log, NULL);
795             build_log[log_size] = '\0';
796             printf("%s\n",build_log);
797             printf("编译失败!");
798             delete []build_log;
799             return -1;
800         }
801 
802         // OpenCL部分 =============== 7 创建Kernel
803         //cl_kernel kernel = NULL;
804         kernels[i] = clCreateKernel(programs[i],"Co_occurrenceMatrixKernel",&ret);
805         if(ret != CL_SUCCESS || NULL == kernels[i]){
806             printf("clCreateProgramWithSource Error\n");
807             return -1;
808         }
809 
810         // OpenCL部分 =============== 8 设置Kernel参数
811         ret = clSetKernelArg(kernels[i],0,sizeof(cl_mem),&mem_read);
812         ret |= clSetKernelArg(kernels[i],1,sizeof(cl_mem),&mem_write);
813         ret |= clSetKernelArg(kernels[i],2,sizeof(cl_int),&m_iHeigth);
814         ret |= clSetKernelArg(kernels[i],3,sizeof(cl_int),&m_iWidth);
815         ret |= clSetKernelArg(kernels[i],4,sizeof(cl_int),&m_iBandCount);
816         ret |= clSetKernelArg(kernels[i],5,sizeof(cl_int),&tmpRowNum);
817         ret |= clSetKernelArg(kernels[i],6,sizeof(cl_int),&overloadPix);
818         ret |= clSetKernelArg(kernels[i],7,sizeof(cl_int),&m_grayLevel);
819         ret |= clSetKernelArg(kernels[i],8,sizeof(cl_int),&m_winHeigth);
820         ret |= clSetKernelArg(kernels[i],9,sizeof(cl_int),&m_winWidth);
821         ret |= clSetKernelArg(kernels[i],10,sizeof(cl_int),&m_dis);
822         ret |= clSetKernelArg(kernels[i],11,sizeof(cl_int),&m_AngleMode);
823         ret |= clSetKernelArg(kernels[i],12,sizeof(cl_int),&m_TextureMode);
824         if(ret != CL_SUCCESS){
825             printf("clSetKernelArg Error\n");
826             return -1;
827         }
828 
829         // OpenCL部分 =============== 9 设置Group Size
830         cl_uint work_dim = 2;
831         size_t global_work_size[] = {m_iWidth,tmpRowNum};
832         size_t *local_work_size = NULL;
833 
834         // OpenCL部分 =============== 10 执行内核
835         ret = clEnqueueNDRangeKernel(command_queue,kernels[i],work_dim,NULL,global_work_size,
836             local_work_size,0,NULL,NULL);
837         ret |= clFinish(command_queue);
838         if(ret != CL_SUCCESS){
839             printf("clEnqueueNDRangeKernel Error\n");
840             return -1;
841         }
842 
843         writeBuff = (float*)clEnqueueMapBuffer(command_queue,mem_write,CL_TRUE,CL_MAP_READ | CL_MAP_WRITE,
844             0,sizeof(cl_float)*m_iWidth*m_iBandCount*(tmpRowNum - overloadPix),0,NULL,NULL,&ret);
845         //ret = clEnqueueReadBuffer(command_queue,mem_resample,CL_TRUE,0,
846         //    sizeof(cl_float)*tmpRowNum*resampleWidth*mssBandCount,(void*)resampleBuf,0,NULL,NULL);
847         if(ret != CL_SUCCESS){
848             printf("clEnqueueMapBuffer Error\n");
849             return -1;
850         }
851         poOutDS->RasterIO(GF_Write,0,startR+overloadPix/2,m_iWidth,tmpRowNum - overloadPix,writeBuff,
852             m_iWidth,tmpRowNum - overloadPix,GDT_Float32,m_iBandCount,m_BandMap,nPixSpace,sizeof(float)*m_iWidth*m_iBandCount,
853             sizeof(float)*m_iWidth);
854         delete []readBuff;readBuff = NULL;
855         delete []writeBuff;writeBuff = NULL;
856 
857         // OpenCL部分 =============== 12 释放资源
858         if(NULL != mem_read) clReleaseMemObject(mem_read);
859         if(NULL != mem_write) clReleaseMemObject(mem_write);
860         std::cout<<i<<"  "<<loopNum<<std::endl;
861         startR = startR + tmpRowNum - overloadPix;
862     }
863     // OpenCL部分 =============== 12 释放资源
864     int i = 0;
865     while(i < loopNum){
866         if(NULL != kernels[i]) clReleaseKernel(kernels[i]);
867         if(NULL != programs[i]) clReleaseProgram(programs[i]);
868         i++;
869     }
870     if(NULL != command_queue) clReleaseCommandQueue(command_queue);
871     if(NULL != context) clReleaseContext(context);
872     GDALClose((GDALDatasetH)poDS);
873     GDALClose((GDALDatasetH)poOutDS);
874     return 0;
875 }

GPU核函数：

  1 #pragma OPENCL EXTENSION cl_amd_printf:enable
  2 
  3 int GetGrayCoocurrenceMatrix(const int *grayValue,
  4                               float *grayMatrix,
  5                               int winHeight,
  6                               int winWidth,
  7                               int grayLevel,
  8                               int dis,
  9                               int AngleMode)
 10 {
 11     int i,j,r,c;
 12     int count = 0;
 13     switch(AngleMode)
 14     {
 15     case 0:
 16         i = 0;
 17         while(i<winHeight){
 18             j = 0;
 19             while(j<winWidth){
 20                 int endR = i;
 21                 int endC = j + dis;
 22                 if(endC < winWidth){
 23                     r = grayValue[i*winHeight + j];
 24                     c = grayValue[endR*winHeight + endC];
 25                     grayMatrix[r*grayLevel + c] += 1;
 26                     grayMatrix[c*grayLevel + r] += 1;
 27                     count = count+2;
 28                 }
 29                 j++;
 30             }
 31             i++;
 32         }
 33         break;
 34     case 1:
 35         i = 0;
 36         while(i<winHeight){
 37             j = 0;
 38             while(j<winWidth){
 39                 int endR = i - dis;
 40                 int endC = j + dis;
 41                 if(endR>=0 && endR < winHeight && endC >=0 && endC < winWidth){
 42                     r = grayValue[i*winHeight + j];
 43                     c = grayValue[endR*winHeight + endC];
 44                     grayMatrix[r*grayLevel + c] += 1;
 45                     grayMatrix[c*grayLevel + r] += 1;
 46                     count = count+2;
 47                 }
 48                 j++;
 49             }
 50             i++;
 51         }
 52         break;
 53     case 2:
 54         i= 0;
 55         while(i<winHeight){
 56             j = 0;
 57             while(j < winWidth){
 58                 int endR = i - dis;
 59                 int endC = j;
 60                 if(endR >= 0){
 61                     r = grayValue[i*winHeight + j];
 62                     c = grayValue[endR*winHeight + endC];
 63                     grayMatrix[r*grayLevel + c] += 1;
 64                     grayMatrix[c*grayLevel + r] += 1;
 65                     count = count+2;
 66                 }
 67                 j++;
 68             }
 69             i++;
 70         }
 71         break;
 72     case 3:
 73         i = 0;
 74         while(i<winHeight){
 75             j = 0;
 76             while(j<winWidth){
 77                 int endR = i - dis;
 78                 int endC = j - dis;
 79                 if(endR >= 0 && endC >= 0){
 80                     r = grayValue[i*winHeight + j];
 81                     c = grayValue[endR*winHeight + endC];
 82                     grayMatrix[r*grayLevel + c] += 1;
 83                     grayMatrix[c*grayLevel + r] += 1;
 84                     count = count+2;
 85                 }
 86                 j++;
 87             }
 88             i++;
 89         }
 90         break;
 91     default:
 92         return 0;
 93     }
 94     return count;
 95 }
 96 
 97 float CalGLCMStatistics(float *grayMatrix,
 98                         const int grayLevel,
 99                         const int count,
100                         int TextureMode)
101 {
102     int i,j,k;
103     float imean = 0.0f;
104     float jmean = 0.0f;
105     float ivar = 0.0f;
106     float jvar = 0.0f;;
107     float res = 0.0f;
108     switch(TextureMode)
109     {
110     case 0:  //同质性
111         i = 0;
112         while(i < grayLevel){
113             j = 0;
114             while(j < grayLevel){
115                 if(grayMatrix[i*grayLevel + j] != 0){
116                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
117                     res += Pij/(1 + (i-j)*(i-j));
118                 }
119                 j++;
120             }
121             i++;
122         }
123         break;
124     case 1:  //对比度
125         i = 0;
126         while(i < grayLevel){
127             j = 0;
128             while(j < grayLevel){
129                 if(grayMatrix[i*grayLevel + j] != 0){
130                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
131                     res += Pij*(i-j)*(i-j);
132                 }
133                 j++;
134             }
135             i++;
136         }
137         break;
138     case 2:  //非相似性
139         i = 0;
140         while(i < grayLevel){
141             j = 0;
142             while(j < grayLevel){
143                 if(grayMatrix[i*grayLevel + j] != 0){
144                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
145                     res += Pij*abs(i-j);
146                 }
147                 j++;
148             }
149             i++;
150         }
151         break;
152     case 3:  // 相关性
153         i = 0;
154         while(i < grayLevel){
155             j = 0;
156             while(j < grayLevel){
157                 if(grayMatrix[i*grayLevel + j] != 0){
158                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
159                     imean += Pij * i;
160                     jmean += Pij * j;
161                 }
162                 j++;
163             }
164             i++;
165         }
166         //  var
167         i = 0;
168         while(i < grayLevel){
169             j = 0;
170             while(j < grayLevel){
171                 if(grayMatrix[i*grayLevel + j] != 0){
172                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
173                     ivar += Pij * (i - imean)*(i - imean);
174                     jvar += Pij * (j - imean)*(j - imean);
175                 }
176                 j++;
177             }
178             i++;
179         }
180         // xiangguanx
181         i = 0;
182         while(i < grayLevel){
183             j = 0;
184             while(j < grayLevel){
185                 if(grayMatrix[i*grayLevel + j] != 0){
186                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
187                     res += (Pij *(i - imean)*(j - jmean))/sqrt(ivar*jvar);
188                 }
189                 j++;
190             }
191             i++;
192         }
193         break;
194     case 4:  // 均值
195         i = 0;
196         while(i < grayLevel){
197             j = 0;
198             while(j < grayLevel){
199                 if(grayMatrix[i*grayLevel + j] != 0){
200                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
201                     imean += Pij * i;
202                     //jmean += Pij * j;
203                 }
204                 j++;
205             }
206             i++;
207         }
208         res = imean;
209         break;
210     case 5:  // 方差
211         //  mean
212         i = 0;
213         while(i < grayLevel){
214             j = 0;
215             while(j < grayLevel){
216                 if(grayMatrix[i*grayLevel + j] != 0){
217                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
218                     imean += Pij * i;
219                 }
220                 j++;
221             }
222             i++;
223         }
224         //  var
225         i = 0;
226         while(i < grayLevel){
227             j = 0;
228             while(j < grayLevel){
229                 if(grayMatrix[i*grayLevel + j] != 0){
230                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
231                     ivar += Pij * (i - imean)*(i - imean);
232                 }
233                 j++;
234             }
235             i++;
236         }
237         res = ivar;
238         break;
239     case 6:  // 角二阶矩
240         i = 0;
241         while(i < grayLevel){
242             j = 0;
243             while(j < grayLevel){
244                 if(grayMatrix[i*grayLevel + j] != 0){
245                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
246                     res += Pij*Pij;
247                 }
248                 j++;
249             }
250             i++;
251         }
252         break;
253     case 7:  // 熵
254         i = 0;
255         while(i < grayLevel){
256             j = 0;
257             while(j < grayLevel){
258                 if(grayMatrix[i*grayLevel + j] != 0){
259                     float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
260                     res += -1.0f*Pij * log(Pij);
261                 }
262                 j++;
263             }
264             i++;
265         }
266         break;
267     case 8:
268         k = 0;
269         while(k<grayLevel){
270             i = 0;
271             while(i < grayLevel){
272                 j = 0;
273                 while(j < grayLevel){
274                     if(grayMatrix[i*grayLevel + j] != 0 && abs(i-j)== k){
275                         float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
276                         res += Pij*k;
277                     }
278                     j++;
279                 }
280                 i++;
281             }
282             k++;
283         }
284         break;
285     case 9:
286         k = 0;
287         while(k<grayLevel){
288             i = 0;
289             while(i < grayLevel){
290                 j = 0;
291                 while(j < grayLevel){
292                     if(grayMatrix[i*grayLevel + j] != 0 && abs(i-j)== k){
293                         float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
294                         res += Pij*k*k;
295                     }
296                     j++;
297                 }
298                 i++;
299             }
300             k++;
301         }
302         break;
303     case 10:
304         k = 0;
305         while(k<grayLevel){
306             i = 0;
307             while(i < grayLevel){
308                 j = 0;
309                 while(j < grayLevel){
310                     if(grayMatrix[i*grayLevel + j] != 0 && abs(i-j)== k){
311                         float Pij = grayMatrix[i*grayLevel + j]*1.0f/count;
312                         res += Pij*Pij;
313                     }
314                     j++;
315                 }
316                 i++;
317             }
318             k++;
319         }
320         break;
321     }
322     return res;
323 }
324 
325 
326 
327 
328 
329 
330 __kernel void Co_occurrenceMatrixKernel(__global const unsigned char *oriData,
331                                         __global float *writeData,
332                                         int iHeight,
333                                         int iWidth,
334                                         int iBandCount,
335                                         int curRowNum,
336                                         int overloadPix,
337                                         int grayLevel,
338                                         int winHeight,
339                                         int winWidth,
340                                         int dis,
341                                         int AngleMode,
342                                         int TextureMode)
343 {
344     int idx = get_global_id(1);  // 采样行
345     int idy = get_global_id(0);  // 采样列
346     if(idx >= overloadPix/2 && idx < curRowNum - overloadPix/2 && idy < iWidth)
347     {
348         // 获得灰度值
349         int startC = idy - winWidth/2;
350         int endC = idy + winWidth/2;
351         if(startC<0){
352             startC = 0;
353             endC = startC + winWidth -1;
354         }
355         if(endC > iWidth-1){
356             endC = iWidth -1;
357             startC = endC - winWidth +1;
358         }
359         int b = 0;
360         while(b < iBandCount)
361         {
362             int grayValue[625] = {0};
363             int k = 0;
364             while(k < winHeight){
365                 int l = 0 ;
366                 while(l < winWidth){
367                     int tmpP = (int)(idx  - overloadPix/2 + k)*iWidth*iBandCount + b*iWidth + startC+l;
368                     grayValue[k*winHeight + l] = (int)oriData[tmpP];
369                     l++;
370                 }
371                 k++;
372             }
373             float grayMatrix[1024] = {0};
374             int count = GetGrayCoocurrenceMatrix(grayValue,grayMatrix,winHeight,winWidth,grayLevel,dis,AngleMode);
375             float tmp = CalGLCMStatistics(grayMatrix,grayLevel,count,TextureMode);
376             writeData[(idx- overloadPix/2)*iWidth*iBandCount+ b*iWidth + idy] = tmp;
377             b++;
378         }
379     }
380 }

posted on 2016-10-19 22:00 机器学习猪阅读(15136) 评论(0) 编辑收藏举报

刷新页面返回顶部