精华内容
下载资源
问答
  • 今天学习数字图像处理课,学到灰度共生矩阵,兴起,遂敲码解之。 方法摘自百度百科: 灰度直方图是对图像上单个像素具有某个灰度进行统计的结果,而灰度共生矩阵是对图像上保持某距离的两像素分别具有某灰度的状况...
  • 使用opencv以及numpy两个库,使用python写成。网上很多只写了特征值的生成,这个代码增加了滑动窗口,让生成的特征值赋予到像素点,从而生成特征图像。
  • 灰度共生矩阵纹理特征提取代码 matlab实现
  • 实验室师兄写的代码,使用灰度共生矩阵计算图像的纹理复杂度,然后可以用来衡量嵌入水印的多少。
  • 灰度共生矩阵纹理特征提取的Matlab实现
  • 灰度共生矩阵纹理特征提取的Matlab实现.pdf
  • MATLAB灰度共生矩阵纹理特征提取

    热门讨论 2013-07-17 21:13:46
    MATLAB 灰度共生矩阵 纹理特征提取 粗糙度、对比度、方向度等,源代码
  • 纹理复杂度计算-灰度共生矩阵,灰度共生矩阵纹理特征提取,matlab源码.zip
  • 利用灰度共生矩阵,对纹理图像进行分割,里面有代码和测试图像
  • 利用灰度共生矩阵实现图像纹理特征的提取,从而实现对图像纹理特征的研究
  • 灰度共生矩阵提取图像纹理特征,采用matlab实现,包括模糊c均值实现分类。代码完整,正确可运行。
  • 灰度共生矩阵纹理特征提取的Matlab实现
  • 通过灰度共生矩阵获取图像的纹理特征,含有相关函数,可直接运行
  • 基于灰度共生矩阵求的彩色图像的纹理特征,计算求的一些纹理信息。
  • 同时考虑SAR图像局部灰度均值和方差及像素空间分布特征等统计量,在以灰度共生矩阵产生的纹理统计量为特征所生成的图像上,建立多分辨双Markov-GAR模型,采用多分辨MPM的参数估计方法及相应的无监督分割算法,对SAR...
  • 通过MATLAB运用灰度共生矩阵提取熵、能量、对比度、相关等特征
  • Matlab求灰度共生矩阵特征

    热门讨论 2012-05-23 16:40:37
    该程序用于求解数字图像处理的灰度共生矩阵纹理特征值,如熵、对比度,同质性、能量等。只需要将该m文件放在Matlab的安装目录:toolbox/images/images文件夹里,按照参数设定,直接调用即可。
  • 灰度共生矩阵来描述纹理特征,这是由于纹理是由灰度分布在空间位置上反复交替变化而形成的,因而在图像空间中相隔某距离的两个像素间一定存在一定的灰度关系,称为是图像中灰度的空间相关特性,通过研究灰度的空间...
  • 运用MATLAB R2014a来完成灰度共生矩阵特征参数的求解。以纸作为纹理分析的对象。首先需将彩色图像将各颜色分量转化为灰度。所用图像的灰度级为256。为了减少计算量,对原始图像灰度级压缩,将灰度量化成16级。计算...
  • 灰度共生矩阵作为机器视觉检测的方法之一,在最近几年被广泛使用。灰度共生矩阵主要在有三个辨别力最好的特征:对比度、熵和相关性,以上三种方法都是能量的体现
  • 图像纹理特征获取和利用knn进行分类,进行几种纹理特征的计算和数据处理,从0开始

    灰度共生矩阵的纹理特征、颜色特征和形状特征利用KNN进行分类

    纹理特征多种属性

    一、获取合适的纹理特征数据

    灰度共生矩阵图像预处理

    先降分辨率到552*1080,提升程序运行速度,进行灰度共生矩阵的各项特征进行处理,然后进行数据比对,选出对于图片区分率最高的几组数据,并且设定阈值

    # 转换为灰度图像并且降低图片分辨率
    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) 
    img_gray = cv2.resize(img_gray, (552, 1080))    
    

    参考文章:
    求灰度共生矩阵:

    def glcm(arr, d_x, d_y, gray_level=16):
        # 计算并返回归一化后的灰度共生矩阵
        max_gray = arr.max()
        height, width = arr.shape
        arr = arr.astype(np.float64)  # 将uint8类型转换为float64,避免数据失真
        # 若灰度级数大于gray_level,则将图像的灰度级缩小至gray_level,减小灰度共生矩阵的大小。量化后灰度值范围:0 ~ gray_level - 1
        arr = arr * (gray_level - 1) // max_gray  
        ret = np.zeros([gray_level, gray_level])
        for j in range(height - abs(d_y)):
            for i in range(width - abs(d_x)): 
                rows = arr[j][i].astype(int)
                cols = arr[j + d_y][i + d_x].astype(int)
                ret[rows][cols] += 1
        if d_x >= d_y:
            # 归一化, 水平方向或垂直方向
            ret = ret / float(height * (width - 1))  
        else:
            # 归一化, 45度或135度方向
            ret = ret / float((height - 1) * (width - 1))  
    
        return ret
    
    
    # 四个方向的共生矩阵
        glcm_0 = glcm(img_gray, 1, 0)  # 水平方向
        glcm_45 = glcm(img_gray, 1, 1)  # 45度方向
        glcm_90 = glcm(img_gray, 0, 1)  # 垂直方向
        glcm_135 = glcm(img_gray, -1, 1)  # 135度方向
    

    1.灰度共生矩阵的能量

    ​ 计算各个角度的平均能量,公式为:
    ∑ i ∑ j P 2 ( i , j ) \sum_{i}\sum_{j}P^2(i,j) ijP2(i,j)
    ​ 代码为:

    # 计算能量
    def power(array):
        rows, cols = array.shape
        pow = 0
        for i in range(rows):
            for j in range(cols):
                pow += np.square(array[i][j])
    
        return pow
    

    2.对比度

    计算各个方向对比度差异,公式为:
    ∑ i ∑ j ∣ i − j ∣ 2 P ( i , j ) \sum_{i}\sum_{j}\lvert{i-j}\rvert{^2}P(i,j) ijij2P(i,j)

    代码为:

    # 计算对比度
    def duibidu(array):
        rows, cols = array.shape
        duibidu = 0
        for i in range(rows):
            for j in range(cols):
                # duibidu += {np.square{abs(i-j)}}*array[i, j]
                duibidu += np.square(abs(i - j))*array[i, j]
        return duibidu
    

    3.最大概率

    计算四个方向里的最大概率,公式为:
    M a x i j P ( i , j ) Max_{ij}P(i,j) MaxijP(i,j)

    代码块:

    # 计算最大概率值
    def maxp(array):
        a = array[0, 0]
        rows, cols = array.shape
        for i in range(rows):
            for j in range(cols):
                if a <= array[i, j]:
                    a = array[i, j]
                else:
                    a = a
        return a
    

    4.均匀度

    计算四个方向分布的均匀度,公式为:
    ∑ i ∑ j P ( i , j ) 1 + ∣ i − j ∣ \sum_i\sum_j\frac{P(i,j)}{1+\lvert{i-j}\rvert} ij1+ijP(i,j)

    代码块:

    # 计算均匀度
    def junyundu(array):
        rows, cols = array.shape
        junyun = 0
        for i in range(rows):
            for j in range(cols):
                junyun += array[i, j]/(1+abs(i-j))
    
        return junyun
    

    5.逆差分矩

    四个方向的逆差分矩,公式为:
    ∑ i ∑ j P ( i , j ) ∣ i − j ∣ 2 \sum_i\sum_j\frac{P(i,j)}{\lvert{i-j}\rvert^2} ijij2P(i,j)

    代码块:

    # 计算逆差分矩
    def nicha(array):
        rows, cols = array.shape
        x = 0
        for i in range(rows):
            for j in range(cols):
                if i != j:
                    x += array[i, j]/(np.square(abs(i-j)))
                else:
                    continue
        return x
    

    6.差异方差

    计算四个方向的差异方差,公式为:
    ∑ i ∑ j ∣ i − j ∣ P ( i , j ) \sum_i\sum_j\vert{i-j}\rvert{P(i,j)} ijijP(i,j)

    代码块:

    # 计算差异分差
    def chayi(array):
        rows, cols = array.shape
        a = 0
        for i in range(rows):
            for j in range(cols):
                a += abs(i-j)*array[i, j]
        return a
    

    7.熵

    计算四个方向的熵值,公式为:
    − ∑ i ∑ j P ( i , j ) log ⁡ 2 P ( i , j ) -\sum_i\sum_jP(i,j)\log_2P(i,j) ijP(i,j)log2P(i,j)

    def shang(array):
        rows, cols = array.shape
        a = 0
        for i in range(rows):
            for j in range(cols):
                if array[i, j] != 0:
                    a += array[i, j]*(math.log(array[i, j], 2))
                else:
                    continue
        return -a
    

    图像数据选取

    时间原因选择能量值和熵值作为分类的指标

    二、颜色特征

    图片转换到HSV空间进行饱和度和明度的获取H和S

    #   获取图像的色调和饱和度
    def color(image):
        rows, cols, channels = image.shape
        #   对图像进行预处理后再转换空间,减小误差,把烟叶以外区域设置为纯黑色
        for row in range(rows):
            for col in range(cols):
                if image[row, col, 0] < 30:
                    image[row][col][0] = 0
                    image[row][col][1] = 0
                    image[row][col][2] = 0
        # cv2.imshow('black', image)
        # cv2.imwrite('black.jpg', image)
        hsv = cv2.cvtColor(image, cv2.COLOR_BGR2HSV)  # 转换为HSV空间图像,uint8
        # cv2.imshow('hsv', hsv)
        cv2.imwrite('hsv.jpg', hsv)
        a = 0
        b = 0
        c = 0
        color = []
        #   遍历图像求烟叶区域的平均色调和平均饱和度
        for row in range(rows):
            for col in range(cols):
                if hsv[row][col][1] > 0:
                    a += 1
                    b += hsv[row][col][1]
                    c += hsv[row][col][2]
        color.append(round(b/a, 2))
        color.append(round(c/a, 2))
        return color
    

    三、形状特征

    形状特征采用分散度:
    f = L 2 S f = \frac{L^2}{S} f=SL2
    轮廓周长为轮廓像素数量,面积为轮廓内像素数量:
    在这里插入图片描述

    #   计算图像的最大轮廓对应的分散度
    def fensandu(array):
        #   对图像进行二值化,阈值设置为30,阈值为观察灰度值所得
        ret, binaryImg = cv2.threshold(array, 30, 255, cv2.THRESH_BINARY)
        #   获取图像的所有轮廓
        contours, hierarchy = cv2.findContours(binaryImg, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
        #   设定一个面积初值
        a = 0
        area = cv2.contourArea(contours[0])
        #   遍历所有轮廓,获取最外围轮廓
        for c in range(len(contours)):
            #   获取最大轮廓的周长面积
            if cv2.contourArea(contours[c]) >= area:
                a = c
                area = cv2.contourArea(contours[c])
                arclen = cv2.arcLength(contours[c], True)
        #   计算分散度L*L/S,作为函数的返回值
        imgnew = cv2.drawContours(img, contours[a], -1, (0, 255, 0), 3)
        # cv2.imshow('lunkuo', imgnew)
        cv2.imwrite('lunkuo.jpg', imgnew)
        fensandu = round(arclen*arclen/area, 2)
        return fensandu
    

    四、利用KNN分类器进行数据分类

    1.数据预处理

    文件格式如下图所示,5种特征,最后一排为数据类别:
    在这里插入图片描述

    对TXT文件进行处理,把TXT文件分割成特征数组和类别数组,便于后续处理:
    # open_file.py
    import numpy as np
    
    
    def str_float(value):
    # 列表中字符串转化浮点数(读取txt文本可能存在空行,自定义避免报错)
        try:
            return float(value)
        except:
            return None
    
    
    def open_file(filename):
    # 读取txt文件 返回分类特征集和类别集
        fr = open(filename, encoding='utf-8')
        arrayOLines = fr.readlines()
        numberOfLines = len(arrayOLines)
        returnMat = np.zeros((numberOfLines, 5))  # 创建一个行*5的0矩阵
        classLabelVector = []
        index = 0
        # print(arrayOLines)
        for line in arrayOLines:
            line = line.strip()  # 去除每一行的空格或换行符
            listFromLine = line.split("\t")  # 按照回车分割每一行为每一个数组
            returnMat[index, :] = list(map(str_float, listFromLine[0:5]))  # 把每一行数据添加到之前创建的0矩阵之中
            classLabelVector.append((float(listFromLine[-1])))
            index += 1
    	return returnMat, classLabelVector
    

    2.进行数据的监督学习

    KNN进行数据分类

    from numpy import np
    
    # KNN分类算法函数定义
    def kNNClassify(newInput, dataSet, labels, k):
        numSamples = dataSet.shape[0]   # shape[0]表示行数
    
        # # step 1: 计算距离
        diff = tile(newInput, (numSamples, 1)) - dataSet  # 按元素求差值
        # print(diff)
        squaredDiff = diff ** 2  # 将差值平方
        squaredDist = sum(squaredDiff, axis=1)   # 按行累加
        distance = squaredDist ** 0.5  # 将差值平方和求开方,即得距离
        # print(distance)
        # # step 2: 对距离排序
        # argsort() 返回排序后的索引值
        sortedDistIndices = argsort(distance)
        # print(distance)
        # print(sortedDistIndices)
        classCount = {}  # define a dictionary (can be append element)
        for i in range(k):
            # # step 3: 选择k个最近邻
            voteLabel = labels[sortedDistIndices[i]]
            # print(sortedDistIndices[i])
            # # step 4: 计算k个最近邻中各类别出现的次数
            classCount[voteLabel] = classCount.get(voteLabel, 0) + 1
        # print(classCount)
        # # step 5: 返回出现次数最多的类别标签
        maxCount = 0
        for key, value in classCount.items():
            if value > maxCount:
                maxCount = value
                maxIndex = key
    
    
        return maxIndex
    
    
    展开全文
  • 实验室师兄写的代码,使用灰度共生矩阵计算图像的纹理复杂度,然后可以用来衡量嵌入水印的多少。
  • 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(...

    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

    展开全文

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,492
精华内容 596
关键字:

灰度共生矩阵纹理特征