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

1 #include "Co_occurrenceTextureExtration.h"

2 #include

3 #include

4 #include

5 #include

6 #include

7 #include

8 #include

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 boolCCo_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

40 boolCCo_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;iGetRasterBand(m_BandMap[i]);58 GDALRasterBand *poNewBand = poOutDS->GetRasterBand(i+1);59 poBand->ComputeRasterMinMax(FALSE,bandMINMAX);60 for(int j = 0;jRasterIO(GF_Read,0,j,m_iWidth,1,readBuff,m_iWidth,1,GDT_Float64,0,0);64 int k = 0;65 while(kRasterIO(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 returnTRUE;84 }85

86 double* CCo_occurrenceTextureExtration::cdf( int *h,intlength )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 returnp;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 returnNULL;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 returncSourceString;142 }143

144 bool CCo_occurrenceTextureExtration::GetGrayCoocurrenceMatrix(std::vector<:vector>>grayValue,145 std::vector<:vector>> &m_grayMatrix,146 int &m_count)147 {148 inti,j,r,c;149 m_count = 0;150 int ii = 0;151

152 switch(m_AngleMode)153 {154 caseDEGREE_0:155 i = 0;156 while(i=0 && endR < m_winHeigth && endC >=0 && endC = 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 caseDEGREE_135:212 i = 0;213 while(i= 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 returnFALSE;232 }233 returnTRUE;234 }235

236 float CCo_occurrenceTextureExtration::CalGLCMStatistics(std::vector<:vector>>m_grayMatrix,237 const intm_count)238 {239 inti,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

262 i = 0;263 while(i

276 i = 0;277 while(i

290 i = 0;291 while(i

304 i = 0;305 while(i

318 i = 0;319 while(i

332 i = 0;333 while(i

340 }341 j++;342 }343 i++;344 }345 res =imean;346 break;347 case GLCM_VAR: //方差348 //mean

349 i = 0;350 while(i

362 i = 0;363 while(i

375 res =ivar;376 break;377 case GLCM_ASM: //角二阶矩

378 i = 0;379 while(i

392 i = 0;393 while(i

463 boolCCo_occurrenceTextureExtration::ExtraTextureWinEXE()464 {465 //读影像

466 GDALAllRegister();467 GDALDataset *poDS = (GDALDataset*)GDALOpen(m_strReduceLevelFileName.c_str(),GA_ReadOnly);468 if(!poDS){469 returnFALSE;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(bGetRasterBand(b+1);483 GDALRasterBand *poNewBand = poOutDS->GetRasterBand(b+1);484 int i = 0;485 while(i 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 - 1){502 endC = m_iWidth - 1;503 startC = endC - m_winWidth + 1;504 }505 std::vector<:vector>>grayValue;506 grayValue = std::vector<:vector>>(m_winHeigth,std::vector(m_winWidth,0));507 std::vector<:vector>>grayMatrix;508 grayMatrix = std::vector<:vector>>(m_grayLevel,std::vector(m_grayLevel,0.0f));509 int k = 0;510 while(k 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 returnTRUE;546 }547

548 boolCCo_occurrenceTextureExtration::ExtraTexturePixEXE()549 {550 //读影像

551 GDALAllRegister();552 GDALDataset *poDS = (GDALDataset*)GDALOpen(m_strReduceLevelFileName.c_str(),GA_ReadOnly);553 if(!poDS){554 returnFALSE;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 intnLineSpace,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

580 for(int i = 0;i 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

617 std::vector<:vector>>grayValue;618 grayValue = std::vector<:vector>>(m_winHeigth,std::vector(m_winWidth,0));619 std::vector<:vector>>grayMatrix;620 grayMatrix = std::vector<:vector>>(m_grayLevel,std::vector(m_grayLevel,0.0f));621 int k = 0;622 while(k 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 returnTRUE;649 }650

651 intCCo_occurrenceTextureExtration::ExtraTexturePixOpenCLEXE()652 {653 //读影像

654 GDALAllRegister();655 GDALDataset *poDS = (GDALDataset*)GDALOpen(m_strReduceLevelFileName.c_str(),GA_ReadOnly);656 if(!poDS){657 returnFALSE;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 intnLineSpace,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

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 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 创建编译Program776 //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 创建Kernel803 //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<

864 int i = 0;865 while(i